Motion corrected parametric objective function

[rijobro]

Requires on a lot of the other recent pull requests (i.e., reading of multicomponent images and projectors that allow motion.

This will fail all tests, as it hasn't been rebased onto the other pull requests.

[AnderBiguri]

I don't seem to be able to get it working. Compilation errors are confusing (to me) but I think that what they are saying is that this derived class does not match the base class, which is what I was suspecting.

/home/anderbiguri/SIRF/build/sources/STIR/src/include/stir/RegisteredParsingObject.inl:40:23: error: invalid new-expression of abstract class type ‘stir::PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<stir::ParametricDiscretisedDensity<stir::VoxelsOnCartesianGrid<stir::KineticParameters<2, float> > > >’
   Derived * der_ptr = new Derived;
                       ^~~~~~~~~~~
In file included from /home/anderbiguri/SIRF/build/sources/STIR/src/recon_buildblock/recon_buildblock_registries.cxx:60:0:
/home/anderbiguri/SIRF/build/sources/STIR/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData.h:58:7: note:   because the following virtual functions are pure within ‘stir::PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData<stir::ParametricDiscretisedDensity<stir::VoxelsOnCartesianGrid<stir::KineticParameters<2, float> > > >’:
 class PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData:
       ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
In file included from /home/anderbiguri/SIRF/build/sources/STIR/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.h:32:0,
                 from /home/anderbiguri/SIRF/build/sources/STIR/src/recon_buildblock/recon_buildblock_registries.cxx:26:
/home/anderbiguri/SIRF/build/sources/STIR/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMean.h:266:5: note: 	stir::Succeeded stir::PoissonLogLikelihoodWithLinearModelForMean<TargetT>::set_up_before_sensitivity(const std::shared_ptr<const _Tp>&) [with TargetT = stir::ParametricDiscretisedDensity<stir::VoxelsOnCartesianGrid<stir::KineticParameters<2, float> > >]
     set_up_before_sensitivity(shared_ptr<const TargetT > const& target_sptr) = 0;
     ^~~~~~~~~~~~~~~~~~~~~~~~~

So this is in a point where I don't know how to tackle it. I assume this needs to comply with the Base class to be able to work seamlessly within STIR, but I do not know then how to pass the information that is needed there ti that ExamInfo

[rijobro]

Ok some ideas:

1. Change all instances of set_up_before_sensitivity. Probably a lot of changes.
2. Do a naughty const cast (easiest, not good practice)
3. Make the exam info to be a mutable member so that it can be changed for const images
4. Make the time definitions a mutable object so that it can be changed for const exam info
5. Copy the original image, so that it's no longer const. Probably not feasible since the rest of the recon will use the original target_sptr.

[AnderBiguri]

Thanks @rijobro, those are great suggestions!

I guess I'll wait to see what @KrisThielemans prefers :)

[rijobro]

You could probably just comment it out for now. I'm not sure, but I don't think it'll be using the target image's time defs during the recon. It'll be updating it so that when it gets saved to interfile, it shows that the image is a composite of multiple time points. You'd have to check though, but I imagine that the Patlak analysis isn't currently clever enough to be using the time def inside of the target image.

[AnderBiguri]

I was working on Windows for a bit and I moved to Linux now (Windows was acting up a bit) and I get the following error when compiling (after commenting out the problematic line above, as Rich suggested)

[ 14%] Linking CXX executable ctac_to_mu_values
/usr/bin/ld: cannot find -llocal_motion_buildblock
collect2: error: ld returned 1 exit status

I see that -llocal_motion_bildblock points to -llocal_motion_buildblock ../modelling_buildblock/libmodelling_buildblock.a , and that that is indeed correctly at (with the superbuild) ./builds/STIR/build/src/modelling_buildblock/libmodelling_buildblock.a.
any idea?

[KrisThielemans]

@AnderBiguri for your original problem. We can indeed not let this function modify the *target_sptr. There is no guarantee that the user will call set_up() with the same image, as what they will reconstruct for instance. So, this setting should be converted into a check if the time-frame-defs are as expected. I guess I'd do this as

• if it has no time-frame defs, issue warning
• else if it has the wrong time-frame defs, issue error

There should normally be a construct_target member somewhere, which will create the default target (given voxel-sizes etc) for the input data. The set code should be moved there.

[KrisThielemans]

The local_motion_buildblock error could be related to https://github.com/UCL/STIR/pull/240/files#r458742649, but more likely is because you haven't enabled STIR_ENABLE_EXPERIMENTAL=ON.

Obviously, that means we cannot merge this PR until we have moved the relevant motion bits.

[AnderBiguri]

• if it has no time-frame defs, issue warning

Again my C++ lacking. How do I check this? by comparing the getter return to an "empty" (non-set) TimeFrameDefinitions ?

[KrisThielemans]

This doesn't have anything (or at least very little) to do with C++ luckily.

...->get_time_frame_definitions().get_num_time_frames()

or something with a similar name

[AnderBiguri]

Apologies, I know how to get the time frames. But you said "if it has no time-frame defs, issue warning". The class will have time-frame defs as a member, and such member will be created by its default constructor when ExamInfo default constructor is called.

My (C++) question is how do I know ...->get_time_frame_definitions().get_num_time_frames() is "empty" i.e. "it has no time frame defs"

[KrisThielemans]

my wrong terminology, I just meant ...->get_time_frame_definitions().get_num_time_frames()==0

[AnderBiguri]

Thanks! I see now! ah, always the easiest solutions after an hour of reading half of stackoverflow.... Thanks :)

[rijobro]

You could check more than just the number of time frames (start and stop times of each frame). Not sure if time frame definitions class overloads the == and != operators, but that would be pretty easy.

[AnderBiguri]

Yes they do! Checking the number of tf is just to see if they are set, but if there is more then zero, then I'm comparing them properly with ==

[AnderBiguri]

I managed to make it work, but I get images out of this that are completely zero in value as a recon result.

I loaded up my own dynamic images that I am using with indirect patlak for my research, and have forward projected them. I generated multiplicative sinos that are of value 1 and additive sinos of value 0. I have created deformation fields filled with value zero too, to test if the code works without motion.

I've been looking at the pieces of this for a bit now and I have no idea what can make the result all zero. Do you have any idea what to look at that I might have missed?

I have checked that all sinos have "correct values (sino has sinograms, mult has all 1, add has all 0 etc).

Its quite hard in STIR to make a minimal complete example, because there are many files needed for that, so let me know what extra info I can add here:

Code that generates data, in SIRF.

import numpy as np
import scipy
import scipy.io as sio
import scipy.ndimage as snd
import matplotlib.pyplot as plt
import sirf.STIR as stir
import sirf.Reg as reg

from sirf.Utilities import  show_3D_array
import subprocess
import os
import sys


def save_multi(time,fname,prefix,extension=".hs"):
    data=[""]*(len(time)+3)
    data[0]="Multi :=\n"
    data[1]="  total number of data sets := "+str(len(time)-1)+"\n"
    for ii in range(len(time)-1):
        data[ii+2]="  data set["+ str(ii+1) + "] := ./sinos/"+prefix+str(ii)+extension+"\n"
    data[len(time)+2]="end :="
    with open("./sinos/"+fname+".txt","w") as outputfile:
        outputfile.writelines(data)


def save_fdef(filename,time,bed_frame_duration):
    data=[]
    data.append("0 " +str(time[0])+"\n")
    for i in range(len(time)-1):
        data.append("1 "+str(bed_frame_duration)+"\n")
        data.append("0 "+str(time[i+1]-time[i]-bed_frame_duration)+"\n")
    with open("./sinos/"+filename+".fdef","w") as outputfile:
        outputfile.writelines(data)


# dynamic sim parameters. 

filename_prefix="dynamicXcat_WB"
bed_frame_duration=35
wbd_pass_time=6*bed_frame_duration
start_time=3*wbd_pass_time
end_time=wbd_pass_time*8 #60*70
eps=0.0000001
desired_t=np.arange(start_time,end_time,wbd_pass_time) 

save_fdef("test",desired_t,bed_frame_duration)

# load dynamic images
dynamic_image=[]
for i in range(len(desired_t)-1):
    dynamic_image.append(stir.ImageData("./dynamicXcat_WB/"+filename_prefix+"_"+str(desired_t[i])+"_"+str(desired_t[i]+bed_frame_duration)+".hv"))

#create acq model. 
acq_data_template=stir.AcquisitionData('Discovery MI')
acq_model = stir.AcquisitionModelUsingRayTracingMatrix()
acq_model.set_up(acq_data_template,dynamic_image[0])


# make DVFs
# as I donkt know how to create a DVF from an image template, I register nothing at the correct geometry and then fill it with zero
empty=acq_data_template.create_uniform_image(0)
reg_alg = reg.NiftyF3dSym()
reg_alg.set_initial_affine_transformation(reg.AffineTransformation(np.eye(4, dtype=np.float32)))
reg_alg.set_reference_image(empty)
reg_alg.add_floating_image(empty)
reg_alg.process()
dvf=reg_alg.get_deformation_field_forward()
dvf.fill(0)

for i in range(len(desired_t)-1):
    dvf.write("./sinos/dvf_"+str(i)+".nii")
save_multi(desired_t,"dvf_","dvf_",".nii")
sinos=[]

#make sinograms
for i in range(len(desired_t)-1):
    sinos.append(acq_model.forward(dynamic_image[i]))

# save sino, additive, multiplicative 
for i in range(len(desired_t)-1):
    sinos[i].write("./sinos/sino_"+str(i)+".hs")
    additive=sinos[i].get_uniform_copy(0)
    multiplicative=sinos[i].get_uniform_copy(1)
    additive.write("./sinos/add_"+str(i)+".hs")
    multiplicative.write("./sinos/mult_"+str(i)+".hs")
save_multi(desired_t,"sino_","sino_")
save_multi(desired_t,"add_","add_")
save_multi(desired_t,"mult_","mult_")

Parameter file:

OSMAPOSLParameters :=
; example file for running OSEM with interfiltering
; filter specific things are at the end

objective function type:=PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData
PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData Parameters:=

input file := ./sinos/sino_.txt

; if disabled, defaults to maximum segment number in the file
maximum absolute segment number to process := 0
; see Users Guide to see when you need this
zero end planes of segment 0:= 0

;zoom := .15645
;xy output image size (in pixels) := 21
;z output image size (in pixels) := 17
;z offset (in mm) := 0

projector pair type := Matrix
  Projector Pair Using Matrix Parameters :=
  Matrix type := Ray Tracing
  Ray tracing matrix parameters :=
  	number of rays in tangential direction to trace for each bin := 10
  End Ray tracing matrix parameters :=
  End Projector Pair Using Matrix Parameters :=

Dynamic bin normalisation filename := ./sinos/mult_.txt

; we need this for backwards compatibility with the testing script
use subset sensitivities:=1
;sensitivity filename:= sens.img
; if next is set to 1, sensitivity will be recomputed
; and written to file (if "sensitivity filename" is set)
recompute sensitivity := 1
subset sensitivity filenames := subset_sensitivity_%d.hv

; specify additive projection data to handle randoms or so
; see Users Guide for more info
additive sinograms := ./sinos/add_.txt

; patlak related files
Kinetic Model type := Patlak Plot
Patlak Plot Parameters:=

time frame definition filename := ./sinos/test.fdef
starting frame := 1
calibration factor := 1
blood data filename := if_fdg_w_decay.if 
Time Shift := 0
In total counts := 1
In correct scale := 1
end Patlak Plot Parameters:=

Forward displacement field image := ./sinos/dvf_.txt

end PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData Parameters:=

; Number of subsets should be a divisor of num_views/4
number of subsets:= 68
; Use for starting the numbering from something else than 1
start at subiteration number:=1
; Use if you want to start from another subset than 0 (but why?)
start at subset:= 0
number of subiterations:= 7
save estimates at subiteration intervals:= 7
;write update image := 1

;initial estimate := 0
; enable this when you read an initial estimate with negative data
enforce initial positivity condition:=1

inter-update filter subiteration interval:= 0
inter-update filter type := None

inter-iteration filter subiteration interval:= 0

post-filter type := None

output filename prefix := test

END :=

[AnderBiguri]

Update on the above: the sensitivity is all zeros, and that is what is causing the output to be all zeroes.

I am trying to follow who does this subset computing, to see what happens, but I am a bit lost of where things are happening.

[KrisThielemans]

can you try this without def fields first? If that works, can you make sure that the def fields are written as displacement fields? (then they should just contain zeros).

There could be with writing the dvf to file, and STIR reading them properly (x,y,z quickly starts to loose meaning). It should work of course, but...

[KrisThielemans]

Also, does the recon_test_pack/test_modelling.sh still work?

[AnderBiguri]

can you try this without def fields first? If that works, can you make sure that the def fields are written as displacement fields? (then they should just contain zeros).

Without them "it works" (I am running 1 subset with 1 iteration, so "it works" means there is a blob in the image).

With displacements instead of deformations (true, that was my mistake), the result is the same. The only thing I chahnged from the above code is dvf=reg_alg.get_deformation_field_forward() -> dvf=reg_alg.get_displacement_field_forward()

Also, does the recon_test_pack/test_modelling.sh still work?

There is no such file. You mean recon_test_pack/run_tests_modelling.sh?
I get errors/sysntax errors running this, investigating a bit.

(fixed, apparently sh something.sh and ./something.sh are not as equivalent as I thought)

[AnderBiguri]

Also, does the recon_test_pack/run_tests_modelling.sh still work?

Errors!

Test the direct POSMAPOSL Patlak Plot reconstruction
./run_tests_modelling.sh: line 187: 22706 Segmentation fault      (core dumped) INPUT=fwd_dyn_from_p0005-p5.hs INIT=indirect_Patlak.hv ${MPIRUN} P${direct} P${direct}.par > P${direct}.log 2>&1
ERROR

[KrisThielemans]

With displacements instead of deformations (true, that was my mistake), the result is the same

the same as what? Does it now "work" or still zero sensitivity?

You mean recon_test_pack/run_tests_modelling.sh?

yes, sorry

I get errors/sysntax errors running this, investigating a bit.

a bit strange to get syntax errors, as Travis runs it, and it's ok on master apparently. Maybe something went wrong with the merge

[KrisThielemans]

(fixed, apparently sh something.sh and ./something.sh are not as equivalent as I thought)

they aren't if something.sh syarts with #! /bin/bash or anything else that doesn't run sh.

./run_tests_modelling.sh: line 187: 22706 Segmentation fault

oh well. Best thing to do is build with RelWithDebInfo , and change the script to use gdb --args /path/to/.../POSMAPOSL whatever.par. Then run, in gdb do continue, wait till it crashes, and then info stack etc

[AnderBiguri]

With displacements instead of deformations (true, that was my mistake), the result is the same

the same as what? Does it now "work" or still zero sensitivity?

Returns zero images. The subset_sensitivity%d.hv images (when extracted, they are parametric) are zero.
The frame sensitivities (stored by adding write_to_file("dyn_sensitivity_frame0.hv",dyn_sensitivity[this->_patlak_plot_sptr->get_starting_frame()]); inside PoissonLogLikelihoodWithLinearKineticModelAndDynamicProjectionData::add_subset_sensitivity) is also zero (but interestingly has "number of frames =4" inside the .hv header, even if its only 1, I need to remove it to be able to see it in amide).

oh well. Best thing to do is build with RelWithDebInfo , and change the script to use gdb --args /path/to/.../POSMAPOSL whatever.par. Then run, in gdb do continue, wait till it crashes, and then info stack etc

Ok, step by step.

RelWithDebInfo

Im in Linux, how do I do this? -DCMAKE_BUILD_TYPE:STRING="RelWithDebInfo" ?

change the script to use gdb --args /path/to/.../POSMAPOSL whatever.par.

I have no idea what you mean. Change which script, the run_tests_modelling? The line that errors?
After some checking it seems that the line is calling POSMAPSOL in the most obfuscated way possible, you mean just running gdb on that particular call?

[KrisThielemans]

The frame sensitivities ... is also zero

well, that's a nice bug then to look into...

I'll have to let you figure some of this out, but what you say is in the correct direction. See e.g. https://cmake.org/cmake/help/latest/variable/CMAKE_BUILD_TYPE.html. Obviously, you don't have to use gdb, if you prefer another debugger. And maybe you're using a nice IDE. There's some info on http://stir.sourceforge.net/registered/ amazingly enough.

[AnderBiguri]

The frame sensitivities ... is also zero

well, that's a nice bug then to look into...

Yes! But I don't know where in this code are the calls to the frame sensitivities.... Been looking for a while. Thus the comment above (somewhere)

Plus is strange that without motion they are not (I assume, as it does not return all zero for when I call this without DVFs)

[rijobro]

A small disclaimer: this wasn't working before I left, thus why we hadn't merged it. Unfortunately, I can't remember what wasn't working...

Is there a unit test for running the non-rigid displacement (would have to be using the experimental version)? If so, is it currently being tested with the multi file?

[AnderBiguri]

A small disclaimer: this wasn't working before I left, thus why we hadn't merged it. Unfortunately, I can't remember what wasn't working...

Ah! That makes sense and changes a bit how I am testing it, as I keep trying to find bugs in my code (as its 99.9% of the times, my fault, maybe this is the 0.1%!). Its good to know, thanks! I'll check what happens on those DVFs, as they seem to be causing some issues :)

[rijobro]

Apologies, should have flagged that up sooner! But really I can't remember which component(s) were causing issues.

[AnderBiguri]

Apologies, should have flagged that up sooner! But really I can't remember which component(s) were causing issues.

No worries! There were many things that I have being setting wrong until now, so its the right moment to flag that up 👍

[KrisThielemans]

Probably the same problem as #580 (comment)