Need to provide FOV mask as part of the input data

[KrisThielemans]

Dataset as returned by petric.get_data should have another member FOV_mask, which is a sirf.STIR.ImageData which is zero where the solution has to be zero, and 1 where the solution has to be non-negative.

[gschramm]

@KrisThielemans :related to that issue. what is the most efficient way in STIR to generate a binary STIR image that shows whether the binary STIR image is > 0? In numpy it would be sth like y = (x > 0) or y = (x > 0).astype(int)

[KrisThielemans]

Apologies for the delay on this issue. Our plan was to let this FOV mask correspond to wherever the provided OSEM is > 0, but I wanted to check how that looks like first...

Regarding the question, I'm assuming you mean SIRF (for sirf.STIR.ImageData). I think this might be one of our missing operations. It's generally quite hard to pass Python conditions to C++ of course. At present, I think we need to do something like

mask = im.minimum(im.allocate(eps))/eps

that wouldn't be very efficient though. @evgueni-ovtchinnikov @paskino @casperdcl any great ideas?

[gschramm]

(1) ok. just wanted to make sure that I was not missing a method of sirf.STIR.ImageData.
In the meantime, I will create use get_uniform_copy(0) and fill((x.as_array() > 0).astype(float))

(2) I guess using OSEM >0, instead of adjoint_ones > 0 for the fov mask should be fine since the OSEM
should not be 0 anywhere in the FOV.

[KrisThielemans]

Our plan was to let this FOV mask correspond to wherever the provided OSEM is > 0,

I have finally (!) checked that the above assumption isn't quite accurate: for some of the Siemens_mMR data, there are voxels which are non-zero in the reference image but are zero in the OSEM image. I believe that this is because of the penalty (in fact, OSEM can set voxels to zero "by accident" if the data is very noisy in a subset, but I don't think that's happening here).

However, all these voxels are in the end-planes, and are NOT in the evaluation VOIs. So, using OSEM_image.as_array() > 0 should not influence the evaluation results.

[KrisThielemans]

I have now checked that with the settings used for the projector

mask = STIR.TruncateToCylinderProcessor().process(OSEM_image.allocate(1))

gives identical results to adjoint_ones > 0, i.e.

acquisition_model = STIR.AcquisitionModelUsingParallelproj()
sens_mask = acquisition_model.backward(mult_factors)
sens_mask.fill(sens_mask.as_array() > 0)

Also, there are no non-zero voxels in either the OSEM or the reference image outside this mask.

[KrisThielemans]

Summary:

• easiest way to find the FOV mask (i.e. where non-zero voxels can occur) is to use

  mask = STIR.TruncateToCylinderProcessor().process(OSEM_image.allocate(1))

• the projector gives a "automatic" mask, and ignoring the FOV mask should work for most algorithms