Pct feriel

AddTrackerUncertainty.py

@@ -64,9 +64,8 @@ def GetSigmaSc(energy, xOverX0):
     pairs[:,2,0:2] += dYuncertEntry[:,1,:] # Entrance direction X/Y
     pairs[:,1,0:2] += dYuncertExit[:,0,:] # Exit position X/Y
     pairs[:,3,0:2] += dYuncertExit[:,1,:] # Exit direction X/Y
-    
+
     itk.imwrite(itk.GetImageFromArray(pairs, is_vector=True), output)
 
 if __name__ == '__main__':
     AddTrackerUncertainty()
-

BPF.py

@@ -0,0 +1,80 @@
+#!/usr/bin/env python
+
+import click
+import SimpleITK as sitk
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy.signal as signal
+import scipy.special as special
+
+
+@click.command()
+@click.option("-i", "--input", required=True, help="Input file name")
+@click.option("-o", "--output", required=True, help="Output file name")
+@click.option("-m", "--m", required=True, default=500, help="Size of the reconstruction image matrix mxm")
+@click.option("-k", "--k",required=True, default=500, help="Size of the region used for offset correction kxk")
+
+def BPF(input, output, m, k):
+
+    def BPFfilter(N, dx):
+        x = y = np.arange(N) * dx - (N / 2 - 1) *dx #filter sampled such that x=y=0 is a sample otherwise artifacts appear
+        X , Y = np.meshgrid (x, y)
+        k = np.zeros(X.shape)
+        R = np.sqrt(X**2 + Y**2)
+        not_zero = R!=0
+        var_k = np.pi * R[not_zero] / dx
+        k[not_zero] = 1./(4 * np.pi**2 * R[not_zero]**3) * (var_k**2 * special.jn(1,var_k) - np.pi * var_k/2 * (special.jn(1,var_k) * special.struve(0,var_k) - special.jn(0,var_k) * special.struve(1,var_k)))
+        k[R==0] = np.pi / (12 * dx**3)
+        return k
+
+    corrSize = int(k)
+    reconSize = int(m)
+    imgReader = sitk.ImageFileReader()
+    imgReader.SetFileName(input)
+    sino = imgReader.Execute()
+    spacing = sino.GetSpacing()
+    offset = sino.GetOrigin()
+    sino = sitk.GetArrayFromImage(sino)
+    bpSize = sino.shape[1]
+    nb_proj = sino.shape[0]
+    dphi = np.pi / nb_proj
+    proj_axis = 0
+    dx = spacing[-1]
+
+    #Compute backprojection
+    bp = np.sum(sino, axis = proj_axis) * dphi
+    sino = []
+
+    #Compute filter on very large matrix
+    k = BPFfilter(corrSize, dx)
+    a = int((corrSize - bpSize) / 2) #to resize filter to BP size
+
+    #Calculations for offset correction
+    xbp = np.arange(bpSize) * dx - (bpSize - 1) / 2 * dx
+    Xbp, Ybp = np.meshgrid (xbp, xbp)
+    xc = np.arange(corrSize) * dx - (corrSize - 1) / 2 * dx
+    Xc, Yc = np.meshgrid (xc, xc)
+    idx = np.in1d(Xc, Xbp) * np.in1d(Yc, Ybp)
+    idx = np.reshape(idx, Xc.shape)
+    xcind = np.arange(corrSize)
+    Xcind, Ycind = np.meshgrid(xcind, xcind)
+    SecondIntegral = (dx**2 * np.sum(k[corrSize - 1 - Xcind[~idx], corrSize - 1 - Ycind[~idx]] / (np.sqrt(Xc[~idx]**2 + Yc[~idx]**2))))
+    b = int((bpSize - reconSize) / 2) #resize bp to recon size
+
+    reconstruction = np.zeros(bp.shape)
+    for z in np.arange(bp.shape[1]):
+        reconstruction[:,z,:] = signal.fftconvolve(bp[:,z,:], k[a:-a,a:-a], mode='same') * dx**2
+        delta = (dx**2 * np.sum(reconstruction[b:-b,z,b:-b])) * SecondIntegral
+        reconstruction[:,z,:] += delta
+
+    offset_recon = -(reconSize - 1) / 2 * dx
+    reconstruction = reconstruction[b:-b,:,b:-b]
+    reconimg = sitk.GetImageFromArray(reconstruction)
+    reconimg.SetSpacing(spacing[:-1])
+    reconimg.SetOrigin([offset_recon, offset[1],offset_recon])
+    writer = sitk.ImageFileWriter()
+    writer.SetFileName(output)
+    writer.Execute(reconimg)
+
+if __name__ == '__main__':
+    BPF()

pctProtonPairsToBackProjection.h

@@ -29,7 +29,7 @@ class ITK_EXPORT ProtonPairsToBackProjection :
   typedef itk::Image<ProtonPairsPixelType,2>           ProtonPairsImageType;
   typedef ProtonPairsImageType::Pointer                ProtonPairsImagePointer;
 
-  typedef typename itk::Image< unsigned int,
+  typedef typename itk::Image< double,
                                TInputImage::ImageDimension> CountImageType;
   typedef typename CountImageType::Pointer                  CountImagePointer;
 
@@ -90,6 +90,16 @@ class ITK_EXPORT ProtonPairsToBackProjection :
   itkSetMacro(DisableRotation, bool);
   itkBooleanMacro(DisableRotation);
 
+  /** Get / Set the boolean to use Collins-Fekete binning. Default is off. */
+  itkGetMacro(WeightsCF, bool);
+  itkSetMacro(WeightsCF, bool);
+  itkBooleanMacro(WeightsCF);
+
+  /** Get / Set the boolean to output MLP uncertainty map. Default is off. */
+  itkGetMacro(SigmaMap, bool);
+  itkSetMacro(SigmaMap, bool);
+  itkBooleanMacro(SigmaMap);
+
 protected:
   ProtonPairsToBackProjection();
   virtual ~ProtonPairsToBackProjection() {}
@@ -133,6 +143,12 @@ class ITK_EXPORT ProtonPairsToBackProjection :
   bool m_DisableRotation = false;
 
   std::mutex m_Mutex;
+
+  /** Use binning as in Collins-Fekete */
+  bool m_WeightsCF;
+
+  /** Output sigma map */
+  bool m_SigmaMap;
 };
 
 } // end namespace pct

pctProtonPairsToBackProjection.txx

@@ -62,15 +62,16 @@ ProtonPairsToBackProjection<TInputImage, TOutputImage>
       [this, iProj](const ProtonPairsImageType::RegionType & outputRegionForThread)
         {
         // Create MLP depending on type
-        pct::MostLikelyPathFunction<double>::Pointer mlp;
-        if(m_MostLikelyPathType == "polynomial")
-          mlp = pct::ThirdOrderPolynomialMLPFunction<double>::New();
-        else if (m_MostLikelyPathType == "schulte")
-          mlp = pct::SchulteMLPFunction::New();
-        else
-          {
-          itkGenericExceptionMacro("MLP must either be schulte or polynomial, not [" << m_MostLikelyPathType << ']');
-          }
+        //pct::MostLikelyPathFunction<double>::Pointer mlp;
+	pct::SchulteMLPFunction::Pointer mlp = pct::SchulteMLPFunction::New();
+        //if(m_MostLikelyPathType == "polynomial")
+        //  mlp = pct::ThirdOrderPolynomialMLPFunction<double>::New();
+        //else if (m_MostLikelyPathType == "schulte")
+        //  mlp = pct::SchulteMLPFunction::New();
+        //else
+        //  {
+        //  itkGenericExceptionMacro("MLP must either be schulte or polynomial, not [" << m_MostLikelyPathType << ']');
+        //  }
 
         // Get geometry information. We need the rotation matrix alone to rotate the direction
         // and the same matrix combined with mm (physical point) to voxel conversion for the volume.
@@ -96,7 +97,7 @@ ProtonPairsToBackProjection<TInputImage, TOutputImage>
         const typename OutputImageType::SpacingType imgSpacing = this->GetInput()->GetSpacing();
 
         typename OutputImageType::PixelType *imgData = this->GetOutput()->GetBufferPointer();
-        unsigned int *imgCountData = m_Counts->GetBufferPointer();
+        double *imgCountData = m_Counts->GetBufferPointer();
         itk::Vector<float, 3> imgSpacingInv;
         double minSpacing = imgSpacing[0];
         for(unsigned int i=0; i<3; i++)
@@ -179,14 +180,38 @@ ProtonPairsToBackProjection<TInputImage, TOutputImage>
             value = m_ConvFunc->GetValue(eOut, eIn); // convert to WEPL
           ++it;
 
+          VectorType nucInfo(0.);
+
+          if(it.GetIndex()[0] != 0)
+            {
+            nucInfo = it.Get();
+            ++it;
+            }
+
           // Move straight to entrance and exit shapes
           VectorType pSIn = pIn;
           VectorType pSOut = pOut;
           double nearDistIn, nearDistOut, farDistIn, farDistOut;
+
+          VectorType pRotIn (0.);
+          VectorType dRotIn (0.);
+          VectorType pRotOut(0.);
+          VectorType dRotOut (0.);
+          for(unsigned int i=0; i<3; i++)
+            {
+            for(unsigned int j=0; j<3; j++)
+              {
+              pRotIn[i] += rotMat[i][j] * pIn[j];
+              dRotIn[i] += rotMat[i][j] * dIn[j];
+              pRotOut[i] += rotMat[i][j] * pOut[j];
+              dRotOut[i] += rotMat[i][j] * dOut[j];
+              }
+            }
+
           if(m_QuadricIn.GetPointer()!=NULL)
             {
-            if(m_QuadricIn->IsIntersectedByRay(pIn,dIn,nearDistIn,farDistIn) &&
-               m_QuadricOut->IsIntersectedByRay(pOut,dOut,nearDistOut,farDistOut))
+            if(m_QuadricIn->IsIntersectedByRay(pRotIn,dRotIn,nearDistIn,farDistIn) &&
+               m_QuadricOut->IsIntersectedByRay(pRotOut,dRotOut,nearDistOut,farDistOut))
               {
               pSIn  = pIn  + dIn  * nearDistIn;
               if(pSIn[2]<pIn[2]  || pSIn[2]>pOut[2])
@@ -195,6 +220,13 @@ ProtonPairsToBackProjection<TInputImage, TOutputImage>
               if(pSOut[2]<pIn[2] || pSOut[2]>pOut[2])
                 pSOut = pOut + dOut * farDistOut;
               }
+            else
+              {
+            //  pSIn = pOut;
+            //  pSOut = pIn;
+              pSIn[2] = 0;
+              pSOut[2] = 0;
+              }
             }
 
           // Normalize direction with respect to z
@@ -207,9 +239,13 @@ ProtonPairsToBackProjection<TInputImage, TOutputImage>
 
           VectorType dCurr = dIn;
           VectorType pCurr(0.);
+          itk::Matrix<double, 2, 2> mlp_error;
+          mlp_error.Fill(0);
 
           // Init MLP before mm to voxel conversion
           mlp->Init(pSIn, pSOut, dIn, dOut);
+          int totalLength = 0.;
+          std::vector<typename OutputImageType::OffsetValueType> offsets;
 
           for(unsigned int k=0; k<zmm.size(); k++)
             {
@@ -288,10 +324,33 @@ ProtonPairsToBackProjection<TInputImage, TOutputImage>
                idx[1]>=0 && idx[1]<(int)imgSize[1] &&
                idx[2]>=0 && idx[2]<(int)imgSize[2])
               {
-              typename OutputImageType::OffsetValueType offset = this->GetOutput()->ComputeOffset(idx);
+              if(m_WeightsCF)
+                {
+                idx[2] = 0;
+                offsets.push_back(this->GetOutput()->ComputeOffset(idx));
+                totalLength++;
+                }
+              else
+                {
+                typename OutputImageType::OffsetValueType offset = this->GetOutput()->ComputeOffset(idx);
+                m_Mutex.lock();
+                imgData[ offset ] += value;
+                imgCountData[ offset ]++;
+                m_Mutex.unlock();
+                }
+              }
+            }
+          if(m_WeightsCF)
+            {
+            std::vector<typename OutputImageType::OffsetValueType> channels(offsets);
+            std::sort( channels.begin(), channels.end() );
+            channels.erase( std::unique( channels.begin(), channels.end() ), channels.end() );
+            for(unsigned int k=0; k<channels.size(); k++)
+              {
+              double weight =  (double) std::count(offsets.begin(),offsets.end(),channels[k])/totalLength;
               m_Mutex.lock();
-              imgData[ offset ] += value;
-              imgCountData[ offset ]++;
+              imgData[ channels[k] ] += value * weight * weight;
+              imgCountData[ channels[k] ] += weight * weight;
               m_Mutex.unlock();
               }
             }

pctProtonPairsToDistanceDrivenProjection.h

@@ -26,7 +26,7 @@ class ITK_EXPORT ProtonPairsToDistanceDrivenProjection :
   typedef itk::Image<ProtonPairsPixelType,2>                ProtonPairsImageType;
   typedef ProtonPairsImageType::Pointer                     ProtonPairsImagePointer;
 
-  typedef itk::Image<unsigned int, 3>                       CountImageType;
+  typedef itk::Image<float, 3>                       CountImageType;
   typedef CountImageType::Pointer                           CountImagePointer;
 
   typedef itk::Image<float, 3>                              AngleImageType;
@@ -37,7 +37,8 @@ class ITK_EXPORT ProtonPairsToDistanceDrivenProjection :
   typedef typename OutputImageType::RegionType              OutputImageRegionType;
 
   typedef rtk::QuadricShape                                 RQIType;
-
+  typedef rtk::ThreeDCircularProjectionGeometry     GeometryType;
+  typedef typename GeometryType::Pointer            GeometryPointer;
   /** Method for creation through the object factory. */
   itkNewMacro(Self);
 
@@ -68,6 +69,9 @@ class ITK_EXPORT ProtonPairsToDistanceDrivenProjection :
  /** Get/Set the angle of proton pairs per pixel. */
   itkGetMacro(Angle, AngleImagePointer);
 
+ /** Get/Set the MLP uncertainty of proton pairs per pixel. */
+  itkGetMacro(Sigma, OutputImagePointer);
+
   /** Get/Set the ionization potential used in the Bethe-Bloch equation. */
   itkGetMacro(IonizationPotential, double);
   itkSetMacro(IonizationPotential, double);
@@ -84,6 +88,27 @@ class ITK_EXPORT ProtonPairsToDistanceDrivenProjection :
   itkGetConstMacro(ComputeScattering, bool);
   itkBooleanMacro(ComputeScattering);
 
+  /** Get / Set the boolean to use Collins-Fekete binning. Default is off. */
+  itkGetMacro(WeightsCF, bool);
+  itkSetMacro(WeightsCF, bool);
+  itkBooleanMacro(WeightsCF);
+
+  /** Get / Set the boolean to output MLP uncertainty map. Default is off. */
+  itkGetMacro(SigmaMap, bool);
+  itkSetMacro(SigmaMap, bool);
+  itkBooleanMacro(SigmaMap);
+
+  itkGetMacro(Geometry, GeometryPointer);
+  itkSetMacro(Geometry, GeometryPointer);
+
+  itkGetMacro(Index, int);
+  itkSetMacro(Index, int);
+
+  /** Get / Set the boolean to use MLP Krah. Default is off. */
+  itkGetMacro(MLPKrah, bool);
+  itkSetMacro(MLPKrah, bool);
+  itkBooleanMacro(MLPKrah);
+
 protected:
   ProtonPairsToDistanceDrivenProjection();
   virtual ~ProtonPairsToDistanceDrivenProjection() {}
@@ -121,6 +146,8 @@ class ITK_EXPORT ProtonPairsToDistanceDrivenProjection :
   std::vector<OutputImagePointer> m_Outputs;
   std::vector<OutputImagePointer> m_AngleOutputs;
   std::vector<OutputImagePointer> m_AngleSqOutputs;
+  std::vector<OutputImagePointer> m_Sigmas;
+  OutputImagePointer m_Sigma;
 
   /** The two quadric functions defining the object support. */
   RQIType::Pointer m_QuadricIn;
@@ -135,6 +162,21 @@ class ITK_EXPORT ProtonPairsToDistanceDrivenProjection :
   ProtonPairsImageType::Pointer m_ProtonPairs;
   bool                          m_Robust;
   bool                          m_ComputeScattering;
+
+  /** Use binning as in Collins-Fekete */
+  bool m_WeightsCF;
+
+  /** Output sigma map */
+  bool m_SigmaMap;
+
+  /** RTK geometry object */
+  GeometryPointer m_Geometry;
+
+  /** Index of the projection in the geometry object */
+  int m_Index;
+
+  /** Use MLP of Krah */
+  bool m_MLPKrah;
 };
 
 } // end namespace pct