making scatter work with blocks

[danieldeidda]

No description provided.

[KrisThielemans]

some comments as discussed. Also probably drop changes to inverse_SSRB as discussed (although we probably want to do that nearer to the end of this PR)

[danieldeidda]

I am a bit confused by the following:

STIR/src/scatter_buildblock/ScatterSimulation.cxx

Line 93 in 8ef0d10

error("ScatterSimulation: output projection data incompatible with what was used for set_up()");

Are we meant to provide a template which already describes the downsample scanner?

[KrisThielemans]

process_data() requires set_up() first. This has logic to create the template via downsampling (if you didn't give one yet I guess).

[danieldeidda]

process_data() requires set_up() first. This has logic to create the template via downsampling (if you didn't give one yet I guess).

OK so it should work even if we don't give a template in input? In my mind the template was the original scanner that then is downsampled

[KrisThielemans]

for cylindrical, the template is optional (as far as I recall), i.e. it will downsample then according to (optional) parameters.

It'd be nice to end up there in that situation for blocksoncylindrical as well, which I think you nearly had it.

Before @NikEfth created the downsampling functionality, we had to create the template by hand. Now it should be optional.

[danieldeidda]

for cylindrical, the template is optional (as far as I recall), i.e. it will downsample then according to (optional) parameters.

It'd be nice to end up there in that situation for blocksoncylindrical as well, which I think you nearly had it.

Before @NikEfth created the downsampling functionality, we had to create the template by hand. Now it should be optional.

It seems that if I don't give a template name then post_processing() which calls set_output_proj_data(const std::string& filename) throws the error:

set_output_proj_data(const std::string& filename)
{

    if(is_null_ptr(this->proj_data_info_cyl_noarc_cor_sptr))
    {
        error("ScatterSimulation: Template ProjData has not been set. Aborting.");
    }

I do have the downsample ready but it stops before it can even do the downsampling

[KrisThielemans]

Possibly that's a bug.then. I believe we do run ScatterSimulation without template in ScatterEstimation. Might need some debugging. sorry

[danieldeidda]

Possibly that's a bug.then. I believe we do run ScatterSimulation without template in ScatterEstimation. Might need some debugging. sorry

OK after some debugging I think for the ScattereEstimation it makes sense to allow no template because the info about the scanner are provided in that class and post_processing is not called, therefore we don't encounter that error. However, when we run the simulate_scatter utility on its own we need a template otherwise it does not have any clue about scanner geometry. I think the easiest solution would be to re-set output_proj_data after the downsampling in this way everything runs and the simulated scatter looks something reasonably alike to a cylinder

[KrisThielemans]

for the ScattereEstimation it makes sense to allow no template because the info about the scanner are provided in that class and post_processing is not called, therefore we don't encounter that error. However, when we run the simulate_scatter utility on its own we need a template otherwise it does not have any clue about scanner geometry.

yes, makes sense.

re-set output_proj_data after the downsampling in this way everything runs and the simulated scatter looks something reasonably alike to a cylinder

I'm 100% lost. reset where? ScatterEstimation ,but that works already? ScatterSimulation, but you just said it cannot work.

[danieldeidda]

I'm 100% lost. reset where? ScatterEstimation ,but that works already? ScatterSimulation, but you just said it cannot work.

In ScatterSimulation, we give a template which is the original scanner (out_projdata is set equal to that in post_processing()), then we downsample this scanner and change out_projdata info with the downsampled one.

[KrisThielemans]

ok. I think it does make sense to (optionally?) let ScatterSimulation do the downsampling, but then I also think it has to do the upsampling, i.e. output has to be the same size as the template. This would be good, but obviously out of scope for this PR.

[danieldeidda]

This seems to run OK. However, if we look at the downsample scatter and the upsample one in the figure below, we can see that the downsample seem to have some empty part which I am not sure is correct

[KrisThielemans]

who knows how Amide displays output of extract_segments. "image" dimensions might not make sense anyway. A better test is to compute the scatter with full sampling (in m) and compare that to the upsampled sino. Clearly, you will have to use something with asymmetry in z

[danieldeidda]

who knows how Amide displays output of extract_segments. "image" dimensions might not make sense anyway. A better test is to compute the scatter with full sampling (in m) and compare that to the upsampled sino. Clearly, you will have to use something with asymmetry in z

I am running simulate_scatter with original scanner. It has been running for 25h with openMP on. does this seem normal?

[KrisThielemans]

it will take a long time yes, certainly if you have it in full 3D. That's the whole reason we do the downsampling... You should be able to monitor progress in the log files. It'll give an estimate of run-time, which initially is too pessimistic as it ignores the benefit of the caching of "scatter-point <-> detector" lines, but it would still give you some idea.

[danieldeidda]

it will take a long time yes, certainly if you have it in full 3D. That's the whole reason we do the downsampling... You should be able to monitor progress in the log files. It'll give an estimate of run-time, which initially is too pessimistic as it ignores the benefit of the caching of "scatter-point <-> detector" lines, but it would still give you some idea.

I think this look weird, is it possible that the z_shift needs to be applied again?

[danieldeidda]

it will take a long time yes, certainly if you have it in full 3D. That's the whole reason we do the downsampling... You should be able to monitor progress in the log files. It'll give an estimate of run-time, which initially is too pessimistic as it ignores the benefit of the caching of "scatter-point <-> detector" lines, but it would still give you some idea.

I think this look weird, is it possible that the z_shift needs to be applied again?

while the following is the same but using Cylindrical geometry

[KrisThielemans]

I guess it'd be good to see coronal slice of the activity/attenuation images.

Might need a hunt through the ScatterSimulation code how it gets the coordinates, and in particular if there's something checks on cylindrical or not.

[danieldeidda]

I guess it'd be good to see coronal slice of the activity/attenuation images.

Might need a hunt through the ScatterSimulation code how it gets the coordinates, and in particular if there's something checks on cylindrical or not.

I had put this check for Cylindrical only but pheraps needs to be there and in fact detector_coord_A.z() (and B) is different than 0

STIR/src/scatter_buildblock/ScatterSimulation.cxx

Line 385 in bdaa466

if(this->proj_data_info_cyl_noarc_cor_sptr->get_scanner_sptr()->get_scanner_geometry()=="Cylindrical"){

STIR/src/scatter_buildblock/ScatterSimulation.cxx

Line 386 in bdaa466

assert(detector_coord_A.z() == 0);

STIR/src/scatter_buildblock/ScatterSimulation.cxx

Line 387 in bdaa466

assert(detector_coord_B.z() == 0);

before this line find_cartesian_coordinates_of_detection(detector_coord_A, detector_coord_B, Bin(0, 0, 0, 0)); is called which adds a z_shift further away from zero

[danieldeidda]

Might need a hunt through the ScatterSimulation code how it gets the coordinates, and in particular if there's something checks on cylindrical or not.

So I had a look at the code and there is a shift in the following:

STIR/src/scatter_buildblock/ScatterSimulation.cxx

Lines 404 to 405 in bdaa466

	this->shift_detector_coordinates_to_origin =
	CartesianCoordinate3D<float>(this->proj_data_info_cyl_noarc_cor_sptr->get_m(Bin(0, 0, 0, 0)), 0, 0);

when I change the sign of this shift I get something that looks more like the equivalent in the cylindrica geometry

[KrisThielemans]

unfortunately that makes little sense to me. You did a lot of hard work to make sure that find_cartesian* and get_m give the same results for cylindrical and block. That then implies that the shift has to be the same in both cases. So I'm not sure what's going on.

Possibly you could enable the debug code above, removing lines

STIR/src/scatter_buildblock/ScatterSimulation.cxx

Line 385 in bdaa466

if(this->proj_data_info_cyl_noarc_cor_sptr->get_scanner_sptr()->get_scanner_geometry()=="Cylindrical"){

and

STIR/src/scatter_buildblock/ScatterSimulation.cxx

Line 394 in bdaa466

if(this->proj_data_info_cyl_noarc_cor_sptr->get_scanner_sptr()->get_scanner_geometry()=="Cylindrical")

as those were appropriate before we shifted everything back.

[danieldeidda]

unfortunately that makes little sense to me. You did a lot of hard work to make sure that find_cartesian* and get_m give the same results for cylindrical and block. That then implies that the shift has to be the same in both cases. So I'm not sure what's going on.
Possibly you could enable the debug code above, removing lines

STIR/src/scatter_buildblock/ScatterSimulation.cxx

Line 385 in bdaa466

if(this->proj_data_info_cyl_noarc_cor_sptr->get_scanner_sptr()->get_scanner_geometry()=="Cylindrical"){

and

STIR/src/scatter_buildblock/ScatterSimulation.cxx

Line 394 in bdaa466

if(this->proj_data_info_cyl_noarc_cor_sptr->get_scanner_sptr()->get_scanner_geometry()=="Cylindrical")

as those were appropriate before we shifted everything back.
L394 I have already removed and works fine but the 385 will fail. the z position appears to be -axial_length

or pheraps we could have a check like in L385 and if is not zero calculate the difference and force it zero? I know it sounds ugly even while I am writing :)

[KrisThielemans]

yup. ugly 😁. My interpretation is that this check has to be satisfied, as otherwise there is something wrong with the coordinate systems. No?

By the way, I feel there's no need to quote whole messages if the reply is "in sequence".

[danieldeidda]

This now seems to be working. Although I am a bit confused with the sampling in particular the line below:

STIR/src/include/stir/ProjDataInfoGeneric.inl

Line 96 in 0b9c81d

return get_scanner_ptr()->get_ring_spacing()/2; // TODO currently restricted to span=1 get_num_axial_poss_per_ring_inc(segment_num);

to me this should not be divided by 2

[KrisThielemans]

./run_scatter_tests.sh is causing time-outs (run time of more than 5 hours)

[KrisThielemans]

The "2d" sino output is not SSRB'd.

[KrisThielemans]

looks like it's still taking too much time. Something else I guess. Did you run the test locally?

[danieldeidda]

not yet. Still waiting for my computer

[KrisThielemans]

My suggestions cause the tests to fail with

ERROR: ProjDataInfoCylindrical::get_ring_pair_for_segment_axial_pos_num does not work for this type of sampled data

which is thrown in get_ring_pair_for_segment_axial_pos_num

STIR/src/buildblock/ProjDataInfoCylindrical.cxx

Lines 338 to 339 in aaa0f8e

	if (!sampling_corresponds_to_physical_rings)
	error("ProjDataInfoCylindrical::get_ring_pair_for_segment_axial_pos_num does not work for this type of sampled data");

So now I'm confused. This code needs segment_axial_pos_to_ring1_plus_ring2 which previously we did initialise, but potentially wrongly, as it used m_offset etc, which does not necessarily make sense for the block scanner... There seems to be something deeper that needs fixing.

[danieldeidda]

I am kind of lost here :D

[KrisThielemans]

We will fix the initialisation of the ring diff arrays in a future PR.