+ 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++;
+ }
+ }
+ delete 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
+
+
+
+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.isComplex())
+ if (! result.convertRealToComplex ())
+ return false;
+
+ ImageFileArrayConst vLHS = getArray();
+ ImageFileArrayConst vLHSImag = getImaginaryArray();
+ ImageFileArray vRealResult = result.getArray();
+ ImageFileArray vImagResult = result.getImaginaryArray();
+
+ unsigned int ix, iy;
+
+ // alloc output matrix
+ CTSimComplex** complexOut = new CTSimComplex* [m_nx];
+ for (ix = 0; ix < m_nx; ix++)
+ complexOut[ix] = new CTSimComplex [m_ny];
+
+ // fourier each x column
+ CTSimComplex* pY = new CTSimComplex [m_ny];
+ for (ix = 0; ix < m_nx; ix++) {
+ for (iy = 0; iy < m_ny; iy++) {
+ double dImag = 0;
+ if (isComplex())
+ dImag = vLHSImag[ix][iy];
+ pY[iy] = std::complex<double>(vLHS[ix][iy], dImag);
+ }
+ ProcessSignal::finiteFourierTransform (pY, complexOut[ix], m_ny, ProcessSignal::FORWARD);
+ }
+ delete [] pY;
+
+ // fourier each y row
+ CTSimComplex* pX = new CTSimComplex [m_nx];
+ CTSimComplex* complexOutRow = new CTSimComplex [m_nx];
+ for (iy = 0; iy < m_ny; iy++) {
+ for (ix = 0; ix < m_nx; ix++)
+ pX[ix] = complexOut[ix][iy];
+ ProcessSignal::finiteFourierTransform (pX, complexOutRow, m_nx, ProcessSignal::FORWARD);
+ for (ix = 0; ix < m_nx; ix++)
+ complexOut[ix][iy] = complexOutRow[ix];
+ }
+ delete [] pX;
+ delete [] complexOutRow;
+
+ for (ix = 0; ix < m_nx; ix++)
+ for (iy = 0; iy < m_ny; iy++) {
+ vRealResult[ix][iy] = complexOut[ix][iy].real();
+ vImagResult[ix][iy] = complexOut[ix][iy].imag();