From dfa390de2efc04d85b03718a6480f735516df0e8 Mon Sep 17 00:00:00 2001 From: "Kevin M. Rosenberg" Date: Fri, 7 Jul 2000 15:30:59 +0000 Subject: [PATCH] r138: *** empty log message *** --- ChangeLog | 6 ++-- configure.in | 2 +- include/backprojectors.h | 17 +++++++++-- include/filter.h | 12 ++++---- libctsim/backprojectors.cpp | 61 +++++++++++++++++++++++++++++++++++-- libctsim/filter.cpp | 41 +++++++++++++++---------- 6 files changed, 110 insertions(+), 29 deletions(-) diff --git a/ChangeLog b/ChangeLog index 6cd4d1a..d48b1ae 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,10 @@ +2.0.0-b3 - 7/09/00 + 2.0.0-b2 - 7/07/00 - Added zeropadding option to pjrec Cleaned up SignalFilter class + Added zeropad option to pjrec Added zeropad options to html and cgi files - Added fourier_table filter method + Added fourier_table and rfttw filter methods Added FFTW routines to use real/half-complex transformations 2.0.0-b1 - 7/05/00 diff --git a/configure.in b/configure.in index 121694d..776c723 100644 --- a/configure.in +++ b/configure.in @@ -4,7 +4,7 @@ dnl Must reset CDPATH so that bash's cd does not print to stdout dnl CDPATH= AC_INIT(src/pjrec.cpp) -AM_INIT_AUTOMAKE(ctsim,2.0.0-b2) +AM_INIT_AUTOMAKE(ctsim,2.0.0-b3) AM_CONFIG_HEADER(config.h) dnl Checks for programs. diff --git a/include/backprojectors.h b/include/backprojectors.h index 6b8b9c4..0a2a8e1 100644 --- a/include/backprojectors.h +++ b/include/backprojectors.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: backprojectors.h,v 1.5 2000/06/25 17:32:24 kevin Exp $ +** $Id: backprojectors.h,v 1.6 2000/07/07 15:30:59 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 @@ -39,7 +39,8 @@ class Backprojector BPROJ_TABLE, BPROJ_DIFF, BPROJ_DIFF2, - BPROJ_IDIFF2 + BPROJ_IDIFF2, + BPROJ_IDIFF3 } BackprojectID; typedef enum { @@ -59,6 +60,7 @@ class Backprojector static const char BPROJ_DIFF_STR[]= "diff"; static const char BPROJ_DIFF2_STR[]= "diff2"; static const char BPROJ_IDIFF2_STR[]= "idiff2"; + static const char BPROJ_IDIFF3_STR[]= "idiff3"; static const char INTERP_NEAREST_STR[]= "nearest"; static const char INTERP_LINEAR_STR[]= "linear"; @@ -187,3 +189,14 @@ class BackprojectIntDiff2 : public BackprojectDiff }; +class BackprojectIntDiff3 : public BackprojectDiff +{ + public: + BackprojectIntDiff3 (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType) + : BackprojectDiff::BackprojectDiff (proj, im, interpType) + {} + + void BackprojectView (const double* const t, double view_angle); +}; + + diff --git a/include/filter.h b/include/filter.h index e4c4f3f..9ee7390 100644 --- a/include/filter.h +++ b/include/filter.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: filter.h,v 1.11 2000/07/06 18:37:24 kevin Exp $ +** $Id: filter.h,v 1.12 2000/07/07 15:30:59 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 @@ -52,8 +52,10 @@ class SignalFilter { FILTER_METHOD_FOURIER, FILTER_METHOD_FOURIER_TABLE, FILTER_METHOD_FFT, +#if HAVE_FFTW FILTER_METHOD_FFTW, - FILTER_METHOD_RFFTW + FILTER_METHOD_RFFTW, +#endif } FilterMethodID; typedef enum { @@ -77,8 +79,10 @@ class SignalFilter { static const char FILTER_METHOD_FOURIER_STR[]= "fourier"; static const char FILTER_METHOD_FOURIER_TABLE_STR[]="fourier_table"; static const char FILTER_METHOD_FFT_STR[]= "fft"; +#if HAVE_FFTW static const char FILTER_METHOD_FFTW_STR[]= "fftw"; static const char FILTER_METHOD_RFFTW_STR[]= "rfftw"; +#endif static const char DOMAIN_FREQUENCY_STR[]="frequency"; static const char DOMAIN_SPATIAL_STR[]="spatial"; @@ -136,8 +140,6 @@ class SignalFilter { static double spatialResponseAnalytic (FilterID fType, double bw, double x, double param); - static void dotProduct (const double v1[], const complex v2[], complex output[], const int n); - private: double m_bw; int m_nFilterPoints; @@ -155,8 +157,6 @@ class SignalFilter { rfftw_plan m_realPlanForward, m_realPlanBackward; fftw_complex* m_vecComplexFftInput; fftw_plan m_complexPlanForward, m_complexPlanBackward; -#else - complex* m_vecFftInput; #endif bool m_fail; diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp index 1f6e92c..0684c16 100644 --- a/libctsim/backprojectors.cpp +++ b/libctsim/backprojectors.cpp @@ -8,7 +8,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: backprojectors.cpp,v 1.4 2000/06/25 17:32:24 kevin Exp $ +** $Id: backprojectors.cpp,v 1.5 2000/07/07 15:30:59 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 @@ -87,6 +87,8 @@ Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const m_pBackprojectImplem = static_cast(new BackprojectDiff2 (proj, im, m_idInterpolation)); else if (m_idBackproject == BPROJ_IDIFF2) m_pBackprojectImplem = static_cast(new BackprojectIntDiff2 (proj, im, m_idInterpolation)); + else if (m_idBackproject == BPROJ_IDIFF3) + m_pBackprojectImplem = static_cast(new BackprojectIntDiff3 (proj, im, m_idInterpolation)); else { m_fail = true; m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]"; @@ -112,6 +114,8 @@ Backprojector::convertBackprojectNameToID (const char* const backprojName) backprojID = BPROJ_DIFF2; else if (strcasecmp (backprojName, BPROJ_IDIFF2_STR) == 0) backprojID = BPROJ_IDIFF2; + else if (strcasecmp (backprojName, BPROJ_IDIFF3_STR) == 0) + backprojID = BPROJ_IDIFF3; return (backprojID); } @@ -131,6 +135,8 @@ Backprojector::convertBackprojectIDToName (const BackprojectID bprojID) bprojName = BPROJ_DIFF2_STR; else if (bprojID == BPROJ_IDIFF2) bprojName = BPROJ_IDIFF2_STR; + else if (bprojID == BPROJ_IDIFF3) + bprojName = BPROJ_IDIFF3_STR; return (bprojName); } @@ -500,7 +506,58 @@ BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const do if (iDetPos < 0 || iDetPos >= nDet - 1) errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); else - *pImCol++ += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); + *pImCol++ += ((1.-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); + } + } // end for y + } // end for x +} + +// CLASS IDENTICATION +// BackprojectIntDiff3 +// +// PURPOSE +// Highly optimized version of BackprojectIntDiff2 + +void +BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const double view_angle) +{ + double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle + static const int scaleShift = 16; + static const kint32 scale = (1 << scaleShift); + static const double dScale = scale; + static const kint32 halfScale = scale / 2; + + const kint32 det_dx = nearest (xInc * sin (theta) / detInc * scale); + const kint32 det_dy = nearest (yInc * cos (theta) / detInc * scale); + + // calculate L for first point in image (0, 0) + kint32 detPosColStart = nearest ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale); + + for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) { + kint32 curDetPos = detPosColStart; + ImageFileColumn pImCol = v[ix]; + + for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { + if (interpType == Backprojector::INTERP_NEAREST) { + int iDetPos = (curDetPos + halfScale) >> scaleShift; + + if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos + errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); + else + *pImCol++ += filteredProj[iDetPos]; + } else if (interpType == Backprojector::INTERP_LINEAR) { + kint32 detPosFloor = curDetPos / scale; + kint32 detPosRemainder = curDetPos % scale; + if (detPosRemainder < 0) { + detPosFloor--; + detPosRemainder += scale; + } + int iDetPos = iDetCenter + detPosFloor; + double frac = detPosRemainder / dScale; + if (iDetPos < 0 || iDetPos >= nDet - 1) + errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); + else + *pImCol++ += ((1.-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); } } // end for y } // end for x diff --git a/libctsim/filter.cpp b/libctsim/filter.cpp index 6a344d4..65163f3 100644 --- a/libctsim/filter.cpp +++ b/libctsim/filter.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: filter.cpp,v 1.14 2000/07/06 18:37:24 kevin Exp $ +** $Id: filter.cpp,v 1.15 2000/07/07 15:30:59 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 @@ -129,10 +129,21 @@ SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID m_vecFourierSinTable = NULL; m_vecFilter = NULL; - if (m_idFilterMethod == FILTER_METHOD_FFT) - m_idFilterMethod = FILTER_METHOD_FFTW; + if (m_idFilterMethod == FILTER_METHOD_FFT) { +#if HAVE_FFTW + m_idFilterMethod = FILTER_METHOD_RFFTW; +#else + m_fail = true; + m_failMessage = "FFT not yet implemented"; + return; +#endif + } - if (m_idFilterMethod == FILTER_METHOD_FOURIER || FILTER_METHOD_FOURIER_TABLE || m_idFilterMethod == FILTER_METHOD_FFT || m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW) { + if (m_idFilterMethod == FILTER_METHOD_FOURIER || FILTER_METHOD_FOURIER_TABLE || m_idFilterMethod == FILTER_METHOD_FFT +#if HAVE_FFTW + || m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW +#endif + ) { m_nFilterPoints = m_nSignalPoints; if (m_zeropad > 0) { double logBase2 = log(m_nSignalPoints) / log(2); @@ -321,10 +332,12 @@ SignalFilter::convertFilterMethodNameToID (const char* const filterMethodName) fmID = FILTER_METHOD_FOURIER_TABLE; else if (strcasecmp (filterMethodName, FILTER_METHOD_FFT_STR) == 0) fmID = FILTER_METHOD_FFT; +#if HAVE_FFTW else if (strcasecmp (filterMethodName, FILTER_METHOD_FFTW_STR) == 0) fmID = FILTER_METHOD_FFTW; else if (strcasecmp (filterMethodName, FILTER_METHOD_RFFTW_STR) == 0) fmID = FILTER_METHOD_RFFTW; +#endif return (fmID); } @@ -342,10 +355,12 @@ SignalFilter::convertFilterMethodIDToName (const FilterMethodID fmID) return (FILTER_METHOD_FOURIER_TABLE_STR); else if (fmID == FILTER_METHOD_FFT) return (FILTER_METHOD_FFT_STR); +#if HAVE_FFTW else if (fmID == FILTER_METHOD_FFTW) return (FILTER_METHOD_FFTW_STR); else if (fmID == FILTER_METHOD_RFFTW) return (FILTER_METHOD_RFFTW_STR); +#endif return (name); } @@ -391,10 +406,10 @@ SignalFilter::filterSignal (const float input[], double output[]) const inputSignal[i] = 0; // zeropad complex fftSignal[m_nFilterPoints]; finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, -1); - complex filteredSignal[m_nFilterPoints]; - dotProduct (m_vecFilter, fftSignal, filteredSignal, m_nFilterPoints); + for (int i = 0; i < m_nFilterPoints; i++) + fftSignal[i] *= m_vecFilter[i]; double inverseFourier[m_nFilterPoints]; - finiteFourierTransform (filteredSignal, inverseFourier, m_nFilterPoints, 1); + finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, 1); for (int i = 0; i < m_nSignalPoints; i++) output[i] = inverseFourier[i]; } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) { @@ -405,10 +420,10 @@ SignalFilter::filterSignal (const float input[], double output[]) const inputSignal[i] = 0; // zeropad complex fftSignal[m_nFilterPoints]; finiteFourierTransform (inputSignal, fftSignal, -1); - complex filteredSignal[m_nFilterPoints]; - dotProduct (m_vecFilter, fftSignal, filteredSignal, m_nFilterPoints); + for (int i = 0; i < m_nFilterPoints; i++) + fftSignal[i] *= m_vecFilter[i]; double inverseFourier[m_nFilterPoints]; - finiteFourierTransform (filteredSignal, inverseFourier, 1); + finiteFourierTransform (fftSignal, inverseFourier, 1); for (int i = 0; i < m_nSignalPoints; i++) output[i] = inverseFourier[i]; } @@ -930,9 +945,3 @@ SignalFilter::finiteFourierTransform (const complex input[], double outp } -void -SignalFilter::dotProduct (const double v1[], const complex v2[], complex output[], const int n) -{ - for (int i = 0; i < n; i++) - output[i] = v1[i] * v2[i]; -} -- 2.34.1