r182: *** empty log message ***

diff --git a/TODO b/TODO
index 006154230816e56c29446a50d375cb70a2caf967..46eec6e6d08c4d4e24449a2d50f8b34ae4d7d9ed 100644 (file)
--- a/TODO
+++ b/TODO
@@ -2,9 +2,9 @@ Fix BUGS -- see BUGS file
 
 Convert pol to C++ object.
 
 
 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
 
 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 941d905ea76512f78814763a11ce62bd5de2b02a..3a0e1c799aad5d156e402653b0050b933890d51a 100644 (file)
--- 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
 **
 **  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
 **
 **  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_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;
 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"},
 
 const char* SignalFilter::s_aszFilterName[] = {
   {"abs_bandlimit"},

-  {"abs_sinc"},

   {"abs_hamming"},
   {"abs_cosine"},
   {"abs_hamming"},
   {"abs_cosine"},

+  {"abs_sinc"},

   {"shepp"},
   {"bandlimit"},
   {"sinc"},
   {"shepp"},
   {"bandlimit"},
   {"sinc"},


@@ -56,9 +56,9 @@ const char* SignalFilter::s_aszFilterName[] = {
 
 const char* SignalFilter::s_aszFilterTitle[] = {
   {"Abs(w) * Bandlimit"},
 
 const char* SignalFilter::s_aszFilterTitle[] = {
   {"Abs(w) * Bandlimit"},

-  {"Abs(w) * Sinc"},

   {"Abs(w) * Hamming"},
   {"Abs(w) * Cosine"},
   {"Abs(w) * Hamming"},
   {"Abs(w) * Cosine"},

+  {"Abs(w) * Sinc"},

   {"Shepp"},
   {"Bandlimit"},
   {"Sinc"},
   {"Shepp"},
   {"Bandlimit"},
   {"Sinc"},


@@ -401,31 +401,31 @@ SignalFilter::frequencyResponse (int filterID, double bw, double u, double param
 
   switch (filterID) {
   case FILTER_BANDLIMIT:
 
   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:
       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:
       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:
       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:
       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);
       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:
     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:
       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));
       q = 0;
     else
       q = au * (param + (1 - param) * cos(TWOPI * u / bw));

diff --git a/libctsim/imagefile.cpp b/libctsim/imagefile.cpp
index 1b1660dc0375653a4665d990612a11fc89e76140..3f7c0beff717a84bafe3f588816dfc50d29bf0db 100644 (file)
--- 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
 **
 **  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
 **
 **  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);
 
     }
     mean /= (nx * ny);
 

-    static const int nbin = 256;
+    static const int nbin = 1024;

     int hist[ nbin ] = {0};
     double spread = max - min;
     mode = 0;
     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;
     }
 
     mode = (max_binindex * spread / (nbin - 1)) + min;

-      
-    median = 0.;
+
+    int nPixels = nx * ny;
+    slist<double> vecImage;
+    for (int ix = 0; ix < nx; ix++)
+      for (int iy = 0; iy < ny; iy++)
+       vecImage.push_front (v[ix][iy]);
+    vecImage.sort();
+    slist<double>::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 0ccbc09f47e7782559a52c8ae22cec859c78d2c2..4cb6f694fe54de5e67e97cbe5d85573c49d57396 100644 (file)
--- 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
 **
 **  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
 **
 **  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)
 {
 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();
   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;
       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();
     }
   }
       *theApp->getLog() << os.str().c_str();
     }
   }