X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fprojections.cpp;h=2ca1ce38e12af9fbf708595028df1f9dc9c3fbc8;hp=2a80405a43dc1eb0356f1ddaf6d2d9bb6c58ed7e;hb=c0f892798de8f89715266150f7d8e413f2cf29fe;hpb=320062d4d1e53ab9da799c77e9886b1e36cddcf5 diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp index 2a80405..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.70 2001/03/24 05:28:28 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 * @@ -155,7 +156,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 +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()) @@ -739,12 +1011,12 @@ 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 = pProj->m_nDet; int iNumInterpDetWithZeros = ProcessSignal::addZeropadFactor (iInterpDet, iZeropad); double dZeropadRatio = static_cast(iNumInterpDetWithZeros) / static_cast(iInterpDet); @@ -752,28 +1024,38 @@ Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad 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) / SQRT2; double dFFTScale = 1. / static_cast(iInterpDet * iInterpDet); - dFFTScale /= 2; // Not sure why this scaling is necessary + 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 = getDetectorArray(iView).detValues(); - LinearInterpolator projInterp (detval, m_nDet); + DetectorValue* detval = pProj->getDetectorArray(iView).detValues(); + LinearInterpolator projInterp (detval, pProj->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 dInterpPos = (m_nDet / 2.) + (iDet - dMidPoint) * dInterpScale; + pcIn[iDet].re = projInterp.interpolate (dInterpPos) * dInterpScale; 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 (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 [iNumInterpDetWithZeros]; - for (unsigned int iD = 0; iD < iNumInterpDetWithZeros; iD++) + 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); } @@ -785,11 +1067,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]; @@ -804,8 +1089,8 @@ 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 dLength = phmLen(); double xMin = -dLength / 2; double xMax = xMin + dLength; double yMin = -dLength / 2;