alternative MCIR implementation

[KrisThielemans]

This is the experimental version adjusted such that it doesn't conflicts. This verison uses adjoint interpolation as opposed to inverse.

This PR is not really for merging, but for testing.

[KrisThielemans]

@rijobro, I made a PR anyway to see if it compiles

[rijobro]

Thanks. I'm trying to figure out which arguments would need to be exposed to SIRF:

https://github.com/KrisThielemans/STIR-1/blob/420249b3ea3ee394cf1b70c4f11dafd850b6fdc9/src/experimental/motion/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion.cxx#L123-250

From the code snippet, it seems:

• GatedProjData for sinos,
• GatedProjData for additive sinos,
• vector of transformations, and
• vector of BinNormalisations.

How do I add attenuation? As part of the projector pair? I'm slightly surprised that there are multiple BinNormalisations, as I would have expected them to be constant across different gates.

[KrisThielemans]

Seems we need to add set_forward_transformation_sptrs and set_normalisation_sptrs. Currently internal member is VectorWithOffset but could make that std::vector of course.

How do I add attenuation? As part of the projector pair? I'm slightly surprised that there are multiple BinNormalisations, as I would have expected them to be constant across different gates.

This is the same as for your CIL stuff, the stir:::BinNormalisation corresponds to sirf::stir::AcquisitionSensitivityModel. It encodes attenuation as well. Also, we need to stick gate_duration factors in the BinNormalisation.

This would look a bit different from normal SIRF stuff where everything sits in the AcquisitionModel. Here we'd pass a (purely geometric) sirf::stir::AcquisitionModel (which has a ProjectorPair I guess), and a vector of sirf::stir::AcquisitionSensitivityModels, (which have stir:::BinNormalisations).

I don't think you have to expose GatedProjData to SIRF, could just use a vector, and internally convert.

However, It seems you'd need to expose the STIR resampling dataprocessor to SIRF as well. That won't be the same as your NiftyReg resampler. That'd mean converting the NiftyReg def fields to STIR def fields. hmmmm.

[rijobro]

However, It seems you'd need to expose the STIR resampling dataprocessor to SIRF as well. That won't be the same as your NiftyReg resampler. That'd mean converting the NiftyReg def fields to STIR def fields. hmmmm.

I was thinking about that. STIR is capable of reading nifti files with ITK, so I think this conversion would be possible (since in my kinetic MCIR stuff we were doing that, albeit passing via files). However, it all seems like quite a lot of work...

[rijobro]

As well as registry, library CMake changes I had to add/change a few methods in order for the class to not be abstract.

[rijobro]

Two findings:

1. I think cache_transformed_coords is necessary for reasonable reconstruction times
2. As far as I can see, the reading of tensor nifti images (4D) for the displacement, results in reconstructed images containing only zeroes. Splitting those displacements into their x-, y- and z-components seems to solve the problem.

The latter is strange, since I believe I used tensor niftis for POSMAPOSL, e.g., here #240.

[KrisThielemans]

https://travis-ci.org/github/UCL/STIR/jobs/699793304#L4732-L4733

CMake Error at src/experimental/CMakeLists.txt:9 (list):
  list does not recognize sub-command PREPEND

Seems PREPEND was introduced in CMake 3.15

We should really build experimental also in most of the others on Travis and Appveyor.

[KrisThielemans]

Travis/Appveyor are now fine with more runs enabled. Remaining failure should be fixed once merging #591. (I did put that on release_4 though)

[KrisThielemans]

I've put in a number of PRs onrelease_4. Once merged, we should merge again to master`, and then that into here. Then should be good to merge (with an addition to release_5.htm)

[KrisThielemans]

Pending #599 and I guess #587. @rijobro, can you take care of this? (Don't forget to add something in release_5.htm

[rijobro]

Can we replace the New part of the class name with Adjoint?

[KrisThielemans]

Can we replace the New part of the class name with Adjoint?

good idea, but I find PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotionAdjoint confusing, but my alternative is even longer...
PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotionUsingAdjoint?

If you like that, we'd want to rename the source files as well then.

Right now, I'd still keep it in experimental. Would want to get rid if the explicit lists of bin-normalisations etc in the par file first.

let's make the caching default to true.

[rijobro]

At this point, I think keeping class names short has gone out the window, so the latter is fine with me!

let's make the caching default to true.

Yup.

[KrisThielemans]

At this point, I think keeping class names short has gone out the window, so the latter is fine with me!

😄

I merged release_4 to master (#599), but then noticed that you put #587 on release_4. However, unfortunately, it breaks backwards compatibility, and also means we have to merge again.
Would it be easy for you to revert #587 and put it in master instead? If not, we'll live with it.

apologies that I didn't notice.

[rijobro]

Reverted from release_4 and new PR for GatedProjData at #600.

[KrisThielemans]

@rijobro, could you merge master in here? There will be some conflicts.

[rijobro]

Done. Although best to hold off merging for now, still getting all-zero reconstructions for when using motion (even zero motion). Not sure why, I used the same non-rigid class with the kinetic stuff.

[rijobro]

Really lost now. Results are non-zero when no motion is used. Yet, when I use motion (even set to identity rigid matrix, see below), result is zero...

shared_ptr<RigidObject3DTransformation> rigid_sptr = MAKE_SHARED<RigidObject3DTransformation>();
this->_forward_transformations.at(1).reset(new Transform3DObjectImageProcessor<float>(rigid_sptr));

[KrisThielemans]

weird. I did a simple test that uncovered some end-plane problems, see #601, but that cannot cause this.

Are you sure there isn't anything in the rest of your data?

What do you mean "when no motion is used". how did you run that then?

I guess you could create a TrivialDataProcessor (with the adjoint stuff), that does nothing as another check.

[rijobro]

I checked by creating a displacement field image, filled with zeros:

im = reg.ImageData('recon.nii')
disp = reg.NiftiImageData3DDisplacement()
disp.create_from_3D_image(im)
disp.fill(0)
disp.write("no_motion_disp")

I used an example of the output reconstructed image to make sure that metadata matched (just incase there was a problem there).

[rijobro]

I also used:

disp.write_split_xyz_components("no_motion_split_%s")

to check whether the problem was in the reading of 4D displacement fields. But it still doesn't work, implying the error isn't there.

[KrisThielemans]

would it be easy to create a minimum test script (no attenuation, norm), as in recon_test_pack/test_simulate_and_recon.sh but just add the different .par file and no motion. Then I can pick this up a bit easier.

Otherwise, write some images after the transforms I guess to see what happens.

[rijobro]

Ok, I've created a test script that can be run with: sh ./run_test_simulate_and_recon_with_motion_using_adjoint.sh.
If you comment out the motion component in OSMAPOSL_with_motion_test_using_adjoint.par, the resulting image will be non-zero (check for min and max in current estimate... at bottom of reconstruction), whereas with the motion component uncommented, you should see min and max in current estimate 0 0.

To be clear, when i say "comment out the motion", I mean this bit:

transformation type for gate 1 := transformation
Transformation Parameters :=
cache_transformed_coords := 1
Transformation Type := BSplines transformation
	BSplines Transformation Parameters :=
		deformation field from file x-component := my_uniform_cylinder_motion_3d.hv
		deformation field from file y-component := my_uniform_cylinder_motion_3d.hv
		deformation field from file z-component := my_uniform_cylinder_motion_3d.hv
	End BSplines Transformation Parameters :=
End Transformation Parameters :=

[rijobro]

Any reason why I get the following?

stir_math --times-scalar 0 my_uniform_cylinder_motion_3d.hv my_uniform_cylinder.hv
list_image_info my_uniform_cylinder_motion_3d.hv 

Output:

....
Image min: 0
Image max: 1
Image sum: 118422

EDIT: turned out it needed --include-first. Confusing.

[KrisThielemans]

The sensitivity image is zero. Once that happens, the recon will just be zero. So the question is why is the sensitivity zero.

I notice that the code still uses the obsolete PresmoothingForwardProjectorByBin and PostsmoothingBackProjectorByBin. I wonder if there's a problem there. Could you replace them?

[rijobro]

I thought about replacing them, but since each objective function needs its own projector pair, I think we would need the equivalent of a create_shared_copy-esque method. In this way, we could get a copy for each motion state and apply the data processor. However, the copy method doesn't exist, so there's no way of getting a copy of the original projector pair. Or am I wrong?

[rijobro]

i.e., this:

shared_ptr<ForwardProjectorByBin> forward_projector_sptr_this_gate(
		 new PresmoothingForwardProjectorByBin(this->projector_pair_ptr->
						       get_forward_projector_sptr(),
						       this->_forward_transformations[gate_num]) );

Would have to become something like this:

shared_ptr<ForwardProjectorByBin> forward_projector_sptr_this_gate = this->projector_pair_ptr->
						       get_forward_projector_sptr()->get_shared_copy();
forward_projector_sptr_this_gate->set_pre_data_processor(this->_forward_transformations[gate_num]);

[rijobro]

I also notice that the presmoothingforwardprojector doesn't copy the original projector, which could lead to the sort of troubles we've been resolving recently with shared_ptr<const Obj>. Can open an issue if you'd like.

[KrisThielemans]

ah. this might indeed be a case where we need the PresmoothingForwardProjectorByBin and PostsmoothingBackProjectorByBin after all. In fact, I don't see how have create_shared_copy would help, nor do see another way to create a new projector that adds a pre-processor from an existing one. (without creating potentially a lot of overhead, as we do want these to share the cache etc) We'll then have to rename them to something more sensible.

For debugging, maybe you could just use set_pre_data_processor etc for the one gate. At least we'd know if the bug sits in the PostsmoothingBackProjectorByBin (sensitivity doesn't use the forward projector)

[rijobro]

Good spot. Modifying the following:

              shared_ptr<BackProjectorByBin> back_projector_sptr_this_gate
        (new PostsmoothingBackProjectorByBin(this->projector_pair_ptr->
                             get_back_projector_sptr(),
                             transpose_transformer_sptr)
         );

to

              shared_ptr<BackProjectorByBin> back_projector_sptr_this_gate
                      = this->projector_pair_ptr->get_back_projector_sptr();
              back_projector_sptr_this_gate->set_post_data_processor(transpose_transformer_sptr);

results in a non-zero sensitivity image. Guess that suggests the problem is in the PostSmoothing projector.