** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: projections.cpp,v 1.77 2002/05/28 18:43:16 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
******************************************************************************/
#include "ct.h"
+#include <ctime>
+#include "interpolator.h"
const kuint16 Projections::m_signature = ('P'*256 + 'J');
const char* const Projections::s_aszInterpName[] =
{
- {"nearest"},
- {"bilinear"},
+ "nearest",
+ "bilinear",
// {"bicubic"},
};
const char* const Projections::s_aszInterpTitle[] =
{
- {"Nearest"},
- {"Bilinear"},
+ "Nearest",
+ "Bilinear",
// {"Bicubic"},
};
double** ppdDet = adDet.getArray();
std::complex<double>** ppcDetValue = new std::complex<double>* [pProj->m_nView];
- unsigned int iView;
+ int iView;
for (iView = 0; iView < pProj->m_nView; iView++) {
ppcDetValue[iView] = new std::complex<double> [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<double>(detval[iDet], 0);
}
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<int>(static_cast<double>(sqrt(nx*nx+ny*ny)));
int iNumInterpDetWithZeros = ProcessSignal::addZeropadFactor (iInterpDet, iZeropad);
-
+ double dProjScale = iInterpDet / (pProj->viewDiameter() * 0.05);
double dZeropadRatio = static_cast<double>(iNumInterpDetWithZeros) / static_cast<double>(iInterpDet);
- fftw_plan plan = fftw_create_plan (iNumInterpDetWithZeros, FFTW_FORWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
+ fftw_complex* pcIn = static_cast<fftw_complex*> (fftw_malloc (sizeof(fftw_complex) * iNumInterpDetWithZeros));
+ fftw_plan plan = fftw_plan_dft_1d (iNumInterpDetWithZeros, pcIn, pcIn, FFTW_FORWARD, FFTW_ESTIMATE);
- fftw_complex* pcIn = new fftw_complex [iNumInterpDetWithZeros];
std::complex<double>** ppcDetValue = new std::complex<double>* [pProj->m_nView];
- double dInterpScale = (pProj->m_nDet-1) / static_cast<double>(iInterpDet-1) / SQRT2;
+ //double dInterpScale = (pProj->m_nDet-1) / static_cast<double>(iInterpDet-1);
+ double dInterpScale = pProj->m_nDet / static_cast<double>(iInterpDet);
double dFFTScale = 1. / static_cast<double>(iInterpDet * iInterpDet);
int iMidPoint = iInterpDet / 2;
double dMidPoint = static_cast<double>(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<DetectorValue> 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].im = 0;
+ pcIn[iDet][0] = projInterp.interpolate (dInterpPos) * dProjScale;
+ pcIn[iDet][1] = 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;
+ for (int iDet1 = iInterpDet -1; iDet1 >= iMidPoint; iDet1--) {
+ pcIn[iDet1+iZerosAdded][0] = pcIn[iDet1][0];
+ pcIn[iDet1+iZerosAdded][1] = pcIn[iDet1][1];
+ }
+ for (int iDet2 = iMidPoint; iDet2 < iMidPoint + iZerosAdded; iDet2++)
+ pcIn[iDet2][0] = pcIn[iDet2][1] = 0;
}
- fftw_one (plan, pcIn, NULL);
+ fftw_execute (plan);
ppcDetValue[iView] = new std::complex<double> [iNumInterpDetWithZeros];
- for (unsigned int iD = 0; iD < iNumInterpDetWithZeros; iD++) {
- ppcDetValue[iView][iD] = std::complex<double> (pcIn[iD].re * dFFTScale, pcIn[iD].im * dFFTScale);
+ for (int iD = 0; iD < iNumInterpDetWithZeros; iD++) {
+ ppcDetValue[iView][iD] = std::complex<double> (pcIn[iD][0] * dFFTScale, pcIn[iD][1] * dFFTScale);
}
Fourier::shuffleFourierToNaturalOrder (ppcDetValue[iView], iNumInterpDetWithZeros);
}
- delete [] pcIn;
+ fftw_free(pcIn) ;
fftw_destroy_plan (plan);
int iNumDetWithZeros, double dZeropadRatio, double dDetInc)
{
double dLength = viewDiameter();
-// double dLength = phmLen();
double xMin = -dLength / 2;
double xMax = xMin + dLength;
double yMin = -dLength / 2;
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
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;
}
}
}
{
typedef std::complex<double> complexValue;
- BilinearInterpolator<complexValue>* pBilinear;
+ BilinearPolarInterpolator<complexValue>* pBilinear = NULL;
+ BicubicPolyInterpolator<complexValue>* pBicubic = NULL;
if (iInterpolationID == POLAR_INTERP_BILINEAR)
- pBilinear = new BilinearInterpolator<complexValue> (ppcDetValue, nView, nDetWithZeros);
-
- BicubicPolyInterpolator<complexValue>* pBicubic;
- if (iInterpolationID == POLAR_INTERP_BICUBIC)
+ pBilinear = new BilinearPolarInterpolator<complexValue> (ppcDetValue, nView, nDetWithZeros);
+ else if (iInterpolationID == POLAR_INTERP_BICUBIC)
pBicubic = new BicubicPolyInterpolator<complexValue> (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<int> (ppdView[ix][iy]);
unsigned int iDet = nearest<int> (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<double> vInterp = pBilinear->interpolate (ppdView[ix][iy], ppdDet[ix][iy]);
for (int iv = 0; iv < iNViews; iv++) {
unsigned char* pArgBase = pData + lDataPos;
unsigned char* p = pArgBase+0; SwapBytes4IfLittleEndian (p);
- long lProjNumber = *reinterpret_cast<long*>(p);
+ // long lProjNumber = *reinterpret_cast<long*>(p);
p = pArgBase+20; SwapBytes4IfLittleEndian (p);
long lEscale = *reinterpret_cast<long*>(p);
p = pArgBase+28; SwapBytes4IfLittleEndian (p);
- long lTime = *reinterpret_cast<long*>(p);
+ // long lTime = *reinterpret_cast<long*>(p);
p = pArgBase + 4; SwapBytes4IfLittleEndian (p);
double dAlpha = *reinterpret_cast<float*>(p) + HALFPI;
p = pArgBase+12; SwapBytes4IfLittleEndian (p);
- double dAlign = *reinterpret_cast<float*>(p);
+ // double dAlign = *reinterpret_cast<float*>(p);
p = pArgBase + 16; SwapBytes4IfLittleEndian (p);
- double dMaxValue = *reinterpret_cast<float*>(p);
+ // double dMaxValue = *reinterpret_cast<float*>(p);
DetectorArray& detArray = getDetectorArray (iv);
detArray.setViewAngle (dAlpha);
double dViewAngle = m_rotStart;
int iLastFloor = -1;
+ LinearInterpolator<double> 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<double> interp (pdThetaValuesForT, pdRaysumsForT, pProjNew->nView(), false);
detValues[iD] = interp.interpolate (dViewAngle, &iLastFloor);
}
}
///////////////////////////////////////////////////////////////////////////////
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();