projects / ctsim.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
r7939: Initial working version of fftw3 conversion -- tested only with pjrec
[ctsim.git] / src / views.cpp
diff --git a/src/views.cpp b/src/views.cpp
index 2d7874d690f5364f46f81e0a1d832b708c87cc4c..1fe4bbfbbc2c71facb551f464b64bc7dec8a4b36 100644 (file)
--- a/src/views.cpp
+++ b/src/views.cpp
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (c) 1983-2001 Kevin Rosenberg
 **
-**  $Id: views.cpp,v 1.168 2003/01/29 07:30:49 kevin Exp $
+**  $Id$
 **
 **  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
@@ -1088,7 +1088,9 @@ void
 ImageFileView::OnDraw (wxDC* dc)
 {
   if (m_pBitmap && m_pBitmap->Ok()) {
-    *theApp->getLog() << "Drawing bitmap";
+#ifdef DEBUG    
+    *theApp->getLog() << "Drawing bitmap\n";
+#endif
     dc->DrawBitmap(*m_pBitmap, 0, 0, false);
   }
   
@@ -1144,7 +1146,9 @@ ImageFileView::OnUpdate (wxView *WXUNUSED(sender), wxObject *WXUNUSED(hint) )
       delete m_pBitmap;
       m_pBitmap = NULL;
     }
-    *theApp->getLog() << "Making new bitmap bitmap";
+#ifdef DEBUG
+    *theApp->getLog() << "Making new bitmap\n";
+#endif
     m_pBitmap = new wxBitmap (image);
     delete imageData;
     m_pCanvas->SetScrollbars(20, 20, nx/20, ny/20);
@@ -1556,19 +1560,19 @@ ImageFileView::OnPlotFFTRow (wxCommandEvent& event)
   int ny = rIF.ny();
   
   if (v != NULL && yCursor < ny) {
-    fftw_complex* pcIn = new fftw_complex [nx];
+    fftw_complex* pcIn = static_cast<fftw_complex*>(fftw_malloc (sizeof(fftw_complex) * nx));
     
     int i;
     for (i = 0; i < nx; i++) {
-      pcIn[i].re = v[i][yCursor];
+      pcIn[i][0] = v[i][yCursor];
       if (rIF.isComplex())
-        pcIn[i].im = vImag[i][yCursor];
+        pcIn[i][1] = vImag[i][yCursor];
       else
-        pcIn[i].im = 0;
+        pcIn[i][1] = 0;
     }
     
-    fftw_plan plan = fftw_create_plan (nx, FFTW_FORWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
-    fftw_one (plan, pcIn, NULL);
+    fftw_plan plan = fftw_plan_dft_1d (nx, pcIn, pcIn, FFTW_FORWARD, FFTW_ESTIMATE);
+    fftw_execute (plan);
     fftw_destroy_plan (plan);
     
     double* pX = new double [nx];
@@ -1577,9 +1581,9 @@ ImageFileView::OnPlotFFTRow (wxCommandEvent& event)
     double* pYMag = new double [nx];
     for (i = 0; i < nx; i++) {
       pX[i] = i;
-      pYReal[i] = pcIn[i].re / nx;
-      pYImag[i] = pcIn[i].im / nx;
-      pYMag[i] = ::sqrt (pcIn[i].re * pcIn[i].re + pcIn[i].im * pcIn[i].im);
+      pYReal[i] = pcIn[i][0] / nx;
+      pYImag[i] = pcIn[i][1] / nx;
+      pYMag[i] = ::sqrt (pcIn[i][0] * pcIn[i][0] + pcIn[i][1] * pcIn[i][1]);
     }
     Fourier::shuffleFourierToNaturalOrder (pYReal, nx);
     Fourier::shuffleFourierToNaturalOrder (pYImag, nx);
@@ -1624,7 +1628,7 @@ ImageFileView::OnPlotFFTRow (wxCommandEvent& event)
     delete pYReal;
     delete pYImag;
     delete pYMag;
-    delete [] pcIn;
+    fftw_free(pcIn);
     
     if (theApp->getAskDeleteNewDocs())
       pPlotDoc->Modify (true);
@@ -1658,7 +1662,7 @@ ImageFileView::OnPlotFFTCol (wxCommandEvent& event)
       pdTemp[i] = v[xCursor][i];
     Fourier::shuffleNaturalToFourierOrder (pdTemp, ny);
     for (i = 0; i < ny; i++) 
-      pcIn[i].re = pdTemp[i];
+      pcIn[i][0] = pdTemp[i];
     
     for (i = 0; i < ny; i++) {
       if (rIF.isComplex())
@@ -1668,10 +1672,10 @@ ImageFileView::OnPlotFFTCol (wxCommandEvent& event)
     }
     Fourier::shuffleNaturalToFourierOrder (pdTemp, ny);
     for (i = 0; i < ny; i++)
-      pcIn[i].im = pdTemp[i];
+      pcIn[i][1] = pdTemp[i];
     
-    fftw_plan plan = fftw_create_plan (ny, FFTW_BACKWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
-    fftw_one (plan, pcIn, NULL);
+    fftw_plan plan = fftw_plan_dft_1d (ny, pcIn, pcIn, FFTW_BACKWARD, FFTW_ESTIMATE);
+    fftw_execute (plan);
     fftw_destroy_plan (plan);
     
     double* pX = new double [ny];
@@ -1680,9 +1684,9 @@ ImageFileView::OnPlotFFTCol (wxCommandEvent& event)
     double* pYMag = new double [ny];
     for (i = 0; i < ny; i++) {
       pX[i] = i;
-      pYReal[i] = pcIn[i].re / ny;
-      pYImag[i] = pcIn[i].im / ny;
-      pYMag[i] = ::sqrt (pcIn[i].re * pcIn[i].re + pcIn[i].im * pcIn[i].im);
+      pYReal[i] = pcIn[i][0] / ny;
+      pYImag[i] = pcIn[i][1] / ny;
+      pYMag[i] = ::sqrt (pcIn[i][0] * pcIn[i][0] + pcIn[i][1] * pcIn[i][1]);
     }
     
     PlotFileDocument* pPlotDoc = theApp->newPlotDoc();
@@ -2285,7 +2289,6 @@ PhantomFileView::OnRasterize (wxCommandEvent& event)
     *theApp->getLog() << os.str().c_str() << "\n";
     pImageFile->labelAdd (os.str().c_str(), timer.timerEnd());
 
-    OnUpdate(this, NULL);
     pRasterDoc->UpdateAllViews(this);
     pRasterDoc->getView()->setInitialClientSize();
     pRasterDoc->Activate();
@@ -2518,7 +2521,7 @@ ProjectionFileView::ProjectionFileView()
   m_iDefaultFilterMethod = ProcessSignal::FILTER_METHOD_CONVOLUTION;
   m_iDefaultFilterGeneration = ProcessSignal::FILTER_GENERATION_DIRECT;
 #endif
-  m_iDefaultZeropad = 1;
+  m_iDefaultZeropad = 2;
   m_iDefaultBackprojector = Backprojector::BPROJ_IDIFF;
   m_iDefaultInterpolation = Backprojector::INTERP_LINEAR;
   m_iDefaultInterpParam = 1;
@@ -2527,7 +2530,7 @@ ProjectionFileView::ProjectionFileView()
   m_iDefaultPolarNX = 256;
   m_iDefaultPolarNY = 256;
   m_iDefaultPolarInterpolation = Projections::POLAR_INTERP_BILINEAR;
-  m_iDefaultPolarZeropad = 1;
+  m_iDefaultPolarZeropad = 2;
 }
 
 ProjectionFileView::~ProjectionFileView()
@@ -2594,7 +2597,6 @@ ProjectionFileView::OnConvertPolar (wxCommandEvent& event)
     wxString strInterpolation (dialogPolar.getInterpolationName());
     m_iDefaultPolarNX = dialogPolar.getXSize();
     m_iDefaultPolarNY = dialogPolar.getYSize();
-    ImageFileDocument* pPolarDoc = theApp->newImageDoc();
     ImageFile* pIF = new ImageFile (m_iDefaultPolarNX, m_iDefaultPolarNY);
     m_iDefaultPolarInterpolation = Projections::convertInterpNameToID (strInterpolation.c_str());
     
@@ -2604,7 +2606,7 @@ ProjectionFileView::OnConvertPolar (wxCommandEvent& event)
       return;
     }
     
-    pPolarDoc = theApp->newImageDoc ();
+    ImageFileDocument* pPolarDoc = theApp->newImageDoc();
     if (! pPolarDoc) {
       sys_error (ERR_SEVERE, "Unable to create image file");
       return;
@@ -3060,12 +3062,12 @@ ProjectionFileView::CreateChildFrame(wxDocument *doc, wxView *view)
   GetDocumentManager()->FileHistoryAddFilesToMenu(m_pFileMenu);
   GetDocumentManager()->FileHistoryUseMenu(m_pFileMenu);
   
-  wxMenu *convert_menu = new wxMenu;
-  convert_menu->Append (PJMENU_CONVERT_RECTANGULAR, "&Rectangular Image");
-  convert_menu->Append (PJMENU_CONVERT_POLAR, "&Polar Image...\tCtrl-L");
-  convert_menu->Append (PJMENU_CONVERT_FFT_POLAR, "FF&T->Polar Image...\tCtrl-T");
-  convert_menu->AppendSeparator();
-  convert_menu->Append (PJMENU_CONVERT_PARALLEL, "&Interpolate to Parallel");
+  m_pConvertMenu = new wxMenu;
+  m_pConvertMenu->Append (PJMENU_CONVERT_RECTANGULAR, "&Rectangular Image");
+  m_pConvertMenu->Append (PJMENU_CONVERT_POLAR, "&Polar Image...\tCtrl-L");
+  m_pConvertMenu->Append (PJMENU_CONVERT_FFT_POLAR, "FF&T->Polar Image...\tCtrl-T");
+  m_pConvertMenu->AppendSeparator();
+  m_pConvertMenu->Append (PJMENU_CONVERT_PARALLEL, "&Interpolate to Parallel");
   
   //  wxMenu* filter_menu = new wxMenu;
   //  filter_menu->Append (PJMENU_ARTIFACT_REDUCTION, "&Artifact Reduction");
@@ -3074,11 +3076,10 @@ ProjectionFileView::CreateChildFrame(wxDocument *doc, wxView *view)
   analyze_menu->Append (PJMENU_PLOT_HISTOGRAM, "&Plot Histogram");
   analyze_menu->Append (PJMENU_PLOT_TTHETA_SAMPLING, "Plot T-T&heta Sampling...\tCtrl-H");
 
-  wxMenu *reconstruct_menu = new wxMenu;
-  reconstruct_menu->Append (PJMENU_RECONSTRUCT_FBP, "&Filtered Backprojection...\tCtrl-R", "Reconstruct image using filtered backprojection");
-  reconstruct_menu->Append (PJMENU_RECONSTRUCT_FBP_REBIN, "Filtered &Backprojection (Rebin to Parallel)...\tCtrl-B", "Reconstruct image using filtered backprojection");
-  // still buggy
-  //   reconstruct_menu->Append (PJMENU_RECONSTRUCT_FOURIER, "&Fourier...\tCtrl-E", "Reconstruct image using inverse Fourier");
+  m_pReconstructMenu = new wxMenu;
+  m_pReconstructMenu->Append (PJMENU_RECONSTRUCT_FBP, "&Filtered Backprojection...\tCtrl-R", "Reconstruct image using filtered backprojection");
+  m_pReconstructMenu->Append (PJMENU_RECONSTRUCT_FBP_REBIN, "Filtered &Backprojection (Rebin to Parallel)...\tCtrl-B", "Reconstruct image using filtered backprojection");
+  m_pReconstructMenu->Append (PJMENU_RECONSTRUCT_FOURIER, "&Inverse Fourier...\tCtrl-E", "Direct inverse Fourier");
   
   wxMenu *help_menu = new wxMenu;
   help_menu->Append(MAINMENU_HELP_CONTENTS, "&Contents\tF1");
@@ -3089,10 +3090,10 @@ ProjectionFileView::CreateChildFrame(wxDocument *doc, wxView *view)
   wxMenuBar *menu_bar = new wxMenuBar;
   
   menu_bar->Append (m_pFileMenu, "&File");
-  menu_bar->Append (convert_menu, "&Convert");
+  menu_bar->Append (m_pConvertMenu, "&Convert");
   //  menu_bar->Append (filter_menu, "Fi&lter");
   menu_bar->Append (analyze_menu, "&Analyze");
-  menu_bar->Append (reconstruct_menu, "&Reconstruct");
+  menu_bar->Append (m_pReconstructMenu, "&Reconstruct");
   menu_bar->Append (help_menu, "&Help");
   
   subframe->SetMenuBar(menu_bar);  
@@ -3155,6 +3156,14 @@ ProjectionFileView::OnUpdate (wxView *WXUNUSED(sender), wxObject *WXUNUSED(hint)
   const Projections& rProj = GetDocument()->getProjections();
   const int nDet = rProj.nDet();
   const int nView = rProj.nView();
+  if (rProj.geometry() == Scanner::GEOMETRY_PARALLEL) { 
+    m_pReconstructMenu->Enable (PJMENU_RECONSTRUCT_FBP_REBIN, false);
+    m_pConvertMenu->Enable (PJMENU_CONVERT_PARALLEL, false);
+  } else {
+    m_pReconstructMenu->Enable (PJMENU_RECONSTRUCT_FBP_REBIN, true);
+    m_pConvertMenu->Enable (PJMENU_CONVERT_PARALLEL, true);
+  }
+
   if (nDet != 0 && nView != 0) {
     const DetectorArray& detarray = rProj.getDetectorArray(0);
     const DetectorValue* detval = detarray.detValues();
Computed Tomography simulator and reconstructor