X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fprojections.cpp;h=2ca1ce38e12af9fbf708595028df1f9dc9c3fbc8;hp=d6a484335943945d8f4ffa02a6bb77e495412fb7;hb=c0f892798de8f89715266150f7d8e413f2cf29fe;hpb=d850a3477e9ccaecfa85e00bc619848fcc29bdb6 diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp index d6a4843..2ca1ce3 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.56 2001/03/10 23:56:58 kevin Exp $ +** $Id: projections.cpp,v 1.77 2002/05/28 18:43:16 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 @@ -50,6 +50,7 @@ const char* const Projections::s_aszInterpTitle[] = const int Projections::s_iInterpCount = sizeof(s_aszInterpName) / sizeof(char*); + /* NAME * Projections Constructor for projections matrix storage * @@ -150,12 +151,12 @@ Projections::initFromScanner (const Scanner& scanner) m_rotInc = scanner.rotInc(); m_detInc = scanner.detInc(); + m_detStart = scanner.detStart(); m_geometry = scanner.geometry(); m_dFocalLength = scanner.focalLength(); m_dSourceDetectorLength = scanner.sourceDetectorLength(); m_dViewDiameter = scanner.viewDiameter(); - m_rotStart = 0; - m_detStart = -(scanner.detLen() / 2); + m_rotStart = scanner.offsetView()*scanner.rotInc(); m_dFanBeamAngle = scanner.fanBeamAngle(); } @@ -166,6 +167,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 +650,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()) @@ -688,33 +960,36 @@ Projections::convertPolar (ImageFile& rIF, int iInterpolationID) if (! v || nx == 0 || ny == 0) return false; - if (m_geometry != Scanner::GEOMETRY_PARALLEL) { - sys_error (ERR_WARNING, "convertPolar supports Parallel only"); - return false; - } + Projections* pProj = this; + if (m_geometry == Scanner::GEOMETRY_EQUIANGULAR || m_geometry == Scanner::GEOMETRY_EQUILINEAR) + pProj = interpolateToParallel(); Array2d adView (nx, ny); Array2d adDet (nx, ny); double** ppdView = adView.getArray(); double** ppdDet = adDet.getArray(); - if (! calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet)) - return false; - - std::complex** ppcDetValue = new std::complex* [m_nView]; + std::complex** ppcDetValue = new std::complex* [pProj->m_nView]; unsigned int iView; - for (iView = 0; iView < m_nView; iView++) { - ppcDetValue[iView] = new std::complex [m_nDet]; - for (unsigned int iDet = 0; iDet < m_nDet; iDet++) - ppcDetValue[iView][iDet] = std::complex(getDetectorArray (iView).detValues()[iDet], 0); + 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 < pProj->m_nDet; iDet++) + ppcDetValue[iView][iDet] = std::complex(detval[iDet], 0); } - interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, m_nView, m_nDet, iInterpolationID); + 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; + if (m_geometry == Scanner::GEOMETRY_EQUIANGULAR || m_geometry == Scanner::GEOMETRY_EQUILINEAR) + delete pProj; + return true; } @@ -722,6 +997,10 @@ Projections::convertPolar (ImageFile& rIF, int iInterpolationID) bool Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad) { +#ifndef HAVE_FFTW + rIF.arrayDataClear(); + return false; +#else unsigned int nx = rIF.nx(); unsigned int ny = rIF.ny(); ImageFileArray v = rIF.getArray(); @@ -732,67 +1011,105 @@ 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; - } - -#ifndef HAVE_FFT - return false; -#else - Array2d adView (nx, ny); - Array2d adDet (nx, ny); - double** ppdView = adView.getArray(); - double** ppdDet = adDet.getArray(); + Projections* pProj = this; + if (m_geometry == Scanner::GEOMETRY_EQUIANGULAR || m_geometry == Scanner::GEOMETRY_EQUILINEAR) + pProj = interpolateToParallel(); - std::complex** ppcDetValue = new std::complex* [m_nView]; - unsigned int iView; - double* pdDet = new double [m_nDet]; - fftw_complex* pcIn = new fftw_complex [m_nDet]; - fftw_plan plan = fftw_create_plan (m_nDet, FFTW_FORWARD, FFTW_IN_PLACE); - - for (iView = 0; iView < m_nView; iView++) { - unsigned int iDet; - for (iDet = 0; iDet < m_nDet; iDet++) { - pcIn[iDet].re = getDetectorArray(iView).detValues()[iDet]; + int iInterpDet = nx; +// int iInterpDet = pProj->m_nDet; + int iNumInterpDetWithZeros = ProcessSignal::addZeropadFactor (iInterpDet, iZeropad); + + 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* [pProj->m_nView]; + double dInterpScale = (pProj->m_nDet-1) / static_cast(iInterpDet-1) / SQRT2; + + 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 to nx length, shift to center at origin, and FFt transform + for (unsigned int iView = 0; iView < m_nView; iView++) { + DetectorValue* detval = pProj->getDetectorArray(iView).detValues(); + LinearInterpolator projInterp (detval, pProj->m_nDet); + for (unsigned int iDet = 0; iDet < iInterpDet; iDet++) { + double dInterpPos = (m_nDet / 2.) + (iDet - dMidPoint) * dInterpScale; + pcIn[iDet].re = projInterp.interpolate (dInterpPos) * dInterpScale; pcIn[iDet].im = 0; } + + Fourier::shuffleFourierToNaturalOrder (pcIn, iInterpDet); + if (iZerosAdded > 0) { + for (unsigned int iDet1 = iMidPoint; iDet1 < iInterpDet; iDet1++) + pcIn[iDet1+iZerosAdded] = pcIn[iDet1]; + for (unsigned int iDet2 = iMidPoint; iDet2 < iMidPoint + iZerosAdded; iDet2++) + pcIn[iDet2].re = pcIn[iDet2].im = 0; + } + fftw_one (plan, pcIn, NULL); - ppcDetValue[iView] = new std::complex [m_nDet]; - for (iDet = 0; iDet < m_nDet; iDet++) - ppcDetValue[iView][iDet] = std::complex (pcIn[iDet].re, pcIn[iDet].im); - Fourier::shuffleFourierToNaturalOrder (ppcDetValue[iView], m_nDet); + + ppcDetValue[iView] = new std::complex [iNumInterpDetWithZeros]; + for (unsigned int iD = 0; iD < iNumInterpDetWithZeros; iD++) { + ppcDetValue[iView][iD] = std::complex (pcIn[iD].re * dFFTScale, pcIn[iD].im * dFFTScale); + } + + Fourier::shuffleFourierToNaturalOrder (ppcDetValue[iView], iNumInterpDetWithZeros); } + delete [] pcIn; fftw_destroy_plan (plan); - delete [] pcIn; - bool bError = calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet); + Array2d adView (nx, ny); + Array2d adDet (nx, ny); + double** ppdView = adView.getArray(); + double** ppdDet = adDet.getArray(); + pProj->calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet, iNumInterpDetWithZeros, dZeropadRatio, + pProj->m_detInc * dInterpScale); - if (! bError) - interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, m_nView, m_nDet, iInterpolationID); + pProj->interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, pProj->m_nView, pProj->m_nDet, + iNumInterpDetWithZeros, iInterpolationID); - for (iView = 0; iView < m_nView; iView++) - delete [] ppcDetValue[iView]; + if (m_geometry == Scanner::GEOMETRY_EQUIANGULAR || m_geometry == Scanner::GEOMETRY_EQUILINEAR) + delete pProj; + + for (int i = 0; i < m_nView; i++) + delete [] ppcDetValue[i]; delete [] ppcDetValue; - return bError; + return true; #endif } -bool -Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double** ppdView, double** ppdDet) +void +Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double** ppdView, double** ppdDet, + int iNumDetWithZeros, double dZeropadRatio, double dDetInc) { - double xMin = -phmLen() / 2; - double xMax = xMin + phmLen(); - double yMin = -phmLen() / 2; - double yMax = yMin + phmLen(); - + double dLength = viewDiameter(); +// double dLength = phmLen(); + double xMin = -dLength / 2; + double xMax = xMin + dLength; + double yMin = -dLength / 2; + double yMax = yMin + dLength; + double xCent = (xMin + xMax) / 2; + double yCent = (yMin + yMax) / 2; + + xMin = (xMin - xCent) * dZeropadRatio + xCent; + xMax = (xMax - xCent) * dZeropadRatio + xCent; + yMin = (yMin - yCent) * dZeropadRatio + yCent; + yMax = (yMax - yCent) * dZeropadRatio + yCent; + double xInc = (xMax - xMin) / nx; // size of cells double yInc = (yMax - yMin) / ny; - - int iDetCenter = (m_nDet - 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; // Calculates polar coordinates (view#, det#) for each point on phantom grid double x = xMin + xInc / 2; // Rectang coords of center of pixel @@ -802,210 +1119,268 @@ Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double double r = ::sqrt (x * x + y * y); double phi = atan2 (y, x); + if (phi < 0) + phi += TWOPI; if (phi >= PI) { phi -= PI; - } else if (phi < 0) { - phi += PI; - } else r = -r; + } ppdView[ix][iy] = (phi - m_rotStart) / m_rotInc; - ppdDet[ix][iy] = (r / m_detInc) + iDetCenter; + ppdDet[ix][iy] = (r / dDetInc) + iDetCenter; } } - - return true; } void Projections::interpolatePolar (ImageFileArray& v, ImageFileArray& vImag, - unsigned int nx, unsigned int ny, std::complex** ppcDetValue, - double** ppdView, double** ppdDet, unsigned int nView, unsigned int nDet, int iInterpolationID) + unsigned int nx, unsigned int ny, std::complex** ppcDetValue, double** ppdView, + double** ppdDet, unsigned int nView, unsigned int nDet, unsigned int nDetWithZeros, int iInterpolationID) { + typedef std::complex complexValue; + + BilinearInterpolator* pBilinear; + if (iInterpolationID == POLAR_INTERP_BILINEAR) + pBilinear = new BilinearInterpolator (ppcDetValue, nView, nDetWithZeros); + + BicubicPolyInterpolator* pBicubic; + 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; + iDet = m_nDet - iDet; } - if (iDet >= 0 && iDet < nDet && iView >= 0 && iView < nView) { + 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 { - 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]); + } else v[ix][iy] = 0; - } + } else if (iInterpolationID == POLAR_INTERP_BILINEAR) { - unsigned int iFloorView = static_cast(ppdView[ix][iy]); - double dFracView = ppdView[ix][iy] - iFloorView; - unsigned int iFloorDet = static_cast(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 < nDet - 1) - v4 = ppcDetValue[iFloorView][iFloorDet+1]; - else - v4 = v1; - if (iFloorView < nView - 1 && iFloorDet < nDet - 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; - } + std::complex vInterp = pBilinear->interpolate (ppdView[ix][iy], ppdDet[ix][iy]); + v[ix][iy] = vInterp.real(); + if (vImag) + vImag[ix][iy] = vInterp.imag(); } 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(); } } } } - bool Projections::initFromSomatomAR_STAR (int iNViews, int iNDets, unsigned char* pData, unsigned long lDataLength) { init (iNViews, iNDets); m_geometry = Scanner::GEOMETRY_EQUIANGULAR; - m_dFanBeamAngle = iNDets * convertDegreesToRadians (3.06976 / 60); - m_dFocalLength = 51; - m_dSourceDetectorLength = 89; + m_dFocalLength = 510; + m_dSourceDetectorLength = 890; m_detInc = convertDegreesToRadians (3.06976 / 60); - m_detStart = -m_dFanBeamAngle / 2; + m_dFanBeamAngle = iNDets * m_detInc; + m_detStart = -(m_dFanBeamAngle / 2); m_rotInc = TWOPI / static_cast(iNViews); - m_rotStart = HALFPI; + m_rotStart = 0; m_dViewDiameter = sin (m_dFanBeamAngle / 2) * m_dFocalLength * 2; - if (iNDets != 1024) - return false; - bool bValid = (iNViews == 750 && lDataLength == 1560000L) || (iNViews == 950 && lDataLength == 1976000L) || (iNViews == 1500 && lDataLength == 3120000); - if (! bValid) + if (! ((iNViews == 750 && lDataLength == 1560000L) || (iNViews == 950 && lDataLength == 1976000L) + || (iNViews == 1500 && lDataLength == 3120000))) return false; + double dCenter = (iNDets - 1.) / 2.; // change from (Nm+1)/2 because of 0 vs. 1 indexing + double* pdCosScale = new double [iNDets]; + for (int i = 0; i < iNDets; i++) + pdCosScale[i] = 1. / cos ((i - dCenter) * m_detInc); + long lDataPos = 0; for (int iv = 0; iv < iNViews; iv++) { unsigned char* pArgBase = pData + lDataPos; - unsigned char* p = pArgBase+0; -#ifndef WORDS_BIGENDIAN - SwapBytes4 (p); -#endif + unsigned char* p = pArgBase+0; SwapBytes4IfLittleEndian (p); long lProjNumber = *reinterpret_cast(p); - p = pArgBase+20; -#ifndef WORDS_BIGENDIAN - SwapBytes4 (p); -#endif + p = pArgBase+20; SwapBytes4IfLittleEndian (p); long lEscale = *reinterpret_cast(p); - p = pArgBase+28; -#ifndef WORDS_BIGENDIAN - SwapBytes4 (p); -#endif + p = pArgBase+28; SwapBytes4IfLittleEndian (p); long lTime = *reinterpret_cast(p); - p = pArgBase + 4; -#ifndef WORDS_BIGENDIAN - SwapBytes4 (p); -#endif + p = pArgBase + 4; SwapBytes4IfLittleEndian (p); double dAlpha = *reinterpret_cast(p) + HALFPI; - p = pArgBase+12; -#ifndef WORDS_BIGENDIAN - SwapBytes4 (p); -#endif + p = pArgBase+12; SwapBytes4IfLittleEndian (p); double dAlign = *reinterpret_cast(p); - p = pArgBase + 16; -#ifndef WORDS_BIGENDIAN - SwapBytes4 (p); -#endif + p = pArgBase + 16; SwapBytes4IfLittleEndian (p); double dMaxValue = *reinterpret_cast(p); - lDataPos += 32; - double dEScale = pow (2.0, -lEscale); - double dBetaInc = convertDegreesToRadians (3.06976 / 60); - int iCenter = (iNDets + 1) / 2; - - DetectorArray& detArray = getDetectorArray( iv ); + DetectorArray& detArray = getDetectorArray (iv); detArray.setViewAngle (dAlpha); DetectorValue* detval = detArray.detValues(); - double dTempScale = 2294.4871 * dEScale; + double dViewScale = 1. / (2294.4871 * ::pow (2.0, -lEscale)); + lDataPos += 32; for (int id = 0; id < iNDets; id++) { - int iV = pData[lDataPos+1] + 256 * pData[lDataPos]; + int iV = pData[lDataPos+1] + (pData[lDataPos] << 8); if (iV > 32767) // two's complement signed conversion iV = iV - 65536; - double dCosScale = cos ((id + 1 - iCenter) * dBetaInc); - detval[id] = iV / (dTempScale * dCosScale); + detval[id] = iV * dViewScale * pdCosScale[id]; lDataPos += 2; } +#if 1 + for (int k = iNDets - 2; k >= 0; k--) + detval[k+1] = detval[k]; + detval[0] = 0; +#endif } + delete pdCosScale; return true; } +Projections* +Projections::interpolateToParallel () const +{ + if (m_geometry == Scanner::GEOMETRY_PARALLEL) + return const_cast(this); + + int nDet = m_nDet; + int nView = m_nView; + Projections* pProjNew = new Projections (nView, nDet); + pProjNew->m_geometry = Scanner::GEOMETRY_PARALLEL; + pProjNew->m_dFocalLength = m_dFocalLength; + pProjNew->m_dSourceDetectorLength = m_dSourceDetectorLength; + pProjNew->m_dViewDiameter = m_dViewDiameter; + pProjNew->m_dFanBeamAngle = m_dFanBeamAngle; + pProjNew->m_calcTime = 0; + pProjNew->m_remark = m_remark; + pProjNew->m_remark += "; Interpolate to Parallel"; + pProjNew->m_label.setLabelType (Array2dFileLabel::L_HISTORY); + pProjNew->m_label.setLabelString (pProjNew->m_remark); + pProjNew->m_label.setCalcTime (pProjNew->m_calcTime); + pProjNew->m_label.setDateTime (pProjNew->m_year, pProjNew->m_month, pProjNew->m_day, pProjNew->m_hour, pProjNew->m_minute, pProjNew->m_second); + + pProjNew->m_rotStart = 0; +#ifdef CONVERT_PARALLEL_PI + pProjNew->m_rotInc = PI / nView;; +#else + pProjNew->m_rotInc = TWOPI / nView; +#endif + pProjNew->m_detStart = -m_dViewDiameter / 2; + pProjNew->m_detInc = m_dViewDiameter / nDet; + if (isEven (nDet)) // even + pProjNew->m_detInc = m_dViewDiameter / (nDet - 1); + + ParallelRaysums parallel (this, ParallelRaysums::THETA_RANGE_NORMALIZE_TO_TWOPI); + + double* pdThetaValuesForT = new double [pProjNew->nView()]; + double* pdRaysumsForT = new double [pProjNew->nView()]; + + // interpolate to evenly spaced theta (views) + double dDetPos = pProjNew->m_detStart; + for (int iD = 0; iD < pProjNew->nDet(); iD++, dDetPos += pProjNew->m_detInc) { + parallel.getThetaAndRaysumsForT (iD, pdThetaValuesForT, pdRaysumsForT); + + double dViewAngle = m_rotStart; + int iLastFloor = -1; + for (int iV = 0; iV < pProjNew->nView(); iV++, dViewAngle += pProjNew->m_rotInc) { + DetectorValue* detValues = pProjNew->getDetectorArray (iV).detValues(); + LinearInterpolator interp (pdThetaValuesForT, pdRaysumsForT, pProjNew->nView(), false); + detValues[iD] = interp.interpolate (dViewAngle, &iLastFloor); + } + } + delete pdThetaValuesForT; + delete pdRaysumsForT; + + // interpolate to evenly space t (detectors) + double* pdOriginalDetPositions = new double [pProjNew->nDet()]; + parallel.getDetPositions (pdOriginalDetPositions); + + double* pdDetValueCopy = new double [pProjNew->nDet()]; + double dViewAngle = m_rotStart; + for (int iV = 0; iV < pProjNew->nView(); iV++, dViewAngle += pProjNew->m_rotInc) { + DetectorArray& detArray = pProjNew->getDetectorArray (iV); + DetectorValue* detValues = detArray.detValues(); + detArray.setViewAngle (dViewAngle); + + for (int i = 0; i < pProjNew->nDet(); i++) + pdDetValueCopy[i] = detValues[i]; + + double dDetPos = pProjNew->m_detStart; + int iLastFloor = -1; + 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); + } + delete pdDetValueCopy; + delete pdOriginalDetPositions; + + return pProjNew; +} + + +/////////////////////////////////////////////////////////////////////////////// +// +// Class ParallelRaysums +// +// Used for converting divergent beam raysums into Parallel raysums +// +/////////////////////////////////////////////////////////////////////////////// -ParallelRaysums::ParallelRaysums (Projections* pProjections) -: m_iNumCoordinates(0) +ParallelRaysums::ParallelRaysums (const Projections* pProjections, int iThetaRange) +: m_iNumCoordinates(0), m_iNumView(pProjections->nView()), m_iNumDet(pProjections->nDet()), + m_iThetaRange (iThetaRange), m_pCoordinates(NULL) { - int nDet = pProjections->nDet(); - int nView = pProjections->nView(); int iGeometry = pProjections->geometry(); double dDetInc = pProjections->detInc(); double dDetStart = pProjections->detStart(); double dFocalLength = pProjections->focalLength(); - m_iNumCoordinates = nDet * nView; + m_iNumCoordinates = m_iNumView * m_iNumDet; + m_pCoordinates = new ParallelRaysumCoordinate [m_iNumCoordinates]; m_vecpCoordinates.reserve (m_iNumCoordinates); for (int i = 0; i < m_iNumCoordinates; i++) - m_vecpCoordinates[i] = new ParallelRaysumCoordinate; + m_vecpCoordinates[i] = m_pCoordinates + i; int iCoordinate = 0; - for (int iV = 0; iV < nView; iV++) { + for (int iV = 0; iV < m_iNumView; iV++) { double dViewAngle = pProjections->getDetectorArray(iV).viewAngle(); + const DetectorValue* detValues = pProjections->getDetectorArray(iV).detValues(); double dDetPos = dDetStart; - for (int iD = 0; iD < nDet; iD++) { + for (int iD = 0; iD < m_iNumDet; iD++) { ParallelRaysumCoordinate* pC = m_vecpCoordinates[iCoordinate++]; if (iGeometry == Scanner::GEOMETRY_PARALLEL) { - pC->m_dTheta = normalizeAngle (dViewAngle); + pC->m_dTheta = dViewAngle; pC->m_dT = dDetPos; - } else if (iGeometry == Scanner::GEOMETRY_EQUILINEAR) { double dFanAngle = atan (dDetPos / pProjections->sourceDetectorLength()); - pC->m_dTheta = normalizeAngle (dViewAngle + dFanAngle); + pC->m_dTheta = dViewAngle + dFanAngle; pC->m_dT = dFocalLength * sin(dFanAngle); } else if (iGeometry == Scanner::GEOMETRY_EQUIANGULAR) { // fan angle is same as dDetPos - pC->m_dTheta = normalizeAngle (dViewAngle + dDetPos); + pC->m_dTheta = dViewAngle + dDetPos; pC->m_dT = dFocalLength * sin (dDetPos); } + if (m_iThetaRange != THETA_RANGE_UNCONSTRAINED) { + pC->m_dTheta = normalizeAngle (pC->m_dTheta); + if (m_iThetaRange == THETA_RANGE_FOLD_TO_PI && pC->m_dTheta >= PI) { + pC->m_dTheta -= PI; + pC->m_dT = -pC->m_dT; + } + } + pC->m_dRaysum = detValues[iD]; dDetPos += dDetInc; } } @@ -1013,15 +1388,14 @@ ParallelRaysums::ParallelRaysums (Projections* pProjections) ParallelRaysums::~ParallelRaysums() { - for (int i = 0; i < m_iNumCoordinates; i++) - delete m_vecpCoordinates[i]; + delete m_pCoordinates; } ParallelRaysums::CoordinateContainer& ParallelRaysums::getSortedByTheta() { if (m_vecpSortedByTheta.size() == 0) { - m_vecpSortedByTheta.reserve (m_iNumCoordinates); + m_vecpSortedByTheta.resize (m_iNumCoordinates); for (int i = 0; i < m_iNumCoordinates; i++) m_vecpSortedByTheta[i] = m_vecpCoordinates[i]; std::sort (m_vecpSortedByTheta.begin(), m_vecpSortedByTheta.end(), ParallelRaysumCoordinate::compareByTheta); @@ -1034,7 +1408,7 @@ ParallelRaysums::CoordinateContainer& ParallelRaysums::getSortedByT() { if (m_vecpSortedByT.size() == 0) { - m_vecpSortedByT.reserve (m_iNumCoordinates); + m_vecpSortedByT.resize (m_iNumCoordinates); for (int i = 0; i < m_iNumCoordinates; i++) m_vecpSortedByT[i] = m_vecpCoordinates[i]; std::sort (m_vecpSortedByT.begin(), m_vecpSortedByT.end(), ParallelRaysumCoordinate::compareByT); @@ -1068,3 +1442,28 @@ ParallelRaysums::getLimits (double* dMinT, double* dMaxT, double* dMinTheta, dou *dMaxTheta = dTheta; } } + +void +ParallelRaysums::getThetaAndRaysumsForT (int iTheta, double* pTheta, double* pRaysum) +{ + const CoordinateContainer& coordsT = getSortedByT(); + + int iBase = iTheta * m_iNumView; + for (int i = 0; i < m_iNumView; i++) { + int iPos = iBase + i; + pTheta[i] = coordsT[iPos]->m_dTheta; + pRaysum[i] = coordsT[iPos]->m_dRaysum; + } +} + +void +ParallelRaysums::getDetPositions (double* pdDetPos) +{ + const CoordinateContainer& coordsT = getSortedByT(); + + int iPos = 0; + for (int i = 0; i < m_iNumDet; i++) { + pdDetPos[i] = coordsT[iPos]->m_dT; + iPos += m_iNumView; + } +}