+#ifdef HAVE_CTN_DICOM
+#include "ctndicom.h"
+#endif
+#include "interpolator.h"
+
+const double ImageFile::s_dRedGrayscaleFactor = 0.299;
+const double ImageFile::s_dGreenGrayscaleFactor = 0.587;
+const double ImageFile::s_dBlueGrayscaleFactor = 0.114;
+
+
+const int ImageFile::EXPORT_FORMAT_INVALID = -1;
+const int ImageFile::EXPORT_FORMAT_TEXT = 0;
+const int ImageFile::EXPORT_FORMAT_PGM = 1;
+const int ImageFile::EXPORT_FORMAT_PGMASCII = 2;
+#ifdef HAVE_PNG
+const int ImageFile::EXPORT_FORMAT_PNG = 3;
+const int ImageFile::EXPORT_FORMAT_PNG16 = 4;
+#endif
+#ifdef HAVE_CTN_DICOM
+const int ImageFile::EXPORT_FORMAT_DICOM = 5;
+#endif
+const int ImageFile::EXPORT_FORMAT_RAW = 6;
+
+const char* ImageFile::s_aszExportFormatName[] =
+{
+ "text",
+ "pgm",
+ "pgmascii",
+#ifdef HAVE_PNG
+ "png",
+ "png16",
+#endif
+#ifdef HAVE_CTN_DICOM
+ "dicom",
+#endif
+};
+
+const char* ImageFile::s_aszExportFormatTitle[] =
+{
+ "Text",
+ "PGM",
+ "PGM ASCII",
+#ifdef HAVE_PNG
+ "PNG",
+ "PNG 16-bit",
+#endif
+#ifdef HAVE_CTN_DICOM
+ "Dicom",
+#endif
+};
+const int ImageFile::s_iExportFormatCount = sizeof(s_aszExportFormatName) / sizeof(const char*);
+
+
+const int ImageFile::IMPORT_FORMAT_INVALID = -1;
+const int ImageFile::IMPORT_FORMAT_PPM = 0;
+#ifdef HAVE_PNG
+const int ImageFile::IMPORT_FORMAT_PNG = 1;
+#endif
+#ifdef HAVE_CTN_DICOM
+const int ImageFile::IMPORT_FORMAT_DICOM = 2;
+#endif
+
+
+const char* ImageFile::s_aszImportFormatName[] =
+{
+ "ppm",
+#ifdef HAVE_PNG
+ "png",
+#endif
+#ifdef HAVE_CTN_DICOM
+ "dicom",
+#endif
+};
+
+const char* ImageFile::s_aszImportFormatTitle[] =
+{
+ "PPM",
+#ifdef HAVE_PNG
+ "PNG",
+#endif
+#ifdef HAVE_CTN_DICOM
+ "Dicom",
+#endif
+};
+const int ImageFile::s_iImportFormatCount = sizeof(s_aszImportFormatName) / sizeof(const char*);
+
+
+
+F32Image::F32Image (int nx, int ny, int dataType)
+: Array2dFile (nx, ny, sizeof(kfloat32), Array2dFile::PIXEL_FLOAT32, dataType)
+{
+}
+
+F32Image::F32Image (void)
+: Array2dFile()
+{
+ setPixelFormat (Array2dFile::PIXEL_FLOAT32);
+ setPixelSize (sizeof(kfloat32));
+ setDataType (Array2dFile::DATA_TYPE_REAL);
+}
+
+F64Image::F64Image (int nx, int ny, int dataType)
+: Array2dFile (nx, ny, sizeof(kfloat64), Array2dFile::PIXEL_FLOAT64, dataType)
+{
+}
+
+F64Image::F64Image (void)
+: Array2dFile ()
+{
+ setPixelFormat (PIXEL_FLOAT64);
+ setPixelSize (sizeof(kfloat64));
+ setDataType (Array2dFile::DATA_TYPE_REAL);
+}
+
+void
+ImageFile::getCenterCoordinates (unsigned int& iXCenter, unsigned int& iYCenter)
+{
+ if (isEven (m_nx))
+ iXCenter = m_nx / 2;
+ else
+ iXCenter = (m_nx - 1) / 2;
+
+ if (isEven (m_ny))
+ iYCenter = m_ny / 2;
+ else
+ iYCenter = (m_ny - 1) / 2;
+}
+
+
+void
+ImageFile::filterResponse (const char* const domainName, double bw, const char* const filterName,
+ double filt_param, double dInputScale, double dOutputScale)
+{
+ ImageFileArray v = getArray();
+ SignalFilter filter (filterName, domainName, bw, filt_param);
+
+ unsigned int iXCenter, iYCenter;
+ getCenterCoordinates (iXCenter, iYCenter);
+
+ for (unsigned int ix = 0; ix < m_nx; ix++)
+ for (unsigned int iy = 0; iy < m_ny; iy++) {
+ long lD2 = ((ix - iXCenter) * (ix - iXCenter)) + ((iy - iYCenter) * (iy - iYCenter));
+ double r = ::sqrt (static_cast<double>(lD2)) * dInputScale;
+ v[ix][iy] = filter.response (r) * dOutputScale;
+ }
+}
+
+
+// ImageFile::comparativeStatistics Calculate comparative stats
+//
+// OUTPUT
+// d Normalized root mean squared distance measure
+// r Normalized mean absolute distance measure
+// e Worst case distance measure
+//
+// REFERENCES
+// G.T. Herman, Image Reconstruction From Projections, 1980
+
+bool
+ImageFile::comparativeStatistics (const ImageFile& imComp, double& d, double& r, double& e) const
+{
+ if (imComp.nx() != m_nx && imComp.ny() != m_ny) {
+ sys_error (ERR_WARNING, "Image sizes differ [ImageFile::comparativeStatistics]");
+ return false;
+ }
+ ImageFileArrayConst v = getArray();
+ if (v == NULL || m_nx == 0 || m_ny == 0)
+ return false;
+
+ ImageFileArrayConst vComp = imComp.getArray();
+
+ double myMean = 0.;
+ for (unsigned int ix = 0; ix < m_nx; ix++) {
+ for (unsigned int iy = 0; iy < m_ny; iy++) {
+ myMean += v[ix][iy];
+ }
+ }
+ myMean /= (m_nx * m_ny);
+
+ double sqErrorSum = 0.;
+ double absErrorSum = 0.;
+ double sqDiffFromMeanSum = 0.;
+ double absValueSum = 0.;
+ for (unsigned int ix2 = 0; ix2 < m_nx; ix2++) {
+ for (unsigned int iy = 0; iy < m_ny; iy++) {
+ double diff = v[ix2][iy] - vComp[ix2][iy];
+ sqErrorSum += diff * diff;
+ absErrorSum += fabs(diff);
+ double diffFromMean = v[ix2][iy] - myMean;
+ sqDiffFromMeanSum += diffFromMean * diffFromMean;
+ absValueSum += fabs(v[ix2][iy]);
+ }
+ }
+
+ d = ::sqrt (sqErrorSum / sqDiffFromMeanSum);
+ r = absErrorSum / absValueSum;
+
+ int hx = m_nx / 2;
+ int hy = m_ny / 2;
+ double eMax = -1;
+ for (int ix3 = 0; ix3 < hx; ix3++) {
+ for (int iy = 0; iy < hy; iy++) {
+ double avgPixel = 0.25 * (v[2*ix3][2*iy] + v[2*ix3+1][2*iy] + v[2*ix3][2*iy+1] + v[2*ix3+1][2*iy+1]);
+ double avgPixelComp = 0.25 * (vComp[2*ix3][2*iy] + vComp[2*ix3+1][2*iy] + vComp[2*ix3][2*iy+1] + vComp[2*ix3+1][2*iy+1]);
+ double error = fabs (avgPixel - avgPixelComp);
+ if (error > eMax)
+ eMax = error;
+ }
+ }
+
+ e = eMax;
+
+ return true;
+}
+
+
+bool
+ImageFile::printComparativeStatistics (const ImageFile& imComp, std::ostream& os) const
+{
+ double d, r, e;
+
+ if (comparativeStatistics (imComp, d, r, e)) {
+ os << " Normalized root mean squared distance (d): " << d << std::endl;
+ os << " Normalized mean absolute distance (r): " << r << std::endl;
+ os << "Worst case distance (2x2 pixel average) (e): " << e << std::endl;
+ return true;
+ }
+ return false;
+}
+
+
+void
+ImageFile::printStatistics (std::ostream& os) const
+{
+ double min, max, mean, mode, median, stddev;
+
+ statistics (min, max, mean, mode, median, stddev);
+ if (isComplex())
+ 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;
+
+ if (isComplex()) {
+ 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
+{
+ 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;
+ vecImage.resize (nx * ny);
+ for (int ix = 0; ix < nx; ix++) {
+ for (int iy = 0; iy < ny; iy++)
+ vecImage[iVec++] = v[ix][iy];
+ }
+
+ vectorNumericStatistics (vecImage, nx * ny, min, max, mean, mode, median, stddev);
+}
+
+void
+ImageFile::getMinMax (double& min, double& max) const
+{
+ int nx = m_nx;
+ int ny = m_ny;
+ ImageFileArrayConst v = getArray();
+
+ if (v == NULL || nx == 0 || ny == 0)
+ return;
+
+ min = v[0][0];
+ max = v[0][0];
+ for (int ix = 0; ix < nx; ix++) {
+ for (int iy = 0; iy < ny; iy++) {
+ if (v[ix][iy] > max)
+ max = v[ix][iy];
+ if (v[ix][iy] < min)
+ min = v[ix][iy];
+ }
+ }
+}
+
+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++) {
+ CTSimComplex c (*vRealCol, *vImagCol);
+ *vRealCol++ = std::abs (c);
+ vImagCol++;
+ }
+ }
+
+ return reallocComplexToReal();
+}
+
+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;
+ }
+
+ if (isComplex() || rRHS.isComplex() && ! result.isComplex())
+ result.convertRealToComplex();
+
+ ImageFileArrayConst vLHS = getArray();
+ ImageFileArrayConst vLHSImag = getImaginaryArray();
+ ImageFileArrayConst vRHS = rRHS.getArray();
+ ImageFileArrayConst vRHSImag = rRHS.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++) {
+ vResult[ix][iy] = vLHS[ix][iy] - vRHS[ix][iy];
+ if (result.isComplex()) {
+ vResultImag[ix][iy] = 0;
+ if (isComplex())
+ vResultImag[ix][iy] += vLHSImag[ix][iy];
+ if (rRHS.isComplex())
+ vResultImag[ix][iy] -= vRHSImag[ix][iy];
+ }
+ }
+ }
+
+ 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;
+ }
+
+ if (isComplex() || rRHS.isComplex() && ! result.isComplex())
+ result.convertRealToComplex();
+
+ ImageFileArrayConst vLHS = getArray();
+ ImageFileArrayConst vLHSImag = getImaginaryArray();
+ ImageFileArrayConst vRHS = rRHS.getArray();
+ ImageFileArrayConst vRHSImag = rRHS.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++) {
+ vResult[ix][iy] = vLHS[ix][iy] + vRHS[ix][iy];
+ if (result.isComplex()) {
+ vResultImag[ix][iy] = 0;
+ if (isComplex())
+ vResultImag[ix][iy] += vLHSImag[ix][iy];
+ if (rRHS.isComplex())
+ vResultImag[ix][iy] += vRHSImag[ix][iy];
+ }
+ }
+ }
+
+ 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;
+ }
+
+ if (isComplex() || rRHS.isComplex() && ! result.isComplex())
+ result.convertRealToComplex();
+
+ ImageFileArrayConst vLHS = getArray();
+ ImageFileArrayConst vLHSImag = getImaginaryArray();
+ ImageFileArrayConst vRHS = rRHS.getArray();
+ ImageFileArrayConst vRHSImag = rRHS.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);
+ dImag = 0;
+ if (rRHS.isComplex())
+ dImag = vRHSImag[ix][iy];
+ std::complex<double> cRHS (vRHS[ix][iy], dImag);
+ std::complex<double> cResult = cLHS * cRHS;
+ vResult[ix][iy] = cResult.real();
+ vResultImag[ix][iy] = cResult.imag();
+ } else
+ vResult[ix][iy] = vLHS[ix][iy] * vRHS[ix][iy];
+ }
+ }
+
+
+ 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;
+ }
+
+ if (isComplex() || rRHS.isComplex() && ! result.isComplex())
+ result.convertRealToComplex();
+
+ ImageFileArrayConst vLHS = getArray();
+ ImageFileArrayConst vLHSImag = getImaginaryArray();
+ ImageFileArrayConst vRHS = rRHS.getArray();
+ ImageFileArrayConst vRHSImag = rRHS.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);
+ dImag = 0;
+ if (rRHS.isComplex())
+ dImag = vRHSImag[ix][iy];
+ std::complex<double> cRHS (vRHS[ix][iy], dImag);
+ std::complex<double> cResult = cLHS / cRHS;
+ vResult[ix][iy] = cResult.real();
+ vResultImag[ix][iy] = cResult.imag();
+ } else {
+ if (vRHS != 0)
+ vResult[ix][iy] = vLHS[ix][iy] / vRHS[ix][iy];
+ else
+ vResult[ix][iy] = 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;
+ }
+
+ if (isComplex() && ! result.isComplex())
+ result.convertRealToComplex();
+
+ ImageFileArrayConst vLHS = getArray();
+ ImageFileArray vResult = result.getArray();
+
+ for (unsigned int ix = 0; ix < m_nx; ix++) {
+ ImageFileColumnConst in = vLHS[ix];
+ ImageFileColumn out = vResult[ix];
+ for (unsigned 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;
+ }
+
+ if (isComplex() && ! result.isComplex())
+ result.convertRealToComplex();
+
+ 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 && 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]);
+ }
+ }
+
+
+ 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::log]");
+ return false;
+ }
+
+ if (isComplex() && ! result.isComplex())
+ result.convertRealToComplex();
+
+ ImageFileArrayConst vLHS = getArray();
+ 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::log (cLHS);
+ vResult[ix][iy] = cResult.real();
+ vResultImag[ix][iy] = cResult.imag();
+ } else {
+ if (vLHS[ix][iy] > 0)
+ vResult[ix][iy] = ::log (vLHS[ix][iy]);
+ else
+ vResult[ix][iy] = 0;
+ }
+ }
+ }
+
+
+ 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;
+ }
+
+ if (isComplex() && ! result.isComplex())
+ result.convertRealToComplex();
+
+ ImageFileArrayConst vLHS = getArray();
+ 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::exp (cLHS);
+ vResult[ix][iy] = cResult.real();
+ vResultImag[ix][iy] = cResult.imag();
+ } else
+ vResult[ix][iy] = ::exp (vLHS[ix][iy]);
+ }
+ }
+
+
+ return true;
+}
+
+bool
+ImageFile::scaleImage (ImageFile& result) const
+{
+ unsigned int nx = m_nx;
+ unsigned int ny = m_ny;
+ unsigned int newNX = result.nx();
+ unsigned int newNY = result.ny();
+
+ double dXScale = static_cast<double>(newNX) / static_cast<double>(nx);
+ double dYScale = static_cast<double>(newNY) / static_cast<double>(ny);
+
+ if (isComplex() && ! result.isComplex())
+ result.convertRealToComplex();
+
+ ImageFileArrayConst vReal = getArray();
+ ImageFileArrayConst vImag = getImaginaryArray();
+ ImageFileArray vResult = result.getArray();
+ ImageFileArray vResultImag = result.getImaginaryArray();
+
+ BilinearInterpolator<ImageFileValue> realInterp (vReal, nx, ny);
+ BilinearInterpolator<ImageFileValue> imagInterp (vImag, nx, ny);
+
+ for (unsigned int ix = 0; ix < newNX; ix++) {
+ for (unsigned int iy = 0; iy < newNY; iy++) {
+ double dXPos = ix / dXScale;
+ double dYPos = iy / dYScale;
+ vResult[ix][iy] = realInterp.interpolate (dXPos, dYPos);
+ if (result.isComplex())
+ if (isComplex())
+ vResultImag[ix][iy] = imagInterp.interpolate (dXPos, dYPos);
+ else
+ vResultImag[ix][iy] = 0;
+ }
+ }
+
+ return true;
+}
+
+#ifdef HAVE_FFTW
+bool
+ImageFile::fft (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.convertRealToComplex ())
+ return false;
+ }
+
+ fftw_complex* in = static_cast<fftw_complex*> (fftw_malloc (sizeof(fftw_complex) * m_nx * m_ny));
+
+ ImageFileArrayConst vReal = getArray();
+ ImageFileArrayConst vImag = getImaginaryArray();
+
+ unsigned int ix, iy;
+ unsigned int iArray = 0;
+ for (ix = 0; ix < m_nx; ix++) {
+ for (iy = 0; iy < m_ny; iy++) {
+ in[iArray][0] = vReal[ix][iy];
+ if (isComplex())
+ in[iArray][1] = vImag[ix][iy];
+ else
+ in[iArray][1] = 0;
+ iArray++;
+ }
+ }
+
+ fftw_plan plan = fftw_plan_dft_2d (m_nx, m_ny, in, in, FFTW_FORWARD, FFTW_ESTIMATE);
+ fftw_execute (plan);
+
+ ImageFileArray vRealResult = result.getArray();
+ ImageFileArray vImagResult = result.getImaginaryArray();
+ iArray = 0;
+ unsigned int iScale = m_nx * m_ny;
+ for (ix = 0; ix < m_nx; ix++) {
+ for (iy = 0; iy < m_ny; iy++) {
+ vRealResult[ix][iy] = in[iArray][0] / iScale;
+ vImagResult[ix][iy] = in[iArray][1] / iScale;
+ iArray++;
+ }
+ }
+ fftw_free(in);
+ fftw_destroy_plan (plan);
+
+ Fourier::shuffleFourierToNaturalOrder (result);
+
+ return true;
+}
+
+
+bool
+ImageFile::ifft (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.convertRealToComplex ())
+ return false;
+ }
+
+ ImageFileArrayConst vReal = getArray();
+ ImageFileArrayConst vImag = getImaginaryArray();
+ ImageFileArray vRealResult = result.getArray();
+ ImageFileArray vImagResult = result.getImaginaryArray();
+ unsigned int ix, iy;
+ for (ix = 0; ix < m_nx; ix++) {
+ for (iy = 0; iy < m_ny; iy++) {
+ vRealResult[ix][iy] = vReal[ix][iy];
+ if (isComplex())
+ vImagResult[ix][iy] = vImag[ix][iy];
+ else
+ vImagResult[ix][iy] = 0;
+ }
+ }
+
+ Fourier::shuffleNaturalToFourierOrder (result);
+
+ fftw_complex* in = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * m_nx * m_ny));
+
+ unsigned int iArray = 0;
+ for (ix = 0; ix < m_nx; ix++) {
+ for (iy = 0; iy < m_ny; iy++) {
+ in[iArray][0] = vRealResult[ix][iy];
+ in[iArray][1] = vImagResult[ix][iy];
+ iArray++;
+ }
+ }
+
+ fftw_plan plan = fftw_plan_dft_2d (m_nx, m_ny, in, in, FFTW_BACKWARD, FFTW_ESTIMATE);
+
+ fftw_execute (plan);
+
+ iArray = 0;
+ for (ix = 0; ix < m_nx; ix++) {
+ for (iy = 0; iy < m_ny; iy++) {
+ vRealResult[ix][iy] = in[iArray][0];
+ vImagResult[ix][iy] = in[iArray][1];
+ iArray++;
+ }
+ }
+ fftw_destroy_plan (plan);
+ fftw_free(in);
+
+ return true;
+}
+
+bool
+ImageFile::fftRows (ImageFile& result) const
+{
+ if (m_nx != result.nx() || m_ny != result.ny()) {
+ sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::fftRows]");
+ return false;
+ }
+
+ if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
+ if (! result.convertRealToComplex ())
+ return false;
+ }
+
+ ImageFileArrayConst vReal = getArray();
+ ImageFileArrayConst vImag = getImaginaryArray();
+
+ fftw_complex* in = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * m_nx));
+ fftw_plan plan = fftw_plan_dft_1d (m_nx, in, in, FFTW_FORWARD, FFTW_ESTIMATE);
+
+ std::complex<double>* pcRow = new std::complex<double> [m_nx];
+ for (unsigned int iy = 0; iy < m_ny; iy++) {
+ unsigned int ix;
+ for (ix = 0; ix < m_nx; ix++) {
+ in[ix][0] = vReal[ix][iy];
+ if (isComplex())
+ in[ix][1] = vImag[ix][iy];
+ else
+ in[ix][1] = 0;
+ }
+
+ fftw_execute (plan);
+
+ for (ix = 0; ix < m_nx; ix++)
+ pcRow[ix] = std::complex<double>(in[ix][0], in[ix][1]);
+
+ Fourier::shuffleFourierToNaturalOrder (pcRow, m_nx);
+ for (ix = 0; ix < m_nx; ix++) {
+ vReal[ix][iy] = pcRow[ix].real() / m_nx;
+ vImag[ix][iy] = pcRow[ix].imag() / m_nx;
+ }
+ }
+ delete [] pcRow;
+
+ fftw_destroy_plan (plan);
+ fftw_free(in);
+
+ return true;
+}
+
+bool
+ImageFile::ifftRows (ImageFile& result) const
+{
+ if (m_nx != result.nx() || m_ny != result.ny()) {
+ sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::fftRows]");
+ return false;
+ }
+
+ if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
+ if (! result.convertRealToComplex ())
+ return false;
+ }
+
+ ImageFileArrayConst vReal = getArray();
+ ImageFileArrayConst vImag = getImaginaryArray();
+
+ fftw_complex* in = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * m_nx));
+ fftw_plan plan = fftw_plan_dft_1d (m_nx, in, in, FFTW_BACKWARD, FFTW_ESTIMATE);
+ std::complex<double>* pcRow = new std::complex<double> [m_nx];
+
+ unsigned int ix, iy;
+ // unsigned int iArray = 0;
+ for (iy = 0; iy < m_ny; iy++) {
+ for (ix = 0; ix < m_nx; ix++) {
+ double dImag = 0;
+ if (isComplex())
+ dImag = vImag[ix][iy];
+ pcRow[ix] = std::complex<double> (vReal[ix][iy], dImag);
+ }
+
+ Fourier::shuffleNaturalToFourierOrder (pcRow, m_nx);
+
+ for (ix = 0; ix < m_nx; ix++) {
+ in[ix][0] = pcRow[ix].real();
+ in[ix][1] = pcRow[ix].imag();
+ }
+
+ fftw_execute (plan);
+
+ for (ix = 0; ix < m_nx; ix++) {
+ vReal[ix][iy] = in[ix][0];
+ vImag[ix][iy] = in[ix][1];
+ }
+ }
+ delete [] pcRow;
+
+ fftw_destroy_plan (plan);
+ fftw_free(in);
+
+ return true;
+}
+
+bool
+ImageFile::fftCols (ImageFile& result) const
+{
+ if (m_nx != result.nx() || m_ny != result.ny()) {
+ sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::fftRows]");
+ return false;
+ }
+
+ if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
+ if (! result.convertRealToComplex ())
+ return false;
+ }
+
+ ImageFileArrayConst vReal = getArray();
+ ImageFileArrayConst vImag = getImaginaryArray();
+
+ fftw_complex* in = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * m_ny));
+ fftw_plan plan = fftw_plan_dft_1d (m_ny, in, in, FFTW_FORWARD, FFTW_ESTIMATE);
+
+ std::complex<double>* pcCol = new std::complex<double> [m_ny];
+ for (unsigned int ix = 0; ix < m_nx; ix++) {
+ unsigned int iy;
+ for (iy = 0; iy < m_ny; iy++) {
+ in[iy][0] = vReal[ix][iy];
+ if (isComplex())
+ in[iy][1] = vImag[ix][iy];
+ else
+ in[iy][1] = 0;
+ }
+
+ fftw_execute (plan);
+
+ for (iy = 0; iy < m_ny; iy++)
+ pcCol[iy] = std::complex<double>(in[iy][0], in[iy][1]);
+
+ Fourier::shuffleFourierToNaturalOrder (pcCol, m_ny);
+ for (iy = 0; iy < m_ny; iy++) {
+ vReal[ix][iy] = pcCol[iy].real() / m_ny;
+ vImag[ix][iy] = pcCol[iy].imag() / m_ny;
+ }
+ }
+ delete [] pcCol;
+
+ fftw_destroy_plan (plan);
+ fftw_free(in);
+
+ return true;
+}
+
+bool
+ImageFile::ifftCols (ImageFile& result) const
+{
+ if (m_nx != result.nx() || m_ny != result.ny()) {
+ sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::fftRows]");
+ return false;
+ }
+
+ if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
+ if (! result.convertRealToComplex ())
+ return false;
+ }
+
+ ImageFileArrayConst vReal = getArray();
+ ImageFileArrayConst vImag = getImaginaryArray();
+
+ fftw_complex* in = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * m_ny));
+ fftw_plan plan = fftw_plan_dft_1d (m_ny, in, in, FFTW_BACKWARD, FFTW_ESTIMATE);
+ std::complex<double>* pcCol = new std::complex<double> [m_ny];
+
+ unsigned int ix, iy;
+ // unsigned int iArray = 0;
+ for (ix = 0; ix < m_nx; ix++) {
+ for (iy = 0; iy < m_ny; iy++) {
+ double dImag = 0;
+ if (isComplex())
+ dImag = vImag[ix][iy];
+ pcCol[iy] = std::complex<double> (vReal[ix][iy], dImag);
+ }
+
+ Fourier::shuffleNaturalToFourierOrder (pcCol, m_ny);
+
+ for (iy = 0; iy < m_ny; iy++) {
+ in[iy][0] = pcCol[iy].real();
+ in[iy][1] = pcCol[iy].imag();
+ }
+
+ fftw_execute (plan);
+
+ for (iy = 0; iy < m_ny; iy++) {
+ vReal[ix][iy] = in[iy][0];
+ vImag[ix][iy] = in[iy][1];
+ }
+ }
+ delete [] pcCol;
+
+ fftw_destroy_plan (plan);
+ fftw_free(in);
+
+ return true;
+}
+
+#endif // HAVE_FFTW
+