+ 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();
+
+ for (unsigned int ix = 0; ix < newNX; ix++) {
+ for (unsigned int iy = 0; iy < newNY; iy++) {
+ double dXPos = ix / dXScale;
+ double dYPos = iy / dYScale;
+ unsigned int scaleNX = static_cast<unsigned int> (dXPos);
+ unsigned int scaleNY = static_cast<unsigned int> (dYPos);
+ double dXFrac = dXPos - scaleNX;
+ double dYFrac = dYPos - scaleNY;
+ if (scaleNX >= nx - 1 || scaleNY >= ny - 1) {
+ vResult[ix][iy] = vReal[scaleNX][scaleNY];
+ if (result.isComplex()) {
+ if (isComplex())
+ vResultImag[ix][iy] = vImag[scaleNX][scaleNY];
+ else
+ vResultImag[ix][iy] = 0;
+ }
+ } else {
+ vResult[ix][iy] = (1 - dXFrac) * (1 - dYFrac) * vReal[scaleNX][scaleNY] +
+ dXFrac * (1 - dYFrac) * vReal[scaleNX+1][scaleNY] +
+ dYFrac * (1 - dXFrac) * vReal[scaleNX][scaleNY+1] +
+ dXFrac * dYFrac * vReal[scaleNX+1][scaleNY+1];
+ if (result.isComplex()) {
+ if (isComplex())
+ vResultImag[ix][iy] = (1 - dXFrac) * (1 - dYFrac) * vImag[scaleNX][scaleNY] +
+ dXFrac * (1 - dYFrac) * vImag[scaleNX+1][scaleNY] +
+ dYFrac * (1 - dXFrac) * vImag[scaleNX][scaleNY+1] +
+ dXFrac * dYFrac * vImag[scaleNX+1][scaleNY+1];
+ 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 = new 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].re = vReal[ix][iy];
+ if (isComplex())
+ in[iArray].im = vImag[ix][iy];
+ else
+ in[iArray].im = 0;
+ iArray++;
+ }