X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fprojections.cpp;h=e5c79e7364f67a15a0a1ec5ce4b08f01ed0d2df7;hp=c55aa37a5fcc94b6292ac9b912f93e858bbb5b7f;hb=728e9fcbe0b785e56e7dcd2bd305786c647050af;hpb=999a754d1519a49ca062ee87b22bf601c1ee9f21 diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp index c55aa37..e5c79e7 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.57 2001/03/11 06:34:37 kevin Exp $ +** $Id: projections.cpp,v 1.72 2001/03/29 21:25:49 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 * @@ -688,33 +689,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 +726,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(); @@ -736,63 +744,95 @@ Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad 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(); + int iInterpDet = nx; + 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* [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]; + double dInterpScale = (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 (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 = (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(); + calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet, iNumInterpDetWithZeros, dZeropadRatio, + m_detInc * dInterpScale); - if (! bError) - interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, m_nView, m_nDet, iInterpolationID); + interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, m_nView, m_nDet, iNumInterpDetWithZeros, + iInterpolationID); - for (iView = 0; iView < m_nView; iView++) - delete [] ppcDetValue[iView]; + 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,177 +842,139 @@ 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 () +Projections::interpolateToParallel () const { if (m_geometry == Scanner::GEOMETRY_PARALLEL) - return this; + return const_cast(this); int nDet = m_nDet; int nView = m_nView; @@ -991,13 +993,17 @@ Projections::interpolateToParallel () 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 (nDet % 2 == 0) // even + if (isEven (nDet)) // even pProjNew->m_detInc = m_dViewDiameter / (nDet - 1); - ParallelRaysums parallel (this); + ParallelRaysums parallel (this, ParallelRaysums::THETA_RANGE_NORMALIZE_TO_TWOPI); double* pdThetaValuesForT = new double [pProjNew->nView()]; double* pdRaysumsForT = new double [pProjNew->nView()]; @@ -1005,13 +1011,14 @@ Projections::interpolateToParallel () // 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); + 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(); - - detValues[iD] = parallel.interpolate (pdThetaValuesForT, pdRaysumsForT, pProjNew->nView(), dViewAngle); + LinearInterpolator interp (pdThetaValuesForT, pdRaysumsForT, pProjNew->nView(), false); + detValues[iD] = interp.interpolate (dViewAngle, &iLastFloor); } } delete pdThetaValuesForT; @@ -1029,12 +1036,13 @@ Projections::interpolateToParallel () detArray.setViewAngle (dViewAngle); for (int i = 0; i < pProjNew->nDet(); i++) - pdDetValueCopy[i] = detValues[i]; + pdDetValueCopy[i] = detValues[i]; double dDetPos = pProjNew->m_detStart; - for (int iD = 0; iD < pProjNew->nDet(); iD++, dDetPos += pProjNew->m_detInc) { - detValues[iD] = parallel.interpolate (pdOriginalDetPositions, pdDetValueCopy, pProjNew->nDet(), dDetPos); - } + 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; @@ -1051,8 +1059,9 @@ Projections::interpolateToParallel () // /////////////////////////////////////////////////////////////////////////////// -ParallelRaysums::ParallelRaysums (Projections* pProjections) -: m_iNumCoordinates(0), m_iNumView(pProjections->nView()), m_iNumDet(pProjections->nDet()) +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 iGeometry = pProjections->geometry(); double dDetInc = pProjections->detInc(); @@ -1060,9 +1069,10 @@ ParallelRaysums::ParallelRaysums (Projections* pProjections) double dFocalLength = pProjections->focalLength(); 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 < m_iNumView; iV++) { @@ -1074,21 +1084,24 @@ ParallelRaysums::ParallelRaysums (Projections* pProjections) 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 (pC->m_dTheta >= PI) { // convert T/Theta to 0-PI interval - pC->m_dTheta -= PI; - pC->m_dT = -pC->m_dT - pProjections->detInc(); + 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; @@ -1098,8 +1111,7 @@ ParallelRaysums::ParallelRaysums (Projections* pProjections) ParallelRaysums::~ParallelRaysums() { - for (int i = 0; i < m_iNumCoordinates; i++) - delete m_vecpCoordinates[i]; + delete m_pCoordinates; } ParallelRaysums::CoordinateContainer& @@ -1178,31 +1190,3 @@ ParallelRaysums::getDetPositions (double* pdDetPos) iPos += m_iNumView; } } - -// locate by bisection, O(log2(n)) -double -ParallelRaysums::interpolate (double* pdX, double* pdY, int n, double dX) -{ - int iLower = -1; - int iUpper = n; - - while (iUpper - iLower > 1) { - int iMiddle = (iUpper + iLower) >> 1; - if (dX >= pdX[iMiddle]) - iLower = iMiddle; - else - iUpper = iMiddle; - } - if (dX <= pdX[0]) - return pdY[0]; - else if (dX >= pdX[n-1]) - return pdY[1]; - - if (iLower < 0 || iLower >= n) { - sys_error (ERR_SEVERE, "Coordinate out of range [locateThetaBase]"); - return 0; - } - - return pdY[iLower] + (pdY[iUpper] - pdY[iLower]) * ((dX - pdX[iLower]) / (pdX[iUpper] - pdX[iLower])); -} -