X-Git-Url: http://git.kpe.io/?a=blobdiff_plain;f=libctsim%2Fbackprojectors.cpp;h=0d3cd804892ec01ed3b8e5299122fbfd9781b6de;hb=7a50ad8860d015c3e06038475929bc62011b561c;hp=0bcd4077cd527281d11769b4084c3d29e07cce1d;hpb=a05f3cb550877e94aa118cc04b361c0c8fdb3dc3;p=ctsim.git

diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp
index 0bcd407..0d3cd80 100644
--- a/libctsim/backprojectors.cpp
+++ b/libctsim/backprojectors.cpp
@@ -8,7 +8,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: backprojectors.cpp,v 1.13 2000/08/31 08:38:58 kevin Exp $
+**  $Id: backprojectors.cpp,v 1.18 2000/12/06 15:17:51 kevin Exp $
 **
 **  This program is free software; you can redistribute it and/or modify
 **  it under the terms of the GNU General Public License (version 2) as
@@ -256,13 +256,14 @@ Backprojector::convertInterpIDToTitle (const int interpID)
 // PURPOSE
 //   Pure virtual base class for all backprojectors.
 
-Backproject::Backproject (const Projections& proj, ImageFile& im, const int interpType, const int interpFactor)
+Backproject::Backproject (const Projections& proj, ImageFile& im, int interpType, const int interpFactor)
     : proj(proj), im(im), interpType(interpType), m_interpFactor(interpFactor)
 {
   detInc = proj.detInc();
   nDet = proj.nDet();
-  iDetCenter = nDet / 2;	// index refering to L=0 projection 
-  rotInc = proj.rotInc();
+  iDetCenter = (nDet - 1) / 2;	// index refering to L=0 projection 
+  rotScale = proj.rotInc();
+  rotScale /= (proj.nView() * proj.rotInc() / PI); // scale by number of PI rotations
 
   v = im.getArray();
   nx = im.nx();
@@ -288,7 +289,7 @@ Backproject::ScaleImageByRotIncrement ()
 {
   for (int ix = 0; ix < nx; ix++)
     for (int iy = 0; iy < ny; iy++)
-      v[ix][iy] *= rotInc;
+      v[ix][iy] *= rotScale;
 }
 
 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double r, double phi, double L, int iDetPos)
@@ -298,7 +299,8 @@ void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, doubl
 }
 
 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double L, int iDetPos)
-{
+{
+#if 1
   ostringstream os;
   os << "ix=" << ix << ", iy=" << iy << ", theta=" << theta << ", L=" << L << ", detinc=" << detInc << "\n";
   os << "ndet=" << nDet << ", detInc=" << detInc << ", iDetCenter=" << iDetCenter << "\n";
@@ -306,7 +308,8 @@ void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, doubl
   os << "yMin=" << yMin << ", yMax=" << yMax << ", yInc=" << yInc << "\n";
   os << "iDetPos index outside bounds: " << iDetPos << " [backprojector]";;
 
-  sys_error (ERR_WARNING, os.str().c_str());
+  sys_error (ERR_WARNING, os.str().c_str());
+#endif
 }
 
 
@@ -335,7 +338,7 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double
 	if (iDetPos < 0 || iDetPos >= nDet)	// check for impossible: index outside of raysum pos 
 	    errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos);
 	else
-	  v[ix][iy] += rotInc * filteredProj[iDetPos];
+	  v[ix][iy] += rotScale * filteredProj[iDetPos];
       } else if (interpType == Backprojector::INTERP_LINEAR) {
 	  double p = L / detInc;	// position along detector
 	  double pFloor = floor (p);
@@ -344,7 +347,7 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double
 	  if (iDetPos < 0 || iDetPos >= nDet - 1) 	// check for impossible: index outside of raysum pos 
 	    errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos);
 	  else
-	    v[ix][iy] += rotInc * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
+	    v[ix][iy] += rotScale * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
       }
     }
   }
@@ -358,10 +361,10 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double
 //   Precalculates trigometric function value for each point in image for backprojection.
 
 BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, int interpType, const int interpFactor)
-  : Backproject::Backproject (proj, im, interpType, interpFactor)
+  : Backproject (proj, im, interpType, interpFactor)
 {
-  arrayR.initSetSize (nx, ny);
-  arrayPhi.initSetSize (nx, ny);
+  arrayR.initSetSize (im.nx(), im.ny());
+  arrayPhi.initSetSize (im.nx(), im.ny());
   r = arrayR.getArray();
   phi = arrayPhi.getArray();
 
@@ -420,7 +423,7 @@ BackprojectTable::BackprojectView (const double* const filteredProj, const doubl
 //   Iterates in x & y direction by adding difference in L position
 
 BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, int interpType, const int interpFactor)
-  :  Backproject::Backproject (proj, im, interpType, interpFactor)
+  :  Backproject (proj, im, interpType, interpFactor)
 {
   // calculate center of first pixel v[0][0] 
   double x = xMin + xInc / 2;
@@ -493,7 +496,7 @@ BackprojectDiff2::BackprojectView (const double* const filteredProj, const doubl
   double detPosColStart = start_r * cos (theta - start_phi) / detInc;
 	
 #ifdef DEBUG
-  printf ("start_r=%8.5f, start_phi=%8.5f, rotInc=%8.5f\n", start_r, start_phi, rotInc);
+  printf ("start_r=%8.5f, start_phi=%8.5f, rotScale=%8.5f\n", start_r, start_phi, rotScale);
 #endif
   for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
     double curDetPos = detPosColStart;
@@ -598,7 +601,7 @@ BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const do
   kint32 detPosColStart = nearest<kint32> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
 	
   // precalculate scaled difference for linear interpolation
-  double deltaFilteredProj [nDet];
+  double* deltaFilteredProj = new double [nDet];
   if (interpType == Backprojector::INTERP_LINEAR) {
     for (int i = 0; i < nDet - 1; i++)
       deltaFilteredProj[i] = (filteredProj[i+1] - filteredProj[i]) * dInvScale;
@@ -630,7 +633,9 @@ BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const do
 	  *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
       }	// end for iy
     } //end linear
-  } // end for ix
+  } // end for ix
+
+  delete deltaFilteredProj;
 }
 
 
@@ -665,7 +670,7 @@ BackprojectEquiangular::BackprojectView (const double* const filteredProj, const
 	if (iDetPos < 0 || iDetPos >= nDet - 1) {
 	  ; //	    errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], gamma, iDetPos);
 	} else
-	  pImCol[iy] += (((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1])) / dL2;
+	  pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / dL2;
       }
     }	// end for y 
   }	// end for x 
@@ -705,7 +710,7 @@ BackprojectEquilinear::BackprojectView (const double* const filteredProj, const
 	if (iDetPos < 0 || iDetPos >= nDet - 1)
 	  ; //	    errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], dDetPos, iDetPos);
 	else
-	  pImCol[iy] += (((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1])) / (dU * dU);
+	  pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / (dU * dU);
       }
     }	// end for y 
   }	// end for x