Merge pull request #555 from simonCtBern/float_precision_fix_issue_544

float precision fix #544

Common/CUDA/Siddon_projection.cu

@@ -674,100 +674,90 @@ void computeDeltas_Siddon(Geometry geo,int i, Point3D* uvorigin, Point3D* deltaU
     //End point
     Point3D P,Pu0,Pv0;
 
-    P.x  =-(geo.DSD[i]-geo.DSO[i]);   P.y  = geo.dDetecU*(0-((float)geo.nDetecU/2)+0.5);       P.z  = geo.dDetecV*(((float)geo.nDetecV/2)-0.5-0);
-    Pu0.x=-(geo.DSD[i]-geo.DSO[i]);   Pu0.y= geo.dDetecU*(1-((float)geo.nDetecU/2)+0.5);       Pu0.z= geo.dDetecV*(((float)geo.nDetecV/2)-0.5-0);
-    Pv0.x=-(geo.DSD[i]-geo.DSO[i]);   Pv0.y= geo.dDetecU*(0-((float)geo.nDetecU/2)+0.5);       Pv0.z= geo.dDetecV*(((float)geo.nDetecV/2)-0.5-1);
-    // Geomtric trasnformations:
+    P.x  =-(geo.DSD[i]-geo.DSO[i]);   P.y  = geo.dDetecU*(-((double)geo.nDetecU/2.0)+0.5);       P.z  = geo.dDetecV*(((double)geo.nDetecV/2.0)-0.5);
+    Pu0.x=0;                          Pu0.y= geo.dDetecU;                                    Pu0.z= 0;
+    Pv0.x=0;                          Pv0.y= 0;                                              Pv0.z= geo.dDetecV*(-1);
+
+    // Geometric transformations:
     // Now we have the Real world (OXYZ) coordinates of the bottom corner and its two neighbours.
-    // The obkjective is to get a position of the detector in a coordinate system where:
+    // The objective is to get a position of the detector in a coordinate system where:
     // 1-units are voxel size (in each direction can be different)
     // 2-The image has the its first voxel at (0,0,0)
     // 3-The image never rotates
-    
+
     // To do that, we need to compute the "deltas" the detector, or "by how much
     // (in new xyz) does the voxels change when and index is added". To do that
     // several geometric steps needs to be changed
-    
+
     //1.Roll,pitch,jaw
     // The detector can have a small rotation.
     // according to
     //"A geometric calibration method for cone beam CT systems" Yang K1, Kwan AL, Miller DF, Boone JM. Med Phys. 2006 Jun;33(6):1695-706.
     // Only the Z rotation will have a big influence in the image quality when they are small.
     // Still all rotations are supported
-    
+
     // To roll pitch jaw, the detector has to be in centered in OXYZ.
-    P.x=0;Pu0.x=0;Pv0.x=0;
-
+    // NB: do not apply offsets to Pu0 and Pv0: they are directions, and are invariant through translations
+    P.x=0;
+
     // Roll pitch yaw
     rollPitchYaw(geo,i,&P);
     rollPitchYaw(geo,i,&Pu0);
     rollPitchYaw(geo,i,&Pv0);
-    //Now ltes translate the points where they should be:
+    //Now let's translate the points where they should be:
+    // NB: do not apply offsets to Pu0 and Pv0: they are directions, and are invariant through translations
     P.x=P.x-(geo.DSD[i]-geo.DSO[i]);
-    Pu0.x=Pu0.x-(geo.DSD[i]-geo.DSO[i]);
-    Pv0.x=Pv0.x-(geo.DSD[i]-geo.DSO[i]);
-
+
     //1: Offset detector
-    
-    
+
+
     //S doesnt need to chagne
-    
-    
+
+
     //3: Rotate (around z)!
     Point3D Pfinal, Pfinalu0, Pfinalv0;
     Pfinal.x  =P.x;
     Pfinal.y  =P.y  +geo.offDetecU[i]; Pfinal.z  =P.z  +geo.offDetecV[i];
-    Pfinalu0.x=Pu0.x;
-    Pfinalu0.y=Pu0.y  +geo.offDetecU[i]; Pfinalu0.z  =Pu0.z  +geo.offDetecV[i];
-    Pfinalv0.x=Pv0.x;
-    Pfinalv0.y=Pv0.y  +geo.offDetecU[i]; Pfinalv0.z  =Pv0.z  +geo.offDetecV[i];
-
+    Pfinalu0 = Pu0;
+    Pfinalv0 = Pv0;
+
     eulerZYZ(geo,&Pfinal);
     eulerZYZ(geo,&Pfinalu0);
     eulerZYZ(geo,&Pfinalv0);
     eulerZYZ(geo,&S);
-    
+
     //2: Offset image (instead of offseting image, -offset everything else)
-
+    // NB: do not apply offsets to Pfinalu0 and Pfinalv0: they are directions, and are invariant through translations
+
     Pfinal.x  =Pfinal.x-geo.offOrigX[i];     Pfinal.y  =Pfinal.y-geo.offOrigY[i];     Pfinal.z  =Pfinal.z-geo.offOrigZ[i];
-    Pfinalu0.x=Pfinalu0.x-geo.offOrigX[i];   Pfinalu0.y=Pfinalu0.y-geo.offOrigY[i];   Pfinalu0.z=Pfinalu0.z-geo.offOrigZ[i];
-    Pfinalv0.x=Pfinalv0.x-geo.offOrigX[i];   Pfinalv0.y=Pfinalv0.y-geo.offOrigY[i];   Pfinalv0.z=Pfinalv0.z-geo.offOrigZ[i];
     S.x=S.x-geo.offOrigX[i];               S.y=S.y-geo.offOrigY[i];               S.z=S.z-geo.offOrigZ[i];
-    
+
     // As we want the (0,0,0) to be in a corner of the image, we need to translate everything (after rotation);
     Pfinal.x  =Pfinal.x+geo.sVoxelX/2;      Pfinal.y  =Pfinal.y+geo.sVoxelY/2;          Pfinal.z  =Pfinal.z  +geo.sVoxelZ/2;
-    Pfinalu0.x=Pfinalu0.x+geo.sVoxelX/2;    Pfinalu0.y=Pfinalu0.y+geo.sVoxelY/2;        Pfinalu0.z=Pfinalu0.z+geo.sVoxelZ/2;
-    Pfinalv0.x=Pfinalv0.x+geo.sVoxelX/2;    Pfinalv0.y=Pfinalv0.y+geo.sVoxelY/2;        Pfinalv0.z=Pfinalv0.z+geo.sVoxelZ/2;
     S.x      =S.x+geo.sVoxelX/2;          S.y      =S.y+geo.sVoxelY/2;              S.z      =S.z      +geo.sVoxelZ/2;
-    
+
     //4. Scale everything so dVoxel==1
     Pfinal.x  =Pfinal.x/geo.dVoxelX;      Pfinal.y  =Pfinal.y/geo.dVoxelY;        Pfinal.z  =Pfinal.z/geo.dVoxelZ;
     Pfinalu0.x=Pfinalu0.x/geo.dVoxelX;    Pfinalu0.y=Pfinalu0.y/geo.dVoxelY;      Pfinalu0.z=Pfinalu0.z/geo.dVoxelZ;
     Pfinalv0.x=Pfinalv0.x/geo.dVoxelX;    Pfinalv0.y=Pfinalv0.y/geo.dVoxelY;      Pfinalv0.z=Pfinalv0.z/geo.dVoxelZ;
     S.x      =S.x/geo.dVoxelX;          S.y      =S.y/geo.dVoxelY;            S.z      =S.z/geo.dVoxelZ;
-    
-    
+
+
     //mexPrintf("COR: %f \n",geo.COR[i]);
     //5. apply COR. Wherever everything was, now its offesetd by a bit
-    float CORx, CORy;
+    // NB: do not apply offsets to Pfinalu0 and Pfinalv0: they are directions, and are invariant through translations
+    double CORx, CORy;
     CORx=-geo.COR[i]*sin(geo.alpha)/geo.dVoxelX;
     CORy= geo.COR[i]*cos(geo.alpha)/geo.dVoxelY;
     Pfinal.x+=CORx;   Pfinal.y+=CORy;
-    Pfinalu0.x+=CORx;   Pfinalu0.y+=CORy;
-    Pfinalv0.x+=CORx;   Pfinalv0.y+=CORy;
     S.x+=CORx; S.y+=CORy;
-    
+
     // return
-    
+
     *uvorigin=Pfinal;
-
-    deltaU->x=Pfinalu0.x-Pfinal.x;
-    deltaU->y=Pfinalu0.y-Pfinal.y;
-    deltaU->z=Pfinalu0.z-Pfinal.z;
-
-    deltaV->x=Pfinalv0.x-Pfinal.x;
-    deltaV->y=Pfinalv0.y-Pfinal.y;
-    deltaV->z=Pfinalv0.z-Pfinal.z;
+
+    *deltaU=Pfinalu0;
+    *deltaV=Pfinalv0;
 
     *source=S;
 }

Common/CUDA/types_TIGRE.hpp

@@ -89,4 +89,21 @@ struct  Geometry {
     float y;
     float z;
 };
+
+struct Point3Ddouble{
+    double x;
+    double y;
+    double z;
+
+    // cast to float member function for "copying" Point3Ddouble to Point3D
+    Point3D to_float()
+    {
+        Point3D castToFloat;
+        castToFloat.x = (float)x;
+        castToFloat.y = (float)y;
+        castToFloat.z = (float)z;
+        return(castToFloat);
+    }
+};
+
 #endif

Common/CUDA/voxel_backprojection.cu

@@ -752,6 +752,69 @@ void freeGeoArray(unsigned int splits,Geometry* geoArray){
     }
     free(geoArray);
 }
+
+
+//______________________________________________________________________________
+//
+//      double precision functions for rotating Point3Ddouble coordinates
+//______________________________________________________________________________
+
+void eulerZYZT(Geometry geo, Point3Ddouble* point){
+
+    Point3Ddouble auxPoint;
+    auxPoint.x=point->x;
+    auxPoint.y=point->y;
+    auxPoint.z=point->z;
+
+    // calculate sin and cos of 3 angles (used multiple times)
+    double sin_alpha, cos_alpha, sin_theta, cos_theta, sin_psi, cos_psi;
+    sin_alpha = sin((double)geo.alpha);
+    cos_alpha = cos((double)geo.alpha);
+    sin_theta = sin((double)geo.theta);
+    cos_theta = cos((double)geo.theta);
+    sin_psi = sin((double)geo.psi);
+    cos_psi = cos((double)geo.psi);
+
+    point->x = auxPoint.x*(cos_psi*cos_theta*cos_alpha-sin_psi*sin_alpha)
+    +auxPoint.y*(-cos_psi*cos_theta*sin_alpha-sin_psi*cos_alpha)
+    +auxPoint.z*cos_psi*sin_theta;
+    point->y = auxPoint.x*(sin_psi*cos_theta*cos_alpha+cos_psi*sin_alpha)
+    +auxPoint.y*(-sin_psi*cos_theta*sin_alpha+cos_psi*cos_alpha)
+    +auxPoint.z*sin_psi*sin_theta;
+    point->z =-auxPoint.x*sin_theta*cos_alpha
+    +auxPoint.y*sin_theta*sin_alpha
+    +auxPoint.z*cos_theta;
+}
+
+void rollPitchYawT(Geometry geo,int i, Point3Ddouble* point){
+
+    Point3Ddouble auxPoint;
+    auxPoint.x=point->x;
+    auxPoint.y=point->y;
+    auxPoint.z=point->z;
+
+    // calculate sin and cos of 3 angles (used multiple times)
+    double sin_dRoll, cos_dRoll, sin_dPitch, cos_dPitch, sin_dYaw, cos_dYaw;
+    sin_dRoll = sin((double)geo.dRoll[i]);
+    cos_dRoll = cos((double)geo.dRoll[i]);
+    sin_dPitch = sin((double)geo.dPitch[i]);
+    cos_dPitch = cos((double)geo.dPitch[i]);
+    sin_dYaw = sin((double)geo.dYaw[i]);
+    cos_dYaw = cos((double)geo.dYaw[i]);
+
+    point->x=cos_dRoll*cos_dPitch*auxPoint.x
+            +sin_dRoll*cos_dPitch*auxPoint.y
+            -sin_dPitch*auxPoint.z;
+
+    point->y=(cos_dRoll*sin_dPitch*sin_dYaw - sin_dRoll*cos_dYaw)*auxPoint.x
+            +(sin_dRoll*sin_dPitch*sin_dYaw + cos_dRoll*cos_dYaw)*auxPoint.y
+            +cos_dPitch*sin_dYaw*auxPoint.z;
+
+    point->z=(cos_dRoll*sin_dPitch*cos_dYaw + sin_dRoll*sin_dYaw)*auxPoint.x
+            +(sin_dRoll*sin_dPitch*cos_dYaw - cos_dRoll*sin_dYaw)*auxPoint.y
+            +cos_dPitch*cos_dYaw*auxPoint.z;
+}
+
 //______________________________________________________________________________
 //
 //      Function:       computeDeltasCube
@@ -763,28 +826,25 @@ void freeGeoArray(unsigned int splits,Geometry* geoArray){
 void computeDeltasCube(Geometry geo,int i, Point3D* xyzorigin, Point3D* deltaX, Point3D* deltaY, Point3D* deltaZ,Point3D* S)
 {
 
-    Point3D P, Px,Py,Pz;
+    // initialize points with double precision
+    Point3Ddouble P, Px,Py,Pz;
+
     // Get coords of Img(0,0,0)
     P.x=-(geo.sVoxelX/2-geo.dVoxelX/2)+geo.offOrigX[i];
     P.y=-(geo.sVoxelY/2-geo.dVoxelY/2)+geo.offOrigY[i];
     P.z=-(geo.sVoxelZ/2-geo.dVoxelZ/2)+geo.offOrigZ[i];
 
-    // Get coors from next voxel in each direction
-    Px.x=P.x+geo.dVoxelX;      Py.x=P.x;                Pz.x=P.x;
+    // Get coords from next voxel in each direction
+    Px.x=P.x+geo.dVoxelX;       Py.x=P.x;                Pz.x=P.x;
     Px.y=P.y;                   Py.y=P.y+geo.dVoxelY;    Pz.y=P.y;
     Px.z=P.z;                   Py.z=P.z;                Pz.z=P.z+geo.dVoxelZ;
 
-
-
-// Rotate image around X axis (this is equivalent of rotating the source and detector) RZ RY RZ
-
+    // Rotate image around X axis (this is equivalent of rotating the source and detector) RZ RY RZ
     eulerZYZT(geo,&P);
     eulerZYZT(geo,&Px);
     eulerZYZT(geo,&Py);
     eulerZYZT(geo,&Pz);
 
-
-
     //detector offset
     P.z =P.z-geo.offDetecV[i];            P.y =P.y-geo.offDetecU[i];
     Px.z =Px.z-geo.offDetecV[i];          Px.y =Px.y-geo.offDetecU[i];
@@ -793,7 +853,6 @@ void computeDeltasCube(Geometry geo,int i, Point3D* xyzorigin, Point3D* deltaX,
 
     //Detector Roll pitch Yaw
     //
-    //
     // first, we need to offset everything so (0,0,0) is the center of the detector
     // Only X is required for that
     P.x=P.x+(geo.DSD[i]-geo.DSO[i]);
@@ -810,13 +869,12 @@ void computeDeltasCube(Geometry geo,int i, Point3D* xyzorigin, Point3D* deltaX,
     Py.x=Py.x-(geo.DSD[i]-geo.DSO[i]);
     Pz.x=Pz.x-(geo.DSD[i]-geo.DSO[i]);
     //Done for P, now source
-    Point3D source;
+    Point3Ddouble source;
     source.x=geo.DSD[i]; //already offseted for rotation
     source.y=-geo.offDetecU[i];
     source.z=-geo.offDetecV[i];
     rollPitchYawT(geo,i,&source);
 
-
     source.x=source.x-(geo.DSD[i]-geo.DSO[i]);//   source.y=source.y-auxOff.y;    source.z=source.z-auxOff.z;
 
 //       mexPrintf("%f,%f,%f\n",source.x,source.y,source.z);
@@ -834,49 +892,11 @@ void computeDeltasCube(Geometry geo,int i, Point3D* xyzorigin, Point3D* deltaX,
     deltaY->x=Py.x-P.x;   deltaY->y=Py.y-P.y;    deltaY->z=Py.z-P.z;
     deltaZ->x=Pz.x-P.x;   deltaZ->y=Pz.y-P.y;    deltaZ->z=Pz.z-P.z;
 
-
-    *xyzorigin=P;
-    *S=source;
+    // cast the results from the double precision calculations back to float
+    *xyzorigin=P.to_float();
+    *S=source.to_float();
 }
 
-void eulerZYZT(Geometry geo, Point3D* point){
-
-    Point3D auxPoint;
-    auxPoint.x=point->x;
-    auxPoint.y=point->y;
-    auxPoint.z=point->z;
-
-    point->x = auxPoint.x*(cos(geo.psi)*cos(geo.theta)*cos(geo.alpha)-sin(geo.psi)*sin(geo.alpha))
-    +auxPoint.y*(-cos(geo.psi)*cos(geo.theta)*sin(geo.alpha)-sin(geo.psi)*cos(geo.alpha))
-    +auxPoint.z*cos(geo.psi)*sin(geo.theta);
-    point->y = auxPoint.x*(sin(geo.psi)*cos(geo.theta)*cos(geo.alpha)+cos(geo.psi)*sin(geo.alpha))
-    +auxPoint.y*(-sin(geo.psi)*cos(geo.theta)*sin(geo.alpha)+cos(geo.psi)*cos(geo.alpha))
-    +auxPoint.z*sin(geo.psi)*sin(geo.theta);
-    point->z =-auxPoint.x*sin(geo.theta)*cos(geo.alpha)
-    +auxPoint.y*sin(geo.theta)*sin(geo.alpha)
-    +auxPoint.z*cos(geo.theta);
-}
-void rollPitchYawT(Geometry geo,int i, Point3D* point){
-    Point3D auxPoint;
-    auxPoint.x=point->x;
-    auxPoint.y=point->y;
-    auxPoint.z=point->z;
-
-    point->x=cos(geo.dRoll[i])*cos(geo.dPitch[i])*auxPoint.x
-            +sin(geo.dRoll[i])*cos(geo.dPitch[i])*auxPoint.y
-            -sin(geo.dPitch[i])*auxPoint.z;
-
-
-    point->y=(cos(geo.dRoll[i])*sin(geo.dPitch[i])*sin(geo.dYaw[i]) - sin(geo.dRoll[i])*cos(geo.dYaw[i]))*auxPoint.x
-            +(sin(geo.dRoll[i])*sin(geo.dPitch[i])*sin(geo.dYaw[i]) + cos(geo.dRoll[i])*cos(geo.dYaw[i]))*auxPoint.y
-            +cos(geo.dPitch[i])*sin(geo.dYaw[i])*auxPoint.z;
-
-
-    point->z=(cos(geo.dRoll[i])*sin(geo.dPitch[i])*cos(geo.dYaw[i]) + sin(geo.dRoll[i])*sin(geo.dYaw[i]))*auxPoint.x
-            +(sin(geo.dRoll[i])*sin(geo.dPitch[i])*cos(geo.dYaw[i]) - cos(geo.dRoll[i])*sin(geo.dYaw[i]))*auxPoint.y
-            +cos(geo.dPitch[i])*cos(geo.dYaw[i])*auxPoint.z;
-
-}
 void checkFreeMemory(const GpuIds& gpuids,size_t *mem_GPU_global){
     size_t memfree;
     size_t memtotal;

Common/CUDA/voxel_backprojection.hpp

@@ -48,10 +48,10 @@ Codes  : https://github.com/CERN/TIGRE
 
 #ifndef BACKPROJECTION_HPP
 #define BACKPROJECTION_HPP
-void rollPitchYawT(Geometry geo,int i, Point3D* point);
+void rollPitchYawT(Geometry geo,int i, Point3Ddouble* point);
 int  voxel_backprojection(float* projections, Geometry geo, float* result,float const * const alphas,int nalpha, const GpuIds& gpuids);
 void splitCTbackprojection(const GpuIds& gpuids,Geometry geo,int nalpha, unsigned int* split_image, unsigned int * split_projections);
-void eulerZYZT(Geometry geo, Point3D* point);
+void eulerZYZT(Geometry geo, Point3Ddouble* point);
 void computeDeltasCube(Geometry geo,int i, Point3D* xyzorigin, Point3D* deltaX, Point3D* deltaY, Point3D* deltaZ,Point3D* S);
 void createGeoArray(unsigned int image_splits, Geometry geo,Geometry* geoArray, unsigned int nangles);
 void freeGeoArray(unsigned int splits,Geometry* geoArray);