X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fimagefile.cpp;h=b6525e1fda69077f252b2994aede0cb997d01d21;hp=54f38a18550998ee53bfd6b0af12b55bfbeaa55e;hb=7ae47cb0ff0a16d1c36797576155263434cc73ff;hpb=5c6b29ab4885308cc3381af6e0a68f4804956d2e

diff --git a/libctsim/imagefile.cpp b/libctsim/imagefile.cpp
index 54f38a1..b6525e1 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.25 2000/12/29 15:45:06 kevin Exp $
+**  $Id: imagefile.cpp,v 1.27 2000/12/29 20:15:37 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
@@ -28,8 +28,8 @@
 #include "ct.h"
 
 
-F32Image::F32Image (int nx, int ny)
-: Array2dFile (nx, ny, sizeof(kfloat32), Array2dFile::PIXEL_FLOAT32)
+F32Image::F32Image (int nx, int ny, int dataType)
+: Array2dFile (nx, ny, sizeof(kfloat32), Array2dFile::PIXEL_FLOAT32, dataType)
 {
 }
 
@@ -38,10 +38,11 @@ F32Image::F32Image (void)
 {
 	setPixelFormat (Array2dFile::PIXEL_FLOAT32);
 	setPixelSize (sizeof(kfloat32));
+  setDataType (Array2dFile::DATA_TYPE_REAL);
 }
 
-F64Image::F64Image (int nx, int ny)
-: Array2dFile (nx, ny, sizeof(kfloat64), Array2dFile::PIXEL_FLOAT64)
+F64Image::F64Image (int nx, int ny, int dataType)
+: Array2dFile (nx, ny, sizeof(kfloat64), Array2dFile::PIXEL_FLOAT64, dataType)
 {
 }
 
@@ -50,6 +51,7 @@ F64Image::F64Image (void)
 {
 	setPixelFormat (PIXEL_FLOAT64);
 	setPixelSize (sizeof(kfloat64));
+  setDataType (Array2dFile::DATA_TYPE_REAL);
 }
 
 void 
@@ -212,26 +214,46 @@ ImageFile::printStatistics (std::ostream& os) const
 {
     double min, max, mean, mode, median, stddev;
 	
-    statistics (min, max, mean, mode, median, stddev);
-	
+    statistics (min, max, mean, mode, median, stddev);
+    if (dataType() == Array2dFile::DATA_TYPE_COMPLEX)
+      os << "Real Component Statistics" << std::endl;
+
     os << "   min: " << min << std::endl;
     os << "   max: " << max << std::endl;
     os << "  mean: " << mean << std::endl;
     os << "  mode: " << mode << std::endl;
     os << "median: " << median << std::endl;
-    os << "stddev: " << stddev << std::endl;
+    os << "stddev: " << stddev << std::endl;
+
+    if (dataType() == Array2dFile::DATA_TYPE_COMPLEX) {
+      statistics (getImaginaryArray(), min, max, mean, mode, median, stddev);
+      os << std::endl << "Imaginary Component Statistics" << std::endl;
+    os << "   min: " << min << std::endl;
+    os << "   max: " << max << std::endl;
+    os << "  mean: " << mean << std::endl;
+    os << "  mode: " << mode << std::endl;
+    os << "median: " << median << std::endl;
+    os << "stddev: " << stddev << std::endl;
+    }
 }
 
 
 void
 ImageFile::statistics (double& min, double& max, double& mean, double& mode, double& median, double& stddev) const
-{
-    int nx = m_nx;
-    int ny = m_ny;
-    ImageFileArrayConst v = getArray();
-    
-    if (v == NULL || nx == 0 || ny == 0)
-		return;
+{
+  ImageFileArrayConst v = getArray();
+  statistics (v, min, max, mean, mode, median, stddev);
+}
+
+
+void
+ImageFile::statistics (ImageFileArrayConst& v, double& min, double& max, double& mean, double& mode, double& median, double& stddev) const
+{
+    int nx = m_nx;
+    int ny = m_ny;
+    
+    if (v == NULL || nx == 0 || ny == 0)
+		return;
 
 	std::vector<double> vecImage;
 	int iVec = 0;
@@ -242,9 +264,8 @@ ImageFile::statistics (double& min, double& max, double& mean, double& mode, dou
 	}
 
 	vectorNumericStatistics (vecImage, nx * ny, min, max, mean, mode, median, stddev);
-}
-
-
+}
+
 void
 ImageFile::getMinMax (double& min, double& max) const
 {
@@ -267,6 +288,46 @@ ImageFile::getMinMax (double& min, double& max) const
     }
 }
 
+bool
+ImageFile::convertRealToComplex ()
+{
+  if (dataType() != Array2dFile::DATA_TYPE_REAL)
+    return false;
+
+  if (! reallocRealToComplex())
+    return false;
+
+  ImageFileArray vImag = getImaginaryArray();
+  for (unsigned int ix = 0; ix < m_nx; ix++) {
+    ImageFileColumn vCol = vImag[ix];
+    for (unsigned int iy = 0; iy < m_ny; iy++)
+      *vCol++ = 0;
+  }
+
+  return true;
+}
+
+bool
+ImageFile::convertComplexToReal ()
+{
+  if (dataType() != Array2dFile::DATA_TYPE_COMPLEX)
+    return false;
+
+  ImageFileArray vReal = getArray();
+  ImageFileArray vImag = getImaginaryArray();
+  for (unsigned int ix = 0; ix < m_nx; ix++) {
+    ImageFileColumn vRealCol = vReal[ix];
+    ImageFileColumn vImagCol = vImag[ix];
+    for (unsigned int iy = 0; iy < m_ny; iy++) {
+      std::complex<double> c (*vRealCol, *vImagCol);
+      *vRealCol++ = std::abs (c);
+      vImagCol++;
+    }
+  }
+
+  return reallocComplexToReal();
+}
+
 bool
 ImageFile::subtractImages (const ImageFile& rRHS, ImageFile& result) const
 {
@@ -279,11 +340,11 @@ ImageFile::subtractImages (const ImageFile& rRHS, ImageFile& result) const
   ImageFileArrayConst vRHS = rRHS.getArray();
   ImageFileArray vResult = result.getArray();
 
-  for (int ix = 0; ix < m_nx; ix++) {
+  for (unsigned int ix = 0; ix < m_nx; ix++) {
     ImageFileColumnConst in1 = vLHS[ix];
     ImageFileColumnConst in2 = vRHS[ix];
     ImageFileColumn out = vResult[ix];
-    for (int iy = 0; iy < m_ny; iy++)
+    for (unsigned int iy = 0; iy < m_ny; iy++)
         *out++ = *in1++ - *in2++;
   }
 
@@ -302,11 +363,11 @@ ImageFile::addImages (const ImageFile& rRHS, ImageFile& result) const
   ImageFileArrayConst vRHS = rRHS.getArray();
   ImageFileArray vResult = result.getArray();
 
-  for (int ix = 0; ix < m_nx; ix++) {
+  for (unsigned int ix = 0; ix < m_nx; ix++) {
     ImageFileColumnConst in1 = vLHS[ix];
     ImageFileColumnConst in2 = vRHS[ix];
     ImageFileColumn out = vResult[ix];
-    for (int iy = 0; iy < m_ny; iy++)
+    for (unsigned int iy = 0; iy < m_ny; iy++)
         *out++ = *in1++ + *in2++;
   }
 
@@ -325,11 +386,11 @@ ImageFile::multiplyImages (const ImageFile& rRHS, ImageFile& result) const
   ImageFileArrayConst vRHS = rRHS.getArray();
   ImageFileArray vResult = result.getArray();
 
-  for (int ix = 0; ix < m_nx; ix++) {
+  for (unsigned int ix = 0; ix < m_nx; ix++) {
     ImageFileColumnConst in1 = vLHS[ix];
     ImageFileColumnConst in2 = vRHS[ix];
     ImageFileColumn out = vResult[ix];
-    for (int iy = 0; iy < m_ny; iy++)
+    for (unsigned int iy = 0; iy < m_ny; iy++)
         *out++ = *in1++ * *in2++;
   }
 
@@ -348,11 +409,11 @@ ImageFile::divideImages (const ImageFile& rRHS, ImageFile& result) const
   ImageFileArrayConst vRHS = rRHS.getArray();
   ImageFileArray vResult = result.getArray();
 
-  for (int ix = 0; ix < m_nx; ix++) {
+  for (unsigned int ix = 0; ix < m_nx; ix++) {
     ImageFileColumnConst in1 = vLHS[ix];
     ImageFileColumnConst in2 = vRHS[ix];
     ImageFileColumn out = vResult[ix];
-    for (int iy = 0; iy < m_ny; iy++) {
+    for (unsigned int iy = 0; iy < m_ny; iy++) {
       if (*in2 != 0.)
         *out++ = *in1++ / *in2++;
       else
@@ -375,10 +436,10 @@ ImageFile::invertPixelValues (ImageFile& result) const
   ImageFileArrayConst vLHS = getArray();
   ImageFileArray vResult = result.getArray();
 
-  for (int ix = 0; ix < m_nx; ix++) {
+  for (unsigned int ix = 0; ix < m_nx; ix++) {
     ImageFileColumnConst in = vLHS[ix];
     ImageFileColumn out = vResult[ix];
-    for (int iy = 0; iy < m_ny; iy++)
+    for (unsigned int iy = 0; iy < m_ny; iy++)
         *out++ = - *in++;
   }
 
@@ -396,10 +457,10 @@ ImageFile::sqrt (ImageFile& result) const
   ImageFileArrayConst vLHS = getArray();
   ImageFileArray vResult = result.getArray();
 
-  for (int ix = 0; ix < m_nx; ix++) {
+  for (unsigned int ix = 0; ix < m_nx; ix++) {
     ImageFileColumnConst in = vLHS[ix];
     ImageFileColumn out = vResult[ix];
-    for (int iy = 0; iy < m_ny; iy++)
+    for (unsigned int iy = 0; iy < m_ny; iy++)
       if (*in < 0)
         *out++ = -::sqrt(-*in++);
       else
@@ -420,10 +481,10 @@ ImageFile::log (ImageFile& result) const
   ImageFileArrayConst vLHS = getArray();
   ImageFileArray vResult = result.getArray();
 
-  for (int ix = 0; ix < m_nx; ix++) {
+  for (unsigned int ix = 0; ix < m_nx; ix++) {
     ImageFileColumnConst in = vLHS[ix];
     ImageFileColumn out = vResult[ix];
-    for (int iy = 0; iy < m_ny; iy++)
+    for (unsigned int iy = 0; iy < m_ny; iy++)
       if (*in <= 0)
         *out++ = 0;
       else
@@ -444,10 +505,10 @@ ImageFile::exp (ImageFile& result) const
   ImageFileArrayConst vLHS = getArray();
   ImageFileArray vResult = result.getArray();
 
-  for (int ix = 0; ix < m_nx; ix++) {
+  for (unsigned int ix = 0; ix < m_nx; ix++) {
     ImageFileColumnConst in = vLHS[ix];
     ImageFileColumn out = vResult[ix];
-    for (int iy = 0; iy < m_ny; iy++)
+    for (unsigned int iy = 0; iy < m_ny; iy++)
       *out++ = ::exp (*in++);
   }
 
@@ -455,17 +516,34 @@ ImageFile::exp (ImageFile& result) const
 }
 
 bool
-ImageFile::FFTMagnitude (ImageFile& result) const
+ImageFile::inverseFourier (ImageFile& result) const
 {
   if (m_nx != result.nx() || m_ny != result.ny()) {
     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
     return false;
   }
 
+  return true;
+}
+
+bool
+ImageFile::fourier (ImageFile& result) const
+{
+  if (m_nx != result.nx() || m_ny != result.ny()) {
+    sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
+    return false;
+  }
+
+  if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
+    if (! result.reallocRealToComplex ())
+      return false;
+  }
+
   ImageFileArrayConst vLHS = getArray();
-  ImageFileArray vResult = result.getArray();
+  ImageFileArray vRealResult = result.getArray();
+  ImageFileArray vImagResult = result.getImaginaryArray();
 
-  int ix, iy;
+  unsigned int ix, iy;
   double* pY = new double [m_ny];
   std::complex<double>** complexOut = new std::complex<double>* [m_nx];
   for (ix = 0; ix < m_nx; ix++)
@@ -504,8 +582,10 @@ ImageFile::FFTMagnitude (ImageFile& result) const
   delete [] complexOutCol;
 
   for (ix = 0; ix < m_nx; ix++)
-    for (iy = 0; iy < m_ny; iy++)
-      vResult[ix][iy] = std::abs (complexOut[ix][iy]);
+    for (iy = 0; iy < m_ny; iy++) {
+      vRealResult[ix][iy] = complexOut[ix][iy].real();
+      vImagResult[ix][iy] = complexOut[ix][iy].imag();
+    }
 
   for (ix = 0; ix < m_nx; ix++)
     delete [] complexOut[ix];
@@ -516,62 +596,65 @@ ImageFile::FFTMagnitude (ImageFile& result) const
 }
 
 bool
-ImageFile::FFTPhase (ImageFile& result) const
+ImageFile::magnitude (ImageFile& result) const
 {
   if (m_nx != result.nx() || m_ny != result.ny()) {
     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
     return false;
   }
 
-  ImageFileArrayConst vLHS = getArray();
-  ImageFileArray vResult = result.getArray();
-
-  int ix, iy;
-  double* pY = new double [m_ny];
-  std::complex<double>** complexOut = new std::complex<double>* [m_nx];
-  for (ix = 0; ix < m_nx; ix++)
-    complexOut[ix] = new std::complex<double> [m_ny];
-
-  for (ix = 0; ix < m_nx; ix++) {
-    for (iy = 0; iy < m_ny; iy++)
-      pY[iy] = vLHS[ix][iy];
-    ProcessSignal::finiteFourierTransform (pY, complexOut[ix], m_ny,  ProcessSignal::FORWARD);
+  ImageFileArray vReal = getArray();
+  ImageFileArray vImag = NULL;
+  if (dataType() == Array2dFile::DATA_TYPE_COMPLEX)
+    vImag = getImaginaryArray();
+
+  ImageFileArray vRealResult = result.getArray();
+  if (result.dataType() == Array2dFile::DATA_TYPE_COMPLEX) {
+    ImageFileArray vImagResult = result.getImaginaryArray();
+    for (unsigned int ix = 0; ix < m_nx; ix++)
+      for (unsigned int iy = 0; iy < m_ny; iy++)
+        vImagResult[ix][iy] = 0;
   }
-  delete pY;
 
-  std::complex<double>* pX = new std::complex<double> [m_nx];
-  std::complex<double>* complexOutCol = new std::complex<double> [m_nx];
-  for (iy = 0; iy < m_ny; iy++) {
-    for (ix = 0; ix < m_nx; ix++)
-      pX[ix] = complexOut[ix][iy];
-    ProcessSignal::finiteFourierTransform (pX, complexOutCol, m_nx, ProcessSignal::FORWARD);
-    for (ix = 0; ix < m_nx; ix++)
-      complexOut[ix][iy] = complexOutCol[ix];
-  }
-  delete [] pX;
-
-  for (ix = 0; ix < m_nx; ix++)
-    ProcessSignal::shuffleFourierToNaturalOrder (complexOut[ix], m_ny);
+  for (unsigned int ix = 0; ix < m_nx; ix++)
+    for (unsigned int iy = 0; iy < m_ny; iy++) {
+      if (dataType() == Array2dFile::DATA_TYPE_COMPLEX) 
+        vRealResult[ix][iy] = ::sqrt (vReal[ix][iy] * vReal[ix][iy] + vImag[ix][iy] * vImag[ix][iy]);
+      else
+        vRealResult[ix][iy] = vReal[ix][iy];
+    }
 
-  
-  for (iy = 0; iy < m_ny; iy++) {
-    for (ix = 0; ix < m_nx; ix++)
-      complexOutCol[ix] = complexOut[ix][iy];
-    ProcessSignal::shuffleFourierToNaturalOrder (complexOutCol, m_nx);;
-    for (ix = 0; ix < m_nx; ix++)
-      complexOut[ix][iy] = complexOutCol[ix];
+  return true;
+}
 
+bool
+ImageFile::phase (ImageFile& result) const
+{
+  if (m_nx != result.nx() || m_ny != result.ny()) {
+    sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
+    return false;
   }
-  delete [] complexOutCol;
-
-  for (ix = 0; ix < m_nx; ix++)
-    for (iy = 0; iy < m_ny; iy++)
-      vResult[ix][iy] = atan (complexOut[ix][iy].imag / complexOut[ix][iy].real);
 
-  for (ix = 0; ix < m_nx; ix++)
-    delete [] complexOut[ix];
+  ImageFileArray vReal = getArray();
+  ImageFileArray vImag = NULL;
+  if (dataType() == Array2dFile::DATA_TYPE_COMPLEX)
+    vImag = getImaginaryArray();
+
+  ImageFileArray vRealResult = result.getArray();
+  if (result.dataType() == Array2dFile::DATA_TYPE_COMPLEX) {
+    ImageFileArray vImagResult = result.getImaginaryArray();
+    for (unsigned int ix = 0; ix < m_nx; ix++)
+      for (unsigned int iy = 0; iy < m_ny; iy++)
+        vImagResult[ix][iy] = 0;
+  }
 
-  delete [] complexOut;
+  for (unsigned int ix = 0; ix < m_nx; ix++)
+    for (unsigned int iy = 0; iy < m_ny; iy++) {
+      if (dataType() == Array2dFile::DATA_TYPE_COMPLEX) 
+        vRealResult[ix][iy] = ::atan (vImag[ix][iy] / vReal[ix][iy]);
+      else
+        vRealResult[ix][iy] = vReal[ix][iy];
+    }
 
   return true;
 }
@@ -587,10 +670,10 @@ ImageFile::square (ImageFile& result) const
   ImageFileArrayConst vLHS = getArray();
   ImageFileArray vResult = result.getArray();
 
-  for (int ix = 0; ix < m_nx; ix++) {
+  for (unsigned int ix = 0; ix < m_nx; ix++) {
     ImageFileColumnConst in = vLHS[ix];
     ImageFileColumn out = vResult[ix];
-    for (int iy = 0; iy < m_ny; iy++) {
+    for (unsigned int iy = 0; iy < m_ny; iy++) {
         *out++ = *in * *in;
         in++;
     }