projects / ctsim.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
Update copyright date; remove old CVS keyword
[ctsim.git] / libctsim / fourier.cpp
diff --git a/libctsim/fourier.cpp b/libctsim/fourier.cpp
index f7357cdf5363272c496e0634e16b2794b9799deb..d85bfbe5b97d8c332a582ed94eea94469aaca171 100644 (file)
--- a/libctsim/fourier.cpp
+++ b/libctsim/fourier.cpp
@@ -7,9 +7,7 @@
 **   Date Started:  Dec 2000
 **
 **  This is part of the CTSim program
-**  Copyright (c) 1983-2001 Kevin Rosenberg
-**
-**  $Id: fourier.cpp,v 1.5 2001/03/18 18:08:25 kevin Exp $
+**  Copyright (c) 1983-2009 Kevin Rosenberg
 **
 **  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
@@ -99,181 +97,83 @@ Fourier::shuffleNaturalToFourierOrder (ImageFile& im)
   delete [] pRow;
 }
 
-
-// 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
-Fourier::shuffleNaturalToFourierOrder (double* pdVector, const int n)
+#ifdef HAVE_FFTW
+void Fourier::shuffleNaturalToFourierOrder (fftw_complex* pVector, const int n)
 {
-  double* pdTemp = new double [n];
+  fftw_complex* pTemp = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * n));
   int i;
-  if (isOdd(n)) { // Odd
-    int iHalfN = (n - 1) / 2;
 
-    pdTemp[0] = pdVector[iHalfN];
-    for (i = 0; i < iHalfN; i++)
-      pdTemp[i + 1] = pdVector[i + 1 + iHalfN];
-    for (i = 0; i < iHalfN; i++)
-      pdTemp[i + iHalfN + 1] = pdVector[i];
-  } else {     // Even
-    int iHalfN = n / 2;
-    pdTemp[0] = pdVector[iHalfN];
-    for (i = 0; i < iHalfN - 1; i++)
-      pdTemp[i + 1] = pdVector[i + iHalfN + 1];
-    for (i = 0; i < iHalfN; i++)
-      pdTemp[i + iHalfN] = pdVector[i];
-  }
-
-  for (i = 0; i < n; i++)
-    pdVector[i] = pdTemp[i];
-  delete pdTemp;
-}
-
-void
-Fourier::shuffleNaturalToFourierOrder (std::complex<double>* pdVector, const int n)
-{
-  std::complex<double>* pdTemp = new std::complex<double> [n];
-  int i;
   if (isOdd(n)) { // Odd
     int iHalfN = (n - 1) / 2;
 
-    pdTemp[0] = pdVector[iHalfN];
-    for (i = 0; i < iHalfN; i++)
-      pdTemp[i + 1] = pdVector[i + 1 + iHalfN];
-    for (i = 0; i < iHalfN; i++)
-      pdTemp[i + iHalfN + 1] = pdVector[i];
+    pTemp[0][0] = pVector[iHalfN][0];
+    pTemp[0][1] = pVector[iHalfN][1];
+    for (i = 0; i < iHalfN; i++) {
+      pTemp[i + 1][0] = pVector[i + 1 + iHalfN][0];
+      pTemp[i + 1][1] = pVector[i + 1 + iHalfN][1];
+    }
+    for (i = 0; i < iHalfN; i++) {
+      pTemp[i + iHalfN + 1][0] = pVector[i][0];
+      pTemp[i + iHalfN + 1][1] = pVector[i][1];
+    }
   } else {     // Even
     int iHalfN = n / 2;
-    pdTemp[0] = pdVector[iHalfN];
-    for (i = 0; i < iHalfN - 1; i++)
-      pdTemp[i + 1] = pdVector[i + iHalfN + 1];
-    for (i = 0; i < iHalfN; i++)
-      pdTemp[i + iHalfN] = pdVector[i];
+    pTemp[0][0] = pVector[iHalfN][0];
+    pTemp[0][1] = pVector[iHalfN][1];
+    for (i = 0; i < iHalfN - 1; i++) {
+      pTemp[i + 1][0] = pVector[i + iHalfN + 1][0];
+      pTemp[i + 1][1] = pVector[i + iHalfN + 1][1];
+    }
+    for (i = 0; i < iHalfN; i++) {
+      pTemp[i + iHalfN][0] = pVector[i][0];
+      pTemp[i + iHalfN][1] = pVector[i][1];
+    }
   }
 
-  for (i = 0; i < n; i++)
-    pdVector[i] = pdTemp[i];
-  delete [] pdTemp;
-}
-
-
-void
-Fourier::shuffleNaturalToFourierOrder (float* pdVector, const int n)
-{
-  float* pdTemp = new float [n];
-  int i;
-  if (isOdd (n)) { // Odd
-    int iHalfN = (n - 1) / 2;
-
-    pdTemp[0] = pdVector[iHalfN];
-    for (i = 0; i < iHalfN; i++)
-      pdTemp[i + 1] = pdVector[i + 1 + iHalfN];
-    for (i = 0; i < iHalfN; i++)
-      pdTemp[i + iHalfN + 1] = pdVector[i];
-  } else {     // Even
-    int iHalfN = n / 2;
-    pdTemp[0] = pdVector[iHalfN];
-    for (i = 0; i < iHalfN - 1; i++)
-      pdTemp[i + 1] = pdVector[i + iHalfN + 1];
-    for (i = 0; i < iHalfN; i++)
-      pdTemp[i + iHalfN] = pdVector[i];
+  for (i = 0; i < n; i++) {
+    pVector[i][0] = pTemp[i][0];
+    pVector[i][1] = pTemp[i][1];
   }
-
-  for (i = 0; i < n; i++)
-    pdVector[i] = pdTemp[i];
-  delete pdTemp;
+  fftw_free(pTemp);
 }
 
-
-
-void
-Fourier::shuffleFourierToNaturalOrder (double* pdVector, const int n)
+void Fourier::shuffleFourierToNaturalOrder (fftw_complex* pVector, const int n)
 {
-  double* pdTemp = new double [n];
+  fftw_complex* pTemp = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * n));
   int i;
   if (isOdd(n)) { // Odd
     int iHalfN = (n - 1) / 2;
-    
-    pdTemp[iHalfN] = pdVector[0];
-    for (i = 0; i < iHalfN; i++)
-      pdTemp[i + 1 + iHalfN] = pdVector[i + 1];
-    for (i = 0; i < iHalfN; i++)
-      pdTemp[i] = pdVector[i + iHalfN + 1];
-  } else {     // Even
-    int iHalfN = n / 2;
-    pdTemp[iHalfN] = pdVector[0];
-    for (i = 0; i < iHalfN; i++)
-      pdTemp[i] = pdVector[i + iHalfN];
-    for (i = 0; i < iHalfN - 1; i++)
-      pdTemp[i + iHalfN + 1] = pdVector[i+1];
-  }
-  
-  for (i = 0; i < n; i++)
-    pdVector[i] = pdTemp[i];
-  delete pdTemp;
-}
-
 
-void
-Fourier::shuffleFourierToNaturalOrder (std::complex<double>* pdVector, const int n)
-{
-  std::complex<double>* pdTemp = new std::complex<double> [n];
-  int i;
-  if (isOdd(n)) { // Odd
-    int iHalfN = (n - 1) / 2;
-    
-    pdTemp[iHalfN] = pdVector[0];
-    for (i = 0; i < iHalfN; i++)
-      pdTemp[i + 1 + iHalfN] = pdVector[i + 1];
-    for (i = 0; i < iHalfN; i++)
-      pdTemp[i] = pdVector[i + iHalfN + 1];
+    pTemp[iHalfN][0] = pVector[0][0];
+    pTemp[iHalfN][1] = pVector[0][1];
+    for (i = 0; i < iHalfN; i++) {
+      pTemp[i + 1 + iHalfN][0] = pVector[i + 1][0];
+      pTemp[i + 1 + iHalfN][1] = pVector[i + 1][1];
+    }
+    for (i = 0; i < iHalfN; i++) {
+      pTemp[i][0] = pVector[i + iHalfN + 1][0];
+      pTemp[i][1] = pVector[i + iHalfN + 1][1];
+    }
   } else {     // Even
     int iHalfN = n / 2;
-    pdTemp[iHalfN] = pdVector[0];
-    for (i = 0; i < iHalfN; i++)
-      pdTemp[i] = pdVector[i + iHalfN];
-    for (i = 0; i < iHalfN - 1; i++)
-      pdTemp[i + iHalfN + 1] = pdVector[i+1];
+    pTemp[iHalfN][0] = pVector[0][0];
+    pTemp[iHalfN][1] = pVector[0][1];
+    for (i = 0; i < iHalfN; i++) {
+      pTemp[i][0] = pVector[i + iHalfN][0];
+      pTemp[i][1] = pVector[i + iHalfN][1];
+    }
+    for (i = 0; i < iHalfN - 1; i++) {
+      pTemp[i + iHalfN + 1][0] = pVector[i+1][0];
+      pTemp[i + iHalfN + 1][1] = pVector[i+1][1];
+    }
   }
-  
-  for (i = 0; i < n; i++)
-    pdVector[i] = pdTemp[i];
-  delete [] pdTemp;
-}
-
 
-
-
-void
-Fourier::shuffleFourierToNaturalOrder (float* pVector, const int n)
-{
-  float* pTemp = new float [n];
-  int i;
-  if (isOdd (n)) { // Odd
-    int iHalfN = (n - 1) / 2;
-    
-    pTemp[iHalfN] = pVector[0];
-    for (i = 0; i < iHalfN; i++)
-      pTemp[i + 1 + iHalfN] = pVector[i + 1];
-    for (i = 0; i < iHalfN; i++)
-      pTemp[i] = pVector[i + iHalfN + 1];
-  } else {     // Even
-    int iHalfN = n / 2;
-    pTemp[iHalfN] = pVector[0];
-    for (i = 0; i < iHalfN; i++)
-      pTemp[i] = pVector[i + iHalfN];
-    for (i = 0; i < iHalfN - 1; i++)
-      pTemp[i + iHalfN + 1] = pVector[i+1];
+  for (i = 0; i < n; i++) {
+    pVector[i][0] = pTemp[i][0];
+    pVector[i][1] = pTemp[i][1];
   }
-  
-  for (i = 0; i < n; i++)
-    pVector[i] = pTemp[i];
-  delete [] pTemp;
-}
 
+  fftw_free(pTemp);
+}
+#endif
 
Computed Tomography simulator and reconstructor