** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: projections.cpp,v 1.82 2003/03/23 18:37: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
const char* const Projections::s_aszInterpName[] =
{
- {"nearest"},
- {"bilinear"},
+ "nearest",
+ "bilinear",
// {"bicubic"},
};
const char* const Projections::s_aszInterpTitle[] =
{
- {"Nearest"},
- {"Bilinear"},
+ "Nearest",
+ "Bilinear",
// {"Bicubic"},
};
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.5);
+ 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);
double dInterpScale = pProj->m_nDet / static_cast<double>(iInterpDet);
LinearInterpolator<DetectorValue> projInterp (detval, pProj->m_nDet);
for (int iDet = 0; iDet < iInterpDet; iDet++) {
double dInterpPos = (m_nDet / 2.) + (iDet - dMidPoint) * dInterpScale;
- pcIn[iDet].re = projInterp.interpolate (dInterpPos) * dProjScale;
- pcIn[iDet].im = 0;
+ pcIn[iDet][0] = projInterp.interpolate (dInterpPos) * dProjScale;
+ pcIn[iDet][1] = 0;
}
Fourier::shuffleFourierToNaturalOrder (pcIn, iInterpDet);
if (iZerosAdded > 0) {
- for (int iDet1 = iInterpDet -1; iDet1 >= iMidPoint; iDet1--)
- pcIn[iDet1+iZerosAdded] = pcIn[iDet1];
+ 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].re = pcIn[iDet2].im = 0;
+ pcIn[iDet2][0] = pcIn[iDet2][1] = 0;
}
- fftw_one (plan, pcIn, NULL);
+ fftw_execute (plan);
ppcDetValue[iView] = new std::complex<double> [iNumInterpDetWithZeros];
for (int iD = 0; iD < iNumInterpDetWithZeros; iD++) {
- ppcDetValue[iView][iD] = std::complex<double> (pcIn[iD].re * dFFTScale, pcIn[iD].im * dFFTScale);
+ 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);