Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Refactoring for origin shift #260

Open
wants to merge 7 commits into
base: master
Choose a base branch
from

Conversation

ashgillman
Copy link
Collaborator

@ashgillman ashgillman commented Oct 8, 2018

  • Get comments
  • Test scatter
  • Rebase away the first commit - it was just for easier building and testing.

This PR Shouldn't actually change any functionality, just using the correct routines where they should be.

NB: I haven't had time to test changes to scatter.

@ashgillman ashgillman added this to In progress in Image Geometry via automation Oct 8, 2018
@ashgillman ashgillman force-pushed the prepare-for-origin-shift branch 2 times, most recently from cc474c6 to 0365f96 Compare October 11, 2018 05:28
@KrisThielemans
Copy link
Collaborator

KrisThielemans commented Oct 14, 2018

Not sure about DiscretisedDensity know about gantry coordinates (e.g. get_index_coordinates_for_gantry_coordinates). It means a rather tight coupling between images and projection data. I feel that ProjDataInfo should take care of gantry<->bed stuff.

Similarly, ProjDataInfo::get_centre_of_gantry_vector_in_relative_coordinates doesn't make sense, as ProjDataInfo doesn't know what "relative coordinates" are.

Therefore I would not have these functions. Anyone who needs to do this will just have call 2 functions (ProjDataInfo::get_centre_of_gantry_in_bed_coordinates and the relevant DiscretisedDensity::*from_physical_coordinates). What do you think?

Copy link
Collaborator

@KrisThielemans KrisThielemans left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks quite reasonable. Of course I have some suggestions...

Also a minor one: the comments here are obsolete?

->get_index_coordinates_for_gantry_coordinates
(end_of_tor_projected_to_axis_without_offset, proj_data_info_sptr);

const float first_lor_relative_to_tor
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

rename to z_index_start_of_tor? Not sure why the relative_to_tor? (but in the above I don't bother naming it...)

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm still not 100% sure about the add_adjacent_z function, and find some of the variable naming in it a bit confusing. So in this section of code I have been careful to replicate functionality, but must admit I don't 100% understand what is happening, hence more bad variable names.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe

const float z_of_first_voxel_relative_to_start_of_tor                                           
  = z_of_first_voxel_in_index_space       
  - start_of_tor_projected_to_axis_without_offset.z();
= end_of_tconst·float·z_of_end_of_tor_relative_to_start_of_tor
  or_projected_to_axis_without_offset.z()
  - start_of_tor_projected_to_axis_without_offset.z();

?


add_adjacent_z(lor, z_of_first_voxel - left_edge_of_TOR, right_edge_of_TOR -left_edge_of_TOR);
{
CartesianCoordinate3D<float>
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is a bit verbose.... just as well say

const float z_of_first_voxel_in_index_space = static_cast<float>(lor.bin()->coord1());

I'd rename that variable to z_index_of_first_voxel

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe better is z_index_of_traced_lor, with a comment explaining that this is constant because tantheta=0?


// In gantry space
CartesianCoordinate3D<float>
start_of_tor_projected_to_axis_without_offset
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is a somewhat confusing variable. It works because tantheta=0, and therefore x,y are irrelevant, but it's quite hard to understand what's going on here. (it was simpler before). I suggest to introduce a function

// The LOR is parametrised in terms of 'a'
CartesianCoordinate3D<float>
 find_image_index_along_LOR(
                  const float s_in_mm, const float m_in_mm, 
                  const float cphi, const float sphi, 
                  const float tantheta, 
                  const float a,
                  const DiscretisedDensity<3,float> &  density_info,
                 const ProjDataInfo& proj_data_info);

Use this in ray_trace_one_lor and use it here as

// some doc why a=0
const float a=0;
CartesianCoordinate3D<float>image_index_start_of_tor =
 find_image_index_along_LOR(s_in_mm, m_in_mm - sampling_distance_of_adjacent_LORs_z/2...);
...
add_adjacent_z(lor, z_index_of_first_voxel - image_index_start_of_tor.z(), 
                          image_index_end_of_tor.z() - image_index_start_of_tor.z());

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good idea :)

I went a little further and split out a local function

CartesianCoordinate3D<float>
get_point_on_lor_in_index_coordinates
(const float s_in_mm, const float m_in_mm, const float a_in_mm,
 const float cphi, const float sphi, const float tantheta,
 const DiscretisedDensity<3, float>& density_info,
 const ProjDataInfo& proj_data_info)
{
  return density_info.get_index_coordinates_for_relative_coordinates
    (proj_data_info.get_relative_coordinates_for_gantry_coordinates
     (proj_data_info.get_point_on_lor_in_gantry_coordinates
      (s_in_mm, m_in_mm, a_in_mm, cphi, sphi, tantheta)));
}

and method

CartesianCoordinate3D<float>
ProjDataInfo::get_point_on_lor_in_gantry_coordinates
(const float s_in_mm, const float m_in_mm, const float a_in_mm,
 const float cphi, const float sphi, const float tantheta) const
{
  return CartesianCoordinate3D<float>(m_in_mm      - a_in_mm*tantheta, // z
                                      s_in_mm*sphi - a_in_mm*cphi,     // y
                                      s_in_mm*cphi + a_in_mm*sphi);    // x
}

const float offset_in_z,
const float fovrad_in_mm,
const CartesianCoordinate3D<float>& voxel_size,
const bool restrict_to_cylindrical_FOV,
const int num_LORs)
const int num_LORs,
const shared_ptr< DiscretisedDensity<3,float> >&
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

no reason to drag shared_ptr along. just use const DiscretisedDensity<3,float>& (shorter and less dependent on choices made elsewhere). Alternatively, make this a private function, in which case you dont need to pass them along.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sounds good

src/buildblock/ProjDataInfo.cxx Outdated Show resolved Hide resolved
src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx Outdated Show resolved Hide resolved
src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx Outdated Show resolved Hide resolved
src/include/stir/DiscretisedDensity.inl Outdated Show resolved Hide resolved
src/include/stir/ProjDataInfo.h Outdated Show resolved Hide resolved
@KrisThielemans KrisThielemans added this to the v4.0 milestone Oct 14, 2018
@ashgillman ashgillman self-assigned this Oct 15, 2018
@ashgillman
Copy link
Collaborator Author

Not sure about DiscretisedDensity know about gantry coordinates (e.g. get_index_coordinates_for_gantry_coordinates). It means a rather tight coupling between images and projection data. I feel that ProjDataInfo should take care of gantry<->bed stuff.

Similarly, ProjDataInfo::get_centre_of_gantry_vector_in_relative_coordinates doesn't make sense, as ProjDataInfo doesn't know what "relative coordinates" are.

Therefore I would not have these functions. Anyone who needs to do this will just have call 2 functions (ProjDataInfo::get_centre_of_gantry_in_bed_coordinates and the relevant DiscretisedDensity::*from_physical_coordinates). What do you think?

I agree, but had some difficulty drawing a clear line in the separation here. It is quite obviously ugly as it is now, requiring a reference to a ProjDataInfo in the function, so I am happy for suggestions.

The problem is that in the projection space, the FoR is the centre of the gantry. Whereas in Discretised Density, it is from the first ring. So querying ProjDataInfo for the centre in its own coordinates would be (0, 0, 0). So somewhere, there should be a function that essentially adds the vector from the centre of first ring to the center of the gantry. That's what get_relative_coordinates_for_gantry_coordinates() does. I agree that it actually makes more sense to have it in ProjDataInfo now, since it uses nothing from DiscDens. (I had originally thought this function would need to know where the image boundaries were, i.e., what get_index_coordinates_for_gantry_coordinates() does, but this can be replaced by density->get_index_coordinates_for_relative_coordinates(projdatainfo->get_relative_coordinates_for_gantry_coordinates(coords)))

But this would still have ProjDataInfo knowing a little bit about relative coordinates. To avoid ambiguity, I think there should be a function somewhere that converts sinogram-based coordinates (from get_m, etc, which are from center) to relative coordinates.

@ashgillman
Copy link
Collaborator Author

Its a bit confusing because I am using the terminology gantry_coordinates to mean gantry-based but also with FoR at centre. You could define bed_coordinates to have FoR either at centre or at end of first ring.

@ashgillman
Copy link
Collaborator Author

Also a minor one: the comments here are obsolete?

Yes, and I have modified the explanation following slightly

@ashgillman
Copy link
Collaborator Author

Looking at the failed test, it seems these changes break the symmetries test. But I'm not too sure why this would be, nothing that I changed in PRMT seems to be symmetry-specific.

The errors look like:

1: Error. comparing lors
1: Current bin:  segment = 4, axial pos 8, view = 7, tangential_pos_num = 6
1: {19, -3, -8}:0.21517    ||
1: {19, -2, -8}:0.0784836    ||
1: {19, -2, -7}:1.00928    ||
1: {19, -1, -7}:1.08776    ||
1: {19, 0, -7}:0.529053    ||
1: {19, 0, -6}:0.558712    ||
1: {19, 1, -6}:1.08777    ||
1: {19, 2, -6}:0.979621    ||
1: {19, 2, -5}:0.108144    ||
1: {19, 3, -5}:1.08777    ||
1: {19, 4, -5}:1.08777    ||
1: {19, 5, -5}:0.342424    ||
1: {19, 5, -4}:0.745341    ||
1: {19, 6, -4}:1.08777    ||
1: {19, 7, -4}:0.792992    ||
1:                        ||   {21, -3, -8}:0.215167
1:                        ||   {21, -2, -8}:0.0784869
1:                        ||   {21, -2, -7}:1.00928
1:                        ||   {21, -1, -7}:1.08777
1:                        ||   {21, 0, -7}:0.529054
1:                        ||   {21, 0, -6}:0.558711
1:                        ||   {21, 1, -6}:1.08777
1:                        ||   {21, 2, -6}:0.979621
1:                        ||   {21, 2, -5}:0.108144
1:                        ||   {21, 3, -5}:1.08777
1:                        ||   {21, 4, -5}:1.08777
1:                        ||   {21, 5, -5}:0.342423
1:                        ||   {21, 5, -4}:0.745342
1:                        ||   {21, 6, -4}:1.08776
1:                        ||   {21, 7, -4}:0.792991

So looks like the z voxel is different if symmetries are/aren't used.

Copy link
Collaborator

@KrisThielemans KrisThielemans left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

some hopefully helpful comments.

I'm a bit reluctant to merge in all the ORIGINTODO mods. we should have them, but probably not on master.... Easy to remove that commit?

@@ -41,6 +41,7 @@
#include "stir/CartesianCoordinate3D.h"
#include "stir/Array.h"
#include "stir/IO/ExamData.h"
#include "stir/ProjDataInfo.h"
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

shouldn't be needed if it wasn't before (will need to be in cxx, but probably already there)

@@ -344,6 +347,7 @@ set_up(
// test if our 2D code does not have problems
{
// currently 2D code relies on the LOR falling in the middle of a voxel (in z-direction)
// If reactivated, update to use proper index -> physical space conversion
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

not sure what is meant by if reactivated

src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx Outdated Show resolved Hide resolved
(min(max_index.y(), -min_index.y())+.45F)*voxel_size.y());
CartesianCoordinate3D<float> max_pos =
density_info_sptr
->get_relative_coordinates_for_indices(max_index + 0.45F);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think best to use physical. At some point, we will allow x/y-offset shifts. The "relative" coordinates won't see this, and therefore we would use the wrong FOV

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

agree, I'd like to phase out all relative coordinates, either removing the methods or changing them to private.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

let's just make them private (or protected) in a first (2nd?) stage

// within the TOR, leveraging our traced LOR

// since tantheta==0, z is constant
const float z_of_traced_lor
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

best to call this variable (and the next ones) z_index...

src/include/stir/ProjDataInfo.inl Outdated Show resolved Hide resolved
@ashgillman
Copy link
Collaborator Author

@KrisThielemans
I've run into an issue with where to put the gantry <-> physical conversions. In our plan, we wanted this to be in ProjDataInfo. However, this requires (at least in order to be backwards compatible) knowledge of the chosen DiscretisedDensity, because the projectors currently always center the image axially (this makes it ugly). However, this won't compile as it causes circular dependency issues, as DD already depends on PDI.

I've got something like this:

CartesianCoordinate3D<float>
ProjDataInfo::
get_physical_coordinates_for_gantry_coordinates
(const CartesianCoordinate3D<float>& coords, const shared_ptr<DiscretisedDensity<3,float> >& image_sptr) const
{
  // TODO: bed postion
  return coords + get_gantry_origin_in_physical_coordinates(image_sptr);
}

and

CartesianCoordinate3D<float>
ProjDataInfo::
get_gantry_origin_in_physical_coordinates(const shared_ptr<DiscretisedDensity<3,float> >& image_sptr) const
{
  // TODO - STIR currently enforces that an image is centred in the z-plane.
  // In future, this would make more sense to be implemented as:
  // float middle_of_first_ring_to_middle_of_last
  //   = (scanner_ptr->get_num_rings() - 1) * scanner_ptr->get_ring_spacing();
  // return CartesianCoordinate3D<float>(middle_of_first_ring_to_middle_of_last / 2.F, 0, 0);

  // shift by vector to centre of image in z (assume this corresponds to centre of gantry)
  CartesianCoordinate3D<int> min_index;
  CartesianCoordinate3D<int> max_index;
  image_sptr->get_regular_range(min_index, max_index);
  CartesianCoordinate3D<float> middle_of_image = image_sptr->get_physical_coordinates_for_indices(
    (min_index + max_index) / 2.);
  return CartesianCoordinate3D<float>(middle_of_image.z(), 0, 0);
}

@KrisThielemans
Copy link
Collaborator

it's ugly indeed, but I also can't see a way around it for now.

However, this won't compile as it causes circular dependency issues, as DD already depends on PDI.

I'm not sure why that creates problems. as long as we pre-declare DD in PDI.h and PDI in DD.h, then the implementations can sit in the .cxx. Everything should compile and link ok. Maybe you can send the compilation error?

By the way, I'd like to avoid using shared_ptr args as much as possible. Use references whenever possible/appropriate, and that seems definitely the case for get_gantry_origin_in_physical_coordinates etc

@ashgillman
Copy link
Collaborator Author

Okay, so shared pointer should only be used if the object stores a reference as a member?

ProjDataInfo.cxx.o builds okay, but later it gets to DiscretisedDensity.cxx.o, and I get an error then:

_ make
[  1%] Building CXX object src/buildblock/CMakeFiles/buildblock.dir/DiscretisedDensity.cxx.o
In file included from /datasets/scratch/HB_ATLAS_APPS/hb-18-cdc/sirf/master/sources/STIR/src/include/stir/DiscretisedDensity.h:45:0,
                 from /datasets/scratch/HB_ATLAS_APPS/hb-18-cdc/sirf/master/sources/STIR/src/buildblock/DiscretisedDensity.cxx:34:
/datasets/scratch/HB_ATLAS_APPS/hb-18-cdc/sirf/master/sources/STIR/src/include/stir/ProjDataInfo.h:386:54: error: 'DiscretisedDensity' does not name a type
   (const CartesianCoordinate3D<float>& coords, const DiscretisedDensity<3,float>& image) const;
                                                      ^~~~~~~~~~~~~~~~~~
/datasets/scratch/HB_ATLAS_APPS/hb-18-cdc/sirf/master/sources/STIR/src/include/stir/ProjDataInfo.h:386:72: error: expected ',' or '...' before '<' token
   (const CartesianCoordinate3D<float>& coords, const DiscretisedDensity<3,float>& image) const;
                                                                        ^
/datasets/scratch/HB_ATLAS_APPS/hb-18-cdc/sirf/master/sources/STIR/src/include/stir/ProjDataInfo.h:404:51: error: 'DiscretisedDensity' does not name a type
   get_gantry_origin_in_physical_coordinates(const DiscretisedDensity<3,float>& image) const;
                                                   ^~~~~~~~~~~~~~~~~~
/datasets/scratch/HB_ATLAS_APPS/hb-18-cdc/sirf/master/sources/STIR/src/include/stir/ProjDataInfo.h:404:69: error: expected ',' or '...' before '<' token
   get_gantry_origin_in_physical_coordinates(const DiscretisedDensity<3,float>& image) const;
                                                                     ^
In file included from /datasets/scratch/HB_ATLAS_APPS/hb-18-cdc/sirf/master/sources/STIR/src/include/stir/ProjDataInfo.h:409:0,
                 from /datasets/scratch/HB_ATLAS_APPS/hb-18-cdc/sirf/master/sources/STIR/src/include/stir/DiscretisedDensity.h:45,
                 from /datasets/scratch/HB_ATLAS_APPS/hb-18-cdc/sirf/master/sources/STIR/src/buildblock/DiscretisedDensity.cxx:34:
/datasets/scratch/HB_ATLAS_APPS/hb-18-cdc/sirf/master/sources/STIR/src/include/stir/ProjDataInfo.inl:123:52: error: 'DiscretisedDensity' does not name a type (const CartesianCoordinate3D<float>& coords, const DiscretisedDensity<3,float>& image) const
                                                    ^~~~~~~~~~~~~~~~~~
/datasets/scratch/HB_ATLAS_APPS/hb-18-cdc/sirf/master/sources/STIR/src/include/stir/ProjDataInfo.inl:123:70: error: expected ',' or '...' before '<' token
 (const CartesianCoordinate3D<float>& coords, const DiscretisedDensity<3,float>& image) const
                                                                      ^
/datasets/scratch/HB_ATLAS_APPS/hb-18-cdc/sirf/master/sources/STIR/src/include/stir/ProjDataInfo.inl: In member function 'stir::CartesianCoordinate3D<float> stir::ProjDataInfo::get_physical_coordinates_for_gantry_coordinates(const stir::CartesianCoordinate3D<float>&, int) const':
/datasets/scratch/HB_ATLAS_APPS/hb-18-cdc/sirf/master/sources/STIR/src/include/stir/ProjDataInfo.inl:126:61: error: 'image' was not declared in this scope
   return coords + get_gantry_origin_in_physical_coordinates(image);
                                                             ^~~~~
/datasets/scratch/HB_ATLAS_APPS/hb-18-cdc/sirf/master/sources/STIR/src/include/stir/ProjDataInfo.inl:126:61: note: suggested alternative: 'cimagq'
   return coords + get_gantry_origin_in_physical_coordinates(image);
                                                             ^~~~~
                                                             cimagq
src/buildblock/CMakeFiles/buildblock.dir/build.make:662: recipe for target 'src/buildblock/CMakeFiles/buildblock.dir/DiscretisedDensity.cxx.o' failed
make[2]: *** [src/buildblock/CMakeFiles/buildblock.dir/DiscretisedDensity.cxx.o] Error 1
CMakeFiles/Makefile2:707: recipe for target 'src/buildblock/CMakeFiles/buildblock.dir/all' failed
make[1]: *** [src/buildblock/CMakeFiles/buildblock.dir/all] Error 2
Makefile:140: recipe for target 'all' failed
make: *** [all] Error 2

So it seems to be some weird circular thing maybe?

@ashgillman
Copy link
Collaborator Author

I #include DD in PDI.cxx and PDI.h. I guess I should pre-declare instead? I'll have to Google how to do that.

@KrisThielemans
Copy link
Collaborator

I #include DD in PDI.cxx and PDI.h. I guess I should pre-declare instead? I'll have to Google how to do that.

in .h, you should use

template <blabla> class DD;

In the .cxx, you include DD.h

@KrisThielemans
Copy link
Collaborator

so shared pointer should only be used if the object stores a reference as a member?

that's confusingly phrased :-) but yes, I feel more and more it should only be used if we really are going to share ownership (there's too many places in STIR that do...). It makes it semantically clear, it allows calling it with actual objects (e.g. of type VoxelsOnCartesianGrid) and unique_ptr or bare pointers.

@KrisThielemans
Copy link
Collaborator

should add this to the developer's guide really...

@KrisThielemans KrisThielemans modified the milestones: v4.1, v6.0 May 23, 2020
@KrisThielemans KrisThielemans modified the milestones: v6.0, v7.0 Jan 20, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
Image Geometry
  
In progress
Development

Successfully merging this pull request may close these issues.

None yet

2 participants