Helical4cylindrical #418

lesaintjerome · 2021-05-19T12:21:13Z

First attempt for an implementation of Katsevich algorithm as described in the paper "Exact helical reconstruction using native cone-beam geometries", Noo et al., PMB, 2003.
Both flat panel and curved detector geometries are implemented though the curved geometry is still under debugging.
As for the flat panel implementation, the whole thing more or less works. Some strange "tiling" artifacts remain.

I thought it was worth publishing the PR right now, even if a lot of work is still to be done, in order to get some feedback on the technical aspects of the coding.

The algorithm runs in 4 steps :

Compute a derivative of the projections
Cosine-weight
Forward rebinning + Hilbert transform (from ACoussat PR) + backward rebinning
Backprojection.

A new geometry class is created, to account for the specific helical trajectory.
Also, a few applications are added, which basically run each of the steps independently. There might not be necessary in the final version.

Future work :

Debug curved detector
Implement a complete reconstruction image filter (instead of an application that links all the steps).

SimonRit

Thanks for the huge work. I have provided a few comments. Please check CONTRIBUTING.md for formatting your commit messages. It can be prefixed with WIP until it's ready to be merged (in your opinion).

SimonRit · 2021-05-25T07:29:04Z

.gitattributes

@@ -1,6 +1,6 @@
 .git*                                   export-ignore
 # Custom attribute to mark sources as using our C++/C code style.
-[attr]our-c-style  whitespace=tab-in-indent,no-lf-at-eof  hooks.style=KWStyle,clangformat,uncrustify
+[attr]our-c-style  whitespace=tab-in-indent,no-lf-at-eof  hooks.style=KWStyle,clangformat


This change can be kept our of this PR but you can propose it in a separate PR.

I rolled back to the original version (with the uncrustify option) and compilation is OK.
Note though that all this was commited with -n option. Some style guidelines may not be complied with.

SimonRit · 2021-05-25T07:29:34Z

.gitignore

@@ -33,3 +33,6 @@ CMakeLists.txt.user*

 # Mac System File
 .DS_Store
+
+
+.project


What is a .project file?

Eclipse file to create project from git dir. Removed.

SimonRit · 2021-05-25T07:30:15Z

applications/CMakeLists.txt

+#All executables below are specific to the tomo reconstruction from projections aquired along a helical trajectory of the gantry
+add_subdirectory(rtkhelicalgeometry)
+add_subdirectory(rtkkatsevichderivative)
+add_subdirectory(rtkkatsevichforwardbinning)
+add_subdirectory(rtkhilbertonkappalines)
+add_subdirectory(rtkpilines)
+add_subdirectory(rtkkatsevich)
+add_subdirectory(rtkkatsevichtrash)
+


I agree that those can be removed and be replace by, e.g., a single rtkkatsevitch command line tool. I'm skipping the code of these applications in my review. Maybe rtkhelicalgeometry can be merged in rtksimulatedgeometry with a new --pitch option (0 by default)?

The final katsevichrecons application is now submitted for review.
I agree that rtkhelicalgeometry could be merged in rtksimulatedgeometry. But we probably need to decide what to do with the ThreeDHelicalProjectionGeometry class : merge into ThreeDCircular or not ? Probably best.

SimonRit · 2021-05-25T07:35:08Z

include/rtkFFTHilbertImageFilter.h

@@ -0,0 +1,112 @@
+/*=========================================================================


If those files are directly copied from #371, I should probably merge it first and you rebase your PR on it.

SimonRit · 2021-05-25T07:36:06Z

include/rtkHilbertTransformOnKappaLinesImageFilter.hxx

+  // Mettre a jour la taille de l'image de sortie. Ainsi que l'origine.
+  // Sur le modele de itk::ShrinkImageFilter


Avoid French comments...

Fixed. The filter will probably disappear in final version.

SimonRit · 2021-05-25T07:41:16Z

include/rtkKatsevichBackwardBinningImageFilter.hxx

+  //  // Input 0 is Output ie result of forward binning
+  //  typename Superclass::InputImagePointer inputPtr0 = const_cast<TInputImage *>(this->GetInput());
+  //  std::cout << "Input largest : " << this->GetInput()->GetLargestPossibleRegion() << std::endl;
+  //  std::cout << "inputPtr0 : " << inputPtr0->GetRequestedRegion() << std::endl;
+  //  if (!inputPtr0)
+  //    return;
+  //  inputPtr0->SetRequestedRegion(this->GetInput()->GetLargestPossibleRegion());
+  //
+  //  //TO DELETE : SPACING AND ORIGIN ALWAYS (1,1,1) AND (0,0,0) -> NOT OKAY !!!!
+  //  //std::cout<<"Input size: "<<this->GetInput()->GetLargestPossibleRegion().GetSize()<<std::endl;
+  //  //std::cout<<"Input spacing : "<<this->GetInput()->GetSpacing()<<std::endl;
+  //  //std::cout<<"Input origin : "<<this->GetInput()->GetOrigin()<<std::endl;
+  //
+  //


To be deleted or implemented

SimonRit · 2021-05-25T07:42:46Z

include/rtkKatsevichDerivativeImageFilter.h

+/** \class KatsevichDerivativeImageFilter
+ * \brief Computes the derivatives of the projections wrt the projection index (the rotation angle).
+ *
+ * This filter implements the derivative as described in Noo et al., PMB, 2003


Is there only one implementation in it? I would specify the relevant equations.

Agreed. Reference to eq in Noo et al. added.

SimonRit · 2021-05-25T07:43:28Z

include/rtkKatsevichForwardBinningImageFilter.h

+ * Backprojects a stack of projection images (input 1) in a 3D volume (input 0)
+ * using linear interpolation according to a specified geometry. The operation
+ * is voxel-based, meaning that the center of each voxel is projected in the
+ * projection images to determine the interpolation location.
+ *
+ * \test rtkfovtest.cxx
+ *
+ * \author Simon Rit
+ *
+ * \ingroup RTK Projector


To be updated

SimonRit · 2021-05-25T07:45:00Z

include/rtkThreeDCircularProjectionGeometry.h

@@ -106,7 +106,7 @@ class RTK_EXPORT ThreeDCircularProjectionGeometry : public ProjectionGeometry<3>
   * @param detectorRowVector absolute direction vector indicating the
   * orientation of the detector's rows r (sometimes referred to as v1)
   * @param detectorColumnVector absolute direction vector indicating the
-   * orientation of the detector's columns c (sometimes referred to as v2)
+   * orientation of the detector's columns c (sometimes refered to as v2)


This one is weird ! This file should not have been touched. Should I correct it or can I rollback to the original version ?

SimonRit · 2021-05-25T07:47:48Z

include/rtkThreeDHelicalProjectionGeometry.h

+ * The source and the detector rotate around a helix paremeterized
+ * with the SourceToDetectorDistance and the SourceToIsocenterDistance and a helical pitch.
+ * The position of each projection along the helix is parameterized
+ * by the HelicalAngle.
+ * The detector is assumed to have no horizontal offset (ProjectionOffsetX = 0). It can be also rotated with
+ * InPlaneAngles and OutOfPlaneAngles. All angles are in radians except for the function AddProjection that takes angles
+ * in degrees. The source is assumed to have no horizontal offset (SourceOffsetX = 0).


We can probably improve a bit the definition and the implementation but this can be left for later at this stage.

lesaintjerome · 2021-05-27T10:33:48Z

Thanks for the feedback. Will dig into it shortly.
Meanwhile, I added a reconstruction filter and the corresponding application. The latter will probably the only remaining application in the final PR.

First attempt for an implementation of Katsevich algorithm as described in the paper "Exact helical reconstruction using native cone-beam geometries", Noo et al., PMB, 2003. Both flat panel and curved detector geometries are implemented though the curved geometry is still under debugging. As for the flat panel implementation, the whole thing more or less works. Some strange "tiling" artifacts remain. The algorithm runs in 4 steps : Compute a derivative of the projections Cosine-weight Forward rebinning + Hilbert transform (from ACoussat PR) + backward rebinning Backprojection. First commit - 2 First commit - nth Add post-cosine weight for curved detector filtering Various debug stuff and computation of FOV radius Misc Update Hilbert transform with ACoussant PR + minor changes requiring future cleanup Clean code from debugging couts Add a Katsevich Reconstruction filter and corresponding application

lesaintjerome · 2021-05-28T15:58:34Z

Squashed all commits into 1, with WIP prefixing.

Point by point response to PR review.

To avoid redundancy in the code, the rtk::KatsevichBackProjectionImageFilter inherits from rtk::BackProjectionImageFilter. WARNING : Compilation threw a couple of warning due to Geometry getter/setter (implemented in macros) whcih overrides the mother function, without saying it. But don't know how to deal with that.

acoussat · 2021-06-03T07:49:48Z

Impressive PR! Out of curiosity, did you find any issue regarding the Hilbert filter? Did you have to fix anything?

lesaintjerome · 2021-06-03T08:35:44Z

Impressive PR! Out of curiosity, did you find any issue regarding the Hilbert filter? Did you have to fix anything?

The initial version of your PR - without PixelShift - produced severe oscillations on the transformed signal. But the newer one fixed it. Further testing will confirm.

SimonRit reviewed May 25, 2021

View reviewed changes

lesaintjerome force-pushed the helical4cylindrical branch from 5ef1f55 to 19970ed Compare May 28, 2021 15:55

lesaintjerome added 2 commits May 29, 2021 13:16

WIP: Response to PR review

5fe19c2

Point by point response to PR review.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Helical4cylindrical #418

Helical4cylindrical #418

lesaintjerome commented May 19, 2021

SimonRit left a comment

SimonRit May 25, 2021

lesaintjerome May 28, 2021

SimonRit May 25, 2021

lesaintjerome May 28, 2021

SimonRit May 25, 2021

lesaintjerome May 29, 2021

SimonRit May 25, 2021

SimonRit May 25, 2021

lesaintjerome May 28, 2021

SimonRit May 25, 2021

lesaintjerome May 28, 2021

SimonRit May 25, 2021

lesaintjerome May 29, 2021

SimonRit May 25, 2021

lesaintjerome May 29, 2021

SimonRit May 25, 2021

lesaintjerome May 29, 2021

SimonRit May 25, 2021

lesaintjerome commented May 27, 2021

lesaintjerome commented May 28, 2021

acoussat commented Jun 3, 2021

lesaintjerome commented Jun 3, 2021

		@@ -0,0 +1,112 @@
		/*=========================================================================

		// Mettre a jour la taille de l'image de sortie. Ainsi que l'origine.
		// Sur le modele de itk::ShrinkImageFilter

Helical4cylindrical #418

Are you sure you want to change the base?

Helical4cylindrical #418

Conversation

lesaintjerome commented May 19, 2021

SimonRit left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

lesaintjerome commented May 27, 2021

lesaintjerome commented May 28, 2021

acoussat commented Jun 3, 2021

lesaintjerome commented Jun 3, 2021