X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fprojections.cpp;h=871f62e6e52a4543df9262355872d3b22957f392;hp=86976c7235a8a371755660dc92001f17bde3a635;hb=8a7697ce57b56cdc43698cd1241ad98d49f9b5ac;hpb=bd1d464294e037da19ccc80d8cc60475768eb2ca diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp index 86976c7..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.74 2001/09/24 09:40:42 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,15 +37,15 @@ 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"}, }; @@ -238,30 +240,30 @@ Projections::Helical180LI_Equiangular(int interpView) // is there enough data in the data set? Should have 2(Pi+fanAngle) // coverage minimum - if ( m_nView < static_cast((2*( M_PI + fanAngle ) ) / dbeta) -1 ){ + 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 M_PI+fanAngle + if (interpView < 0) // use default position at PI+fanAngle { - interpView = static_cast ((M_PI+fanAngle)/dbeta); + interpView = static_cast ((PI+fanAngle)/dbeta); } else { - // check if there is M_PI+fanAngle data on either side of the + // check if there is PI+fanAngle data on either side of the // of the specified image plane - if ( interpView*dbeta < M_PI+fanAngle || - interpView*dbeta + M_PI + fanAngle > m_nView*dbeta) + 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((M_PI+fanAngle)/dbeta); + offsetView = interpView - static_cast((PI+fanAngle)/dbeta); } - int last_interp_view = static_cast ((M_PI+fanAngle)/dbeta); + int last_interp_view = static_cast ((PI+fanAngle)/dbeta); // make a new array for data... @@ -282,20 +284,20 @@ Projections::Helical180LI_Equiangular(int interpView) for ( int iDet = 0; iDet < m_nDet; iDet++) { double gamma = (iDet -(m_nDet-1)/2)* dgamma ; int newiView, newiDet; - if (beta < M_PI+fanAngle) { //if (M_PI +fanAngle - beta > dbeta ) + if (beta < PI+fanAngle) { //if (PI +fanAngle - beta > dbeta ) //newbeta = beta; //newgamma = gamma; newiDet = iDet; newiView = iView; } - else // (beta > M_PI+fanAngle) + else // (beta > PI+fanAngle) { //newbeta = beta +2*gamma - 180; //newgamma = -gamma; newiDet = -iDet + (m_nDet -1); - // newiView = nearest((beta + 2*gamma - M_PI)/dbeta); - //newiView = static_cast(( (iView*dbeta) + 2*(iDet-(m_nDet-1)/2)*dgamma - M_PI)/dbeta); - newiView = nearest(( (iView*dbeta) + 2*(iDet-(m_nDet-1)/2)*dgamma - M_PI)/dbeta); + // 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 @@ -304,43 +306,43 @@ Projections::Helical180LI_Equiangular(int interpView) #endif if ( ( beta > fanAngle - 2*gamma) - && ( beta < 2*M_PI + 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)/(M_PI+2*gamma) + (beta +2*gamma - fanAngle)/(PI+2*gamma) * detval[iDet]; } else if ( beta > 2*fanAngle - && beta <= M_PI - 2*gamma) { // in region 3 + && beta <= PI - 2*gamma) { // in region 3 newdetval[newiDet] += - (beta +2*gamma - fanAngle)/(M_PI+2*gamma) + (beta +2*gamma - fanAngle)/(PI+2*gamma) * detval[iDet]; } - else if ( beta > M_PI -2*gamma - && beta <= M_PI + fanAngle ) { // in region 4 + else if ( beta > PI -2*gamma + && beta <= PI + fanAngle ) { // in region 4 newdetval[newiDet] += - (beta +2*gamma - fanAngle)/(M_PI+2*gamma) + (beta +2*gamma - fanAngle)/(PI+2*gamma) * detval[iDet]; } - else if ( beta > M_PI + fanAngle - && beta <= M_PI +2*fanAngle -2*gamma) { // in region 5 + else if ( beta > PI + fanAngle + && beta <= PI +2*fanAngle -2*gamma) { // in region 5 newdetval[newiDet] += - (2*M_PI - beta - 2*gamma + fanAngle)/(M_PI-2*gamma) + (2*PI - beta - 2*gamma + fanAngle)/(PI-2*gamma) *detval[iDet]; } - else if ( beta > M_PI +2*fanAngle -2*gamma - && beta <= 2*M_PI) { // in region 6 + else if ( beta > PI +2*fanAngle -2*gamma + && beta <= 2*PI) { // in region 6 newdetval[newiDet] += - (2*M_PI - beta - 2*gamma + fanAngle)/(M_PI-2*gamma) + (2*PI - beta - 2*gamma + fanAngle)/(PI-2*gamma) *detval[iDet]; } - else if ( beta > 2*M_PI - && beta <= 2*M_PI + fanAngle -2*gamma){ // in region 7 + else if ( beta > 2*PI + && beta <= 2*PI + fanAngle -2*gamma){ // in region 7 newdetval[newiDet] += - (2*M_PI - beta -2*gamma + fanAngle)/(M_PI-2*gamma) + (2*PI - beta -2*gamma + fanAngle)/(PI-2*gamma) *detval[iDet]; } else @@ -374,7 +376,7 @@ Projections::HalfScanFeather(void) double fanAngle = m_dFanBeamAngle; // is there enough data? - if ( m_nView != static_cast(( M_PI+fanAngle ) / dbeta) +1 ){ + 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); } @@ -392,12 +394,12 @@ Projections::HalfScanFeather(void) double beta2 = iView2 * dbeta; for ( int iDet2 = 0; iDet2 < m_nDet; iDet2++) { double gamma2 = (iDet2 -(m_nDet-1)/2)* dgamma ; - if ( ( beta2 >= M_PI - 2*gamma2) ) { // in redundant data region + if ( ( beta2 >= PI - 2*gamma2) ) { // in redundant data region int iView1, iDet1; iDet1 = (m_nDet -1) - iDet2; - //iView1 = nearest((beta2 + 2*gamma2 - M_PI)/dbeta); + //iView1 = nearest((beta2 + 2*gamma2 - PI)/dbeta); iView1 = nearest(( (iView2*dbeta) - + 2*(iDet2-(m_nDet-1)/2)*dgamma - M_PI)/dbeta); + + 2*(iDet2-(m_nDet-1)/2)*dgamma - PI)/dbeta); DetectorValue* detval2 = (m_projData[iView2])->detValues(); @@ -410,10 +412,10 @@ Projections::HalfScanFeather(void) gamma1 = -gamma2; if ( beta1 <= (fanAngle - 2*gamma1) ) x = beta1 / ( fanAngle - 2*gamma1); - else if ( (fanAngle - 2*gamma1 <= beta1 ) && beta1 <= M_PI - 2*gamma1) + else if ( (fanAngle - 2*gamma1 <= beta1 ) && beta1 <= PI - 2*gamma1) x = 1; - else if ( (M_PI - 2*gamma1 <= beta1 ) && ( beta1 <=M_PI + fanAngle) ) - x = (M_PI +fanAngle - beta1)/(fanAngle + 2*gamma1); + 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); @@ -427,7 +429,7 @@ Projections::HalfScanFeather(void) } } // heuristic scaling, why this factor? - double scalefactor = m_nView * m_rotInc / M_PI; + 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++) { @@ -650,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()) @@ -970,11 +972,11 @@ Projections::convertPolar (ImageFile& rIF, int iInterpolationID) double** ppdDet = adDet.getArray(); std::complex** ppcDetValue = new std::complex* [pProj->m_nView]; - unsigned int iView; + 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 < pProj->m_nDet; iDet++) + for (int iDet = 0; iDet < pProj->m_nDet; iDet++) ppcDetValue[iView][iDet] = std::complex(detval[iDet], 0); } @@ -1015,45 +1017,45 @@ Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad if (m_geometry == Scanner::GEOMETRY_EQUIANGULAR || m_geometry == Scanner::GEOMETRY_EQUILINEAR) pProj = interpolateToParallel(); - int iInterpDet = nx; -// int iInterpDet = pProj->m_nDet; + 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* [pProj->m_nView]; - double dInterpScale = (pProj->m_nDet-1) / static_cast(iInterpDet-1) / SQRT2; + //double dInterpScale = (pProj->m_nDet-1) / static_cast(iInterpDet-1); + double dInterpScale = pProj->m_nDet / static_cast(iInterpDet); 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++) { + // 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 (unsigned int iDet = 0; iDet < iInterpDet; iDet++) { + for (int iDet = 0; iDet < iInterpDet; iDet++) { double dInterpPos = (m_nDet / 2.) + (iDet - dMidPoint) * dInterpScale; - pcIn[iDet].re = projInterp.interpolate (dInterpPos) * dInterpScale; + pcIn[iDet].re = projInterp.interpolate (dInterpPos) * dProjScale; pcIn[iDet].im = 0; } Fourier::shuffleFourierToNaturalOrder (pcIn, iInterpDet); if (iZerosAdded > 0) { - for (unsigned int iDet1 = iMidPoint; iDet1 < iInterpDet; iDet1++) + for (int iDet1 = iInterpDet -1; iDet1 >= iMidPoint; iDet1--) pcIn[iDet1+iZerosAdded] = pcIn[iDet1]; - for (unsigned int iDet2 = iMidPoint; iDet2 < iMidPoint + iZerosAdded; iDet2++) + 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++) { + for (int iD = 0; iD < iNumInterpDetWithZeros; iD++) { ppcDetValue[iView][iD] = std::complex (pcIn[iD].re * dFFTScale, pcIn[iD].im * dFFTScale); } @@ -1089,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; @@ -1106,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 @@ -1119,15 +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; } } } @@ -1139,30 +1140,29 @@ Projections::interpolatePolar (ImageFileArray& v, ImageFileArray& vImag, { typedef std::complex complexValue; - BilinearInterpolator* pBilinear; + BilinearPolarInterpolator* pBilinear = NULL; + BicubicPolyInterpolator* pBicubic = NULL; if (iInterpolationID == POLAR_INTERP_BILINEAR) - pBilinear = new BilinearInterpolator (ppcDetValue, nView, nDetWithZeros); - - BicubicPolyInterpolator* pBicubic; - if (iInterpolationID == POLAR_INTERP_BICUBIC) + 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) { std::complex vInterp = pBilinear->interpolate (ppdView[ix][iy], ppdDet[ix][iy]); @@ -1206,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); @@ -1292,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(), false); detValues[iD] = interp.interpolate (dViewAngle, &iLastFloor); } } @@ -1337,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();