From 2a39ee3b125e3e2e68bbba2ac15a65039456ff7e Mon Sep 17 00:00:00 2001 From: "Kevin M. Rosenberg" Date: Tue, 22 Aug 2000 07:02:48 +0000 Subject: [PATCH] r181: *** empty log message *** --- ChangeLog | 10 +- include/procsignal.h | 8 +- libctgraphics/ezplot.cpp | 4 +- libctsim/filter.cpp | 16 +-- libctsim/procsignal.cpp | 227 ++++++++++++++++++++++++++++++--------- libctsim/projections.cpp | 5 +- src/dialogs.cpp | 15 ++- src/dialogs.h | 7 +- src/views.cpp | 7 +- tools/pjrec.cpp | 8 +- 10 files changed, 222 insertions(+), 85 deletions(-) diff --git a/ChangeLog b/ChangeLog index 72fcec0..50bb4f2 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,10 +1,12 @@ -2.0.0-b9 - 8/15/00 +2.0.0-b9 - 8/22/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 Filter-Generation option to reconstruction - Fixed compilation for non-SGP architectures - Decomposed SignalFilter class into ProcessSignal and SignalFilter classes + Fixed compilation for non-SGP architectures + Decomposed SignalFilter class into ProcessSignal and SignalFilter classes + Added Filter-Generation option to reconstruction to allow direct or + inverse_fourier construction of filters + 2.0.0-b8 - 8/1/00 Added line color support to SGP Fixed lineAbs bug diff --git a/include/procsignal.h b/include/procsignal.h index 1b08aea..8995dd6 100644 --- a/include/procsignal.h +++ b/include/procsignal.h @@ -9,7 +9,7 @@ ** 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 $ +** $Id: procsignal.h,v 1.2 2000/08/22 07:02:48 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 @@ -58,9 +58,9 @@ class ProcessSignal { 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 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, const int iTraceLevel = TRACE_NONE); - 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 (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, const int iTraceLevel = TRACE_NONE); ~ProcessSignal(); @@ -150,7 +150,7 @@ class ProcessSignal { 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); + 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, const int iTraceLevel); // transforms that use precalculated trig tables, therefore don't // require number of data points (n) as an argument diff --git a/libctgraphics/ezplot.cpp b/libctgraphics/ezplot.cpp index 348dd3f..7038132 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.10 2000/08/19 22:59:06 kevin Exp $ +** $Id: ezplot.cpp,v 1.11 2000/08/22 07:02:48 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 @@ -827,7 +827,7 @@ int EZPlot::axis_scale (double min, double max, int nint, double *minp, double *maxp, int *nintp) { if (min >= max || nint < 1) { - sys_error (ERR_WARNING, "Invalid params: min=%lf, min=%lf, num intervals=%d [axis_scale]", min, max, nint); + sys_error (ERR_WARNING, "Invalid params: min=%lf, max=%lf, num intervals=%d [axis_scale]", min, max, nint); return (FALSE); } diff --git a/libctsim/filter.cpp b/libctsim/filter.cpp index a88fad5..941d905 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.25 2000/08/19 22:59:06 kevin Exp $ +** $Id: filter.cpp,v 1.26 2000/08/22 07:02:48 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 @@ -401,31 +401,31 @@ SignalFilter::frequencyResponse (int filterID, double bw, double u, double param switch (filterID) { case FILTER_BANDLIMIT: - if (au >= bw / 2) + if (fabs(au) > fabs(bw / 2) + F_EPSILON) q = 0.; else q = 1; break; case FILTER_ABS_BANDLIMIT: - if (au >= bw / 2) + if (fabs(au) > fabs(bw / 2) + F_EPSILON) q = 0.; else q = au; break; case FILTER_TRIANGLE: - if (au >= bw) + if (fabs(au) > fabs(bw / 2) + F_EPSILON) q = 0; else q = 1 - au / bw; break; case FILTER_COSINE: - if (au >= bw / 2) + if (fabs(au) > fabs(bw / 2) + F_EPSILON) q = 0; else q = cos(PI * u / bw); break; case FILTER_ABS_COSINE: - if (au >= bw / 2) + if (fabs(au) > fabs(bw / 2) + F_EPSILON) q = 0; else q = au * cos(PI * u / bw); @@ -437,13 +437,13 @@ SignalFilter::frequencyResponse (int filterID, double bw, double u, double param q = au * bw * sinc (PI * bw * u, 1.); break; case FILTER_G_HAMMING: - if (au >= bw / 2) + if (fabs(au) > fabs(bw / 2) + F_EPSILON) q = 0; else q = param + (1 - param) * cos (TWOPI * u / bw); break; case FILTER_ABS_G_HAMMING: - if (au >= bw / 2) + if (fabs(au) > fabs(bw / 2) + F_EPSILON) q = 0; else q = au * (param + (1 - param) * cos(TWOPI * u / bw)); diff --git a/libctsim/procsignal.cpp b/libctsim/procsignal.cpp index 9fdbb68..70d759c 100644 --- a/libctsim/procsignal.cpp +++ b/libctsim/procsignal.cpp @@ -9,7 +9,7 @@ ** 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 $ +** $Id: procsignal.cpp,v 1.2 2000/08/22 07:02:48 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 @@ -49,8 +49,8 @@ const char* ProcessSignal::s_aszFilterMethodName[] = { }; const char* ProcessSignal::s_aszFilterMethodTitle[] = { {"Convolution"}, - {"Direct Fourier"}, - {"Fouier Trigometric Table Lookout"}, + {"Fourier"}, + {"Fouier Trigometric Table"}, {"FFT"}, #if HAVE_FFTW {"FFTW"}, @@ -77,7 +77,7 @@ const int ProcessSignal::s_iFilterGenerationCount = sizeof(s_aszFilterGeneration // 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) +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, int iTraceLevel = TRACE_NONE) : m_adFourierCosTable(NULL), m_adFourierSinTable(NULL), m_adFilter(NULL), m_fail(false) { m_idFilterMethod = convertFilterMethodNameToID (szFilterMethodName); @@ -109,12 +109,12 @@ ProcessSignal::ProcessSignal (const char* szFilterName, const char* szFilterMeth return; } - init (m_idFilter, m_idFilterMethod, dBandwidth, dSignalIncrement, nSignalPoints, dFilterParam, m_idDomain, m_idFilterGeneration, iZeropad, iPreinterpolationFactor); + init (m_idFilter, m_idFilterMethod, dBandwidth, dSignalIncrement, nSignalPoints, dFilterParam, m_idDomain, m_idFilterGeneration, iZeropad, iPreinterpolationFactor, iTraceLevel); } 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) +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, int iTraceLevel) { m_idFilter = idFilter; m_idDomain = idDomain; @@ -124,7 +124,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw m_fail = true; return; } - m_traceLevel = TRACE_NONE; + m_traceLevel = iTraceLevel; m_nameFilterMethod = convertFilterMethodIDToName (m_idFilterMethod); m_nameFilterGeneration = convertFilterGenerationIDToName (m_idFilterGeneration); m_dBandwidth = dBandwidth; @@ -152,7 +152,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw if (! m_bFrequencyFiltering) { if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) { - m_nFilterPoints = 2 * m_nSignalPoints - 1; + m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1; m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1); m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1); m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1); @@ -160,63 +160,177 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw 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_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1; 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); + if (m_traceLevel >= TRACE_PLOT) { + SGPDriver sgpDriver ("Frequency Filter: Natural Order"); + SGP sgp (sgpDriver); + EZPlot ezplot (sgp); + + ezplot.ezset ("title Filter Response: Natural Order"); + ezplot.addCurve (adFrequencyFilter, m_nFilterPoints); + ezplot.plot(); + cio_put_str ("Press any key to continue"); + cio_kb_getc (); + } + shuffleNaturalToFourierOrder (adFrequencyFilter, m_nFilterPoints); - ProcessSignal::finiteFourierTransform (adFrequencyFilter, adInverseFilter, m_nFilterPoints, 1); - for (int i = 0; i < m_nFilterPoints; i++) - m_adFilter [i] = adInverseFilter[i]; + if (m_traceLevel >= TRACE_PLOT) { + SGPDriver sgpDriver ("Frequency Filter: Fourier Order"); + SGP sgp (sgpDriver); + EZPlot ezplot (sgp); + + ezplot.ezset ("title Filter Response: Fourier Order"); + ezplot.addCurve (adFrequencyFilter, m_nFilterPoints); + ezplot.plot(); + cio_put_str ("Press any key to continue"); + cio_kb_getc (); + } + ProcessSignal::finiteFourierTransform (adFrequencyFilter, m_adFilter, m_nFilterPoints, -1); + if (m_traceLevel >= TRACE_PLOT) { + SGPDriver sgpDriver ("Inverse Fourier Frequency: Fourier Order"); + SGP sgp (sgpDriver); + EZPlot ezplot (sgp); + + ezplot.ezset ("title Inverse Fourier Frequency: Fourier Order"); + ezplot.addCurve (m_adFilter, m_nFilterPoints); + ezplot.plot(); + cio_put_str ("Press any key to continue"); + cio_kb_getc (); + } + shuffleFourierToNaturalOrder (m_adFilter, m_nFilterPoints); + if (m_traceLevel >= TRACE_PLOT) { + SGPDriver sgpDriver ("Inverse Fourier Frequency: Natural Order"); + SGP sgp (sgpDriver); + EZPlot ezplot (sgp); + + ezplot.ezset ("title Inverse Fourier Frequency: Natural Order"); + ezplot.addCurve (m_adFilter, m_nFilterPoints); + ezplot.plot(); + cio_put_str ("Press any key to continue"); + cio_kb_getc (); + } + for (int i = 0; i < m_nFilterPoints; i++) { + m_adFilter[i] /= m_dSignalInc; + } } } // 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) + if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) { + // calculate number of filter points with zeropadding + m_nFilterPoints = m_nSignalPoints; + if (m_iZeropad > 0) { + double logBase2 = log(m_nFilterPoints) / 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; + } + m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor; - if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) { + if (m_nFilterPoints % 2) { // Odd 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 { // Even + m_dFilterMin = -1. / (2 * m_dSignalInc); + m_dFilterMax = 1. / (2 * m_dSignalInc); + m_dFilterInc = (m_dFilterMax - m_dFilterMin) / m_nFilterPoints; + m_dFilterMax -= m_dFilterInc; + } + + 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); + if (m_traceLevel >= TRACE_PLOT) { + SGPDriver sgpDriver ("Frequency Filter: Natural Order"); + SGP sgp (sgpDriver); + EZPlot ezplot (sgp); + + ezplot.ezset ("title Filter Filter: Natural Order"); + ezplot.addCurve (m_adFilter, m_nFilterPoints); + ezplot.plot(); + cio_put_str ("Press any key to continue"); + cio_kb_getc (); + } + shuffleNaturalToFourierOrder (m_adFilter, m_nFilterPoints); + if (m_traceLevel >= TRACE_PLOT) { + SGPDriver sgpDriver ("Frequency Filter: Fourier Order"); + SGP sgp (sgpDriver); + EZPlot ezplot (sgp); + + ezplot.ezset ("title Filter Filter: Fourier Order"); + ezplot.addCurve (m_adFilter, m_nFilterPoints); + ezplot.plot(); + cio_put_str ("Press any key to continue"); + cio_kb_getc (); + } } 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]; + // calculate number of filter points with zeropadding + int nSpatialPoints = 2 * (m_nSignalPoints - 1) + 1; + m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1); + m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1); + m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (nSpatialPoints - 1); + m_nFilterPoints = nSpatialPoints; + if (m_iZeropad > 0) { + double logBase2 = log(nSpatialPoints) / log(2); + int nextPowerOf2 = static_cast(floor(logBase2)); + if (logBase2 != floor(logBase2)) + nextPowerOf2++; + nextPowerOf2 += (m_iZeropad - 1); + m_nFilterPoints = 1 << nextPowerOf2; } - } + m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor; + if (m_traceLevel >= TRACE_TEXT) + cout << "nFilterPoints = " << m_nFilterPoints << endl; + double adSpatialFilter [m_nFilterPoints]; + SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, nSpatialPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL); + filter.copyFilterData (adSpatialFilter, 0, nSpatialPoints); + if (m_traceLevel >= TRACE_PLOT) { + SGPDriver sgpDriver ("Spatial Filter: Natural Order"); + SGP sgp (sgpDriver); + EZPlot ezplot (sgp); + + ezplot.ezset ("title Spatial Filter: Natural Order"); + ezplot.addCurve (adSpatialFilter, nSpatialPoints); + ezplot.plot(); + cio_put_str ("Press any key to continue"); + cio_kb_getc (); + } + for (int i = nSpatialPoints; i < m_nFilterPoints; i++) + adSpatialFilter[i] = 0; + m_adFilter = new double [m_nFilterPoints]; + complex acInverseFilter [m_nFilterPoints]; + finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, 1); + for (int i = 0; i < m_nFilterPoints; i++) + m_adFilter[i] = abs(acInverseFilter[i]) * m_dSignalInc; + if (m_traceLevel >= TRACE_PLOT) { + SGPDriver sgpDriver ("Spatial Filter: Inverse"); + SGP sgp (sgpDriver); + EZPlot ezplot (sgp); + + ezplot.ezset ("title Spatial Filter: Inverse"); + ezplot.addCurve (m_adFilter, m_nFilterPoints); + ezplot.plot(); + cio_put_str ("Press any key to continue"); + cio_kb_getc (); + } + } + } + // 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; @@ -262,6 +376,7 @@ ProcessSignal::~ProcessSignal (void) { delete [] m_adFourierSinTable; delete [] m_adFourierCosTable; + delete [] m_adFilter; #if HAVE_FFTW if (m_idFilterMethod == FILTER_METHOD_FFTW) { @@ -491,7 +606,7 @@ ProcessSignal::finiteFourierTransform (const double input[], double output[], co finiteFourierTransform (input, complexOutput, n, direction); for (int i = 0; i < n; i++) - output[i] = abs(complexOutput[n]); + output[i] = complexOutput[i].real(); } void @@ -670,13 +785,17 @@ ProcessSignal::shuffleNaturalToFourierOrder (double* pdVector, const int n) 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]; + for (int i = 0; i < iHalfN; i++) + pdTemp[i + 1] = pdVector[i + 1 + iHalfN]; + for (int i = 0; i < iHalfN; i++) + pdTemp[i + iHalfN + 1] = pdVector[i]; } else { // Even int iHalfN = n / 2; pdTemp[0] = pdVector[iHalfN]; + for (int i = 0; i < iHalfN; i++) + pdTemp[i + 1] = pdVector[i + iHalfN]; + for (int i = 0; i < iHalfN - 1; i++) + pdTemp[i + iHalfN + 1] = pdVector[i]; } for (int i = 0; i < n; i++) @@ -693,13 +812,17 @@ ProcessSignal::shuffleFourierToNaturalOrder (double* pdVector, const int n) 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]; + for (int i = 0; i < iHalfN; i++) + pdTemp[i + 1 + iHalfN] = pdVector[i + 1]; + for (int i = 0; i < iHalfN; i++) + pdTemp[i] = pdVector[i + iHalfN + 1]; } else { // Even int iHalfN = n / 2; pdTemp[iHalfN] = pdVector[0]; + for (int i = 0; i < iHalfN; i++) + pdTemp[i] = pdVector[i + iHalfN]; + for (int i = 0; i < iHalfN - 1; i++) + pdTemp[i + iHalfN + 1] = pdVector[i+1]; } for (int i = 0; i < n; i++) diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp index 813d7bc..d6b496d 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.20 2000/08/19 22:59:06 kevin Exp $ +** $Id: projections.cpp,v 1.21 2000/08/22 07:02:48 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 @@ -503,8 +503,7 @@ Projections::reconstruct (ImageFile& im, const char* const filterName, double fi #endif double filterBW = 1. / detInc; - ProcessSignal processSignal (filterName, filterMethodName, filterBW, m_detInc, m_nDet, filt_param, "spatial", filterGenerationName, zeropad, interpFactor); - processSignal.setTraceLevel(trace); + ProcessSignal processSignal (filterName, filterMethodName, filterBW, m_detInc, m_nDet, filt_param, "spatial", filterGenerationName, zeropad, interpFactor, trace); if (processSignal.fail()) { sys_error (ERR_SEVERE, "%s [Projections::reconstruct]", processSignal.failMessage().c_str()); diff --git a/src/dialogs.cpp b/src/dialogs.cpp index 79913b4..4b2738f 100644 --- a/src/dialogs.cpp +++ b/src/dialogs.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: dialogs.cpp,v 1.7 2000/07/23 01:49:03 kevin Exp $ +** $Id: dialogs.cpp,v 1.8 2000/08/22 07:02:48 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 @@ -398,7 +398,7 @@ DialogGetProjectionParameters::getGeometry (void) ///////////////////////////////////////////////////////////////////// -DialogGetReconstructionParameters::DialogGetReconstructionParameters (wxFrame* pParent, int iDefaultXSize = 0, int iDefaultYSize = 0, int iDefaultFilterID = SignalFilter::FILTER_ABS_BANDLIMIT, double dDefaultFilterParam = 1., int iDefaultFilterMethodID = SignalFilter::FILTER_METHOD_CONVOLUTION, int iDefaultZeropad = 3, int iDefaultInterpID = Backprojector::INTERP_LINEAR, int iDefaultInterpParam = 1, int iDefaultBackprojectID = Backprojector::BPROJ_IDIFF3) +DialogGetReconstructionParameters::DialogGetReconstructionParameters (wxFrame* pParent, int iDefaultXSize = 0, int iDefaultYSize = 0, int iDefaultFilterID = SignalFilter::FILTER_ABS_BANDLIMIT, double dDefaultFilterParam = 1., int iDefaultFilterMethodID = ProcessSignal::FILTER_METHOD_CONVOLUTION, int iDefaultFilterGenerationID = ProcessSignal::FILTER_GENERATION_INVALID, int iDefaultZeropad = 3, int iDefaultInterpID = Backprojector::INTERP_LINEAR, int iDefaultInterpParam = 1, int iDefaultBackprojectID = Backprojector::BPROJ_IDIFF3) : wxDialog (pParent, -1, "Set Reconstruction Parameters", wxDefaultPosition, wxDefaultSize, wxDEFAULT_DIALOG_STYLE | wxCAPTION) { wxBoxSizer* pTopSizer = new wxBoxSizer (wxVERTICAL); @@ -411,9 +411,12 @@ DialogGetReconstructionParameters::DialogGetReconstructionParameters (wxFrame* p m_pListBoxFilter->SetSelection (iDefaultFilterID); pTopSizer->Add (m_pListBoxFilter, 0, wxALL | wxALIGN_CENTER | wxEXPAND); - m_pListBoxFilterMethod = new StringValueAndTitleListBox (this, SignalFilter::getFilterMethodCount(), SignalFilter::getFilterMethodTitleArray(), SignalFilter::getFilterMethodNameArray()); + m_pListBoxFilterMethod = new StringValueAndTitleListBox (this, ProcessSignal::getFilterMethodCount(), ProcessSignal::getFilterMethodTitleArray(), ProcessSignal::getFilterMethodNameArray()); m_pListBoxFilterMethod->SetSelection (iDefaultFilterMethodID); pTopSizer->Add (m_pListBoxFilterMethod, 0, wxALL | wxALIGN_CENTER | wxEXPAND); + m_pListBoxFilterGeneration = new StringValueAndTitleListBox (this, ProcessSignal::getFilterGenerationCount(), ProcessSignal::getFilterGenerationTitleArray(), ProcessSignal::getFilterGenerationNameArray()); + m_pListBoxFilterGeneration->SetSelection (iDefaultFilterGenerationID); + pTopSizer->Add (m_pListBoxFilterGeneration, 0, wxALL | wxALIGN_CENTER | wxEXPAND); m_pListBoxBackproject = new StringValueAndTitleListBox (this, Backprojector::getBackprojectCount(), Backprojector::getBackprojectTitleArray(), Backprojector::getBackprojectNameArray()); m_pListBoxBackproject->SetSelection (iDefaultBackprojectID); @@ -552,3 +555,9 @@ DialogGetReconstructionParameters::getBackprojectName (void) { return m_pListBoxBackproject->getSelectionStringValue(); } + +const char* +DialogGetReconstructionParameters::getFilterGenerationName (void) +{ + return m_pListBoxFilterGeneration->getSelectionStringValue(); +} diff --git a/src/dialogs.h b/src/dialogs.h index 6dcc298..e4f890c 100644 --- a/src/dialogs.h +++ b/src/dialogs.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: dialogs.h,v 1.8 2000/07/28 08:28:08 kevin Exp $ +** $Id: dialogs.h,v 1.9 2000/08/22 07:02:48 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 @@ -33,6 +33,7 @@ #include #include "scanner.h" #include "phantom.h" +#include "procsignal.h" #include "filter.h" // CLASS StringValueAndTitleListBox @@ -136,7 +137,7 @@ class DialogGetProjectionParameters : public wxDialog class DialogGetReconstructionParameters : public wxDialog { public: - DialogGetReconstructionParameters (wxFrame* pParent, int iDefaultXSize = 0, int iDefaultYSize = 0, int iDefaultFilterID = SignalFilter::FILTER_ABS_BANDLIMIT, double dDefaultFilterParam = 1., int iDefaultFilterMethodID = SignalFilter::FILTER_METHOD_CONVOLUTION, int iDefaultZeropad = 3, int iDefaultInterpID = Backprojector::INTERP_LINEAR, int iDefaultInterpParam = 1, int iDefaultBackprojectID = Backprojector::BPROJ_IDIFF3); + DialogGetReconstructionParameters (wxFrame* pParent, int iDefaultXSize = 0, int iDefaultYSize = 0, int iDefaultFilterID = SignalFilter::FILTER_ABS_BANDLIMIT, double dDefaultFilterParam = 1., int iDefaultFilterMethodID = ProcessSignal::FILTER_METHOD_CONVOLUTION, int iDefaultFilterGeneration = ProcessSignal::FILTER_GENERATION_INVALID, int iDefaultZeropad = 3, int iDefaultInterpID = Backprojector::INTERP_LINEAR, int iDefaultInterpParam = 1, int iDefaultBackprojectID = Backprojector::BPROJ_IDIFF3); virtual ~DialogGetReconstructionParameters (void); unsigned int getXSize(void); @@ -145,6 +146,7 @@ class DialogGetReconstructionParameters : public wxDialog double getFilterParam(void); const char* getFilterMethodName(void); unsigned int getZeropad(void); + const char* getFilterGenerationName(void); const char* getInterpName(void); unsigned int getInterpParam(void); const char* getBackprojectName(void); @@ -158,6 +160,7 @@ class DialogGetReconstructionParameters : public wxDialog StringValueAndTitleListBox* m_pListBoxFilter; StringValueAndTitleListBox* m_pListBoxFilterMethod; + StringValueAndTitleListBox* m_pListBoxFilterGeneration; StringValueAndTitleListBox* m_pListBoxInterp; StringValueAndTitleListBox* m_pListBoxBackproject; diff --git a/src/views.cpp b/src/views.cpp index 1ca02eb..0ccbc09 100644 --- a/src/views.cpp +++ b/src/views.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: views.cpp,v 1.13 2000/08/02 18:06:00 kevin Exp $ +** $Id: views.cpp,v 1.14 2000/08/22 07:02:48 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 @@ -625,7 +625,7 @@ ProjectionFileView::OnProperties (wxCommandEvent& event) void ProjectionFileView::OnReconstruct (wxCommandEvent& event) { - DialogGetReconstructionParameters dialogReconstruction (m_frame, 256, 256, SignalFilter::FILTER_ABS_BANDLIMIT, 1., SignalFilter::FILTER_METHOD_CONVOLUTION, 3, Backprojector::INTERP_LINEAR, 1, Backprojector::BPROJ_IDIFF3); + DialogGetReconstructionParameters dialogReconstruction (m_frame, 256, 256, SignalFilter::FILTER_ABS_BANDLIMIT, 1., ProcessSignal::FILTER_METHOD_CONVOLUTION, ProcessSignal::FILTER_GENERATION_INVALID, 3, Backprojector::INTERP_LINEAR, 1, Backprojector::BPROJ_IDIFF3); int retVal = dialogReconstruction.ShowModal(); if (retVal == wxID_OK) { int xSize = dialogReconstruction.getXSize(); @@ -634,6 +634,7 @@ ProjectionFileView::OnReconstruct (wxCommandEvent& event) double optFilterParam = dialogReconstruction.getFilterParam(); wxString optFilterMethodName = dialogReconstruction.getFilterMethodName(); int optZeropad = dialogReconstruction.getZeropad(); + wxString optFilterGenerationName = dialogReconstruction.getFilterGenerationName(); wxString optInterpName = dialogReconstruction.getInterpName(); int optInterpParam = dialogReconstruction.getInterpParam(); wxString optBackprojectName = dialogReconstruction.getBackprojectName(); @@ -642,7 +643,7 @@ ProjectionFileView::OnReconstruct (wxCommandEvent& event) ImageFile& imageFile = pReconDoc->getImageFile(); const Projections& rProj = GetDocument()->getProjections(); imageFile.setArraySize (xSize, ySize); - rProj.reconstruct (imageFile, optFilterName.c_str(), optFilterParam, optFilterMethodName.c_str(), optZeropad, optInterpName.c_str(), optInterpParam, optBackprojectName.c_str(), TRACE_NONE); + rProj.reconstruct (imageFile, optFilterName.c_str(), optFilterParam, optFilterMethodName.c_str(), optZeropad, optFilterGenerationName.c_str(), optInterpName.c_str(), optInterpParam, optBackprojectName.c_str(), TRACE_NONE); pReconDoc->Modify(true); pReconDoc->UpdateAllViews(this); ostringstream os; diff --git a/tools/pjrec.cpp b/tools/pjrec.cpp index 1638dbe..10adbff 100644 --- a/tools/pjrec.cpp +++ b/tools/pjrec.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: pjrec.cpp,v 1.12 2000/08/19 22:59:06 kevin Exp $ +** $Id: pjrec.cpp,v 1.13 2000/08/22 07:02:48 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 @@ -49,7 +49,7 @@ static struct option my_options[] = {0, 0, 0, 0} }; -static const char* g_szIdStr = "$Id: pjrec.cpp,v 1.12 2000/08/19 22:59:06 kevin Exp $"; +static const char* g_szIdStr = "$Id: pjrec.cpp,v 1.13 2000/08/22 07:02:48 kevin Exp $"; void pjrec_usage (const char *program) @@ -72,12 +72,12 @@ pjrec_usage (const char *program) cout << " --filter Filter name" << endl; cout << " abs_bandlimit Abs * Bandlimiting (default)" << endl; cout << " abs_sinc Abs * Sinc" << endl; - cout << " abs_cos Abs * Cosine" << endl; + cout << " abs_cosine Abs * Cosine" << endl; cout << " abs_hamming Abs * Hamming" << endl; cout << " shepp Shepp-Logan" << endl; cout << " bandlimit Bandlimiting" << endl; cout << " sinc Sinc" << endl; - cout << " cos Cosine" << endl; + cout << " cosine Cosine" << endl; cout << " triangle Triangle" << endl; cout << " hamming Hamming" << endl; cout << " --filter-method Filter method before backprojections\n";; -- 2.34.1