X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fprojections.cpp;h=871f62e6e52a4543df9262355872d3b22957f392;hp=4d4a02d2ed560990419e565490fe8cfd0b3e142c;hb=8a7697ce57b56cdc43698cd1241ad98d49f9b5ac;hpb=3ea498d51ce4597e9649cd21f155b51175ea0bea diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp index 4d4a02d..871f62e 100644 --- a/libctsim/projections.cpp +++ b/libctsim/projections.cpp @@ -8,7 +8,7 @@ ** This is part of the CTSim program ** Copyright (c) 1983-2001 Kevin Rosenberg ** -** $Id: projections.cpp,v 1.68 2001/03/21 21:45:31 kevin Exp $ +** $Id$ ** ** 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 @@ -25,6 +25,8 @@ ******************************************************************************/ #include "ct.h" +#include +#include "interpolator.h" const kuint16 Projections::m_signature = ('P'*256 + 'J'); @@ -35,21 +37,22 @@ const int Projections::POLAR_INTERP_BICUBIC = 2; const char* const Projections::s_aszInterpName[] = { - {"nearest"}, - {"bilinear"}, + "nearest", + "bilinear", // {"bicubic"}, }; const char* const Projections::s_aszInterpTitle[] = { - {"Nearest"}, - {"Bilinear"}, + "Nearest", + "Bilinear", // {"Bicubic"}, }; const int Projections::s_iInterpCount = sizeof(s_aszInterpName) / sizeof(char*); + /* NAME * Projections Constructor for projections matrix storage * @@ -155,7 +158,7 @@ Projections::initFromScanner (const Scanner& scanner) m_dFocalLength = scanner.focalLength(); m_dSourceDetectorLength = scanner.sourceDetectorLength(); m_dViewDiameter = scanner.viewDiameter(); - m_rotStart = 0; + m_rotStart = scanner.offsetView()*scanner.rotInc(); m_dFanBeamAngle = scanner.fanBeamAngle(); } @@ -166,6 +169,277 @@ Projections::setNView (int nView) // used by MPI to reduce # of views init (nView, m_nDet); } +// Helical 180 Linear Interpolation. +// This member function takes a set of helical scan projections and +// performs a linear interpolation between pairs of complementary rays +// to produce a single projection data set approximating what would be +// measured at a single axial plane. +// Complementary rays are rays which traverse the same path through the +// phantom in opposite directions. +// +// For parallel beam geometry, a ray with a given gantry angle beta and a +// detector iDet will have a complementary ray at beta + pi and nDet-iDet +// +// For equiangular or equilinear beam geometry the complementary ray to +// gantry angle beta and fan-beam angle gamma is at +// beta-hat = beta +2*gamma + pi, and gamma-hat = -gamma. +// Note that beta-hat - beta depends on gamma and is not constant. +// +// The algorithm used here is from Crawford and King, Med. Phys. 17(6) +// 1990 p967; what they called method "C", CSH-HH. It uses interpolation only +// between pairs of complementary rays on either side of an image plane. +// Input data must sample gantry angles from zero to +// (2*pi + 2* fan-beam-angle). The data set produced contains gantry +// angles from 0 to Pi+fan-beam-angle. This is a "halfscan" data set, +// which still contains redundant data, and can be used with a half scan +// reconstruction to produce an image. +// In this particular implementation a lower triangle from (beta,gamma) = +// (0,-fanAngle/2)->(2*fanAngle,-fanAngle/2)->(0,fanAngle/2) contains +// zeros, but is actually redundant with data contained in the region +// (pi+fanAngle,-fanAngle/2)->(pi+fanAngle, fanAngle/2) ->(pi-fanAngle, +// fanAngle/2). +// +int +Projections::Helical180LI(int interpolation_view) +{ + if (m_geometry == Scanner::GEOMETRY_INVALID) + { + std::cerr << "Invalid geometry " << m_geometry << std::endl; + return (2); + } + else if (m_geometry == Scanner::GEOMETRY_PARALLEL) + { + std::cerr << "Helical 180LI not yet implemented for PARALLEL geometry" + << std::endl; + return (2); + } + else if (m_geometry == Scanner::GEOMETRY_EQUILINEAR) + { + std::cerr << "Helical 180LI not yet implemented for EQUILINEAR geometry" + << std::endl; + return (2); + } + else if (m_geometry == Scanner::GEOMETRY_EQUIANGULAR) + { + return Helical180LI_Equiangular(interpolation_view); + } + else + { + std::cerr << "Invalid geometry in projection data file" << m_geometry + << std::endl; + return (2); + } +} +int +Projections::Helical180LI_Equiangular(int interpView) +{ + double dbeta = m_rotInc; + double dgamma = m_detInc; + double fanAngle = m_dFanBeamAngle; + int offsetView=0; + + // is there enough data in the data set? Should have 2(Pi+fanAngle) + // coverage minimum + if ( m_nView < static_cast((2*( PI + fanAngle ) ) / dbeta) -1 ){ + std::cerr << "Data set does not include 360 +2*FanBeamAngle views" + << std::endl; + return (1); + } + + if (interpView < 0) // use default position at PI+fanAngle + { + interpView = static_cast ((PI+fanAngle)/dbeta); + } + else + { + // check if there is PI+fanAngle data on either side of the + // of the specified image plane + if ( interpView*dbeta < PI+fanAngle || + interpView*dbeta + PI + fanAngle > m_nView*dbeta) + { + std::cerr << "There isn't PI+fanAngle of data on either side of the requested interpolation view" << std::endl; + return(1); + } + offsetView = interpView - static_cast((PI+fanAngle)/dbeta); + + } + int last_interp_view = static_cast ((PI+fanAngle)/dbeta); + + +// make a new array for data... + class DetectorArray ** newdetarray = new DetectorArray * [last_interp_view+1]; + for ( int i=0 ; i <= last_interp_view ; i++ ){ + newdetarray[i] = new DetectorArray (m_nDet); + newdetarray[i]->setViewAngle((i+offsetView)*dbeta); + DetectorValue* newdetval = (newdetarray[i])->detValues(); + // and initialize the data to zero + for (int j=0; j < m_nDet; j++) + newdetval[j] = 0.; + } + + int last_acq_view = 2*last_interp_view; + for ( int iView = 0 ; iView <= last_acq_view; iView++) { + double beta = iView * dbeta; + + for ( int iDet = 0; iDet < m_nDet; iDet++) { + double gamma = (iDet -(m_nDet-1)/2)* dgamma ; + int newiView, newiDet; + if (beta < PI+fanAngle) { //if (PI +fanAngle - beta > dbeta ) + //newbeta = beta; + //newgamma = gamma; + newiDet = iDet; + newiView = iView; + } + else // (beta > PI+fanAngle) + { + //newbeta = beta +2*gamma - 180; + //newgamma = -gamma; + newiDet = -iDet + (m_nDet -1); + // newiView = nearest((beta + 2*gamma - PI)/dbeta); + //newiView = static_cast(( (iView*dbeta) + 2*(iDet-(m_nDet-1)/2)*dgamma - PI)/dbeta); + newiView = nearest(( (iView*dbeta) + 2*(iDet-(m_nDet-1)/2)*dgamma - PI)/dbeta); + } + +#ifdef DEBUG +//std::cout << beta << " "<< gamma << " " << newbeta << " " << newgamma <<" " << iView-offsetView << " " << iDet << " " << newiView << " " << newiDet << std::endl; +//std::cout << iView-offsetView << " " << iDet << " " << newiView << " " << newiDet << std::endl; +#endif + + if ( ( beta > fanAngle - 2*gamma) + && ( beta < 2*PI + fanAngle -2*gamma) ) + { // not in region 1 or 8 + DetectorValue* detval = (m_projData[iView+offsetView])->detValues(); + DetectorValue* newdetval = (newdetarray[newiView])->detValues(); + if ( beta > fanAngle - 2*gamma + && beta <= 2*fanAngle ) { // in region 2 + newdetval[newiDet] += + (beta +2*gamma - fanAngle)/(PI+2*gamma) + * detval[iDet]; + } else if ( beta > 2*fanAngle + && beta <= PI - 2*gamma) { // in region 3 + newdetval[newiDet] += + (beta +2*gamma - fanAngle)/(PI+2*gamma) + * detval[iDet]; + } + else if ( beta > PI -2*gamma + && beta <= PI + fanAngle ) { // in region 4 + newdetval[newiDet] += + (beta +2*gamma - fanAngle)/(PI+2*gamma) + * detval[iDet]; + } + else if ( beta > PI + fanAngle + && beta <= PI +2*fanAngle -2*gamma) { // in region 5 + newdetval[newiDet] += + (2*PI - beta - 2*gamma + fanAngle)/(PI-2*gamma) + *detval[iDet]; + } + else if ( beta > PI +2*fanAngle -2*gamma + && beta <= 2*PI) { // in region 6 + newdetval[newiDet] += + (2*PI - beta - 2*gamma + fanAngle)/(PI-2*gamma) + *detval[iDet]; + } + else if ( beta > 2*PI + && beta <= 2*PI + fanAngle -2*gamma){ // in region 7 + newdetval[newiDet] += + (2*PI - beta -2*gamma + fanAngle)/(PI-2*gamma) + *detval[iDet]; + } + else + { + ; // outside region of interest + } + } + } + } + deleteProjData(); + m_projData = newdetarray; + m_nView = last_interp_view+1; + + return (0); +} +// HalfScanFeather: +// A HalfScan Projection Data Set for equiangular geometry, +// covering gantry angles from 0 to pi+fanBeamAngle +// and fan angle gamma from -fanBeamAngle/2 to fanBeamAngle/2 +// contains redundant information. If one copy of this data is left as +// zero, (as in the Helical180LI routine above) overweighting is avoided, +// but the discontinuity in the data introduces ringing in the image. +// This routine makes a copy of the data and applies a weighting to avoid +// over-representation, as given in Appendix C of Crawford and King, Med +// Phys 17 1990, p967. +int +Projections::HalfScanFeather(void) +{ + double dbeta = m_rotInc; + double dgamma = m_detInc; + double fanAngle = m_dFanBeamAngle; + +// is there enough data? + if ( m_nView != static_cast(( PI+fanAngle ) / dbeta) +1 ){ + std::cerr << "Data set does seem have enough data to be a halfscan data set" << std::endl; + return (1); + } + if (m_geometry == Scanner::GEOMETRY_INVALID) { + std::cerr << "Invalid geometry " << m_geometry << std::endl; + return (2); + } + + if (m_geometry == Scanner::GEOMETRY_PARALLEL) { + std::cerr << "HalfScanFeather not yet implemented for PARALLEL geometry"<< std::endl; + return (2); + } + + for ( int iView2 = 0 ; iView2 < m_nView; iView2++) { + double beta2 = iView2 * dbeta; + for ( int iDet2 = 0; iDet2 < m_nDet; iDet2++) { + double gamma2 = (iDet2 -(m_nDet-1)/2)* dgamma ; + if ( ( beta2 >= PI - 2*gamma2) ) { // in redundant data region + int iView1, iDet1; + iDet1 = (m_nDet -1) - iDet2; + //iView1 = nearest((beta2 + 2*gamma2 - PI)/dbeta); + iView1 = nearest(( (iView2*dbeta) + + 2*(iDet2-(m_nDet-1)/2)*dgamma - PI)/dbeta); + + + DetectorValue* detval2 = (m_projData[iView2])->detValues(); + DetectorValue* detval1 = (m_projData[iView1])->detValues(); + + detval1[iDet1] = detval2[iDet2] ; + + double x, w1,w2,beta1, gamma1; + beta1= iView1*dbeta; + gamma1 = -gamma2; + if ( beta1 <= (fanAngle - 2*gamma1) ) + x = beta1 / ( fanAngle - 2*gamma1); + else if ( (fanAngle - 2*gamma1 <= beta1 ) && beta1 <= PI - 2*gamma1) + x = 1; + else if ( (PI - 2*gamma1 <= beta1 ) && ( beta1 <=PI + fanAngle) ) + x = (PI +fanAngle - beta1)/(fanAngle + 2*gamma1); + else { + std::cerr << "Shouldn't be here!"<< std::endl; + return(4); + } + w1 = (3*x - 2*x*x)*x; + w2 = 1-w1; + detval1[iDet1] *= w1; + detval2[iDet2] *= w2; + + } + } + } + // heuristic scaling, why this factor? + double scalefactor = m_nView * m_rotInc / PI; + for ( int iView = 0 ; iView < m_nView; iView++) { + DetectorValue* detval = (m_projData[iView])->detValues(); + for ( int iDet = 0; iDet < m_nDet; iDet++) { + detval[iDet] *= scalefactor; + } + } + + return (0); +} + // NAME // newProjData @@ -378,7 +652,7 @@ Projections::read (const char* filename) #ifdef MSVC frnetorderstream fileRead (m_filename.c_str(), std::ios::in | std::ios::binary); #else - frnetorderstream fileRead (m_filename.c_str(), std::ios::in | std::ios::binary | std::ios::nocreate); + frnetorderstream fileRead (m_filename.c_str(), std::ios::in | std::ios::binary); // | std::ios::nocreate); #endif if (fileRead.fail()) @@ -697,21 +971,21 @@ Projections::convertPolar (ImageFile& rIF, int iInterpolationID) double** ppdView = adView.getArray(); double** ppdDet = adDet.getArray(); - std::complex** ppcDetValue = new std::complex* [m_nView]; - unsigned int iView; - for (iView = 0; iView < m_nView; iView++) { - ppcDetValue[iView] = new std::complex [m_nDet]; + std::complex** ppcDetValue = new std::complex* [pProj->m_nView]; + int iView; + for (iView = 0; iView < pProj->m_nView; iView++) { + ppcDetValue[iView] = new std::complex [pProj->m_nDet]; DetectorValue* detval = pProj->getDetectorArray (iView).detValues(); - for (unsigned int iDet = 0; iDet < m_nDet; iDet++) + for (int iDet = 0; iDet < pProj->m_nDet; iDet++) ppcDetValue[iView][iDet] = std::complex(detval[iDet], 0); } - pProj->calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet, m_nDet, 1., m_detInc); + pProj->calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet, pProj->m_nDet, 1., pProj->m_detInc); pProj->interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, pProj->m_nView, pProj->m_nDet, pProj->m_nDet, iInterpolationID); - for (iView = 0; iView < m_nView; iView++) + for (iView = 0; iView < pProj->m_nView; iView++) delete [] ppcDetValue[iView]; delete [] ppcDetValue; @@ -739,39 +1013,51 @@ Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad if (! v || nx == 0 || ny == 0) return false; - if (m_geometry != Scanner::GEOMETRY_PARALLEL) { - sys_error (ERR_WARNING, "convertFFTPolar supports Parallel only"); - return false; - } + Projections* pProj = this; + if (m_geometry == Scanner::GEOMETRY_EQUIANGULAR || m_geometry == Scanner::GEOMETRY_EQUILINEAR) + pProj = interpolateToParallel(); - int iInterpDet = nx; + int iInterpDet = static_cast(static_cast(sqrt(nx*nx+ny*ny))); int iNumInterpDetWithZeros = ProcessSignal::addZeropadFactor (iInterpDet, iZeropad); - + double dProjScale = iInterpDet / (pProj->viewDiameter() * 0.5); double dZeropadRatio = static_cast(iNumInterpDetWithZeros) / static_cast(iInterpDet); fftw_plan plan = fftw_create_plan (iNumInterpDetWithZeros, FFTW_FORWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM); fftw_complex* pcIn = new fftw_complex [iNumInterpDetWithZeros]; - std::complex** ppcDetValue = new std::complex* [m_nView]; - double dInterpScale = (m_nDet-1) / static_cast(iInterpDet-1) / SQRT2; + std::complex** ppcDetValue = new std::complex* [pProj->m_nView]; + //double dInterpScale = (pProj->m_nDet-1) / static_cast(iInterpDet-1); + double dInterpScale = pProj->m_nDet / static_cast(iInterpDet); - for (unsigned int iView = 0; iView < m_nView; iView++) { - DetectorValue* detval = getDetectorArray(iView).detValues(); - LinearInterpolator projInterp (detval, m_nDet); - for (unsigned int iDet = 0; iDet < iInterpDet; iDet++) { -// double dInterpPos = iInterpDet * dInterpScale; - double dInterpPos = (m_nDet / 2.) + (iDet - iInterpDet/2.) * dInterpScale; - pcIn[iDet].re = projInterp.interpolate (dInterpPos); + double dFFTScale = 1. / static_cast(iInterpDet * iInterpDet); + int iMidPoint = iInterpDet / 2; + double dMidPoint = static_cast(iInterpDet) / 2.; + int iZerosAdded = iNumInterpDetWithZeros - iInterpDet; + + // For each view, interpolate, shift to center at origin, and FFT + for (int iView = 0; iView < m_nView; iView++) { + DetectorValue* detval = pProj->getDetectorArray(iView).detValues(); + LinearInterpolator projInterp (detval, pProj->m_nDet); + for (int iDet = 0; iDet < iInterpDet; iDet++) { + double dInterpPos = (m_nDet / 2.) + (iDet - dMidPoint) * dInterpScale; + pcIn[iDet].re = projInterp.interpolate (dInterpPos) * dProjScale; pcIn[iDet].im = 0; } - for (unsigned int iDet2 = iInterpDet; iDet2 < iNumInterpDetWithZeros; iDet2++) - pcIn[iDet2].re = pcIn[iDet2].im = 0; + + Fourier::shuffleFourierToNaturalOrder (pcIn, iInterpDet); + if (iZerosAdded > 0) { + for (int iDet1 = iInterpDet -1; iDet1 >= iMidPoint; iDet1--) + pcIn[iDet1+iZerosAdded] = pcIn[iDet1]; + for (int iDet2 = iMidPoint; iDet2 < iMidPoint + iZerosAdded; iDet2++) + pcIn[iDet2].re = pcIn[iDet2].im = 0; + } fftw_one (plan, pcIn, NULL); ppcDetValue[iView] = new std::complex [iNumInterpDetWithZeros]; - for (unsigned int iD = 0; iD < iNumInterpDetWithZeros; iD++) - ppcDetValue[iView][iD] = std::complex (pcIn[iD].re / iInterpDet / (iInterpDet/2), pcIn[iD].im / iInterpDet / (iInterpDet/2)); + for (int iD = 0; iD < iNumInterpDetWithZeros; iD++) { + ppcDetValue[iView][iD] = std::complex (pcIn[iD].re * dFFTScale, pcIn[iD].im * dFFTScale); + } Fourier::shuffleFourierToNaturalOrder (ppcDetValue[iView], iNumInterpDetWithZeros); } @@ -783,11 +1069,14 @@ Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad Array2d adDet (nx, ny); double** ppdView = adView.getArray(); double** ppdDet = adDet.getArray(); - calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet, iNumInterpDetWithZeros, dZeropadRatio, - m_detInc * dInterpScale); + pProj->calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet, iNumInterpDetWithZeros, dZeropadRatio, + pProj->m_detInc * dInterpScale); - interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, m_nView, m_nDet, iNumInterpDetWithZeros, - iInterpolationID); + pProj->interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, pProj->m_nView, pProj->m_nDet, + iNumInterpDetWithZeros, iInterpolationID); + + if (m_geometry == Scanner::GEOMETRY_EQUIANGULAR || m_geometry == Scanner::GEOMETRY_EQUILINEAR) + delete pProj; for (int i = 0; i < m_nView; i++) delete [] ppcDetValue[i]; @@ -802,8 +1091,7 @@ void Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double** ppdView, double** ppdDet, int iNumDetWithZeros, double dZeropadRatio, double dDetInc) { -// double dLength = viewDiameter(); - double dLength = phmLen(); + double dLength = viewDiameter(); double xMin = -dLength / 2; double xMax = xMin + dLength; double yMin = -dLength / 2; @@ -819,10 +1107,10 @@ Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double double xInc = (xMax - xMin) / nx; // size of cells double yInc = (yMax - yMin) / ny; + double dDetCenter = (iNumDetWithZeros - 1) / 2.; // index refering to L=0 projection // +1 is correct for frequency data, ndet-1 is correct for projections - int iDetCenter = (iNumDetWithZeros - 1) / 2; // index refering to L=0 projection - if (isEven (iNumDetWithZeros)) - iDetCenter = (iNumDetWithZeros + 1) / 2; + // if (isEven (iNumDetWithZeros)) + // dDetCenter = (iNumDetWithZeros + 0) / 2; // Calculates polar coordinates (view#, det#) for each point on phantom grid double x = xMin + xInc / 2; // Rectang coords of center of pixel @@ -832,16 +1120,15 @@ Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double double r = ::sqrt (x * x + y * y); double phi = atan2 (y, x); - if (phi < 0) + if (phi <= -m_rotInc / 2) phi += TWOPI; - - if (phi >= PI) { + if (phi >= PI - (m_rotInc / 2)) { phi -= PI; r = -r; } ppdView[ix][iy] = (phi - m_rotStart) / m_rotInc; - ppdDet[ix][iy] = (r / dDetInc) + iDetCenter; + ppdDet[ix][iy] = (r / dDetInc) + dDetCenter; } } } @@ -852,72 +1139,41 @@ Projections::interpolatePolar (ImageFileArray& v, ImageFileArray& vImag, double** ppdDet, unsigned int nView, unsigned int nDet, unsigned int nDetWithZeros, int iInterpolationID) { typedef std::complex complexValue; - BilinearInterpolator bilinear (ppcDetValue, nView, nDetWithZeros); + BilinearPolarInterpolator* pBilinear = NULL; + BicubicPolyInterpolator* pBicubic = NULL; + if (iInterpolationID == POLAR_INTERP_BILINEAR) + pBilinear = new BilinearPolarInterpolator (ppcDetValue, nView, nDetWithZeros); + else if (iInterpolationID == POLAR_INTERP_BICUBIC) + pBicubic = new BicubicPolyInterpolator (ppcDetValue, nView, nDetWithZeros); + for (unsigned int ix = 0; ix < ny; ix++) { for (unsigned int iy = 0; iy < ny; iy++) { - if (iInterpolationID == POLAR_INTERP_NEAREST) { unsigned int iView = nearest (ppdView[ix][iy]); unsigned int iDet = nearest (ppdDet[ix][iy]); - if (iView == nView) { - iView = 0; - iDet = m_nDet - iDet; - } + if (iView == nView) + iView = 0; if (iDet >= 0 && iDet < nDetWithZeros && iView >= 0 && iView < nView) { v[ix][iy] = ppcDetValue[iView][iDet].real(); if (vImag) vImag[ix][iy] = ppcDetValue[iView][iDet].imag(); - } else + } else { v[ix][iy] = 0; + if (vImag) + vImag[ix][iy] = 0; + } } else if (iInterpolationID == POLAR_INTERP_BILINEAR) { -#if 1 - std::complex vInterp = bilinear.interpolate (ppdView[ix][iy], ppdDet[ix][iy]); + std::complex vInterp = pBilinear->interpolate (ppdView[ix][iy], ppdDet[ix][iy]); v[ix][iy] = vInterp.real(); if (vImag) vImag[ix][iy] = vInterp.imag(); -#else - int iFloorView = ::floor (ppdView[ix][iy]); - double dFracView = ppdView[ix][iy] - iFloorView; - int iFloorDet = ::floor (ppdDet[ix][iy]); - double dFracDet = ppdDet[ix][iy] - iFloorDet; - - if (iFloorDet >= 0 && iFloorView >= 0) { - std::complex v1 = ppcDetValue[iFloorView][iFloorDet]; - std::complex v2, v3, v4; - if (iFloorView < nView - 1) - v2 = ppcDetValue[iFloorView + 1][iFloorDet]; - else - v2 = ppcDetValue[0][iFloorDet]; - if (iFloorDet < nDetWithZeros - 1) - v4 = ppcDetValue[iFloorView][iFloorDet+1]; - else - v4 = v1; - if (iFloorView < nView - 1 && iFloorDet < nDetWithZeros - 1) - v3 = ppcDetValue [iFloorView+1][iFloorDet+1]; - else if (iFloorView < nView - 1) - v3 = v2; - else - v3 = ppcDetValue[0][iFloorDet+1]; - - std::complex vInterp = (1 - dFracView) * (1 - dFracDet) * v1 + - dFracView * (1 - dFracDet) * v2 + dFracView * dFracDet * v3 + - dFracDet * (1 - dFracView) * v4; - v[ix][iy] = vInterp.real(); - if (vImag) - vImag[ix][iy] = vInterp.imag(); - } else { - // sys_error (ERR_SEVERE, "Can't find projection data for ix=%d,iy=%d with radView=%f and radDet=%f", ix, iy, ppdView[ix][iy], ppdDet[ix][iy]); - v[ix][iy] = 0; - if (vImag) - vImag[ix][iy] = 0; - } -#endif } else if (iInterpolationID == POLAR_INTERP_BICUBIC) { - v[ix][iy] =0; - if (vImag) - vImag[ix][iy] = 0; + std::complex vInterp = pBicubic->interpolate (ppdView[ix][iy], ppdDet[ix][iy]); + v[ix][iy] = vInterp.real(); + if (vImag) + vImag[ix][iy] = vInterp.imag(); } } } @@ -950,22 +1206,22 @@ Projections::initFromSomatomAR_STAR (int iNViews, int iNDets, unsigned char* pDa for (int iv = 0; iv < iNViews; iv++) { unsigned char* pArgBase = pData + lDataPos; unsigned char* p = pArgBase+0; SwapBytes4IfLittleEndian (p); - long lProjNumber = *reinterpret_cast(p); + // long lProjNumber = *reinterpret_cast(p); p = pArgBase+20; SwapBytes4IfLittleEndian (p); long lEscale = *reinterpret_cast(p); p = pArgBase+28; SwapBytes4IfLittleEndian (p); - long lTime = *reinterpret_cast(p); + // long lTime = *reinterpret_cast(p); p = pArgBase + 4; SwapBytes4IfLittleEndian (p); double dAlpha = *reinterpret_cast(p) + HALFPI; p = pArgBase+12; SwapBytes4IfLittleEndian (p); - double dAlign = *reinterpret_cast(p); + // double dAlign = *reinterpret_cast(p); p = pArgBase + 16; SwapBytes4IfLittleEndian (p); - double dMaxValue = *reinterpret_cast(p); + // double dMaxValue = *reinterpret_cast(p); DetectorArray& detArray = getDetectorArray (iv); detArray.setViewAngle (dAlpha); @@ -1036,9 +1292,9 @@ Projections::interpolateToParallel () const double dViewAngle = m_rotStart; int iLastFloor = -1; + LinearInterpolator interp (pdThetaValuesForT, pdRaysumsForT, pProjNew->nView(), false); for (int iV = 0; iV < pProjNew->nView(); iV++, dViewAngle += pProjNew->m_rotInc) { DetectorValue* detValues = pProjNew->getDetectorArray (iV).detValues(); - LinearInterpolator interp (pdThetaValuesForT, pdRaysumsForT, pProjNew->nView()); detValues[iD] = interp.interpolate (dViewAngle, &iLastFloor); } } @@ -1057,11 +1313,11 @@ Projections::interpolateToParallel () const detArray.setViewAngle (dViewAngle); for (int i = 0; i < pProjNew->nDet(); i++) - pdDetValueCopy[i] = detValues[i]; + pdDetValueCopy[i] = detValues[i]; double dDetPos = pProjNew->m_detStart; int iLastFloor = -1; - LinearInterpolator interp (pdOriginalDetPositions, pdDetValueCopy, pProjNew->nDet()); + LinearInterpolator interp (pdOriginalDetPositions, pdDetValueCopy, pProjNew->nDet(), false); for (int iD = 0; iD < pProjNew->nDet(); iD++, dDetPos += pProjNew->m_detInc) detValues[iD] = interp.interpolate (dDetPos, &iLastFloor); } @@ -1081,8 +1337,8 @@ Projections::interpolateToParallel () const /////////////////////////////////////////////////////////////////////////////// ParallelRaysums::ParallelRaysums (const Projections* pProjections, int iThetaRange) -: m_iNumCoordinates(0), m_iNumView(pProjections->nView()), m_iNumDet(pProjections->nDet()), - m_iThetaRange (iThetaRange), m_pCoordinates(NULL) +: m_pCoordinates(NULL), m_iNumCoordinates(0), m_iNumView(pProjections->nView()), m_iNumDet(pProjections->nDet()), + m_iThetaRange (iThetaRange) { int iGeometry = pProjections->geometry(); double dDetInc = pProjections->detInc();