From a5139613f6b92aeb8db9daa87f021c78d17d82a5 Mon Sep 17 00:00:00 2001 From: "Kevin M. Rosenberg" Date: Tue, 22 Aug 2000 16:49:56 +0000 Subject: [PATCH] r182: *** empty log message *** --- TODO | 4 ++-- libctsim/filter.cpp | 26 +++++++++++++------------- libctsim/imagefile.cpp | 22 ++++++++++++++++++---- src/views.cpp | 10 +++++++--- 4 files changed, 40 insertions(+), 22 deletions(-) diff --git a/TODO b/TODO index 0061542..46eec6e 100644 --- a/TODO +++ b/TODO @@ -2,9 +2,9 @@ Fix BUGS -- see BUGS file Convert pol to C++ object. -Add loading / saving of arbitrary phantom objections in ASCII format. - Adding animation of X-ray projection collection Add animation of filtering & backprojections +FFT filtering requires FFTW library. Add a public domain FFT routine. + diff --git a/libctsim/filter.cpp b/libctsim/filter.cpp index 941d905..3a0e1c7 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.26 2000/08/22 07:02:48 kevin Exp $ +** $Id: filter.cpp,v 1.27 2000/08/22 16:49:56 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 @@ -31,9 +31,9 @@ int SignalFilter::N_INTEGRAL=500; //static member const int SignalFilter::FILTER_INVALID = -1 ; const int SignalFilter::FILTER_ABS_BANDLIMIT = 0; // filter times |x = | -const int SignalFilter::FILTER_ABS_SINC = 1; -const int SignalFilter::FILTER_ABS_G_HAMMING = 2; -const int SignalFilter::FILTER_ABS_COSINE = 3; +const int SignalFilter::FILTER_ABS_G_HAMMING = 1; +const int SignalFilter::FILTER_ABS_COSINE = 2; +const int SignalFilter::FILTER_ABS_SINC = 3; const int SignalFilter::FILTER_SHEPP = 4; const int SignalFilter::FILTER_BANDLIMIT = 5; const int SignalFilter::FILTER_SINC = 6; @@ -43,9 +43,9 @@ const int SignalFilter::FILTER_TRIANGLE = 9; const char* SignalFilter::s_aszFilterName[] = { {"abs_bandlimit"}, - {"abs_sinc"}, {"abs_hamming"}, {"abs_cosine"}, + {"abs_sinc"}, {"shepp"}, {"bandlimit"}, {"sinc"}, @@ -56,9 +56,9 @@ const char* SignalFilter::s_aszFilterName[] = { const char* SignalFilter::s_aszFilterTitle[] = { {"Abs(w) * Bandlimit"}, - {"Abs(w) * Sinc"}, {"Abs(w) * Hamming"}, {"Abs(w) * Cosine"}, + {"Abs(w) * Sinc"}, {"Shepp"}, {"Bandlimit"}, {"Sinc"}, @@ -401,31 +401,31 @@ SignalFilter::frequencyResponse (int filterID, double bw, double u, double param switch (filterID) { case FILTER_BANDLIMIT: - if (fabs(au) > fabs(bw / 2) + F_EPSILON) + if (fabs(au) >= fabs(bw / 2) + F_EPSILON) q = 0.; else q = 1; break; case FILTER_ABS_BANDLIMIT: - if (fabs(au) > fabs(bw / 2) + F_EPSILON) + if (fabs(au) >= fabs(bw / 2) + F_EPSILON) q = 0.; else q = au; break; case FILTER_TRIANGLE: - if (fabs(au) > fabs(bw / 2) + F_EPSILON) + if (fabs(au) >= fabs(bw / 2) + F_EPSILON) q = 0; else q = 1 - au / bw; break; case FILTER_COSINE: - if (fabs(au) > fabs(bw / 2) + F_EPSILON) + if (fabs(au) >= fabs(bw / 2) + F_EPSILON) q = 0; else q = cos(PI * u / bw); break; case FILTER_ABS_COSINE: - if (fabs(au) > fabs(bw / 2) + F_EPSILON) + 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 (fabs(au) > fabs(bw / 2) + F_EPSILON) + 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 (fabs(au) > fabs(bw / 2) + F_EPSILON) + if (fabs(au) >= fabs(bw / 2) + F_EPSILON) q = 0; else q = au * (param + (1 - param) * cos(TWOPI * u / bw)); diff --git a/libctsim/imagefile.cpp b/libctsim/imagefile.cpp index 1b1660d..3f7c0be 100644 --- a/libctsim/imagefile.cpp +++ b/libctsim/imagefile.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: imagefile.cpp,v 1.9 2000/07/15 08:36:13 kevin Exp $ +** $Id: imagefile.cpp,v 1.10 2000/08/22 16:49:56 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 @@ -219,7 +219,7 @@ ImageFile::statistics (double& min, double& max, double& mean, double& mode, dou } mean /= (nx * ny); - static const int nbin = 256; + static const int nbin = 1024; int hist[ nbin ] = {0}; double spread = max - min; mode = 0; @@ -244,8 +244,22 @@ ImageFile::statistics (double& min, double& max, double& mean, double& mode, dou } mode = (max_binindex * spread / (nbin - 1)) + min; - - median = 0.; + + int nPixels = nx * ny; + slist vecImage; + for (int ix = 0; ix < nx; ix++) + for (int iy = 0; iy < ny; iy++) + vecImage.push_front (v[ix][iy]); + vecImage.sort(); + slist::const_iterator iter = vecImage.begin(); + for (int i = 0; i < (nPixels / 2) - 1; i++) + iter++; // Advance iterator to (nPixels / 2) - 1; + + if (nPixels % 2) { // Odd + iter++; + median = *iter; + } else // Even + median = (*iter++ + *iter) / 2; } diff --git a/src/views.cpp b/src/views.cpp index 0ccbc09..4cb6f69 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.14 2000/08/22 07:02:48 kevin Exp $ +** $Id: views.cpp,v 1.15 2000/08/22 16:49:56 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,11 @@ ProjectionFileView::OnProperties (wxCommandEvent& event) void ProjectionFileView::OnReconstruct (wxCommandEvent& event) { - 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); +#if HAVE_FFTW + DialogGetReconstructionParameters dialogReconstruction (m_frame, 256, 256, SignalFilter::FILTER_ABS_BANDLIMIT, 1., ProcessSignal::FILTER_METHOD_RFFTW, ProcessSignal::FILTER_GENERATION_INVERSE_FOURIER, 1, Backprojector::INTERP_LINEAR, 1, Backprojector::BPROJ_IDIFF3); +#else + DialogGetReconstructionParameters dialogReconstruction (m_frame, 256, 256, SignalFilter::FILTER_ABS_BANDLIMIT, 1., ProcessSignal::FILTER_METHOD_CONVOLUTION, ProcessSignal::FILTER_GENERATION_DIRECT, 1, Backprojector::INTERP_LINEAR, 1, Backprojector::BPROJ_IDIFF3); +#endif int retVal = dialogReconstruction.ShowModal(); if (retVal == wxID_OK) { int xSize = dialogReconstruction.getXSize(); @@ -647,7 +651,7 @@ ProjectionFileView::OnReconstruct (wxCommandEvent& event) pReconDoc->Modify(true); pReconDoc->UpdateAllViews(this); ostringstream os; - os << "Reconstruct " << rProj.getFilename() << ": xSize=" << xSize << ", ySize=" << ySize << ", Filter=" << optFilterName.c_str() << ", FilterParam=" << optFilterParam << ", FilterMethod=" << optFilterMethodName.c_str() << ", Zeropad=" << optZeropad << ", Interpolation=" << optInterpName.c_str() << ", InterpolationParam=" << optInterpParam << ", Backprojection=" << optBackprojectName.c_str() << "\n"; + os << "Reconstruct " << rProj.getFilename() << ": xSize=" << xSize << ", ySize=" << ySize << ", Filter=" << optFilterName.c_str() << ", FilterParam=" << optFilterParam << ", FilterMethod=" << optFilterMethodName.c_str() << ", FilterGeneration=" << optFilterGenerationName.c_str() << ", Zeropad=" << optZeropad << ", Interpolation=" << optInterpName.c_str() << ", InterpolationParam=" << optInterpParam << ", Backprojection=" << optBackprojectName.c_str() << "\n"; *theApp->getLog() << os.str().c_str(); } } -- 2.34.1