From: Kevin M. Rosenberg Date: Sat, 19 Aug 2000 23:00:05 +0000 (+0000) Subject: r180: *** empty log message *** X-Git-Tag: debian-4.5.3-3~837 X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=commitdiff_plain;h=e36dfad3f0818b4c3457fbe7277faa6f4ca28dfe r180: *** empty log message *** --- diff --git a/ChangeLog b/ChangeLog index c1641d6..72fcec0 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,10 +1,10 @@ -2.0.0-b9 - 8/5/00 +2.0.0-b9 - 8/15/00 Added RCS Id strings to executable files Added RPM Spec file for RPM package creation Added loading of ASCII phanthom definitions from files - Added frequency_filter option to reconstruction + Added Filter-Generation option to reconstruction Fixed compilation for non-SGP architectures - Decomposed SignalFilter class into SignalProcess and CreateFilter classes + Decomposed SignalFilter class into ProcessSignal and SignalFilter classes 2.0.0-b8 - 8/1/00 Added line color support to SGP Fixed lineAbs bug diff --git a/configure b/configure index 6ce145f..79e9cb6 100755 --- a/configure +++ b/configure @@ -4224,7 +4224,7 @@ ctlibs="$ctlibs_base -lctsim $ctlibs_graphics -lctsupport $ctlibs_tools" if test -n "$lamdir" ; then lamprograms="pjrec-lam phm2if-lam phm2pj-lam" - lamdefs="$CLFAGS" + lamdefs="$CFLAGS" fi diff --git a/include/Makefile.am b/include/Makefile.am index 28950d3..a7f4ddf 100644 --- a/include/Makefile.am +++ b/include/Makefile.am @@ -1,4 +1,5 @@ -noinst_HEADERS=ct.h ezplot.h pol.h sgp.h array2d.h imagefile.h backprojectors.h mpiworld.h fnetorderstream.h phantom.h timer.h sstream scanner.h projections.h ctsupport.h filter.h array2dfile.h trace.h transformmatrix.h +noinst_HEADERS=ct.h ezplot.h pol.h sgp.h array2d.h imagefile.h backprojectors.h mpiworld.h fnetorderstream.h phantom.h timer.h sstream scanner.h projections.h ctsupport.h filter.h array2dfile.h trace.h transformmatrix.h procsignal.h + diff --git a/include/ct.h b/include/ct.h index 2cc4b55..a0383af 100644 --- a/include/ct.h +++ b/include/ct.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: ct.h,v 1.30 2000/08/03 09:57:33 kevin Exp $ +** $Id: ct.h,v 1.31 2000/08/19 22:59:06 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 @@ -155,6 +155,7 @@ using namespace std; #include "scanner.h" #include "backprojectors.h" #include "filter.h" +#include "procsignal.h" #include "projections.h" #include "trace.h" diff --git a/include/filter.h b/include/filter.h index ad2fc26..492f93d 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.19 2000/08/09 22:52:52 kevin Exp $ +** $Id: filter.h,v 1.20 2000/08/19 22:59:06 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 @@ -40,9 +40,17 @@ #include +// CLASS IDENTIFICATION +// SignalFilter A filter used to process signals +// +// CONTAINS +// signal vector +// +// Can create either a time/spatial waveform or a frequency signal +// Waveforms can be created either by direct calculation or by inverse fourier transform + class SignalFilter { public: - static const int FILTER_INVALID; static const int FILTER_ABS_BANDLIMIT; // filter times |x| static const int FILTER_ABS_SINC; @@ -55,54 +63,20 @@ class SignalFilter { static const int FILTER_COSINE; static const int FILTER_TRIANGLE; - static const int FILTER_METHOD_INVALID; - static const int FILTER_METHOD_CONVOLUTION; - static const int FILTER_METHOD_FOURIER; - static const int FILTER_METHOD_FOURIER_TABLE; - static const int FILTER_METHOD_FFT; -#if HAVE_FFTW - static const int FILTER_METHOD_FFTW; - static const int FILTER_METHOD_RFFTW; -#endif - static const int DOMAIN_INVALID; static const int DOMAIN_FREQUENCY; static const int DOMAIN_SPATIAL; - static const int FREQUENCY_FILTER_INVALID; - static const int FREQUENCY_FILTER_DIRECT_FREQUENCY; - static const int FREQUENCY_FILTER_INVERSE_SPATIAL; + SignalFilter (const char* szFilterName, double dFilterMinimum, double dFilterMaximum, int nFilterPoints, double dBandwidth, double dFilterParam, const char* szDomainName); - SignalFilter (const char* filterName, const char* filterMethodName,double bw, double signalIncrement, int n, double param, const char* domainName, const char* frequencyFilterName, const int zeropad = 0, const int preinterpolationFactor = 1); + SignalFilter (const int idFilter, double dFilterMinimum, double dFilterMaximum, int nFilterPoints, double dBandwidth, double dFilterParam, const int idDomain); - SignalFilter (const int filt_type, int filterMethodID, double bw, double signalIncrement, int n, double param, const int domain, int filterFilterID, const int zeropad = 0, const int preinterpolationFactor = 1); - - SignalFilter (const char* filterName, const char* domainName, double bw, double param); + SignalFilter (const char* szFilterName, const char* szDomainName, double dBandwidth, double dFilterParam); ~SignalFilter (void); double* getFilter (void) const - { return m_vecFilter; } - - int getNFilterPoints (void) const - { return m_nFilterPoints; } - - double convolve (const double f[], const double dx, const int n, const int np) const; - - double convolve (const float f[], const double dx, const int n, const int np) const; - - void filterSignal (const double input[], double output[]) const; - void filterSignal (const float input[], double output[]) const; - - static void finiteFourierTransform (const double input[], complex output[], const int n, const int direction); - static void finiteFourierTransform (const complex input[], complex output[], const int n, const int direction); - static void finiteFourierTransform (const complex input[], double output[], const int n, const int direction); - - void finiteFourierTransform (const double input[], complex output[], const int direction) const; - void finiteFourierTransform (const complex input[], complex output[], const int direction) const; - void finiteFourierTransform (const complex input[], double output[], const int direction) const; - - void setTraceLevel (int traceLevel) {m_traceLevel = traceLevel; } + { return m_adFilter; } bool fail(void) const {return m_fail;} const string& failMessage(void) const {return m_failMessage;} @@ -111,10 +85,12 @@ class SignalFilter { const string& nameDomain(void) const { return m_nameDomain;} const int idFilter(void) const { return m_idFilter;} const int idDomain(void) const { return m_idDomain;} - const int idFrequencyFilter() const { return m_idFrequencyFilter;} - const double getFilterMin(void) const {return m_filterMin;} - const double getFilterMax(void) const {return m_filterMax;} - const double getFilterIncrement(void) const {return m_filterInc;} + + int getNFilterPoints (void) const { return m_nFilterPoints; } + const double getFilterMin(void) const {return m_dFilterMin;} + const double getFilterMax(void) const {return m_dFilterMax;} + const double getFilterIncrement(void) const {return m_dFilterInc;} + void copyFilterData(double *pdFilter, const int iStart, const int nPoints) const; double response (double x); @@ -135,13 +111,6 @@ class SignalFilter { static const char* convertFilterIDToName (const int idFilter); static const char* convertFilterIDToTitle (const int idFilter); - static const int getFilterMethodCount() {return s_iFilterMethodCount;} - static const char** getFilterMethodNameArray() {return s_aszFilterMethodName;} - static const char** getFilterMethodTitleArray() {return s_aszFilterMethodTitle;} - static int convertFilterMethodNameToID (const char* const filterMethodName); - static const char* convertFilterMethodIDToName (const int idFilterMethod); - static const char* convertFilterMethodIDToTitle (const int idFilterMethod); - static const int getDomainCount() {return s_iDomainCount;} static const char** getDomainNameArray() {return s_aszDomainName;} static const char** getDomainTitleArray() {return s_aszDomainTitle;} @@ -149,78 +118,46 @@ class SignalFilter { static const char* convertDomainIDToName (const int idDomain); static const char* convertDomainIDToTitle (const int idDomain); - static const int getFrequencyFilterCount() {return s_iFrequencyFilterCount;} - static const char** getFrequencyFilterNameArray() {return s_aszFrequencyFilterName;} - static const char** getFrequencyFilterTitleArray() {return s_aszFrequencyFilterTitle;} - static int convertFrequencyFilterNameToID (const char* const ffName); - static const char* convertFrequencyFilterIDToName (const int idFF); - static const char* convertFrequencyFilterIDToTitle (const int idFF); - private: - double m_bw; int m_nFilterPoints; - int m_nSignalPoints; - double m_signalInc; - double m_filterMin; - double m_filterMax; - double m_filterInc; - double* m_vecFilter; - double* m_vecFourierCosTable; - double* m_vecFourierSinTable; - complex* m_complexVecFilter; -#ifdef HAVE_FFTW - fftw_real* m_vecRealFftInput, *m_vecRealFftSignal; - rfftw_plan m_realPlanForward, m_realPlanBackward; - fftw_complex* m_vecComplexFftInput, *m_vecComplexFftSignal; - fftw_plan m_complexPlanForward, m_complexPlanBackward; -#endif + double m_dBandwidth; + double m_dFilterParam; + double m_dFilterInc; + double m_dFilterMin; + double m_dFilterMax; + double* m_adFilter; - bool m_fail; - string m_failMessage; string m_nameFilter; - string m_nameFilterMethod; string m_nameDomain; - string m_nameFrequencyFilter; int m_idFilter; - int m_idFilterMethod; - int m_idFrequencyFilter; int m_idDomain; - double m_filterParam; - int m_traceLevel; - int m_zeropad; - int m_nOutputPoints; - int m_preinterpolationFactor; + + bool m_fail; + string m_failMessage; static const char* s_aszFilterName[]; static const char* s_aszFilterTitle[]; static const int s_iFilterCount; - static const char* s_aszFilterMethodName[]; - static const char* s_aszFilterMethodTitle[]; - static const int s_iFilterMethodCount; static const char* s_aszDomainName[]; static const char* s_aszDomainTitle[]; static const int s_iDomainCount; - static const char* s_aszFrequencyFilterName[]; - static const char* s_aszFrequencyFilterTitle[]; - static const int s_iFrequencyFilterCount; - static int N_INTEGRAL; static const bool haveAnalyticSpatial (const int filterID); - void init (const int filt_type, const int filterMethod, double bw, double signalIncrement, int n, double param, const int domain, const int frequencyFilter, const int zeropad, const int preInterpScale); + void init (const int idFilter, double dFilterMin, double dFilterMax, int nFilterPoints, double dBandwidth, double dFilterParam, const int idDomain); - double spatialResponseCalc (double x, double param) const; + void createFrequencyFilter (double* x) const; + void createSpatialFilter (double* x) const; - double spatialResponseAnalytic (double x, double param) const; - - double frequencyResponse (double u, double param) const; + double spatialResponseCalc (double x) const; + double spatialResponseAnalytic (double x) const; + double frequencyResponse (double u) const; static double sinc (double x, double mult) { return (fabs(x) > F_EPSILON ? (sin (x * mult) / x) : 1.0); } - - static double integral_abscos (double u, double w); - + static double integral_abscos (double u, double w) + { return (fabs (u) > F_EPSILON ? (cos (u * w) - 1) / (u * u) + w / u * sin (u * w) : (w * w / 2)); } }; #endif diff --git a/include/procsignal.h b/include/procsignal.h new file mode 100644 index 0000000..1b08aea --- /dev/null +++ b/include/procsignal.h @@ -0,0 +1,167 @@ +/***************************************************************************** +** FILE IDENTIFICATION +** +** Name: filter.h +** Purpose: Signal filter header file +** Programmer: Kevin Rosenberg +** Date Started: June 2000 +** +** This is part of the CTSim program +** Copyright (C) 1983-2000 Kevin Rosenberg +** +** $Id: procsignal.h,v 1.1 2000/08/19 23:00:05 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 +** published by the Free Software Foundation. +** +** This program is distributed in the hope that it will be useful, +** but WITHOUT ANY WARRANTY; without even the implied warranty of +** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +** GNU General Public License for more details. +** +** You should have received a copy of the GNU General Public License +** along with this program; if not, write to the Free Software +** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +******************************************************************************/ + +#ifndef PROCSIGNAL_H +#define PROCSIGNAL_H + + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif +#ifdef HAVE_FFTW +#include +#include +#endif + +#include + + +class SignalFilter; + +class ProcessSignal { + public: + static const int FILTER_METHOD_INVALID; + static const int FILTER_METHOD_CONVOLUTION; + static const int FILTER_METHOD_FOURIER; + static const int FILTER_METHOD_FOURIER_TABLE; + static const int FILTER_METHOD_FFT; +#if HAVE_FFTW + static const int FILTER_METHOD_FFTW; + static const int FILTER_METHOD_RFFTW; +#endif + + static const int FILTER_GENERATION_INVALID; + static const int FILTER_GENERATION_DIRECT; + static const int FILTER_GENERATION_INVERSE_FOURIER; + + ProcessSignal (const char* szFilterName, const char* szFilterMethodName,double bw, double signalIncrement, int n, double param, const char* szDomainName, const char* szFilterGenerationName, const int zeropad = 0, const int preinterpolationFactor = 1); + + ProcessSignal (const int idFilter, int idFilterMethod, double dBandwidth, double dSignalIncrement, int nSignalPoints, double dFilterParam, const int idDomain, int idFilterGeneration, const int iZeropad = 0, const int iPreinterpolationFactor = 1); + + ~ProcessSignal(); + + void filterSignal (const double input[], double output[]) const; + void filterSignal (const float input[], double output[]) const; + + bool fail(void) const {return m_fail;} + const string& failMessage(void) const {return m_failMessage;} + + void setTraceLevel (int traceLevel) {m_traceLevel = traceLevel; } + + int getNFilterPoints (void) const { return m_nFilterPoints; } + const double getFilterMin(void) const {return m_dFilterMin;} + const double getFilterMax(void) const {return m_dFilterMax;} + const double getFilterIncrement(void) const {return m_dFilterInc;} + double* getFilter(void) {return m_adFilter;} + const double* getFilter(void) const {return m_adFilter;} + + const int idFilterGeneration() const { return m_idFilterGeneration;} + + static const int getFilterGenerationCount() {return s_iFilterGenerationCount;} + static const char** getFilterGenerationNameArray() {return s_aszFilterGenerationName;} + static const char** getFilterGenerationTitleArray() {return s_aszFilterGenerationTitle;} + static int convertFilterGenerationNameToID (const char* const fgName); + static const char* convertFilterGenerationIDToName (const int idFG); + static const char* convertFilterGenerationIDToTitle (const int idFG); + + static const int getFilterMethodCount() {return s_iFilterMethodCount;} + static const char** getFilterMethodNameArray() {return s_aszFilterMethodName;} + static const char** getFilterMethodTitleArray() {return s_aszFilterMethodTitle;} + static int convertFilterMethodNameToID (const char* const filterMethodName); + static const char* convertFilterMethodIDToName (const int idFilterMethod); + static const char* convertFilterMethodIDToTitle (const int idFilterMethod); + + // transforms using direct trigometric calculation + static void finiteFourierTransform (const double input[], double output[], const int n, const int direction); + static void finiteFourierTransform (const double input[], complex output[], const int n, const int direction); + static void finiteFourierTransform (const complex input[], complex output[], const int n, const int direction); + static void finiteFourierTransform (const complex input[], double output[], const int n, const int direction); + + + static void shuffleNaturalToFourierOrder (double* pdVector, const int n); + + static void shuffleFourierToNaturalOrder (double* pdVector, const int n); + + private: + string m_nameFilterMethod; + string m_nameFilterGeneration; + int m_idFilterMethod; + int m_idFilterGeneration; + int m_nSignalPoints; + double* m_adFourierCosTable; + double* m_adFourierSinTable; + int m_nFilterPoints; + double m_dSignalInc; + double m_dFilterInc; + double m_dFilterMin; + double m_dFilterMax; + double* m_adFilter; + bool m_bFrequencyFiltering; + + // Variables also kept in SignalFilter class + int m_idFilter; + int m_idDomain; + + int m_traceLevel; + double m_dBandwidth; + double m_dFilterParam; + int m_iZeropad; + int m_nOutputPoints; + int m_iPreinterpolationFactor; + + bool m_fail; + string m_failMessage; + + static const char* s_aszFilterMethodName[]; + static const char* s_aszFilterMethodTitle[]; + static const int s_iFilterMethodCount; + static const char* s_aszFilterGenerationName[]; + static const char* s_aszFilterGenerationTitle[]; + static const int s_iFilterGenerationCount; + +#ifdef HAVE_FFTW + fftw_real* m_adRealFftInput, *m_adRealFftSignal; + rfftw_plan m_realPlanForward, m_realPlanBackward; + fftw_complex* m_adComplexFftInput, *m_adComplexFftSignal; + fftw_plan m_complexPlanForward, m_complexPlanBackward; +#endif + + void init (const int idFilter, int idFilterMethod, double dBandwidth, double dSignalIncrement, int nSignalPoints, double dFilterParam, const int idDomain, int idFilterGeneration, const int iZeropad, const int iPreinterpolationFactor); + + // transforms that use precalculated trig tables, therefore don't + // require number of data points (n) as an argument + void finiteFourierTransform (const double input[], complex output[], const int direction) const; + void finiteFourierTransform (const complex input[], complex output[], const int direction) const; + void finiteFourierTransform (const complex input[], double output[], const int direction) const; + + double convolve (const double f[], const double dx, const int n, const int np) const; + double convolve (const float f[], const double dx, const int n, const int np) const; + +}; + + +#endif diff --git a/libctgraphics/ezplot.cpp b/libctgraphics/ezplot.cpp index 4579668..348dd3f 100644 --- a/libctgraphics/ezplot.cpp +++ b/libctgraphics/ezplot.cpp @@ -6,7 +6,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: ezplot.cpp,v 1.9 2000/07/29 19:50:08 kevin Exp $ +** $Id: ezplot.cpp,v 1.10 2000/08/19 22:59:06 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 @@ -253,8 +253,8 @@ EZPlot::plot () } /* find nice endpoints for axes */ - axis_scale (xmin, xmax, o_xmajortick - 1, &xgw_min, &xgw_max, &x_nint); - axis_scale (ymin, ymax, o_ymajortick - 1, &ygw_min, &ygw_max, &y_nint); + if (! axis_scale (xmin, xmax, o_xmajortick - 1, &xgw_min, &xgw_max, &x_nint) || ! axis_scale (ymin, ymax, o_ymajortick - 1, &ygw_min, &ygw_max, &y_nint)) + return; /* check if user set x-axis extents */ if (s_xmin == TRUE) { diff --git a/libctsim/Makefile.am b/libctsim/Makefile.am index 96037cc..3138c01 100644 --- a/libctsim/Makefile.am +++ b/libctsim/Makefile.am @@ -1,5 +1,5 @@ noinst_LIBRARIES = libctsim.a -libctsim_a_SOURCES = filter.cpp scanner.cpp projections.cpp phantom.cpp imagefile.cpp backprojectors.cpp array2dfile.cpp trace.cpp +libctsim_a_SOURCES = filter.cpp scanner.cpp projections.cpp phantom.cpp imagefile.cpp backprojectors.cpp array2dfile.cpp trace.cpp procsignal.cpp INCLUDES=@my_includes@ EXTRA_DIST=Makefile.nt diff --git a/libctsim/filter.cpp b/libctsim/filter.cpp index d7d8b97..a88fad5 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.24 2000/08/03 09:57:33 kevin Exp $ +** $Id: filter.cpp,v 1.25 2000/08/19 22:59:06 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 @@ -69,40 +69,6 @@ const char* SignalFilter::s_aszFilterTitle[] = { const int SignalFilter::s_iFilterCount = sizeof(s_aszFilterName) / sizeof(const char*); -const int SignalFilter::FILTER_METHOD_INVALID = -1; -const int SignalFilter::FILTER_METHOD_CONVOLUTION = 0; -const int SignalFilter::FILTER_METHOD_FOURIER = 1; -const int SignalFilter::FILTER_METHOD_FOURIER_TABLE = 2; -const int SignalFilter::FILTER_METHOD_FFT = 3; -#if HAVE_FFTW -const int SignalFilter::FILTER_METHOD_FFTW = 4; -const int SignalFilter::FILTER_METHOD_RFFTW =5 ; -#endif - -const char* SignalFilter::s_aszFilterMethodName[] = { - {"convolution"}, - {"fourier"}, - {"fouier_table"}, - {"fft"}, -#if HAVE_FFTW - {"fftw"}, - {"rfftw"}, -#endif -}; - -const char* SignalFilter::s_aszFilterMethodTitle[] = { - {"Convolution"}, - {"Direct Fourier"}, - {"Fouier Trigometric Table Lookout"}, - {"FFT"}, -#if HAVE_FFTW - {"FFTW"}, - {"Real/Half-Complex FFTW"}, -#endif -}; - -const int SignalFilter::s_iFilterMethodCount = sizeof(s_aszFilterMethodName) / sizeof(const char*); - const int SignalFilter::DOMAIN_INVALID = -1; const int SignalFilter::DOMAIN_FREQUENCY = 0; @@ -121,23 +87,6 @@ const char* SignalFilter::s_aszDomainTitle[] = { const int SignalFilter::s_iDomainCount = sizeof(s_aszDomainName) / sizeof(const char*); -const int SignalFilter::FREQUENCY_FILTER_INVALID = -1; -const int SignalFilter::FREQUENCY_FILTER_DIRECT_FREQUENCY = 0; -const int SignalFilter::FREQUENCY_FILTER_INVERSE_SPATIAL = 1; - -const char* SignalFilter::s_aszFrequencyFilterName[] = { - {"direct_frequency"}, - {"inverse_spatial"}, -}; - -const char* SignalFilter::s_aszFrequencyFilterTitle[] = { - {"Direct Frequency"}, - {"Inverse Spatial"}, -}; - -const int SignalFilter::s_iFrequencyFilterCount = sizeof(s_aszFrequencyFilterName) / sizeof(const char*); - - /* NAME * SignalFilter::SignalFilter Construct a signal * @@ -147,243 +96,132 @@ const int SignalFilter::s_iFrequencyFilterCount = sizeof(s_aszFrequencyFilterNam * int filt_type Type of filter wanted * double bw Bandwidth of filter * double filterMin, filterMax Filter limits - * int nSignalPoints Number of points in signal + * int nFilterPoints Number of points in signal * double param General input parameter to filters * int domain FREQUENCY or SPATIAL domain wanted */ -SignalFilter::SignalFilter (const char* filterName, const char* filterMethodName, double bw, double signalIncrement, int nSignalPoints, double param, const char* domainName, const char* frequencyFilterName, int zeropad = 0, int preinterpolationFactor = 1) - : m_vecFilter(NULL), m_vecFourierCosTable(NULL), m_vecFourierSinTable(NULL), m_fail(false) +SignalFilter::SignalFilter (const char* szFilterName, double dFilterMinimum, double dFilterMaximum, int nFilterPoints, double dBandwidth, double dFilterParam, const char* szDomainName) + : m_adFilter(NULL), m_fail(false) { - m_idFilter = convertFilterNameToID (filterName); + m_idFilter = convertFilterNameToID (szFilterName); if (m_idFilter == FILTER_INVALID) { m_fail = true; m_failMessage = "Invalid Filter name "; - m_failMessage += filterName; - return; - } - m_idFilterMethod = convertFilterMethodNameToID (filterMethodName); - if (m_idFilterMethod == FILTER_METHOD_INVALID) { - m_fail = true; - m_failMessage = "Invalid filter method name "; - m_failMessage += filterMethodName; + m_failMessage += szFilterName; return; } - m_idDomain = convertDomainNameToID (domainName); + m_idDomain = convertDomainNameToID (szDomainName); if (m_idDomain == DOMAIN_INVALID) { m_fail = true; m_failMessage = "Invalid domain name "; - m_failMessage += domainName; + m_failMessage += szDomainName; return; } - m_idFrequencyFilter = convertFrequencyFilterNameToID (frequencyFilterName); - if (m_idFrequencyFilter == FREQUENCY_FILTER_INVALID) { - m_fail = true; - m_failMessage = "Invalid frequency filter name "; - m_failMessage += frequencyFilterName; - return; - } - init (m_idFilter, m_idFilterMethod, bw, signalIncrement, nSignalPoints, param, m_idDomain, m_idFrequencyFilter, zeropad, preinterpolationFactor); + init (m_idFilter, dFilterMinimum, dFilterMaximum, nFilterPoints, dBandwidth, dFilterParam, m_idDomain); } -SignalFilter::SignalFilter (const int filterID, const int filterMethodID, double bw, double signalIncrement, int nSignalPoints, double param, const int domainID, int frequencyFilterID, int zeropad = 0, int preinterpolationFactor = 1) - : m_vecFilter(NULL), m_vecFourierCosTable(NULL), m_vecFourierSinTable(NULL), m_fail(false) +SignalFilter::SignalFilter (const int idFilter, double dFilterMinimum, double dFilterMaximum, int nFilterPoints, double dBandwidth, double dFilterParam, const int idDomain) + : m_adFilter(NULL), m_fail(false) { - init (filterID, filterMethodID, bw, signalIncrement, nSignalPoints, param, domainID, frequencyFilterID, zeropad, preinterpolationFactor); + init (idFilter, dFilterMinimum, dFilterMaximum, nFilterPoints, dBandwidth, dFilterParam, idDomain); } -SignalFilter::SignalFilter (const char* filterName, const char* domainName, double bw, double param) - : m_vecFilter(NULL), m_vecFourierCosTable(NULL), m_vecFourierSinTable(NULL), m_fail(false) +SignalFilter::SignalFilter (const char* szFilterName, const char* szDomainName, double dBandwidth, double dFilterParam) + : m_adFilter(NULL), m_fail(false) { - m_bw = bw; - m_nSignalPoints = 0; m_nFilterPoints = 0; - m_filterParam = param; - m_idFilter = convertFilterNameToID (filterName); + m_dBandwidth = dBandwidth; + m_dFilterParam = dFilterParam; + m_idFilter = convertFilterNameToID (szFilterName); if (m_idFilter == FILTER_INVALID) { m_fail = true; m_failMessage = "Invalid Filter name "; - m_failMessage += filterName; + m_failMessage += szFilterName; return; } - m_idDomain = convertDomainNameToID (domainName); + m_idDomain = convertDomainNameToID (szDomainName); if (m_idDomain == DOMAIN_INVALID) { m_fail = true; m_failMessage = "Invalid domain name "; - m_failMessage += domainName; + m_failMessage += szDomainName; return; } } void -SignalFilter::init (const int filterID, const int filterMethodID, double bw, double signalIncrement, int nSignalPoints, double filterParam, const int domainID, const int frequencyFilterID, int zeropad, int preinterpolationFactor) +SignalFilter::init (const int idFilter, double dFilterMinimum, double dFilterMaximum, int nFilterPoints, double dBandwidth, double dFilterParam, const int idDomain) { - m_bw = bw; - m_idFilter = filterID; - m_idDomain = domainID; - m_idFilterMethod = filterMethodID; - m_idFrequencyFilter = frequencyFilterID; - if (m_idFilter == FILTER_INVALID || m_idDomain == DOMAIN_INVALID || m_idFilterMethod == FILTER_METHOD_INVALID || m_idFrequencyFilter == FREQUENCY_FILTER_INVALID) { + m_idFilter = idFilter; + m_idDomain = idDomain; + if (m_idFilter == FILTER_INVALID || m_idDomain == DOMAIN_INVALID) { m_fail = true; return; } - m_traceLevel = TRACE_NONE; - m_nameFilter = convertFilterIDToName (m_idFilter); - m_nameDomain = convertDomainIDToName (m_idDomain); - m_nameFilterMethod = convertFilterMethodIDToName (m_idFilterMethod); - m_nameFrequencyFilter = convertFrequencyFilterIDToName (m_idFrequencyFilter); - m_nSignalPoints = nSignalPoints; - m_signalInc = signalIncrement; - m_filterParam = filterParam; - m_zeropad = zeropad; - m_preinterpolationFactor = preinterpolationFactor; - - m_vecFourierCosTable = NULL; - m_vecFourierSinTable = NULL; - m_vecFilter = NULL; - - if (m_idFilterMethod == FILTER_METHOD_FFT) { -#if HAVE_FFTW - m_idFilterMethod = FILTER_METHOD_RFFTW; -#else + if (nFilterPoints < 2) { m_fail = true; - m_failMessage = "FFT not yet implemented"; + m_failMessage = "Number of filter points "; + m_failMessage += nFilterPoints; + m_failMessage = " less than 2"; return; -#endif - } - - if (m_idFilterMethod == FILTER_METHOD_FOURIER || m_idFilterMethod == 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); - int nextPowerOf2 = static_cast(floor(logBase2)); - if (logBase2 != floor(logBase2)) - nextPowerOf2++; - nextPowerOf2 += (m_zeropad - 1); - m_nFilterPoints = 1 << nextPowerOf2; - if (m_traceLevel >= TRACE_TEXT) - cout << "nFilterPoints = " << m_nFilterPoints << endl; - } - m_nOutputPoints = m_nFilterPoints * m_preinterpolationFactor; - m_filterMin = -1. / (2 * m_signalInc); - m_filterMax = 1. / (2 * m_signalInc); - m_filterInc = (m_filterMax - m_filterMin) / m_nFilterPoints; - m_vecFilter = new double [m_nFilterPoints]; - int halfFilter = m_nFilterPoints / 2; - for (int i = 0; i <= halfFilter; i++) - m_vecFilter[i] = static_cast(i) / halfFilter/ (2. * m_signalInc); - for (int i = 1; i <= halfFilter; i++) - m_vecFilter[m_nFilterPoints - i] = static_cast(i) / halfFilter / (2. * m_signalInc); } - // precalculate sin and cosine tables for fourier transform - if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) { - int nFourier = max(m_nFilterPoints,m_nOutputPoints) * max(m_nFilterPoints, m_nOutputPoints) + 1; - double angleIncrement = (2. * PI) / m_nFilterPoints; - m_vecFourierCosTable = new double[ nFourier ]; - m_vecFourierSinTable = new double[ nFourier ]; - double angle = 0; - for (int i = 0; i < nFourier; i++) { - m_vecFourierCosTable[i] = cos (angle); - m_vecFourierSinTable[i] = sin (angle); - angle += angleIncrement; - } - } + m_nameFilter = convertFilterIDToName (m_idFilter); + m_nameDomain = convertDomainIDToName (m_idDomain); + m_nFilterPoints = nFilterPoints; + m_dFilterParam = dFilterParam; + m_dBandwidth = dBandwidth; + m_dFilterMin = dFilterMinimum; + m_dFilterMax = dFilterMaximum; -#if HAVE_FFTW - if (m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW) { - for (int i = 0; i < m_nFilterPoints; i++) //fftw uses unnormalized fft - m_vecFilter[i] /= m_nFilterPoints; - } + m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1); + m_adFilter = new double [m_nFilterPoints]; - if (m_idFilterMethod == FILTER_METHOD_RFFTW) { - m_realPlanForward = rfftw_create_plan (m_nFilterPoints, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE); - m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE); - m_vecRealFftInput = new fftw_real [ m_nFilterPoints ]; - m_vecRealFftSignal = new fftw_real [ m_nOutputPoints ]; - for (int i = 0; i < m_nFilterPoints; i++) - m_vecRealFftInput[i] = 0; - } else if (m_idFilterMethod == FILTER_METHOD_FFTW) { - m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE); - m_complexPlanBackward = fftw_create_plan (m_nOutputPoints, FFTW_BACKWARD, FFTW_ESTIMATE); - m_vecComplexFftInput = new fftw_complex [ m_nFilterPoints ]; - m_vecComplexFftSignal = new fftw_complex [ m_nOutputPoints ]; - for (int i = 0; i < m_nFilterPoints; i++) - m_vecComplexFftInput[i].re = m_vecComplexFftInput[i].im = 0; - for (int i = 0; i < m_nOutputPoints; i++) - m_vecComplexFftSignal[i].re = m_vecComplexFftSignal[i].im = 0; - } -#endif - - if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) { - m_nFilterPoints = 2 * m_nSignalPoints - 1; - m_filterMin = -m_signalInc * (m_nSignalPoints - 1); - m_filterMax = m_signalInc * (m_nSignalPoints - 1); - m_filterInc = (m_filterMax - m_filterMin) / (m_nFilterPoints - 1); - m_vecFilter = new double[ m_nFilterPoints ]; - - if (m_idFilter == FILTER_SHEPP) { - double a = 2 * m_bw; - double c = - 4. / (a * a); - int center = (m_nFilterPoints - 1) / 2; - int sidelen = center; - m_vecFilter[center] = 4. / (a * a); - - for (int i = 1; i <= sidelen; i++ ) - m_vecFilter [center + i] = m_vecFilter [center - i] = c / (4 * (i * i) - 1); - } else if (m_idDomain == DOMAIN_FREQUENCY) { - double x; - int i; - for (x = m_filterMin, i = 0; i < m_nFilterPoints; x += m_filterInc, i++) - m_vecFilter[i] = frequencyResponse (x, m_filterParam); - } else if (m_idDomain == DOMAIN_SPATIAL) { - double x; - int i; - for (x = m_filterMin, i = 0; i < m_nFilterPoints; x += m_filterInc, i++) { - if (haveAnalyticSpatial(m_idFilter)) - m_vecFilter[i] = spatialResponseAnalytic (x, m_filterParam); - else - m_vecFilter[i] = spatialResponseCalc (x, m_filterParam); -#if LIMIT_BANDWIDTH_TRIAL - if (i < m_nFilterPoints / 4 || i > (m_nFilterPoints * 3) / 4) - m_vecFilter[i] = 0; -#endif - } - } else { - m_failMessage = "Illegal domain name "; - m_failMessage += m_idDomain; - m_fail = true; - } - } + if (m_idDomain == DOMAIN_FREQUENCY) + createFrequencyFilter (m_adFilter); + else if (m_idDomain == DOMAIN_SPATIAL) + createSpatialFilter (m_adFilter); } + SignalFilter::~SignalFilter (void) { - delete [] m_vecFilter; - delete [] m_vecFourierSinTable; - delete [] m_vecFourierCosTable; - -#if HAVE_FFTW - if (m_idFilterMethod == FILTER_METHOD_FFTW) { - fftw_destroy_plan(m_complexPlanForward); - fftw_destroy_plan(m_complexPlanBackward); - delete [] m_vecComplexFftInput; - delete [] m_vecComplexFftSignal; - } - if (m_idFilterMethod == FILTER_METHOD_RFFTW) { - rfftw_destroy_plan(m_realPlanForward); - rfftw_destroy_plan(m_realPlanBackward); - delete [] m_vecRealFftInput; - delete [] m_vecRealFftSignal; - } -#endif + delete [] m_adFilter; +} + +void +SignalFilter::createFrequencyFilter (double* adFilter) const +{ + double x; + int i; + for (x = m_dFilterMin, i = 0; i < m_nFilterPoints; x += m_dFilterInc, i++) + adFilter[i] = frequencyResponse (x); } +void +SignalFilter::createSpatialFilter (double* adFilter) const +{ + if (m_idFilter == FILTER_SHEPP) { + double a = 2 * m_dBandwidth; + double c = - 4. / (a * a); + int center = (m_nFilterPoints - 1) / 2; + int sidelen = center; + m_adFilter[center] = 4. / (a * a); + + for (int i = 1; i <= sidelen; i++ ) + m_adFilter [center + i] = m_adFilter [center - i] = c / (4 * (i * i) - 1); + } else { + double x = m_dFilterMin; + for (int i = 0; i < m_nFilterPoints; i++, x += m_dFilterInc) { + if (haveAnalyticSpatial(m_idFilter)) + m_adFilter[i] = spatialResponseAnalytic (x); + else + m_adFilter[i] = spatialResponseCalc (x); + } + } +} + int SignalFilter::convertFilterNameToID (const char *filterName) { @@ -420,42 +258,6 @@ SignalFilter::convertFilterIDToTitle (const int filterID) return (title); } -int -SignalFilter::convertFilterMethodNameToID (const char* const filterMethodName) -{ - int fmID = FILTER_METHOD_INVALID; - - for (int i = 0; i < s_iFilterMethodCount; i++) - if (strcasecmp (filterMethodName, s_aszFilterMethodName[i]) == 0) { - fmID = i; - break; - } - - return (fmID); -} - -const char * -SignalFilter::convertFilterMethodIDToName (const int fmID) -{ - static const char *name = ""; - - if (fmID >= 0 && fmID < s_iFilterMethodCount) - return (s_aszFilterMethodName [fmID]); - - return (name); -} - -const char * -SignalFilter::convertFilterMethodIDToTitle (const int fmID) -{ - static const char *title = ""; - - if (fmID >= 0 && fmID < s_iFilterMethodCount) - return (s_aszFilterTitle [fmID]); - - return (title); -} - int SignalFilter::convertDomainNameToID (const char* const domainName) { @@ -492,110 +294,6 @@ SignalFilter::convertDomainIDToTitle (const int domainID) return (title); } -int -SignalFilter::convertFrequencyFilterNameToID (const char* const ffName) -{ - int ffID = FREQUENCY_FILTER_INVALID; - - for (int i = 0; i < s_iFrequencyFilterCount; i++) - if (strcasecmp (ffName, s_aszFrequencyFilterName[i]) == 0) { - ffID = i; - break; - } - - return (ffID); -} - -const char * -SignalFilter::convertFrequencyFilterIDToName (const int ffID) -{ - static const char *name = ""; - - if (ffID >= 0 && ffID < s_iFrequencyFilterCount) - return (s_aszFrequencyFilterName [ffID]); - - return (name); -} - -const char * -SignalFilter::convertFrequencyFilterIDToTitle (const int ffID) -{ - static const char *name = ""; - - if (ffID >= 0 && ffID < s_iFrequencyFilterCount) - return (s_aszFrequencyFilterTitle [ffID]); - - return (name); -} - -void -SignalFilter::filterSignal (const float input[], double output[]) const -{ - if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) { - for (int i = 0; i < m_nSignalPoints; i++) - output[i] = convolve (input, m_signalInc, i, m_nSignalPoints); - } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) { - double inputSignal[m_nFilterPoints]; - for (int i = 0; i < m_nSignalPoints; i++) - inputSignal[i] = input[i]; - for (int i = m_nSignalPoints; i < m_nFilterPoints; i++) - inputSignal[i] = 0; // zeropad - complex fftSignal[m_nFilterPoints]; - finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, -1); - for (int i = 0; i < m_nFilterPoints; i++) - fftSignal[i] *= m_vecFilter[i]; - double inverseFourier[m_nFilterPoints]; - 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) { - double inputSignal[m_nFilterPoints]; - for (int i = 0; i < m_nSignalPoints; i++) - inputSignal[i] = input[i]; - for (int i = m_nSignalPoints; i < m_nFilterPoints; i++) - inputSignal[i] = 0; // zeropad - complex fftSignal[m_nFilterPoints]; - finiteFourierTransform (inputSignal, fftSignal, -1); - for (int i = 0; i < m_nFilterPoints; i++) - fftSignal[i] *= m_vecFilter[i]; - double inverseFourier[m_nFilterPoints]; - finiteFourierTransform (fftSignal, inverseFourier, 1); - for (int i = 0; i < m_nSignalPoints; i++) - output[i] = inverseFourier[i]; - } -#if HAVE_FFTW - else if (m_idFilterMethod == FILTER_METHOD_RFFTW) { - for (int i = 0; i < m_nSignalPoints; i++) - m_vecRealFftInput[i] = input[i]; - - fftw_real fftOutput [ m_nFilterPoints ]; - rfftw_one (m_realPlanForward, m_vecRealFftInput, fftOutput); - for (int i = 0; i < m_nFilterPoints; i++) - m_vecRealFftSignal[i] = m_vecFilter[i] * fftOutput[i]; - for (int i = m_nFilterPoints; i < m_nOutputPoints; i++) - m_vecRealFftSignal[i] = 0; - - fftw_real ifftOutput [ m_nOutputPoints ]; - rfftw_one(m_realPlanBackward, m_vecRealFftSignal, ifftOutput); - for (int i = 0; i < m_nSignalPoints * m_preinterpolationFactor; i++) - output[i] = ifftOutput[i]; - } else if (m_idFilterMethod == FILTER_METHOD_FFTW) { - for (int i = 0; i < m_nSignalPoints; i++) - m_vecComplexFftInput[i].re = input[i]; - - fftw_complex fftOutput [ m_nFilterPoints ]; - fftw_one(m_complexPlanForward, m_vecComplexFftInput, fftOutput); - for (int i = 0; i < m_nFilterPoints; i++) { - m_vecComplexFftSignal[i].re = m_vecFilter[i] * fftOutput[i].re; - m_vecComplexFftSignal[i].im = m_vecFilter[i] * fftOutput[i].im; - } - fftw_complex ifftOutput [ m_nOutputPoints ]; - fftw_one(m_complexPlanBackward, m_vecComplexFftSignal, ifftOutput); - for (int i = 0; i < m_nSignalPoints * m_preinterpolationFactor; i++) - output[i] = ifftOutput[i].re; - } -#endif -} double SignalFilter::response (double x) @@ -603,9 +301,9 @@ SignalFilter::response (double x) double response = 0; if (m_idDomain == DOMAIN_SPATIAL) - response = spatialResponse (m_idFilter, m_bw, x, m_filterParam); + response = spatialResponse (m_idFilter, m_dBandwidth, x, m_dFilterParam); else if (m_idDomain == DOMAIN_FREQUENCY) - response = frequencyResponse (m_idFilter, m_bw, x, m_filterParam); + response = frequencyResponse (m_idFilter, m_dBandwidth, x, m_dFilterParam); return (response); } @@ -620,6 +318,16 @@ SignalFilter::spatialResponse (int filterID, double bw, double x, double param) return spatialResponseCalc (filterID, bw, x, param, N_INTEGRAL); } +void +SignalFilter::copyFilterData (double* pdFilter, const int iStart, const int nPoints) const +{ + int iFirst = clamp (iStart, 0, m_nFilterPoints - 1); + int iLast = clamp (iFirst + nPoints - 1, 0, m_nFilterPoints - 1); + + for (int i = iFirst; i <= iLast; i++) + pdFilter[i - iFirst] = m_adFilter[i]; +} + /* NAME * filter_spatial_response_calc Calculate filter by discrete inverse fourier * transform of filters's frequency @@ -636,9 +344,9 @@ SignalFilter::spatialResponse (int filterID, double bw, double x, double param) */ double -SignalFilter::spatialResponseCalc (double x, double param) const +SignalFilter::spatialResponseCalc (double x) const { - return (spatialResponseCalc (m_idFilter, m_bw, x, param, N_INTEGRAL)); + return (spatialResponseCalc (m_idFilter, m_dBandwidth, x, m_dFilterParam, N_INTEGRAL)); } double @@ -679,9 +387,9 @@ SignalFilter::spatialResponseCalc (int filterID, double bw, double x, double par */ double -SignalFilter::frequencyResponse (double u, double param) const +SignalFilter::frequencyResponse (double u) const { - return frequencyResponse (m_idFilter, m_bw, u, param); + return frequencyResponse (m_idFilter, m_dBandwidth, u, m_dFilterParam); } @@ -765,9 +473,9 @@ SignalFilter::frequencyResponse (int filterID, double bw, double u, double param */ double -SignalFilter::spatialResponseAnalytic (double x, double param) const +SignalFilter::spatialResponseAnalytic (double x) const { - return spatialResponseAnalytic (m_idFilter, m_bw, x, param); + return spatialResponseAnalytic (m_idFilter, m_dBandwidth, x, m_dFilterParam); } const bool @@ -876,234 +584,4 @@ SignalFilter::spatialResponseAnalytic (int filterID, double bw, double x, double * Returns the value of integral of u*cos(u)*dV for V = 0 to w */ -double -SignalFilter::integral_abscos (double u, double w) -{ - return (fabs (u) > F_EPSILON - ? (cos(u * w) - 1) / (u * u) + w / u * sin (u * w) - : (w * w / 2)); -} - - -/* NAME - * convolve Discrete convolution of two functions - * - * SYNOPSIS - * r = convolve (f1, f2, dx, n, np, func_type) - * double r Convolved result - * double f1[], f2[] Functions to be convolved - * double dx Difference between successive x values - * int n Array index to center convolution about - * int np Number of points in f1 array - * int func_type EVEN or ODD or EVEN_AND_ODD function f2 - * - * NOTES - * f1 is the projection data, its indices range from 0 to np - 1. - * The index for f2, the filter, ranges from -(np-1) to (np-1). - * There are 3 ways to handle the negative vertices of f2: - * 1. If we know f2 is an EVEN function, then f2[-n] = f2[n]. - * All filters used in reconstruction are even. - * 2. If we know f2 is an ODD function, then f2[-n] = -f2[n] - * 3. If f2 is both ODD AND EVEN, then we must store the value of f2 - * for negative indices. Since f2 must range from -(np-1) to (np-1), - * if we add (np - 1) to f2's array index, then f2's index will - * range from 0 to 2 * (np - 1), and the origin, x = 0, will be - * stored at f2[np-1]. - */ - -double -SignalFilter::convolve (const double func[], const double dx, const int n, const int np) const -{ - double sum = 0.0; - -#if UNOPTIMIZED_CONVOLUTION - for (int i = 0; i < np; i++) - sum += func[i] * m_vecFilter[n - i + (np - 1)]; -#else - double* f2 = m_vecFilter + n + (np - 1); - for (int i = 0; i < np; i++) - sum += *func++ * *f2--; -#endif - - return (sum * dx); -} - - -double -SignalFilter::convolve (const float func[], const double dx, const int n, const int np) const -{ - double sum = 0.0; - -#if UNOPTIMIZED_CONVOLUTION -for (int i = 0; i < np; i++) - sum += func[i] * m_vecFilter[n - i + (np - 1)]; -#else -double* f2 = m_vecFilter + n + (np - 1); -for (int i = 0; i < np; i++) - sum += *func++ * *f2--; -#endif - - return (sum * dx); -} - - -void -SignalFilter::finiteFourierTransform (const double input[], complex output[], const int n, int direction) -{ - if (direction < 0) - direction = -1; - else - direction = 1; - - double angleIncrement = direction * 2 * PI / n; - for (int i = 0; i < n; i++) { - double sumReal = 0; - double sumImag = 0; - for (int j = 0; j < n; j++) { - double angle = i * j * angleIncrement; - sumReal += input[j] * cos(angle); - sumImag += input[j] * sin(angle); - } - if (direction < 0) { - sumReal /= n; - sumImag /= n; - } - output[i] = complex (sumReal, sumImag); - } -} - - -void -SignalFilter::finiteFourierTransform (const complex input[], complex output[], const int n, int direction) -{ - if (direction < 0) - direction = -1; - else - direction = 1; - - double angleIncrement = direction * 2 * PI / n; - for (int i = 0; i < n; i++) { - complex sum (0,0); - for (int j = 0; j < n; j++) { - double angle = i * j * angleIncrement; - complex exponentTerm (cos(angle), sin(angle)); - sum += input[j] * exponentTerm; - } - if (direction < 0) { - sum /= n; - } - output[i] = sum; - } -} - -void -SignalFilter::finiteFourierTransform (const complex input[], double output[], const int n, int direction) -{ - if (direction < 0) - direction = -1; - else - direction = 1; - - double angleIncrement = direction * 2 * PI / n; - for (int i = 0; i < n; i++) { - double sumReal = 0; - for (int j = 0; j < n; j++) { - double angle = i * j * angleIncrement; - sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle); - } - if (direction < 0) { - sumReal /= n; - } - output[i] = sumReal; - } -} - -void -SignalFilter::finiteFourierTransform (const double input[], complex output[], int direction) const -{ - if (direction < 0) - direction = -1; - else - direction = 1; - - for (int i = 0; i < m_nFilterPoints; i++) { - double sumReal = 0, sumImag = 0; - for (int j = 0; j < m_nFilterPoints; j++) { - int tableIndex = i * j; - if (direction > 0) { - sumReal += input[j] * m_vecFourierCosTable[tableIndex]; - sumImag += input[j] * m_vecFourierSinTable[tableIndex]; - } else { - sumReal += input[j] * m_vecFourierCosTable[tableIndex]; - sumImag -= input[j] * m_vecFourierSinTable[tableIndex]; - } - } - if (direction < 0) { - sumReal /= m_nFilterPoints; - sumImag /= m_nFilterPoints; - } - output[i] = complex (sumReal, sumImag); - } -} - -// (a+bi) * (c + di) = (ac - bd) + (ad + bc)i -void -SignalFilter::finiteFourierTransform (const complex input[], complex output[], int direction) const -{ - if (direction < 0) - direction = -1; - else - direction = 1; - - for (int i = 0; i < m_nFilterPoints; i++) { - double sumReal = 0, sumImag = 0; - for (int j = 0; j < m_nFilterPoints; j++) { - int tableIndex = i * j; - if (direction > 0) { - sumReal += input[j].real() * m_vecFourierCosTable[tableIndex] - - input[j].imag() * m_vecFourierSinTable[tableIndex]; - sumImag += input[j].real() * m_vecFourierSinTable[tableIndex] - + input[j].imag() * m_vecFourierCosTable[tableIndex]; - } else { - sumReal += input[j].real() * m_vecFourierCosTable[tableIndex] - - input[j].imag() * -m_vecFourierSinTable[tableIndex]; - sumImag += input[j].real() * -m_vecFourierSinTable[tableIndex] - + input[j].imag() * m_vecFourierCosTable[tableIndex]; - } - } - if (direction < 0) { - sumReal /= m_nFilterPoints; - sumImag /= m_nFilterPoints; - } - output[i] = complex (sumReal, sumImag); - } -} - -void -SignalFilter::finiteFourierTransform (const complex input[], double output[], int direction) const -{ - if (direction < 0) - direction = -1; - else - direction = 1; - - for (int i = 0; i < m_nFilterPoints; i++) { - double sumReal = 0; - for (int j = 0; j < m_nFilterPoints; j++) { - int tableIndex = i * j; - if (direction > 0) { - sumReal += input[j].real() * m_vecFourierCosTable[tableIndex] - - input[j].imag() * m_vecFourierSinTable[tableIndex]; - } else { - sumReal += input[j].real() * m_vecFourierCosTable[tableIndex] - - input[j].imag() * -m_vecFourierSinTable[tableIndex]; - } - } - if (direction < 0) { - sumReal /= m_nFilterPoints; - } - output[i] = sumReal; - } -} - diff --git a/libctsim/procsignal.cpp b/libctsim/procsignal.cpp new file mode 100644 index 0000000..9fdbb68 --- /dev/null +++ b/libctsim/procsignal.cpp @@ -0,0 +1,709 @@ +/***************************************************************************** +** File IDENTIFICATION +** +** Name: filter.cpp +** Purpose: Routines for signal-procesing filters +** Progammer: Kevin Rosenberg +** Date Started: Aug 1984 +** +** This is part of the CTSim program +** Copyright (C) 1983-2000 Kevin Rosenberg +** +** $Id: procsignal.cpp,v 1.1 2000/08/19 23:00:05 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 +** published by the Free Software Foundation. +** +** This program is distributed in the hope that it will be useful, +** but WITHOUT ANY WARRANTY; without even the implied warranty of +** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +** GNU General Public License for more details. +** +** You should have received a copy of the GNU General Public License +** along with this program; if not, write to the Free Software +** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +******************************************************************************/ + +#include "ct.h" + +// FilterMethod ID/Names +const int ProcessSignal::FILTER_METHOD_INVALID = -1; +const int ProcessSignal::FILTER_METHOD_CONVOLUTION = 0; +const int ProcessSignal::FILTER_METHOD_FOURIER = 1; +const int ProcessSignal::FILTER_METHOD_FOURIER_TABLE = 2; +const int ProcessSignal::FILTER_METHOD_FFT = 3; +#if HAVE_FFTW +const int ProcessSignal::FILTER_METHOD_FFTW = 4; +const int ProcessSignal::FILTER_METHOD_RFFTW =5 ; +#endif +const char* ProcessSignal::s_aszFilterMethodName[] = { + {"convolution"}, + {"fourier"}, + {"fouier_table"}, + {"fft"}, +#if HAVE_FFTW + {"fftw"}, + {"rfftw"}, +#endif +}; +const char* ProcessSignal::s_aszFilterMethodTitle[] = { + {"Convolution"}, + {"Direct Fourier"}, + {"Fouier Trigometric Table Lookout"}, + {"FFT"}, +#if HAVE_FFTW + {"FFTW"}, + {"Real/Half-Complex FFTW"}, +#endif +}; +const int ProcessSignal::s_iFilterMethodCount = sizeof(s_aszFilterMethodName) / sizeof(const char*); + +// FilterGeneration ID/Names +const int ProcessSignal::FILTER_GENERATION_INVALID = -1; +const int ProcessSignal::FILTER_GENERATION_DIRECT = 0; +const int ProcessSignal::FILTER_GENERATION_INVERSE_FOURIER = 1; +const char* ProcessSignal::s_aszFilterGenerationName[] = { + {"direct"}, + {"inverse_fourier"}, +}; +const char* ProcessSignal::s_aszFilterGenerationTitle[] = { + {"Direct"}, + {"Inverse Fourier"}, +}; +const int ProcessSignal::s_iFilterGenerationCount = sizeof(s_aszFilterGenerationName) / sizeof(const char*); + + +// CLASS IDENTIFICATION +// ProcessSignal +// +ProcessSignal::ProcessSignal (const char* szFilterName, const char* szFilterMethodName, double dBandwidth, double dSignalIncrement, int nSignalPoints, double dFilterParam, const char* szDomainName, const char* szFilterGenerationName, int iZeropad = 0, int iPreinterpolationFactor = 1) + : m_adFourierCosTable(NULL), m_adFourierSinTable(NULL), m_adFilter(NULL), m_fail(false) +{ + m_idFilterMethod = convertFilterMethodNameToID (szFilterMethodName); + if (m_idFilterMethod == FILTER_METHOD_INVALID) { + m_fail = true; + m_failMessage = "Invalid filter method name "; + m_failMessage += szFilterMethodName; + return; + } + m_idFilterGeneration = convertFilterGenerationNameToID (szFilterGenerationName); + if (m_idFilterGeneration == FILTER_GENERATION_INVALID) { + m_fail = true; + m_failMessage = "Invalid frequency filter name "; + m_failMessage += szFilterGenerationName; + return; + } + m_idFilter = SignalFilter::convertFilterNameToID (szFilterName); + if (m_idFilter == SignalFilter::FILTER_INVALID) { + m_fail = true; + m_failMessage = "Invalid Filter name "; + m_failMessage += szFilterName; + return; + } + m_idDomain = SignalFilter::convertDomainNameToID (szDomainName); + if (m_idDomain == SignalFilter::DOMAIN_INVALID) { + m_fail = true; + m_failMessage = "Invalid domain name "; + m_failMessage += szDomainName; + return; + } + + init (m_idFilter, m_idFilterMethod, dBandwidth, dSignalIncrement, nSignalPoints, dFilterParam, m_idDomain, m_idFilterGeneration, iZeropad, iPreinterpolationFactor); +} + + +void +ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandwidth, double dSignalIncrement, int nSignalPoints, double dFilterParam, const int idDomain, const int idFilterGeneration, const int iZeropad, const int iPreinterpolationFactor) +{ + m_idFilter = idFilter; + m_idDomain = idDomain; + m_idFilterMethod = idFilterMethod; + m_idFilterGeneration = idFilterGeneration; + if (m_idFilter == SignalFilter::FILTER_INVALID || m_idDomain == SignalFilter::DOMAIN_INVALID || m_idFilterMethod == FILTER_METHOD_INVALID || m_idFilterGeneration == FILTER_GENERATION_INVALID) { + m_fail = true; + return; + } + m_traceLevel = TRACE_NONE; + m_nameFilterMethod = convertFilterMethodIDToName (m_idFilterMethod); + m_nameFilterGeneration = convertFilterGenerationIDToName (m_idFilterGeneration); + m_dBandwidth = dBandwidth; + m_nSignalPoints = nSignalPoints; + m_dSignalInc = dSignalIncrement; + m_dFilterParam = dFilterParam; + m_iZeropad = iZeropad; + m_iPreinterpolationFactor = iPreinterpolationFactor; + + 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 + } + + bool m_bFrequencyFiltering = true; + if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) + m_bFrequencyFiltering = false; + + // Spatial-based filtering + if (! m_bFrequencyFiltering) { + + if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) { + m_nFilterPoints = 2 * m_nSignalPoints - 1; + m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1); + m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1); + m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1); + SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL); + m_adFilter = new double[ m_nFilterPoints ]; + filter.copyFilterData (m_adFilter, 0, m_nFilterPoints); + } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) { + m_nFilterPoints = m_nSignalPoints; + m_dFilterMin = -1. / (2 * m_dSignalInc); + m_dFilterMax = 1. / (2 * m_dSignalInc); + m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1); + SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY); + m_adFilter = new double[ m_nFilterPoints ]; + double adFrequencyFilter [m_nFilterPoints]; + double adInverseFilter [m_nFilterPoints]; + filter.copyFilterData (adFrequencyFilter, 0, m_nFilterPoints); + shuffleNaturalToFourierOrder (adFrequencyFilter, m_nFilterPoints); + ProcessSignal::finiteFourierTransform (adFrequencyFilter, adInverseFilter, m_nFilterPoints, 1); + for (int i = 0; i < m_nFilterPoints; i++) + m_adFilter [i] = adInverseFilter[i]; + } + } + + // Frequency-based filtering + else if (m_bFrequencyFiltering) { + + // calculate number of filter points with zeropadding + m_nFilterPoints = m_nSignalPoints; + if (m_iZeropad > 0) { + double logBase2 = log(m_nSignalPoints) / log(2); + int nextPowerOf2 = static_cast(floor(logBase2)); + if (logBase2 != floor(logBase2)) + nextPowerOf2++; + nextPowerOf2 += (m_iZeropad - 1); + m_nFilterPoints = 1 << nextPowerOf2; + if (m_traceLevel >= TRACE_TEXT) + cout << "nFilterPoints = " << m_nFilterPoints << endl; + } + m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor; + + if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) { + m_dFilterMin = -1. / (2 * m_dSignalInc); + m_dFilterMax = 1. / (2 * m_dSignalInc); + m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1); + SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY); + m_adFilter = new double [m_nFilterPoints]; + filter.copyFilterData (m_adFilter, 0, m_nFilterPoints); + shuffleNaturalToFourierOrder (m_adFilter, m_nFilterPoints); + } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) { + m_nFilterPoints = 2 * m_nSignalPoints - 1; + m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1); + m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1); + m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1); + double adSpatialFilter [m_nFilterPoints]; + double adInverseFilter [m_nFilterPoints]; + SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL); + filter.copyFilterData (adSpatialFilter, 0, m_nFilterPoints); + m_adFilter = new double [m_nFilterPoints]; + finiteFourierTransform (adSpatialFilter, adInverseFilter, m_nFilterPoints, -1); + for (int i = 0; i < m_nFilterPoints; i++) + m_adFilter [i] = adInverseFilter[i]; + } + } + + // precalculate sin and cosine tables for fourier transform + if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) { + int nFourier = max(m_nFilterPoints,m_nOutputPoints) * max(m_nFilterPoints, m_nOutputPoints) + 1; + double angleIncrement = (2. * PI) / m_nFilterPoints; + m_adFourierCosTable = new double[ nFourier ]; + m_adFourierSinTable = new double[ nFourier ]; + double angle = 0; + for (int i = 0; i < nFourier; i++) { + m_adFourierCosTable[i] = cos (angle); + m_adFourierSinTable[i] = sin (angle); + angle += angleIncrement; + } + } + +#if HAVE_FFTW + if (m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW) { + for (int i = 0; i < m_nFilterPoints; i++) //fftw uses unnormalized fft + m_adFilter[i] /= m_nFilterPoints; + } + + if (m_idFilterMethod == FILTER_METHOD_RFFTW) { + m_realPlanForward = rfftw_create_plan (m_nFilterPoints, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE); + m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE); + m_adRealFftInput = new fftw_real [ m_nFilterPoints ]; + m_adRealFftSignal = new fftw_real [ m_nOutputPoints ]; + for (int i = 0; i < m_nFilterPoints; i++) + m_adRealFftInput[i] = 0; + } else if (m_idFilterMethod == FILTER_METHOD_FFTW) { + m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE); + m_complexPlanBackward = fftw_create_plan (m_nOutputPoints, FFTW_BACKWARD, FFTW_ESTIMATE); + m_adComplexFftInput = new fftw_complex [ m_nFilterPoints ]; + m_adComplexFftSignal = new fftw_complex [ m_nOutputPoints ]; + for (int i = 0; i < m_nFilterPoints; i++) + m_adComplexFftInput[i].re = m_adComplexFftInput[i].im = 0; + for (int i = 0; i < m_nOutputPoints; i++) + m_adComplexFftSignal[i].re = m_adComplexFftSignal[i].im = 0; + } +#endif + +} + +ProcessSignal::~ProcessSignal (void) +{ + delete [] m_adFourierSinTable; + delete [] m_adFourierCosTable; + +#if HAVE_FFTW + if (m_idFilterMethod == FILTER_METHOD_FFTW) { + fftw_destroy_plan(m_complexPlanForward); + fftw_destroy_plan(m_complexPlanBackward); + delete [] m_adComplexFftInput; + delete [] m_adComplexFftSignal; + } + if (m_idFilterMethod == FILTER_METHOD_RFFTW) { + rfftw_destroy_plan(m_realPlanForward); + rfftw_destroy_plan(m_realPlanBackward); + delete [] m_adRealFftInput; + delete [] m_adRealFftSignal; + } +#endif +} + +int +ProcessSignal::convertFilterMethodNameToID (const char* const filterMethodName) +{ + int fmID = FILTER_METHOD_INVALID; + + for (int i = 0; i < s_iFilterMethodCount; i++) + if (strcasecmp (filterMethodName, s_aszFilterMethodName[i]) == 0) { + fmID = i; + break; + } + + return (fmID); +} + +const char * +ProcessSignal::convertFilterMethodIDToName (const int fmID) +{ + static const char *name = ""; + + if (fmID >= 0 && fmID < s_iFilterMethodCount) + return (s_aszFilterMethodName [fmID]); + + return (name); +} + +const char * +ProcessSignal::convertFilterMethodIDToTitle (const int fmID) +{ + static const char *title = ""; + + if (fmID >= 0 && fmID < s_iFilterMethodCount) + return (s_aszFilterMethodTitle [fmID]); + + return (title); +} + + +int +ProcessSignal::convertFilterGenerationNameToID (const char* const fgName) +{ + int fgID = FILTER_GENERATION_INVALID; + + for (int i = 0; i < s_iFilterGenerationCount; i++) + if (strcasecmp (fgName, s_aszFilterGenerationName[i]) == 0) { + fgID = i; + break; + } + + return (fgID); +} + +const char * +ProcessSignal::convertFilterGenerationIDToName (const int fgID) +{ + static const char *name = ""; + + if (fgID >= 0 && fgID < s_iFilterGenerationCount) + return (s_aszFilterGenerationName [fgID]); + + return (name); +} + +const char * +ProcessSignal::convertFilterGenerationIDToTitle (const int fgID) +{ + static const char *name = ""; + + if (fgID >= 0 && fgID < s_iFilterGenerationCount) + return (s_aszFilterGenerationTitle [fgID]); + + return (name); +} + +void +ProcessSignal::filterSignal (const float input[], double output[]) const +{ + if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) { + for (int i = 0; i < m_nSignalPoints; i++) + output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints); + } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) { + double inputSignal[m_nFilterPoints]; + for (int i = 0; i < m_nSignalPoints; i++) + inputSignal[i] = input[i]; + for (int i = m_nSignalPoints; i < m_nFilterPoints; i++) + inputSignal[i] = 0; // zeropad + complex fftSignal[m_nFilterPoints]; + finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, -1); + for (int i = 0; i < m_nFilterPoints; i++) + fftSignal[i] *= m_adFilter[i]; + double inverseFourier[m_nFilterPoints]; + 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) { + double inputSignal[m_nFilterPoints]; + for (int i = 0; i < m_nSignalPoints; i++) + inputSignal[i] = input[i]; + for (int i = m_nSignalPoints; i < m_nFilterPoints; i++) + inputSignal[i] = 0; // zeropad + complex fftSignal[m_nFilterPoints]; + finiteFourierTransform (inputSignal, fftSignal, -1); + for (int i = 0; i < m_nFilterPoints; i++) + fftSignal[i] *= m_adFilter[i]; + double inverseFourier[m_nFilterPoints]; + finiteFourierTransform (fftSignal, inverseFourier, 1); + for (int i = 0; i < m_nSignalPoints; i++) + output[i] = inverseFourier[i]; + } +#if HAVE_FFTW + else if (m_idFilterMethod == FILTER_METHOD_RFFTW) { + for (int i = 0; i < m_nSignalPoints; i++) + m_adRealFftInput[i] = input[i]; + + fftw_real fftOutput [ m_nFilterPoints ]; + rfftw_one (m_realPlanForward, m_adRealFftInput, fftOutput); + for (int i = 0; i < m_nFilterPoints; i++) + m_adRealFftSignal[i] = m_adFilter[i] * fftOutput[i]; + for (int i = m_nFilterPoints; i < m_nOutputPoints; i++) + m_adRealFftSignal[i] = 0; + + fftw_real ifftOutput [ m_nOutputPoints ]; + rfftw_one (m_realPlanBackward, m_adRealFftSignal, ifftOutput); + for (int i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++) + output[i] = ifftOutput[i]; + } else if (m_idFilterMethod == FILTER_METHOD_FFTW) { + for (int i = 0; i < m_nSignalPoints; i++) + m_adComplexFftInput[i].re = input[i]; + + fftw_complex fftOutput [ m_nFilterPoints ]; + fftw_one (m_complexPlanForward, m_adComplexFftInput, fftOutput); + for (int i = 0; i < m_nFilterPoints; i++) { + m_adComplexFftSignal[i].re = m_adFilter[i] * fftOutput[i].re; + m_adComplexFftSignal[i].im = m_adFilter[i] * fftOutput[i].im; + } + fftw_complex ifftOutput [ m_nOutputPoints ]; + fftw_one (m_complexPlanBackward, m_adComplexFftSignal, ifftOutput); + for (int i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++) + output[i] = ifftOutput[i].re; + } +#endif +} + + +/* NAME + * convolve Discrete convolution of two functions + * + * SYNOPSIS + * r = convolve (f1, f2, dx, n, np, func_type) + * double r Convolved result + * double f1[], f2[] Functions to be convolved + * double dx Difference between successive x values + * int n Array index to center convolution about + * int np Number of points in f1 array + * int func_type EVEN or ODD or EVEN_AND_ODD function f2 + * + * NOTES + * f1 is the projection data, its indices range from 0 to np - 1. + * The index for f2, the filter, ranges from -(np-1) to (np-1). + * There are 3 ways to handle the negative vertices of f2: + * 1. If we know f2 is an EVEN function, then f2[-n] = f2[n]. + * All filters used in reconstruction are even. + * 2. If we know f2 is an ODD function, then f2[-n] = -f2[n] + * 3. If f2 is both ODD AND EVEN, then we must store the value of f2 + * for negative indices. Since f2 must range from -(np-1) to (np-1), + * if we add (np - 1) to f2's array index, then f2's index will + * range from 0 to 2 * (np - 1), and the origin, x = 0, will be + * stored at f2[np-1]. + */ + +double +ProcessSignal::convolve (const double func[], const double dx, const int n, const int np) const +{ + double sum = 0.0; + +#if UNOPTIMIZED_CONVOLUTION + for (int i = 0; i < np; i++) + sum += func[i] * m_adFilter[n - i + (np - 1)]; +#else + double* f2 = m_adFilter + n + (np - 1); + for (int i = 0; i < np; i++) + sum += *func++ * *f2--; +#endif + + return (sum * dx); +} + + +double +ProcessSignal::convolve (const float func[], const double dx, const int n, const int np) const +{ + double sum = 0.0; + +#if UNOPTIMIZED_CONVOLUTION +for (int i = 0; i < np; i++) + sum += func[i] * m_adFilter[n - i + (np - 1)]; +#else +double* f2 = m_adFilter + n + (np - 1); +for (int i = 0; i < np; i++) + sum += *func++ * *f2--; +#endif + + return (sum * dx); +} + + +void +ProcessSignal::finiteFourierTransform (const double input[], double output[], const int n, int direction) +{ + complex complexOutput[n]; + + finiteFourierTransform (input, complexOutput, n, direction); + for (int i = 0; i < n; i++) + output[i] = abs(complexOutput[n]); +} + +void +ProcessSignal::finiteFourierTransform (const double input[], complex output[], const int n, int direction) +{ + if (direction < 0) + direction = -1; + else + direction = 1; + + double angleIncrement = direction * 2 * PI / n; + for (int i = 0; i < n; i++) { + double sumReal = 0; + double sumImag = 0; + for (int j = 0; j < n; j++) { + double angle = i * j * angleIncrement; + sumReal += input[j] * cos(angle); + sumImag += input[j] * sin(angle); + } + if (direction < 0) { + sumReal /= n; + sumImag /= n; + } + output[i] = complex (sumReal, sumImag); + } +} + + +void +ProcessSignal::finiteFourierTransform (const complex input[], complex output[], const int n, int direction) +{ + if (direction < 0) + direction = -1; + else + direction = 1; + + double angleIncrement = direction * 2 * PI / n; + for (int i = 0; i < n; i++) { + complex sum (0,0); + for (int j = 0; j < n; j++) { + double angle = i * j * angleIncrement; + complex exponentTerm (cos(angle), sin(angle)); + sum += input[j] * exponentTerm; + } + if (direction < 0) { + sum /= n; + } + output[i] = sum; + } +} + +void +ProcessSignal::finiteFourierTransform (const complex input[], double output[], const int n, int direction) +{ + if (direction < 0) + direction = -1; + else + direction = 1; + + double angleIncrement = direction * 2 * PI / n; + for (int i = 0; i < n; i++) { + double sumReal = 0; + for (int j = 0; j < n; j++) { + double angle = i * j * angleIncrement; + sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle); + } + if (direction < 0) { + sumReal /= n; + } + output[i] = sumReal; + } +} + +// Table-based routines + +void +ProcessSignal::finiteFourierTransform (const double input[], complex output[], int direction) const +{ + if (direction < 0) + direction = -1; + else + direction = 1; + + for (int i = 0; i < m_nFilterPoints; i++) { + double sumReal = 0, sumImag = 0; + for (int j = 0; j < m_nFilterPoints; j++) { + int tableIndex = i * j; + if (direction > 0) { + sumReal += input[j] * m_adFourierCosTable[tableIndex]; + sumImag += input[j] * m_adFourierSinTable[tableIndex]; + } else { + sumReal += input[j] * m_adFourierCosTable[tableIndex]; + sumImag -= input[j] * m_adFourierSinTable[tableIndex]; + } + } + if (direction < 0) { + sumReal /= m_nFilterPoints; + sumImag /= m_nFilterPoints; + } + output[i] = complex (sumReal, sumImag); + } +} + +// (a+bi) * (c + di) = (ac - bd) + (ad + bc)i +void +ProcessSignal::finiteFourierTransform (const complex input[], complex output[], int direction) const +{ + if (direction < 0) + direction = -1; + else + direction = 1; + + for (int i = 0; i < m_nFilterPoints; i++) { + double sumReal = 0, sumImag = 0; + for (int j = 0; j < m_nFilterPoints; j++) { + int tableIndex = i * j; + if (direction > 0) { + sumReal += input[j].real() * m_adFourierCosTable[tableIndex] + - input[j].imag() * m_adFourierSinTable[tableIndex]; + sumImag += input[j].real() * m_adFourierSinTable[tableIndex] + + input[j].imag() * m_adFourierCosTable[tableIndex]; + } else { + sumReal += input[j].real() * m_adFourierCosTable[tableIndex] + - input[j].imag() * -m_adFourierSinTable[tableIndex]; + sumImag += input[j].real() * -m_adFourierSinTable[tableIndex] + + input[j].imag() * m_adFourierCosTable[tableIndex]; + } + } + if (direction < 0) { + sumReal /= m_nFilterPoints; + sumImag /= m_nFilterPoints; + } + output[i] = complex (sumReal, sumImag); + } +} + +void +ProcessSignal::finiteFourierTransform (const complex input[], double output[], int direction) const +{ + if (direction < 0) + direction = -1; + else + direction = 1; + + for (int i = 0; i < m_nFilterPoints; i++) { + double sumReal = 0; + for (int j = 0; j < m_nFilterPoints; j++) { + int tableIndex = i * j; + if (direction > 0) { + sumReal += input[j].real() * m_adFourierCosTable[tableIndex] + - input[j].imag() * m_adFourierSinTable[tableIndex]; + } else { + sumReal += input[j].real() * m_adFourierCosTable[tableIndex] + - input[j].imag() * -m_adFourierSinTable[tableIndex]; + } + } + if (direction < 0) { + sumReal /= m_nFilterPoints; + } + output[i] = sumReal; + } +} + +// Odd Number of Points +// Natural Frequency Order: -(n-1)/2...-1,0,1...(n-1)/2 +// Fourier Frequency Order: 0, 1..(n-1)/2,-(n-1)/2...-1 +// Even Number of Points +// Natural Frequency Order: -n/2...-1,0,1...((n/2)-1) +// Fourier Frequency Order: 0,1...((n/2)-1),-n/2...-1 + +void +ProcessSignal::shuffleNaturalToFourierOrder (double* pdVector, const int n) +{ + double* pdTemp = new double [n]; + if (n % 2) { // Odd + int iHalfN = (n - 1) / 2; + + pdTemp[0] = pdVector[iHalfN]; + for (int i = 1; i <= iHalfN; i++) + pdTemp[i] = pdVector[i+iHalfN]; + for (int i = iHalfN+1; i < n; i++) + pdTemp[i] = pdVector[i-iHalfN]; + } else { // Even + int iHalfN = n / 2; + pdTemp[0] = pdVector[iHalfN]; + } + + for (int i = 0; i < n; i++) + pdVector[i] = pdTemp[i]; + delete pdTemp; +} + + +void +ProcessSignal::shuffleFourierToNaturalOrder (double* pdVector, const int n) +{ + double* pdTemp = new double [n]; + if (n % 2) { // Odd + int iHalfN = (n - 1) / 2; + + pdTemp[iHalfN] = pdVector[0]; + for (int i = 1; i <= iHalfN; i++) + pdTemp[i] = pdVector[i+iHalfN]; + for (int i = iHalfN+1; i < n; i++) + pdTemp[i] = pdVector[i-iHalfN]; + } else { // Even + int iHalfN = n / 2; + pdTemp[iHalfN] = pdVector[0]; + } + + for (int i = 0; i < n; i++) + pdVector[i] = pdTemp[i]; + delete pdTemp; +} + diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp index 73e471e..813d7bc 100644 --- a/libctsim/projections.cpp +++ b/libctsim/projections.cpp @@ -8,7 +8,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: projections.cpp,v 1.19 2000/08/03 09:57:33 kevin Exp $ +** $Id: projections.cpp,v 1.20 2000/08/19 22:59:06 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 @@ -486,7 +486,7 @@ Projections::printScanInfo (void) const */ bool -Projections::reconstruct (ImageFile& im, const char* const filterName, double filt_param, const char* const filterMethodName, const int zeropad, const char* frequencyFilterName, const char* const interpName, int interpFactor, const char* const backprojectName, const int trace) const +Projections::reconstruct (ImageFile& im, const char* const filterName, double filt_param, const char* const filterMethodName, const int zeropad, const char* filterGenerationName, const char* const interpName, int interpFactor, const char* const backprojectName, const int trace) const { double detInc = m_detInc; int n_filteredProj = m_nDet * interpFactor; @@ -503,11 +503,11 @@ Projections::reconstruct (ImageFile& im, const char* const filterName, double fi #endif double filterBW = 1. / detInc; - SignalFilter filter (filterName, filterMethodName, filterBW, m_detInc, m_nDet, filt_param, "spatial", frequencyFilterName, zeropad, interpFactor); - filter.setTraceLevel(trace); + ProcessSignal processSignal (filterName, filterMethodName, filterBW, m_detInc, m_nDet, filt_param, "spatial", filterGenerationName, zeropad, interpFactor); + processSignal.setTraceLevel(trace); - if (filter.fail()) { - sys_error (ERR_SEVERE, "%s [Projections::reconstruct]", filter.failMessage().c_str()); + if (processSignal.fail()) { + sys_error (ERR_SEVERE, "%s [Projections::reconstruct]", processSignal.failMessage().c_str()); return false; } @@ -515,23 +515,23 @@ Projections::reconstruct (ImageFile& im, const char* const filterName, double fi cout << "Reconstruct: filter="<= 0) - filterDesc << optFilterName << ": alpha=" << optFilterParam; + if (dOptFilterParam >= 0) + filterDesc << sOptFilterName << ": alpha=" << dOptFilterParam; else - filterDesc << optFilterName; + filterDesc << sOptFilterName; ostringstream label; - label << "pjrec: " << nx << "x" << ny << ", " << filterDesc.str() << ", " << optInterpName << ", preinterpolationFactor=" << optPreinterpolationFactor << ", " << optBackprojName; - remark = label.str(); + label << "pjrec: " << nx << "x" << ny << ", " << filterDesc.str() << ", " << sOptInterpName << ", preinterpolationFactor=" << iOptPreinterpolationFactor << ", " << sOptBackprojectName; + sRemark = label.str(); - if (optVerbose) - cout << "Remark: " << remark << endl; + if (bOptVerbose) + cout << "SRemark: " << sRemark << endl; #ifdef HAVE_MPI } #endif #ifdef HAVE_MPI if (mpiWorld.getRank() == 0) { - projGlobal.read (filenameProj); - if (optVerbose) + projGlobal.read (pszFilenameProj); + if (bOptVerbose) projGlobal.printScanInfo(); mpi_ndet = projGlobal.nDet(); @@ -271,16 +271,16 @@ pjrec_main (int argc, char * argv[]) } TimerCollectiveMPI timerBcast (mpiWorld.getComm()); - mpiWorld.BcastString (optBackprojName); - mpiWorld.BcastString (optFilterName); - mpiWorld.BcastString (optFilterMethodName); - mpiWorld.BcastString (optInterpName); - mpiWorld.getComm().Bcast (&optVerbose, 1, MPI::INT, 0); - mpiWorld.getComm().Bcast (&optDebug, 1, MPI::INT, 0); + mpiWorld.BcastString (sOptBackprojectName); + mpiWorld.BcastString (sOptFilterName); + mpiWorld.BcastString (sOptFilterMethodName); + mpiWorld.BcastString (sOptInterpName); + mpiWorld.getComm().Bcast (&bOptVerbose, 1, MPI::INT, 0); + mpiWorld.getComm().Bcast (&bOptDebug, 1, MPI::INT, 0); mpiWorld.getComm().Bcast (&optTrace, 1, MPI::INT, 0); - mpiWorld.getComm().Bcast (&optFilterParam, 1, MPI::DOUBLE, 0); - mpiWorld.getComm().Bcast (&optZeroPad, 1, MPI::INT, 0); - mpiWorld.getComm().Bcast (&optPreinterpolationFactor, 1, MPI::INT, 0); + mpiWorld.getComm().Bcast (&dOptFilterParam, 1, MPI::DOUBLE, 0); + mpiWorld.getComm().Bcast (&iOptZeropad, 1, MPI::INT, 0); + mpiWorld.getComm().Bcast (&iOptPreinterpolationFactor, 1, MPI::INT, 0); mpiWorld.getComm().Bcast (&mpi_ndet, 1, MPI::INT, 0); mpiWorld.getComm().Bcast (&mpi_nview, 1, MPI::INT, 0); mpiWorld.getComm().Bcast (&mpi_detinc, 1, MPI::DOUBLE, 0); @@ -288,7 +288,7 @@ pjrec_main (int argc, char * argv[]) mpiWorld.getComm().Bcast (&mpi_rotinc, 1, MPI::DOUBLE, 0); mpiWorld.getComm().Bcast (&nx, 1, MPI::INT, 0); mpiWorld.getComm().Bcast (&ny, 1, MPI::INT, 0); - if (optVerbose) + if (bOptVerbose) timerBcast.timerEndAndReport ("Time to broadcast variables"); mpiWorld.setTotalWorkUnits (mpi_nview); @@ -299,8 +299,8 @@ pjrec_main (int argc, char * argv[]) projLocal.setRotInc (mpi_rotinc); TimerCollectiveMPI timerScatter (mpiWorld.getComm()); - ScatterProjectionsMPI (mpiWorld, projGlobal, projLocal, optDebug); - if (optVerbose) + ScatterProjectionsMPI (mpiWorld, projGlobal, projLocal, bOptDebug); + if (bOptVerbose) timerScatter.timerEndAndReport ("Time to scatter projections"); if (mpiWorld.getRank() == 0) { @@ -309,8 +309,8 @@ pjrec_main (int argc, char * argv[]) imLocal = new ImageFile (nx, ny); #else - projGlobal.read (filenameProj); - if (optVerbose) + projGlobal.read (pszFilenameProj); + if (bOptVerbose) projGlobal.printScanInfo(); imGlobal = new ImageFile (nx, ny); @@ -318,28 +318,28 @@ pjrec_main (int argc, char * argv[]) #ifdef HAVE_MPI TimerCollectiveMPI timerReconstruct (mpiWorld.getComm()); - projLocal.reconstruct (*imLocal, optFilterName.c_str(), optFilterParam, optFilterMethodName.c_str(), optZeroPad, optFrequencyFilterName.c_str(), optInterpName.c_str(), optPreinterpolationFactor, optBackprojName.c_str(), optTrace); - if (optVerbose) + projLocal.reconstruct (*imLocal, sOptFilterName.c_str(), dOptFilterParam, sOptFilterMethodName.c_str(), iOptZeropad, sOptFilterGenerationName.c_str(), sOptInterpName.c_str(), iOptPreinterpolationFactor, sOptBackprojectName.c_str(), optTrace); + if (bOptVerbose) timerReconstruct.timerEndAndReport ("Time to reconstruct"); TimerCollectiveMPI timerReduce (mpiWorld.getComm()); ReduceImageMPI (mpiWorld, imLocal, imGlobal); - if (optVerbose) + if (bOptVerbose) timerReduce.timerEndAndReport ("Time to reduce image"); #else - projGlobal.reconstruct (*imGlobal, optFilterName.c_str(), optFilterParam, optFilterMethodName.c_str(), optZeroPad, optFrequencyFilterName.c_str(), optInterpName.c_str(), optPreinterpolationFactor, optBackprojName.c_str(), optTrace); + projGlobal.reconstruct (*imGlobal, sOptFilterName.c_str(), dOptFilterParam, sOptFilterMethodName.c_str(), iOptZeropad, sOptFilterGenerationName.c_str(), sOptInterpName.c_str(), iOptPreinterpolationFactor, sOptBackprojectName.c_str(), optTrace); #endif #ifdef HAVE_MPI if (mpiWorld.getRank() == 0) #endif { - double calcTime = timerProgram.timerEnd(); + double dCalcTime = timerProgram.timerEnd(); imGlobal->labelAdd (projGlobal.getLabel()); - imGlobal->labelAdd (Array2dFileLabel::L_HISTORY, remark.c_str(), calcTime); - imGlobal->fileWrite (filenameImage); - if (optVerbose) - cout << "Run time: " << calcTime << " seconds" << endl; + imGlobal->labelAdd (Array2dFileLabel::L_HISTORY, sRemark.c_str(), dCalcTime); + imGlobal->fileWrite (pszFilenameImage); + if (bOptVerbose) + cout << "Run time: " << dCalcTime << " seconds" << endl; } #ifdef HAVE_MPI MPI::Finalize(); @@ -355,7 +355,7 @@ pjrec_main (int argc, char * argv[]) ////////////////////////////////////////////////////////////////////////////////////// #ifdef HAVE_MPI -static void ScatterProjectionsMPI (MPIWorld& mpiWorld, Projections& projGlobal, Projections& projLocal, const int optDebug) +static void ScatterProjectionsMPI (MPIWorld& mpiWorld, Projections& projGlobal, Projections& projLocal, const bool bOptDebug) { if (mpiWorld.getRank() == 0) { for (int iProc = 0; iProc < mpiWorld.getNumProcessors(); iProc++) {