Update to wx3.0, add SSE optimizations based on target_cpu, fix compile warnings

[ctsim.git] / libctsim / imagefile.cpp

diff --git a/libctsim/imagefile.cpp b/libctsim/imagefile.cpp
index 51a8dea82e94f2924671f0a9998119b055963990..eeda46d37d4300f4662d15ac2924fa88ebdaa311 100644 (file)
--- a/libctsim/imagefile.cpp
+++ b/libctsim/imagefile.cpp
@@ -565,36 +565,39 @@ ImageFile::sqrt (ImageFile& result) const
 
   bool bComplexOutput = result.isComplex();
   ImageFileArrayConst vLHS = getArray();
-  if (! bComplexOutput)   // check if should convert to complex output
-    for (unsigned int ix = 0; ix < m_nx; ix++)
-      for (unsigned int iy = 0; iy < m_ny; iy++)
+  if (! bComplexOutput) {  // check if should convert to complex output
+    for (unsigned int ix = 0; ix < m_nx; ix++) {
+      for (unsigned int iy = 0; iy < m_ny; iy++) {
         if (! bComplexOutput && vLHS[ix][iy] < 0) {
           result.convertRealToComplex();
           bComplexOutput = true;
           break;
         }
+      }
+    }
+  }
 
-        ImageFileArrayConst vLHSImag = getImaginaryArray();
-        ImageFileArray vResult = result.getArray();
-        ImageFileArray vResultImag = result.getImaginaryArray();
-
-        for (unsigned int ix = 0; ix < m_nx; ix++) {
-          for (unsigned int iy = 0; iy < m_ny; iy++) {
-            if (result.isComplex()) {
-              double dImag = 0;
-              if (isComplex())
-                dImag = vLHSImag[ix][iy];
-              std::complex<double> cLHS (vLHS[ix][iy], dImag);
-              std::complex<double> cResult = std::sqrt(cLHS);
-              vResult[ix][iy] = cResult.real();
-              vResultImag[ix][iy] = cResult.imag();
-            } else
-              vResult[ix][iy] = ::sqrt (vLHS[ix][iy]);
-          }
-        }
-
+  ImageFileArrayConst vLHSImag = getImaginaryArray();
+  ImageFileArray vResult = result.getArray();
+  ImageFileArray vResultImag = result.getImaginaryArray();
 
-        return true;
+  for (unsigned int ix = 0; ix < m_nx; ix++) {
+    for (unsigned int iy = 0; iy < m_ny; iy++) {
+      if (result.isComplex()) {
+        double dImag = 0;
+        if (isComplex())
+          dImag = vLHSImag[ix][iy];
+        std::complex<double> cLHS (vLHS[ix][iy], dImag);
+        std::complex<double> cResult = std::sqrt(cLHS);
+        vResult[ix][iy] = cResult.real();
+        vResultImag[ix][iy] = cResult.imag();
+      } else
+        vResult[ix][iy] = ::sqrt (vLHS[ix][iy]);
+    }
+  }
+  
+  
+  return true;
 }
 
 bool
@@ -1039,62 +1042,63 @@ ImageFile::fourier (ImageFile& result) const
     return false;
   }
 
-  if (! result.isComplex())
-    if (! result.convertRealToComplex ())
+  if (! result.isComplex()) {
+    if (! result.convertRealToComplex ()) {
       return false;
-
-    ImageFileArrayConst vLHS = getArray();
-    ImageFileArrayConst vLHSImag = getImaginaryArray();
-    ImageFileArray vRealResult = result.getArray();
-    ImageFileArray vImagResult = result.getImaginaryArray();
-
-    unsigned int ix, iy;
-
-    // alloc output matrix
-    CTSimComplex** complexOut = new CTSimComplex* [m_nx];
-    for (ix = 0; ix < m_nx; ix++)
-      complexOut[ix] = new CTSimComplex [m_ny];
-
-    // fourier each x column
-    CTSimComplex* pY = new CTSimComplex [m_ny];
-    for (ix = 0; ix < m_nx; ix++) {
-      for (iy = 0; iy < m_ny; iy++) {
-        double dImag = 0;
-        if (isComplex())
-          dImag = vLHSImag[ix][iy];
-        pY[iy] = std::complex<double>(vLHS[ix][iy], dImag);
-      }
-      ProcessSignal::finiteFourierTransform (pY, complexOut[ix], m_ny,  ProcessSignal::FORWARD);
     }
-    delete [] pY;
-
-    // fourier each y row
-    CTSimComplex* pX = new CTSimComplex [m_nx];
-    CTSimComplex* complexOutRow = new CTSimComplex [m_nx];
+  }
+  ImageFileArrayConst vLHS = getArray();
+  ImageFileArrayConst vLHSImag = getImaginaryArray();
+  ImageFileArray vRealResult = result.getArray();
+  ImageFileArray vImagResult = result.getImaginaryArray();
+  
+  unsigned int ix, iy;
+  
+  // alloc output matrix
+  CTSimComplex** complexOut = new CTSimComplex* [m_nx];
+  for (ix = 0; ix < m_nx; ix++)
+    complexOut[ix] = new CTSimComplex [m_ny];
+  
+  // fourier each x column
+  CTSimComplex* pY = new CTSimComplex [m_ny];
+  for (ix = 0; ix < m_nx; ix++) {
     for (iy = 0; iy < m_ny; iy++) {
-      for (ix = 0; ix < m_nx; ix++)
-        pX[ix] = complexOut[ix][iy];
-      ProcessSignal::finiteFourierTransform (pX, complexOutRow, m_nx, ProcessSignal::FORWARD);
-      for (ix = 0; ix < m_nx; ix++)
-        complexOut[ix][iy] = complexOutRow[ix];
+      double dImag = 0;
+      if (isComplex())
+        dImag = vLHSImag[ix][iy];
+      pY[iy] = std::complex<double>(vLHS[ix][iy], dImag);
     }
-    delete [] pX;
-    delete [] complexOutRow;
-
+    ProcessSignal::finiteFourierTransform (pY, complexOut[ix], m_ny,  ProcessSignal::FORWARD);
+  }
+  delete [] pY;
+  
+  // fourier each y row
+  CTSimComplex* pX = new CTSimComplex [m_nx];
+  CTSimComplex* complexOutRow = new CTSimComplex [m_nx];
+  for (iy = 0; iy < m_ny; iy++) {
     for (ix = 0; ix < m_nx; ix++)
-      for (iy = 0; iy < m_ny; iy++) {
-        vRealResult[ix][iy] = complexOut[ix][iy].real();
-        vImagResult[ix][iy] = complexOut[ix][iy].imag();
-      }
-
-      Fourier::shuffleFourierToNaturalOrder (result);
-
-      // delete complexOut matrix
-      for (ix = 0; ix < m_nx; ix++)
-        delete [] complexOut[ix];
-      delete [] complexOut;
-
-      return true;
+      pX[ix] = complexOut[ix][iy];
+    ProcessSignal::finiteFourierTransform (pX, complexOutRow, m_nx, ProcessSignal::FORWARD);
+    for (ix = 0; ix < m_nx; ix++)
+      complexOut[ix][iy] = complexOutRow[ix];
+  }
+  delete [] pX;
+  delete [] complexOutRow;
+  
+  for (ix = 0; ix < m_nx; ix++)
+    for (iy = 0; iy < m_ny; iy++) {
+      vRealResult[ix][iy] = complexOut[ix][iy].real();
+      vImagResult[ix][iy] = complexOut[ix][iy].imag();
+    }
+  
+  Fourier::shuffleFourierToNaturalOrder (result);
+  
+  // delete complexOut matrix
+  for (ix = 0; ix < m_nx; ix++)
+    delete [] complexOut[ix];
+  delete [] complexOut;
+  
+  return true;
 }
 
 bool
@@ -1122,7 +1126,7 @@ ImageFile::inverseFourier (ImageFile& result) const
     complexOut[ix] = new CTSimComplex [m_ny];
 
   // put input image into result
-  for (ix = 0; ix < m_nx; ix++)
+  for (ix = 0; ix < m_nx; ix++) {
     for (iy = 0; iy < m_ny; iy++) {
       vRealResult[ix][iy] = vLHSReal[ix][iy];
       if (isComplex())
@@ -1130,44 +1134,45 @@ ImageFile::inverseFourier (ImageFile& result) const
       else
         vImagResult[ix][iy] = 0;
     }
+  }
 
-    Fourier::shuffleNaturalToFourierOrder (result);
-
-    // ifourier each x column
-    CTSimComplex* pCol = new CTSimComplex [m_ny];
-    for (ix = 0; ix < m_nx; ix++) {
-      for (iy = 0; iy < m_ny; iy++) {
-        pCol[iy] = std::complex<double> (vRealResult[ix][iy], vImagResult[ix][iy]);
-      }
-      ProcessSignal::finiteFourierTransform (pCol, complexOut[ix], m_ny,  ProcessSignal::BACKWARD);
-    }
-    delete [] pCol;
+  Fourier::shuffleNaturalToFourierOrder (result);
 
-    // ifourier each y row
-    CTSimComplex* complexInRow = new CTSimComplex [m_nx];
-    CTSimComplex* complexOutRow = new CTSimComplex [m_nx];
+  // ifourier each x column
+  CTSimComplex* pCol = new CTSimComplex [m_ny];
+  for (ix = 0; ix < m_nx; ix++) {
     for (iy = 0; iy < m_ny; iy++) {
-      for (ix = 0; ix < m_nx; ix++)
-        complexInRow[ix] = complexOut[ix][iy];
-      ProcessSignal::finiteFourierTransform (complexInRow, complexOutRow, m_nx, ProcessSignal::BACKWARD);
-      for (ix = 0; ix < m_nx; ix++)
-        complexOut[ix][iy] = complexOutRow[ix];
+      pCol[iy] = std::complex<double> (vRealResult[ix][iy], vImagResult[ix][iy]);
     }
-    delete [] complexInRow;
-    delete [] complexOutRow;
+    ProcessSignal::finiteFourierTransform (pCol, complexOut[ix], m_ny,  ProcessSignal::BACKWARD);
+  }
+  delete [] pCol;
 
+  // ifourier each y row
+  CTSimComplex* complexInRow = new CTSimComplex [m_nx];
+  CTSimComplex* complexOutRow = new CTSimComplex [m_nx];
+  for (iy = 0; iy < m_ny; iy++) {
     for (ix = 0; ix < m_nx; ix++)
-      for (iy = 0; iy < m_ny; iy++) {
-        vRealResult[ix][iy] = complexOut[ix][iy].real();
-        vImagResult[ix][iy] = complexOut[ix][iy].imag();
-      }
-
-      // delete complexOut matrix
-      for (ix = 0; ix < m_nx; ix++)
-        delete [] complexOut[ix];
-      delete [] complexOut;
-
-      return true;
+      complexInRow[ix] = complexOut[ix][iy];
+    ProcessSignal::finiteFourierTransform (complexInRow, complexOutRow, m_nx, ProcessSignal::BACKWARD);
+    for (ix = 0; ix < m_nx; ix++)
+      complexOut[ix][iy] = complexOutRow[ix];
+  }
+  delete [] complexInRow;
+  delete [] complexOutRow;
+  
+  for (ix = 0; ix < m_nx; ix++)
+    for (iy = 0; iy < m_ny; iy++) {
+      vRealResult[ix][iy] = complexOut[ix][iy].real();
+      vImagResult[ix][iy] = complexOut[ix][iy].imag();
+    }
+  
+  // delete complexOut matrix
+  for (ix = 0; ix < m_nx; ix++)
+    delete [] complexOut[ix];
+  delete [] complexOut;
+  
+  return true;
 }
 
 
@@ -1183,18 +1188,18 @@ ImageFile::magnitude (ImageFile& result) const
   ImageFileArray vImag = getImaginaryArray();
   ImageFileArray vRealResult = result.getArray();
 
-  for (unsigned int ix = 0; ix < m_nx; ix++)
+  for (unsigned int ix = 0; ix < m_nx; ix++) {
     for (unsigned int iy = 0; iy < m_ny; iy++) {
       if (isComplex())
         vRealResult[ix][iy] = ::sqrt (vReal[ix][iy] * vReal[ix][iy] + vImag[ix][iy] * vImag[ix][iy]);
       else
         vRealResult[ix][iy] = ::fabs(vReal[ix][iy]);
     }
+  }
+  if (result.isComplex())
+    result.reallocComplexToReal();
 
-    if (result.isComplex())
-      result.reallocComplexToReal();
-
-    return true;
+  return true;
 }
 
 bool
@@ -1311,13 +1316,13 @@ ImageFile::convertExportFormatNameToID (const char* const formatName)
 {
   int formatID = EXPORT_FORMAT_INVALID;
 
-  for (int i = 0; i < s_iExportFormatCount; i++)
+  for (int i = 0; i < s_iExportFormatCount; i++) {
     if (strcasecmp (formatName, s_aszExportFormatName[i]) == 0) {
       formatID = i;
       break;
     }
-
-    return (formatID);
+  }
+  return (formatID);
 }
 
 const char*
@@ -1325,9 +1330,9 @@ ImageFile::convertExportFormatIDToName (int formatID)
 {
   static const char *formatName = "";
 
-  if (formatID >= 0 && formatID < s_iExportFormatCount)
+  if (formatID >= 0 && formatID < s_iExportFormatCount) {
     return (s_aszExportFormatName[formatID]);
-
+  }
   return (formatName);
 }
 
@@ -1347,13 +1352,14 @@ ImageFile::convertImportFormatNameToID (const char* const formatName)
 {
   int formatID = IMPORT_FORMAT_INVALID;
 
-  for (int i = 0; i < s_iImportFormatCount; i++)
+  for (int i = 0; i < s_iImportFormatCount; i++) {
     if (strcasecmp (formatName, s_aszImportFormatName[i]) == 0) {
       formatID = i;
       break;
     }
+  }
 
-    return (formatID);
+  return (formatID);
 }
 
 const char*
@@ -1565,7 +1571,7 @@ ImageFile::readImagePNG (const char* const pszFile)
     return false;
   }
 
-  if (setjmp(png_ptr->jmpbuf)) {
+  if (setjmp(png_jmpbuf(png_ptr))) {
     png_destroy_read_struct(&png_ptr, &info_ptr, &end_info);
     fclose(fp);
     return false;
@@ -1815,7 +1821,7 @@ ImageFile::writeImagePNG (const char* const outfile, int bitdepth, int nxcell, i
     return false;
   }
 
-  if (setjmp (png_ptr->jmpbuf)) {
+  if (setjmp(png_jmpbuf(png_ptr))) {
     png_destroy_write_struct (&png_ptr, &info_ptr);
     fclose (fp);
     return false;