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

diff --git a/libctsim/imagefile.cpp b/libctsim/imagefile.cpp
index b86ff79..54f38a1 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.23 2000/12/21 03:40:58 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;
@@ -266,7 +266,340 @@ ImageFile::getMinMax (double& min, double& max) const
 		}
     }
 }
+
+bool
+ImageFile::subtractImages (const ImageFile& rRHS, ImageFile& result) const
+{
+  if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) {
+    sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");
+    return false;
+  }
+
+  ImageFileArrayConst vLHS = getArray();
+  ImageFileArrayConst vRHS = rRHS.getArray();
+  ImageFileArray vResult = result.getArray();
+
+  for (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++)
+        *out++ = *in1++ - *in2++;
+  }
+
+    return true;
+}
 
+bool
+ImageFile::addImages (const ImageFile& rRHS, ImageFile& result) const
+{
+  if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) {
+    sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");
+    return false;
+  }
+
+  ImageFileArrayConst vLHS = getArray();
+  ImageFileArrayConst vRHS = rRHS.getArray();
+  ImageFileArray vResult = result.getArray();
+
+  for (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++)
+        *out++ = *in1++ + *in2++;
+  }
+
+    return true;
+}
+
+bool
+ImageFile::multiplyImages (const ImageFile& rRHS, ImageFile& result) const
+{
+  if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) {
+    sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");
+    return false;
+  }
+
+  ImageFileArrayConst vLHS = getArray();
+  ImageFileArrayConst vRHS = rRHS.getArray();
+  ImageFileArray vResult = result.getArray();
+
+  for (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++)
+        *out++ = *in1++ * *in2++;
+  }
+
+    return true;
+}
+
+bool
+ImageFile::divideImages (const ImageFile& rRHS, ImageFile& result) const
+{
+  if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) {
+    sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");
+    return false;
+  }
+
+  ImageFileArrayConst vLHS = getArray();
+  ImageFileArrayConst vRHS = rRHS.getArray();
+  ImageFileArray vResult = result.getArray();
+
+  for (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++) {
+      if (*in2 != 0.)
+        *out++ = *in1++ / *in2++;
+      else
+        *out++ = 0;
+    }
+  }
+
+    return true;
+}
+
+
+bool
+ImageFile::invertPixelValues (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();
+
+  for (int ix = 0; ix < m_nx; ix++) {
+    ImageFileColumnConst in = vLHS[ix];
+    ImageFileColumn out = vResult[ix];
+    for (int iy = 0; iy < m_ny; iy++)
+        *out++ = - *in++;
+  }
+
+  return true;
+}
+
+bool
+ImageFile::sqrt (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();
+
+  for (int ix = 0; ix < m_nx; ix++) {
+    ImageFileColumnConst in = vLHS[ix];
+    ImageFileColumn out = vResult[ix];
+    for (int iy = 0; iy < m_ny; iy++)
+      if (*in < 0)
+        *out++ = -::sqrt(-*in++);
+      else
+        *out++ = ::sqrt(*in++);
+  }
+
+  return true;
+}
+
+bool
+ImageFile::log (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();
+
+  for (int ix = 0; ix < m_nx; ix++) {
+    ImageFileColumnConst in = vLHS[ix];
+    ImageFileColumn out = vResult[ix];
+    for (int iy = 0; iy < m_ny; iy++)
+      if (*in <= 0)
+        *out++ = 0;
+      else
+        *out++ = ::log(*in++);
+  }
+
+  return true;
+}
+
+bool
+ImageFile::exp (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();
+
+  for (int ix = 0; ix < m_nx; ix++) {
+    ImageFileColumnConst in = vLHS[ix];
+    ImageFileColumn out = vResult[ix];
+    for (int iy = 0; iy < m_ny; iy++)
+      *out++ = ::exp (*in++);
+  }
+
+  return true;
+}
+
+bool
+ImageFile::FFTMagnitude (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);
+  }
+  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 (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];
+
+  }
+  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 (ix = 0; ix < m_nx; ix++)
+    delete [] complexOut[ix];
+
+  delete [] complexOut;
+
+  return true;
+}
+
+bool
+ImageFile::FFTPhase (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);
+  }
+  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 (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];
+
+  }
+  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];
+
+  delete [] complexOut;
+
+  return true;
+}
+
+bool
+ImageFile::square (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();
+
+  for (int ix = 0; ix < m_nx; ix++) {
+    ImageFileColumnConst in = vLHS[ix];
+    ImageFileColumn out = vResult[ix];
+    for (int iy = 0; iy < m_ny; iy++) {
+        *out++ = *in * *in;
+        in++;
+    }
+  }
+
+  return true;
+}
+
+
 void 
 ImageFile::writeImagePGM (const char *outfile, int nxcell, int nycell, double densmin, double densmax)
 {