projects / ctsim.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
r318: *** empty log message ***
[ctsim.git] / libctsim / imagefile.cpp
diff --git a/libctsim/imagefile.cpp b/libctsim/imagefile.cpp
index a5e9bec75e34ea2442152be631e318c1a93b897f..54f38a18550998ee53bfd6b0af12b55bfbeaa55e 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.24 2000/12/22 04:18:00 kevin Exp $
+**  $Id: imagefile.cpp,v 1.25 2000/12/29 15:45:06 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
@@ -62,7 +62,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 +170,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;
@@ -364,6 +364,242 @@ ImageFile::divideImages (const ImageFile& rRHS, ImageFile& result) const
 }\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 (int ix = 0; ix < m_nx; ix++) {\r
+    ImageFileColumnConst in = vLHS[ix];\r
+    ImageFileColumn out = vResult[ix];\r
+    for (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 (int ix = 0; ix < m_nx; ix++) {\r
+    ImageFileColumnConst in = vLHS[ix];\r
+    ImageFileColumn out = vResult[ix];\r
+    for (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 (int ix = 0; ix < m_nx; ix++) {\r
+    ImageFileColumnConst in = vLHS[ix];\r
+    ImageFileColumn out = vResult[ix];\r
+    for (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 (int ix = 0; ix < m_nx; ix++) {\r
+    ImageFileColumnConst in = vLHS[ix];\r
+    ImageFileColumn out = vResult[ix];\r
+    for (int iy = 0; iy < m_ny; iy++)\r
+      *out++ = ::exp (*in++);\r
+  }\r
+\r
+  return true;\r
+}\r
+\r
+bool\r
+ImageFile::FFTMagnitude (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
+  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
+      vResult[ix][iy] = std::abs (complexOut[ix][iy]);\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::FFTPhase (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
+  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
+      vResult[ix][iy] = atan (complexOut[ix][iy].imag / complexOut[ix][iy].real);\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::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 (int ix = 0; ix < m_nx; ix++) {\r
+    ImageFileColumnConst in = vLHS[ix];\r
+    ImageFileColumn out = vResult[ix];\r
+    for (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