projects / ctsim.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
r322: *** empty log message ***
[ctsim.git] / libctsim / imagefile.cpp
diff --git a/libctsim/imagefile.cpp b/libctsim/imagefile.cpp
index b86ff794ae9cd0242201c21e33184bca9c513e3f..b6525e1fda69077f252b2994aede0cb997d01d21 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
 **
-**  $Id: imagefile.cpp,v 1.23 2000/12/21 03:40:58 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)\r
-: Array2dFile (nx, ny, sizeof(kfloat32), Array2dFile::PIXEL_FLOAT32)\r
+F32Image::F32Image (int nx, int ny, int dataType)\r
+: Array2dFile (nx, ny, sizeof(kfloat32), Array2dFile::PIXEL_FLOAT32, dataType)\r
 {\r
 }\r
 \r
@@ -38,10 +38,11 @@ F32Image::F32Image (void)
 {\r
        setPixelFormat (Array2dFile::PIXEL_FLOAT32);\r
        setPixelSize (sizeof(kfloat32));\r
+  setDataType (Array2dFile::DATA_TYPE_REAL);\r
 }\r
 \r
-F64Image::F64Image (int nx, int ny)\r
-: Array2dFile (nx, ny, sizeof(kfloat64), Array2dFile::PIXEL_FLOAT64)\r
+F64Image::F64Image (int nx, int ny, int dataType)\r
+: Array2dFile (nx, ny, sizeof(kfloat64), Array2dFile::PIXEL_FLOAT64, dataType)\r
 {\r
 }\r
 \r
@@ -50,6 +51,7 @@ F64Image::F64Image (void)
 {\r
        setPixelFormat (PIXEL_FLOAT64);\r
        setPixelSize (sizeof(kfloat64));\r
+  setDataType (Array2dFile::DATA_TYPE_REAL);\r
 }\r
 
 void 
@@ -62,7 +64,7 @@ ImageFile::filterResponse (const char* const domainName, double bw, const char*
        
        for (int i = -hx; i <= hx; i++) {
                for (int j = -hy; j <= hy; j++) {
-                       double r = sqrt (i * i + j * j);
+      double r = ::sqrt (i * i + j * j);
                        
                        v[i+hx][j+hy] = filter.response (r);
                }
@@ -170,7 +172,7 @@ ImageFile::comparativeStatistics (const ImageFile& imComp, double& d, double& r,
                }
     }
        
-    d = sqrt (sqErrorSum / sqDiffFromMeanSum);
+    d = ::sqrt (sqErrorSum / sqDiffFromMeanSum);
     r = absErrorSum / absValueSum;
        
     int hx = m_nx / 2;
@@ -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);\r
+    if (dataType() == Array2dFile::DATA_TYPE_COMPLEX)\r
+      os << "Real Component Statistics" << std::endl;\r
+
     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;\r
+\r
+    if (dataType() == Array2dFile::DATA_TYPE_COMPLEX) {\r
+      statistics (getImaginaryArray(), min, max, mean, mode, median, stddev);\r
+      os << std::endl << "Imaginary Component Statistics" << std::endl;
+    os << "   min: " << min << std::endl;\r
+    os << "   max: " << max << std::endl;\r
+    os << "  mean: " << mean << std::endl;\r
+    os << "  mode: " << mode << std::endl;\r
+    os << "median: " << median << std::endl;\r
+    os << "stddev: " << stddev << std::endl;\r
+    }\r
 }
 
 
 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;
+{\r
+  ImageFileArrayConst v = getArray();\r
+  statistics (v, min, max, mean, mode, median, stddev);\r
+}
+
+
+void\r
+ImageFile::statistics (ImageFileArrayConst& v, double& min, double& max, double& mean, double& mode, double& median, double& stddev) const\r
+{\r
+    int nx = m_nx;\r
+    int ny = m_ny;\r
+    \r
+    if (v == NULL || nx == 0 || ny == 0)\r
+               return;\r
 \r
        std::vector<double> vecImage;\r
        int iVec = 0;\r
@@ -242,9 +264,8 @@ ImageFile::statistics (double& min, double& max, double& mean, double& mode, dou
        }\r
 \r
        vectorNumericStatistics (vecImage, nx * ny, min, max, mean, mode, median, stddev);\r
-}
-
-
+}\r
+\r
 void
 ImageFile::getMinMax (double& min, double& max) const
 {
@@ -266,7 +287,402 @@ ImageFile::getMinMax (double& min, double& max) const
                }
     }
 }
+\r
+bool\r
+ImageFile::convertRealToComplex ()\r
+{\r
+  if (dataType() != Array2dFile::DATA_TYPE_REAL)\r
+    return false;\r
+\r
+  if (! reallocRealToComplex())\r
+    return false;\r
+\r
+  ImageFileArray vImag = getImaginaryArray();\r
+  for (unsigned int ix = 0; ix < m_nx; ix++) {\r
+    ImageFileColumn vCol = vImag[ix];\r
+    for (unsigned int iy = 0; iy < m_ny; iy++)\r
+      *vCol++ = 0;\r
+  }\r
+\r
+  return true;\r
+}\r
+\r
+bool\r
+ImageFile::convertComplexToReal ()\r
+{\r
+  if (dataType() != Array2dFile::DATA_TYPE_COMPLEX)\r
+    return false;\r
+\r
+  ImageFileArray vReal = getArray();\r
+  ImageFileArray vImag = getImaginaryArray();\r
+  for (unsigned int ix = 0; ix < m_nx; ix++) {\r
+    ImageFileColumn vRealCol = vReal[ix];\r
+    ImageFileColumn vImagCol = vImag[ix];\r
+    for (unsigned int iy = 0; iy < m_ny; iy++) {\r
+      std::complex<double> c (*vRealCol, *vImagCol);\r
+      *vRealCol++ = std::abs (c);\r
+      vImagCol++;\r
+    }\r
+  }\r
+\r
+  return reallocComplexToReal();\r
+}\r
+\r
+bool\r
+ImageFile::subtractImages (const ImageFile& rRHS, ImageFile& result) const\r
+{\r
+  if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) {\r
+    sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");\r
+    return false;\r
+  }\r
+\r
+  ImageFileArrayConst vLHS = getArray();\r
+  ImageFileArrayConst vRHS = rRHS.getArray();\r
+  ImageFileArray vResult = result.getArray();\r
+\r
+  for (unsigned int ix = 0; ix < m_nx; ix++) {\r
+    ImageFileColumnConst in1 = vLHS[ix];\r
+    ImageFileColumnConst in2 = vRHS[ix];\r
+    ImageFileColumn out = vResult[ix];\r
+    for (unsigned int iy = 0; iy < m_ny; iy++)\r
+        *out++ = *in1++ - *in2++;\r
+  }\r
+\r
+    return true;\r
+}\r
 
+bool\r
+ImageFile::addImages (const ImageFile& rRHS, ImageFile& result) const\r
+{\r
+  if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) {\r
+    sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");\r
+    return false;\r
+  }\r
+\r
+  ImageFileArrayConst vLHS = getArray();\r
+  ImageFileArrayConst vRHS = rRHS.getArray();\r
+  ImageFileArray vResult = result.getArray();\r
+\r
+  for (unsigned int ix = 0; ix < m_nx; ix++) {\r
+    ImageFileColumnConst in1 = vLHS[ix];\r
+    ImageFileColumnConst in2 = vRHS[ix];\r
+    ImageFileColumn out = vResult[ix];\r
+    for (unsigned int iy = 0; iy < m_ny; iy++)\r
+        *out++ = *in1++ + *in2++;\r
+  }\r
+\r
+    return true;\r
+}\r
+\r
+bool\r
+ImageFile::multiplyImages (const ImageFile& rRHS, ImageFile& result) const\r
+{\r
+  if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) {\r
+    sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");\r
+    return false;\r
+  }\r
+\r
+  ImageFileArrayConst vLHS = getArray();\r
+  ImageFileArrayConst vRHS = rRHS.getArray();\r
+  ImageFileArray vResult = result.getArray();\r
+\r
+  for (unsigned int ix = 0; ix < m_nx; ix++) {\r
+    ImageFileColumnConst in1 = vLHS[ix];\r
+    ImageFileColumnConst in2 = vRHS[ix];\r
+    ImageFileColumn out = vResult[ix];\r
+    for (unsigned int iy = 0; iy < m_ny; iy++)\r
+        *out++ = *in1++ * *in2++;\r
+  }\r
+\r
+    return true;\r
+}\r
+\r
+bool\r
+ImageFile::divideImages (const ImageFile& rRHS, ImageFile& result) const\r
+{\r
+  if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) {\r
+    sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");\r
+    return false;\r
+  }\r
+\r
+  ImageFileArrayConst vLHS = getArray();\r
+  ImageFileArrayConst vRHS = rRHS.getArray();\r
+  ImageFileArray vResult = result.getArray();\r
+\r
+  for (unsigned int ix = 0; ix < m_nx; ix++) {\r
+    ImageFileColumnConst in1 = vLHS[ix];\r
+    ImageFileColumnConst in2 = vRHS[ix];\r
+    ImageFileColumn out = vResult[ix];\r
+    for (unsigned int iy = 0; iy < m_ny; iy++) {\r
+      if (*in2 != 0.)\r
+        *out++ = *in1++ / *in2++;\r
+      else\r
+        *out++ = 0;\r
+    }\r
+  }\r
+\r
+    return true;\r
+}\r
+\r
+\r
+bool\r
+ImageFile::invertPixelValues (ImageFile& result) const\r
+{\r
+  if (m_nx != result.nx() || m_ny != result.ny()) {\r
+    sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
+    return false;\r
+  }\r
+\r
+  ImageFileArrayConst vLHS = getArray();\r
+  ImageFileArray vResult = result.getArray();\r
+\r
+  for (unsigned int ix = 0; ix < m_nx; ix++) {\r
+    ImageFileColumnConst in = vLHS[ix];\r
+    ImageFileColumn out = vResult[ix];\r
+    for (unsigned int iy = 0; iy < m_ny; iy++)\r
+        *out++ = - *in++;\r
+  }\r
+\r
+  return true;\r
+}\r
+\r
+bool\r
+ImageFile::sqrt (ImageFile& result) const\r
+{\r
+  if (m_nx != result.nx() || m_ny != result.ny()) {\r
+    sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
+    return false;\r
+  }\r
+\r
+  ImageFileArrayConst vLHS = getArray();\r
+  ImageFileArray vResult = result.getArray();\r
+\r
+  for (unsigned int ix = 0; ix < m_nx; ix++) {\r
+    ImageFileColumnConst in = vLHS[ix];\r
+    ImageFileColumn out = vResult[ix];\r
+    for (unsigned int iy = 0; iy < m_ny; iy++)\r
+      if (*in < 0)\r
+        *out++ = -::sqrt(-*in++);\r
+      else\r
+        *out++ = ::sqrt(*in++);\r
+  }\r
+\r
+  return true;\r
+}\r
+\r
+bool\r
+ImageFile::log (ImageFile& result) const\r
+{\r
+  if (m_nx != result.nx() || m_ny != result.ny()) {\r
+    sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
+    return false;\r
+  }\r
+\r
+  ImageFileArrayConst vLHS = getArray();\r
+  ImageFileArray vResult = result.getArray();\r
+\r
+  for (unsigned int ix = 0; ix < m_nx; ix++) {\r
+    ImageFileColumnConst in = vLHS[ix];\r
+    ImageFileColumn out = vResult[ix];\r
+    for (unsigned int iy = 0; iy < m_ny; iy++)\r
+      if (*in <= 0)\r
+        *out++ = 0;\r
+      else\r
+        *out++ = ::log(*in++);\r
+  }\r
+\r
+  return true;\r
+}\r
+\r
+bool\r
+ImageFile::exp (ImageFile& result) const\r
+{\r
+  if (m_nx != result.nx() || m_ny != result.ny()) {\r
+    sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
+    return false;\r
+  }\r
+\r
+  ImageFileArrayConst vLHS = getArray();\r
+  ImageFileArray vResult = result.getArray();\r
+\r
+  for (unsigned int ix = 0; ix < m_nx; ix++) {\r
+    ImageFileColumnConst in = vLHS[ix];\r
+    ImageFileColumn out = vResult[ix];\r
+    for (unsigned int iy = 0; iy < m_ny; iy++)\r
+      *out++ = ::exp (*in++);\r
+  }\r
+\r
+  return true;\r
+}\r
+\r
+bool\r
+ImageFile::inverseFourier (ImageFile& result) const\r
+{\r
+  if (m_nx != result.nx() || m_ny != result.ny()) {\r
+    sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
+    return false;\r
+  }\r
+\r
+  return true;\r
+}\r
+\r
+bool\r
+ImageFile::fourier (ImageFile& result) const\r
+{\r
+  if (m_nx != result.nx() || m_ny != result.ny()) {\r
+    sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
+    return false;\r
+  }\r
+\r
+  if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {\r
+    if (! result.reallocRealToComplex ())\r
+      return false;\r
+  }\r
+\r
+  ImageFileArrayConst vLHS = getArray();\r
+  ImageFileArray vRealResult = result.getArray();\r
+  ImageFileArray vImagResult = result.getImaginaryArray();\r
+\r
+  unsigned int ix, iy;\r
+  double* pY = new double [m_ny];\r
+  std::complex<double>** complexOut = new std::complex<double>* [m_nx];\r
+  for (ix = 0; ix < m_nx; ix++)\r
+    complexOut[ix] = new std::complex<double> [m_ny];\r
+\r
+  for (ix = 0; ix < m_nx; ix++) {\r
+    for (iy = 0; iy < m_ny; iy++)\r
+      pY[iy] = vLHS[ix][iy];\r
+    ProcessSignal::finiteFourierTransform (pY, complexOut[ix], m_ny,  ProcessSignal::FORWARD);\r
+  }\r
+  delete pY;\r
+\r
+  std::complex<double>* pX = new std::complex<double> [m_nx];\r
+  std::complex<double>* complexOutCol = new std::complex<double> [m_nx];\r
+  for (iy = 0; iy < m_ny; iy++) {\r
+    for (ix = 0; ix < m_nx; ix++)\r
+      pX[ix] = complexOut[ix][iy];\r
+    ProcessSignal::finiteFourierTransform (pX, complexOutCol, m_nx, ProcessSignal::FORWARD);\r
+    for (ix = 0; ix < m_nx; ix++)\r
+      complexOut[ix][iy] = complexOutCol[ix];\r
+  }\r
+  delete [] pX;\r
+\r
+  for (ix = 0; ix < m_nx; ix++)\r
+    ProcessSignal::shuffleFourierToNaturalOrder (complexOut[ix], m_ny);\r
+\r
+  \r
+  for (iy = 0; iy < m_ny; iy++) {\r
+    for (ix = 0; ix < m_nx; ix++)\r
+      complexOutCol[ix] = complexOut[ix][iy];\r
+    ProcessSignal::shuffleFourierToNaturalOrder (complexOutCol, m_nx);;\r
+    for (ix = 0; ix < m_nx; ix++)\r
+      complexOut[ix][iy] = complexOutCol[ix];\r
+\r
+  }\r
+  delete [] complexOutCol;\r
+\r
+  for (ix = 0; ix < m_nx; ix++)\r
+    for (iy = 0; iy < m_ny; iy++) {\r
+      vRealResult[ix][iy] = complexOut[ix][iy].real();\r
+      vImagResult[ix][iy] = complexOut[ix][iy].imag();\r
+    }\r
+\r
+  for (ix = 0; ix < m_nx; ix++)\r
+    delete [] complexOut[ix];\r
+\r
+  delete [] complexOut;\r
+\r
+  return true;\r
+}\r
+\r
+bool\r
+ImageFile::magnitude (ImageFile& result) const\r
+{\r
+  if (m_nx != result.nx() || m_ny != result.ny()) {\r
+    sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
+    return false;\r
+  }\r
+\r
+  ImageFileArray vReal = getArray();\r
+  ImageFileArray vImag = NULL;\r
+  if (dataType() == Array2dFile::DATA_TYPE_COMPLEX)\r
+    vImag = getImaginaryArray();\r
+\r
+  ImageFileArray vRealResult = result.getArray();\r
+  if (result.dataType() == Array2dFile::DATA_TYPE_COMPLEX) {\r
+    ImageFileArray vImagResult = result.getImaginaryArray();\r
+    for (unsigned int ix = 0; ix < m_nx; ix++)\r
+      for (unsigned int iy = 0; iy < m_ny; iy++)\r
+        vImagResult[ix][iy] = 0;\r
+  }\r
+\r
+  for (unsigned int ix = 0; ix < m_nx; ix++)\r
+    for (unsigned int iy = 0; iy < m_ny; iy++) {\r
+      if (dataType() == Array2dFile::DATA_TYPE_COMPLEX) \r
+        vRealResult[ix][iy] = ::sqrt (vReal[ix][iy] * vReal[ix][iy] + vImag[ix][iy] * vImag[ix][iy]);\r
+      else\r
+        vRealResult[ix][iy] = vReal[ix][iy];\r
+    }\r
+\r
+  return true;\r
+}\r
+\r
+bool\r
+ImageFile::phase (ImageFile& result) const\r
+{\r
+  if (m_nx != result.nx() || m_ny != result.ny()) {\r
+    sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
+    return false;\r
+  }\r
+\r
+  ImageFileArray vReal = getArray();\r
+  ImageFileArray vImag = NULL;\r
+  if (dataType() == Array2dFile::DATA_TYPE_COMPLEX)\r
+    vImag = getImaginaryArray();\r
+\r
+  ImageFileArray vRealResult = result.getArray();\r
+  if (result.dataType() == Array2dFile::DATA_TYPE_COMPLEX) {\r
+    ImageFileArray vImagResult = result.getImaginaryArray();\r
+    for (unsigned int ix = 0; ix < m_nx; ix++)\r
+      for (unsigned int iy = 0; iy < m_ny; iy++)\r
+        vImagResult[ix][iy] = 0;\r
+  }\r
+\r
+  for (unsigned int ix = 0; ix < m_nx; ix++)\r
+    for (unsigned int iy = 0; iy < m_ny; iy++) {\r
+      if (dataType() == Array2dFile::DATA_TYPE_COMPLEX) \r
+        vRealResult[ix][iy] = ::atan (vImag[ix][iy] / vReal[ix][iy]);\r
+      else\r
+        vRealResult[ix][iy] = vReal[ix][iy];\r
+    }\r
+\r
+  return true;\r
+}\r
+\r
+bool\r
+ImageFile::square (ImageFile& result) const\r
+{\r
+  if (m_nx != result.nx() || m_ny != result.ny()) {\r
+    sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
+    return false;\r
+  }\r
+\r
+  ImageFileArrayConst vLHS = getArray();\r
+  ImageFileArray vResult = result.getArray();\r
+\r
+  for (unsigned int ix = 0; ix < m_nx; ix++) {\r
+    ImageFileColumnConst in = vLHS[ix];\r
+    ImageFileColumn out = vResult[ix];\r
+    for (unsigned int iy = 0; iy < m_ny; iy++) {\r
+        *out++ = *in * *in;\r
+        in++;\r
+    }\r
+  }\r
+\r
+  return true;\r
+}\r
+\r
+\r
 void 
 ImageFile::writeImagePGM (const char *outfile, int nxcell, int nycell, double densmin, double densmax)
 {
Computed Tomography simulator and reconstructor