+
+// CLASS IDENTICATION
+// BackprojectIntDiff3
+//
+// PURPOSE
+// Highly optimized version of BackprojectIntDiff2
+
+void
+BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const double view_angle)
+{
+ double theta = view_angle; // add half PI to view angle to get perpendicular theta angle
+ static const int scaleShift = 16;
+ static const kint32 scale = (1 << scaleShift);
+ static const kint32 scaleBitmask = scale - 1;
+ static const kint32 halfScale = scale / 2;
+ static const double dInvScale = 1. / scale;
+
+ const kint32 det_dx = nearest<kint32> (xInc * cos (theta) / detInc * scale);
+ const kint32 det_dy = nearest<kint32> (yInc * sin (theta) / detInc * scale);
+
+ // calculate L for first point in image (0, 0)
+ kint32 detPosColStart = nearest<kint32> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
+
+ // precalculate scaled difference for linear interpolation
+ 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;
+ }
+ deltaFilteredProj[nDet - 1] = 0; // last detector
+
+ int iLastDet = nDet - 1;
+ for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
+ kint32 curDetPos = detPosColStart;
+ ImageFileColumn pImCol = v[ix];
+
+ if (interpType == Backprojector::INTERP_NEAREST) {
+ for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
+ const int iDetPos = (curDetPos + halfScale) >> 16;
+ if (iDetPos >= 0 && iDetPos <= iLastDet)
+ *pImCol++ += filteredProj[iDetPos];
+ } // end for iy
+ } else if (interpType == Backprojector::INTERP_FREQ_PREINTERPOLATION) {
+ for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
+ const int iDetPos = ((curDetPos + halfScale) >> 16) * m_interpFactor;
+ if (iDetPos >= 0 && iDetPos <= iLastDet)
+ *pImCol++ += filteredProj[iDetPos];
+ } // end for iy
+ } else if (interpType == Backprojector::INTERP_LINEAR) {
+ for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
+ const kint32 iDetPos = curDetPos >> scaleShift;
+ const kint32 detRemainder = curDetPos & scaleBitmask;
+ if (iDetPos >= 0 && iDetPos <= iLastDet)
+ *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
+ } // end for iy
+ } //end linear
+ } // end for ix
+
+ delete deltaFilteredProj;
+}
+
+
+void
+BackprojectEquiangular::BackprojectView (const double* const filteredProj, const double view_angle)
+{
+ double beta = view_angle;
+
+ for (int ix = 0; ix < nx; ix++) {
+ ImageFileColumn pImCol = v[ix];
+
+ for (int iy = 0; iy < ny; iy++) {
+ double dAngleDiff = beta - phi[ix][iy];
+ double rcos_t = r[ix][iy] * cos (dAngleDiff);
+ double rsin_t = r[ix][iy] * sin (dAngleDiff);
+ double dFLPlusSin = m_dFocalLength + rsin_t;
+ double gamma = atan (rcos_t / dFLPlusSin);
+ double dL2 = dFLPlusSin * dFLPlusSin + (rcos_t * rcos_t);
+
+ if (interpType == Backprojector::INTERP_NEAREST) {
+ int iDetPos =iDetCenter + nearest<int>(gamma / detInc); // calc index in the filtered raysum vector
+
+ if (iDetPos < 0 || iDetPos >= nDet) { // check for impossible: index outside of raysum pos
+ ; // errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], gamma, iDetPos);
+ } else
+ pImCol[iy] += filteredProj[iDetPos] / dL2;
+ } else if (interpType == Backprojector::INTERP_LINEAR) {
+ double dPos = gamma / detInc; // position along detector
+ double dPosFloor = floor (dPos);
+ int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
+ double frac = dPos - dPosFloor; // fraction distance from det
+ if (iDetPos < 0 || iDetPos >= nDet - 1) {
+ ; // errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], gamma, iDetPos);
+ } else
+ pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / dL2;
+ }
+ } // end for y
+ } // end for x
+}
+
+void
+BackprojectEquilinear::BackprojectView (const double* const filteredProj, const double view_angle)
+{
+ double beta = view_angle;
+
+ for (int ix = 0; ix < nx; ix++) {
+ ImageFileColumn pImCol = v[ix];
+
+ for (int iy = 0; iy < ny; iy++) {
+ double dAngleDiff = beta - phi[ix][iy];
+ double rcos_t = r[ix][iy] * cos (dAngleDiff);
+ double rsin_t = r[ix][iy] * sin (dAngleDiff);
+
+ double dU = (m_dFocalLength + rsin_t) / m_dFocalLength;
+ double dDetPos = rcos_t / dU;
+ // double to scale for imaginary detector that passes through origin
+ // of phantom, see Kak-Slaney Figure 3.22
+ dDetPos *= 2;
+
+ if (interpType == Backprojector::INTERP_NEAREST) {
+ int iDetPos = iDetCenter + nearest<int>(dDetPos / detInc); // calc index in the filtered raysum vector
+
+ if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
+ ; /// errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], dDetPos, iDetPos);
+ else
+ pImCol[iy] += (filteredProj[iDetPos] / (dU * dU));
+ } else if (interpType == Backprojector::INTERP_LINEAR) {
+ double dPos = dDetPos / detInc; // position along detector
+ double dPosFloor = floor (dPos);
+ int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
+ double frac = dPos - dPosFloor; // fraction distance from det
+ if (iDetPos < 0 || iDetPos >= nDet - 1)
+ ; // errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], dDetPos, iDetPos);
+ else
+ pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / (dU * dU);
+ }
+ } // end for y
+ } // end for x
+}
+