X-Git-Url: http://git.kpe.io/?a=blobdiff_plain;f=libctsim%2Fimagefile.cpp;h=c21c17f6858facaa7114cd9a92e90b9667dbc35b;hb=32f6cb0e279bb25a06f025676b5513ab792e734c;hp=0cbd67826ae194750b64cb8df9a70f3dbbdaeeac;hpb=9b2bb510160bdb56f04847f5b55ab61dd8a47976;p=ctsim.git diff --git a/libctsim/imagefile.cpp b/libctsim/imagefile.cpp index 0cbd678..c21c17f 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.29 2001/01/02 05:34:57 kevin Exp $ +** $Id: imagefile.cpp,v 1.31 2001/01/02 08:15:02 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 @@ -88,7 +88,7 @@ ImageFile::filterResponse (const char* const domainName, double bw, const char* { ImageFileArray v = getArray(); SignalFilter filter (filterName, domainName, bw, filt_param); - + #if 1 int iXCenter, iYCenter; if (isEven (m_nx)) @@ -99,24 +99,24 @@ ImageFile::filterResponse (const char* const domainName, double bw, const char* iYCenter = m_ny / 2; else iYCenter = (m_ny - 1) / 2; - - for (int ix = 0; ix < m_nx; ix++) - for (int iy = 0; iy < m_ny; iy++) { + + 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(lD2)) * dInputScale; v[ix][iy] = filter.response (r) * dOutputScale; } #else - int hx = (m_nx - 1) / 2; - int hy = (m_ny - 1) / 2; - - for (int i = -hx; i <= hx; i++) { - for (int j = -hy; j <= hy; j++) { - double r = ::sqrt (i * i + j * j); - - v[i+hx][j+hy] = filter.response (r); - } - } + int hx = (m_nx - 1) / 2; + int hy = (m_ny - 1) / 2; + + for (int i = -hx; i <= hx; i++) { + for (int j = -hy; j <= hy; j++) { + double r = ::sqrt (i * i + j * j); + + v[i+hx][j+hy] = filter.response (r); + } + } #endif } @@ -387,7 +387,7 @@ ImageFile::subtractImages (const ImageFile& rRHS, ImageFile& result) const if (isComplex() || rRHS.isComplex() && ! result.isComplex()) result.convertRealToComplex(); - + ImageFileArrayConst vLHS = getArray(); ImageFileArrayConst vLHSImag = getImaginaryArray(); ImageFileArrayConst vRHS = rRHS.getArray(); @@ -407,7 +407,7 @@ ImageFile::subtractImages (const ImageFile& rRHS, ImageFile& result) const } } } - + return true; } @@ -421,7 +421,7 @@ ImageFile::addImages (const ImageFile& rRHS, ImageFile& result) const if (isComplex() || rRHS.isComplex() && ! result.isComplex()) result.convertRealToComplex(); - + ImageFileArrayConst vLHS = getArray(); ImageFileArrayConst vLHSImag = getImaginaryArray(); ImageFileArrayConst vRHS = rRHS.getArray(); @@ -455,7 +455,7 @@ ImageFile::multiplyImages (const ImageFile& rRHS, ImageFile& result) const if (isComplex() || rRHS.isComplex() && ! result.isComplex()) result.convertRealToComplex(); - + ImageFileArrayConst vLHS = getArray(); ImageFileArrayConst vLHSImag = getImaginaryArray(); ImageFileArrayConst vRHS = rRHS.getArray(); @@ -466,12 +466,14 @@ ImageFile::multiplyImages (const ImageFile& rRHS, ImageFile& result) const for (unsigned int ix = 0; ix < m_nx; ix++) { for (unsigned int iy = 0; iy < m_ny; iy++) { if (result.isComplex()) { - std::complex cLHS (vLHS[ix][iy], 0); + double dImag = 0; if (isComplex()) - cLHS.imag (vLHSImag[ix][iy]); - std::complex cRHS (vRHS[ix][iy], 0); + dImag = vLHSImag[ix][iy]; + std::complex cLHS (vLHS[ix][iy], dImag); + dImag = 0; if (rRHS.isComplex()) - cRHS.imag (vRHSImag[ix][iy]); + dImag = vRHSImag[ix][iy]; + std::complex cRHS (vRHS[ix][iy], dImag); std::complex cResult = cLHS * cRHS; vResult[ix][iy] = cResult.real(); vResultImag[ix][iy] = cResult.imag(); @@ -494,7 +496,7 @@ ImageFile::divideImages (const ImageFile& rRHS, ImageFile& result) const if (isComplex() || rRHS.isComplex() && ! result.isComplex()) result.convertRealToComplex(); - + ImageFileArrayConst vLHS = getArray(); ImageFileArrayConst vLHSImag = getImaginaryArray(); ImageFileArrayConst vRHS = rRHS.getArray(); @@ -505,12 +507,14 @@ ImageFile::divideImages (const ImageFile& rRHS, ImageFile& result) const for (unsigned int ix = 0; ix < m_nx; ix++) { for (unsigned int iy = 0; iy < m_ny; iy++) { if (result.isComplex()) { - std::complex cLHS (vLHS[ix][iy], 0); + double dImag = 0; if (isComplex()) - cLHS.imag (vLHSImag[ix][iy]); - std::complex cRHS (vRHS[ix][iy], 0); + dImag = vLHSImag[ix][iy]; + std::complex cLHS (vLHS[ix][iy], dImag); + dImag = 0; if (rRHS.isComplex()) - cRHS.imag (vRHSImag[ix][iy]); + dImag = vRHSImag[ix][iy]; + std::complex cRHS (vRHS[ix][iy], dImag); std::complex cResult = cLHS / cRHS; vResult[ix][iy] = cResult.real(); vResultImag[ix][iy] = cResult.imag(); @@ -523,7 +527,6 @@ ImageFile::divideImages (const ImageFile& rRHS, ImageFile& result) const } } - return true; } @@ -538,7 +541,7 @@ ImageFile::invertPixelValues (ImageFile& result) const if (isComplex() && ! result.isComplex()) result.convertRealToComplex(); - + ImageFileArrayConst vLHS = getArray(); ImageFileArray vResult = result.getArray(); @@ -559,7 +562,7 @@ ImageFile::sqrt (ImageFile& result) const sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]"); return false; } - + if (isComplex() && ! result.isComplex()) result.convertRealToComplex(); @@ -573,27 +576,28 @@ ImageFile::sqrt (ImageFile& result) const 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()) { - std::complex cLHS (vLHS[ix][iy], 0); - if (isComplex()) - cLHS.imag (vLHSImag[ix][iy]); - std::complex cResult = std::sqrt(cLHS); - vResult[ix][iy] = cResult.real(); - vResultImag[ix][iy] = cResult.imag(); - } else - vResult[ix][iy] = ::sqrt (vLHS[ix][iy]); - } - } - - - return true; + + 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 cLHS (vLHS[ix][iy], dImag); + std::complex 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 @@ -606,7 +610,7 @@ ImageFile::log (ImageFile& result) const if (isComplex() && ! result.isComplex()) result.convertRealToComplex(); - + ImageFileArrayConst vLHS = getArray(); ImageFileArrayConst vLHSImag = getImaginaryArray(); ImageFileArray vResult = result.getArray(); @@ -615,9 +619,10 @@ ImageFile::log (ImageFile& result) const for (unsigned int ix = 0; ix < m_nx; ix++) { for (unsigned int iy = 0; iy < m_ny; iy++) { if (result.isComplex()) { - std::complex cLHS (vLHS[ix][iy], 0); + double dImag = 0; if (isComplex()) - cLHS.imag (vLHSImag[ix][iy]); + dImag = vLHSImag[ix][iy]; + std::complex cLHS (vLHS[ix][iy], dImag); std::complex cResult = std::log (cLHS); vResult[ix][iy] = cResult.real(); vResultImag[ix][iy] = cResult.imag(); @@ -637,10 +642,10 @@ ImageFile::exp (ImageFile& result) const 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(); @@ -649,9 +654,10 @@ ImageFile::exp (ImageFile& result) const for (unsigned int ix = 0; ix < m_nx; ix++) { for (unsigned int iy = 0; iy < m_ny; iy++) { if (result.isComplex()) { - std::complex cLHS (vLHS[ix][iy], 0); + double dImag = 0; if (isComplex()) - cLHS.imag (vLHSImag[ix][iy]); + dImag = vLHSImag[ix][iy]; + std::complex cLHS (vLHS[ix][iy], dImag); std::complex cResult = std::exp (cLHS); vResult[ix][iy] = cResult.real(); vResultImag[ix][iy] = cResult.imag(); @@ -671,13 +677,13 @@ ImageFile::scaleImage (ImageFile& result) const unsigned int ny = m_ny; unsigned int newNX = result.nx(); unsigned int newNY = result.ny(); - + double dXScale = static_cast(newNX) / static_cast(nx); double dYScale = static_cast(newNY) / static_cast(ny); - + if (isComplex() && ! result.isComplex()) result.convertRealToComplex(); - + ImageFileArrayConst vReal = getArray(); ImageFileArrayConst vImag = getImaginaryArray(); ImageFileArray vResult = result.getArray(); @@ -685,14 +691,32 @@ ImageFile::scaleImage (ImageFile& result) const for (unsigned int ix = 0; ix < newNX; ix++) { for (unsigned int iy = 0; iy < newNY; iy++) { - unsigned int scaleNX = static_cast (ix / dXScale); - unsigned int scaleNY = static_cast (iy / dYScale); - vResult[ix][iy] = vReal[scaleNX][scaleNY]; - if (result.isComplex()) { - if (isComplex()) - vResultImag[ix][iy] = vImag[scaleNX][scaleNY]; - else - vResultImag[ix][iy] = 0; + double dXPos = ix / dXScale; + double dYPos = iy / dYScale; + unsigned int scaleNX = static_cast (dXPos); + unsigned int scaleNY = static_cast (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] = vReal[scaleNX][scaleNY] + + dXFrac * (vReal[scaleNX+1][scaleNY] - vReal[scaleNX][scaleNY]) + + dYFrac * (vReal[scaleNX][scaleNY+1] - vReal[scaleNX][scaleNY]); + if (result.isComplex()) { + if (isComplex()) + vResultImag[ix][iy] = vImag[scaleNX][scaleNY] + + dXFrac * (vImag[scaleNX+1][scaleNY] - vImag[scaleNX][scaleNY]) + + dYFrac * (vImag[scaleNX][scaleNY+1] - vImag[scaleNX][scaleNY]); + else + vResultImag[ix][iy] = 0; + } } } } @@ -708,17 +732,17 @@ ImageFile::fft (ImageFile& result) const 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++) @@ -730,29 +754,29 @@ ImageFile::fft (ImageFile& result) const in[iArray].im = 0; iArray++; } - - fftwnd_plan plan = fftw2d_create_plan (m_nx, m_ny, FFTW_FORWARD, FFTW_IN_PLACE); - - fftwnd_one (plan, in, NULL); - - 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].re / iScale; - vImagResult[ix][iy] = in[iArray].im / iScale; - iArray++; - } - - fftwnd_destroy_plan (plan); - delete in; - - - Fourier::shuffleFourierToNaturalOrder (result); - - return true; + + fftwnd_plan plan = fftw2d_create_plan (m_nx, m_ny, FFTW_FORWARD, FFTW_IN_PLACE); + + fftwnd_one (plan, in, NULL); + + 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].re / iScale; + vImagResult[ix][iy] = in[iArray].im / iScale; + iArray++; + } + + fftwnd_destroy_plan (plan); + delete in; + + + Fourier::shuffleFourierToNaturalOrder (result); + + return true; } @@ -763,12 +787,12 @@ ImageFile::ifft (ImageFile& result) const 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(); @@ -782,41 +806,41 @@ ImageFile::ifft (ImageFile& result) const else vImagResult[ix][iy] = 0; } - - Fourier::shuffleNaturalToFourierOrder (result); - - fftw_complex* in = new 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].re = vRealResult[ix][iy]; - in[iArray].im = vImagResult[ix][iy]; - iArray++; - } - - fftwnd_plan plan = fftw2d_create_plan (m_nx, m_ny, FFTW_BACKWARD, FFTW_IN_PLACE); - - fftwnd_one (plan, in, NULL); - - iArray = 0; - for (ix = 0; ix < m_nx; ix++) - for (iy = 0; iy < m_ny; iy++) { - vRealResult[ix][iy] = in[iArray].re; - vImagResult[ix][iy] = in[iArray].im; - iArray++; - } - - fftwnd_destroy_plan (plan); - - delete in; - - return true; + + Fourier::shuffleNaturalToFourierOrder (result); + + fftw_complex* in = new 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].re = vRealResult[ix][iy]; + in[iArray].im = vImagResult[ix][iy]; + iArray++; + } + + fftwnd_plan plan = fftw2d_create_plan (m_nx, m_ny, FFTW_BACKWARD, FFTW_IN_PLACE); + + fftwnd_one (plan, in, NULL); + + iArray = 0; + for (ix = 0; ix < m_nx; ix++) + for (iy = 0; iy < m_ny; iy++) { + vRealResult[ix][iy] = in[iArray].re; + vImagResult[ix][iy] = in[iArray].im; + iArray++; + } + + fftwnd_destroy_plan (plan); + + delete in; + + return true; } #endif // HAVE_FFTW - + bool ImageFile::fourier (ImageFile& result) const { @@ -828,70 +852,59 @@ ImageFile::fourier (ImageFile& result) const 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 -#if 1 - CTSimComplex* pY = new CTSimComplex [m_ny]; - for (ix = 0; ix < m_nx; ix++) { - for (iy = 0; iy < m_ny; iy++) { - pY[iy].real (vLHS[ix][iy]); - if (isComplex()) - pY[iy].imag (vLHSImag[ix][iy]); - else - pY[iy].imag (0); - } - ProcessSignal::finiteFourierTransform (pY, complexOut[ix], m_ny, ProcessSignal::FORWARD); - } -#else - double* pY = new 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); - } -#endif - 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); + + 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][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(); + 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(vLHS[ix][iy], dImag); + } + ProcessSignal::finiteFourierTransform (pY, complexOut[ix], m_ny, ProcessSignal::FORWARD); } + delete [] pY; - Fourier::shuffleFourierToNaturalOrder (result); - - // delete complexOut matrix - for (ix = 0; ix < m_nx; ix++) - delete [] complexOut[ix]; - delete [] complexOut; + // 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; - return true; + 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(); + } + + Fourier::shuffleFourierToNaturalOrder (result); + + // delete complexOut matrix + for (ix = 0; ix < m_nx; ix++) + delete [] complexOut[ix]; + delete [] complexOut; + + return true; } bool @@ -917,7 +930,7 @@ ImageFile::inverseFourier (ImageFile& result) const CTSimComplex** complexOut = new CTSimComplex* [m_nx]; for (ix = 0; ix < m_nx; ix++) complexOut[ix] = new CTSimComplex [m_ny]; - + // put input image into result for (ix = 0; ix < m_nx; ix++) for (iy = 0; iy < m_ny; iy++) { @@ -927,45 +940,44 @@ ImageFile::inverseFourier (ImageFile& result) const else vImagResult[ix][iy] = 0; } - - Fourier::shuffleNaturalToFourierOrder (result); - - // ifourier each x column - CTSimComplex* pCol = new CTSimComplex [m_ny]; - for (ix = 0; ix < m_nx; ix++) { - for (iy = 0; iy < m_ny; iy++) { - pCol[iy].real (vRealResult[ix][iy]); - pCol[iy].imag (vImagResult[ix][iy]); + + Fourier::shuffleNaturalToFourierOrder (result); + + // ifourier each x column + CTSimComplex* pCol = new CTSimComplex [m_ny]; + for (ix = 0; ix < m_nx; ix++) { + for (iy = 0; iy < m_ny; iy++) { + pCol[iy] = std::complex (vRealResult[ix][iy], vImagResult[ix][iy]); + } + ProcessSignal::finiteFourierTransform (pCol, complexOut[ix], m_ny, ProcessSignal::BACKWARD); } - ProcessSignal::finiteFourierTransform (pCol, complexOut[ix], m_ny, ProcessSignal::BACKWARD); - } - delete [] pCol; - - // ifourier each y row - CTSimComplex* complexInRow = new CTSimComplex [m_nx]; - CTSimComplex* complexOutRow = new CTSimComplex [m_nx]; - for (iy = 0; iy < m_ny; iy++) { - for (ix = 0; ix < m_nx; ix++) - complexInRow[ix] = complexOut[ix][iy]; - ProcessSignal::finiteFourierTransform (complexInRow, complexOutRow, m_nx, ProcessSignal::BACKWARD); - for (ix = 0; ix < m_nx; ix++) - complexOut[ix][iy] = complexOutRow[ix]; - } - delete [] complexInRow; - delete [] complexOutRow; - - for (ix = 0; ix < m_nx; ix++) + delete [] pCol; + + // ifourier each y row + CTSimComplex* complexInRow = new CTSimComplex [m_nx]; + CTSimComplex* complexOutRow = new CTSimComplex [m_nx]; for (iy = 0; iy < m_ny; iy++) { - vRealResult[ix][iy] = complexOut[ix][iy].real(); - vImagResult[ix][iy] = complexOut[ix][iy].imag(); + for (ix = 0; ix < m_nx; ix++) + complexInRow[ix] = complexOut[ix][iy]; + ProcessSignal::finiteFourierTransform (complexInRow, complexOutRow, m_nx, ProcessSignal::BACKWARD); + for (ix = 0; ix < m_nx; ix++) + complexOut[ix][iy] = complexOutRow[ix]; } + delete [] complexInRow; + delete [] complexOutRow; - // delete complexOut matrix - for (ix = 0; ix < m_nx; ix++) - delete [] complexOut[ix]; - delete [] complexOut; - - return true; + 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(); + } + + // delete complexOut matrix + for (ix = 0; ix < m_nx; ix++) + delete [] complexOut[ix]; + delete [] complexOut; + + return true; } @@ -980,7 +992,7 @@ ImageFile::magnitude (ImageFile& result) const ImageFileArray vReal = getArray(); ImageFileArray vImag = getImaginaryArray(); ImageFileArray vRealResult = result.getArray(); - + for (unsigned int ix = 0; ix < m_nx; ix++) for (unsigned int iy = 0; iy < m_ny; iy++) { if (isComplex()) @@ -988,10 +1000,10 @@ ImageFile::magnitude (ImageFile& result) const else vRealResult[ix][iy] = vReal[ix][iy]; } - - if (result.isComplex()) - result.convertComplexToReal(); - + + if (result.isComplex()) + result.convertComplexToReal(); + return true; } @@ -1006,7 +1018,7 @@ ImageFile::phase (ImageFile& result) const ImageFileArray vReal = getArray(); ImageFileArray vImag = getImaginaryArray(); ImageFileArray vRealResult = result.getArray(); - + for (unsigned int ix = 0; ix < m_nx; ix++) for (unsigned int iy = 0; iy < m_ny; iy++) { if (isComplex()) @@ -1015,10 +1027,10 @@ ImageFile::phase (ImageFile& result) const vRealResult[ix][iy] = 0; } - if (result.isComplex()) - result.convertComplexToReal(); - - return true; + if (result.isComplex()) + result.convertComplexToReal(); + + return true; } bool @@ -1031,7 +1043,7 @@ ImageFile::square (ImageFile& result) const if (isComplex() && ! result.isComplex()) result.convertRealToComplex(); - + ImageFileArrayConst vLHS = getArray(); ImageFileArrayConst vLHSImag = getImaginaryArray(); ImageFileArray vResult = result.getArray(); @@ -1040,9 +1052,10 @@ ImageFile::square (ImageFile& result) const for (unsigned int ix = 0; ix < m_nx; ix++) { for (unsigned int iy = 0; iy < m_ny; iy++) { if (result.isComplex()) { - std::complex cLHS (vLHS[ix][iy], 0); + double dImag = 0; if (isComplex()) - cLHS.imag (vLHSImag[ix][iy]); + dImag = vLHSImag[ix][iy]; + std::complex cLHS (vLHS[ix][iy], dImag); std::complex cResult = cLHS * cLHS; vResult[ix][iy] = cResult.real(); vResultImag[ix][iy] = cResult.imag(); @@ -1059,24 +1072,24 @@ int ImageFile::convertFormatNameToID (const char* const formatName) { int formatID = FORMAT_INVALID; - + for (int i = 0; i < s_iFormatCount; i++) - if (strcasecmp (formatName, s_aszFormatName[i]) == 0) { - formatID = i; - break; - } - - return (formatID); + if (strcasecmp (formatName, s_aszFormatName[i]) == 0) { + formatID = i; + break; + } + + return (formatID); } const char* ImageFile::convertFormatIDToName (int formatID) { static const char *formatName = ""; - + if (formatID >= 0 && formatID < s_iFormatCount) - return (s_aszFormatName[formatID]); - + return (s_aszFormatName[formatID]); + return (formatName); } @@ -1084,10 +1097,10 @@ const char* ImageFile::convertFormatIDToTitle (const int formatID) { static const char *formatTitle = ""; - + if (formatID >= 0 && formatID < s_iFormatCount) - return (s_aszFormatTitle[formatID]); - + return (s_aszFormatTitle[formatID]); + return (formatTitle); } @@ -1099,7 +1112,7 @@ ImageFile::exportImage (const char* const pszFormat, const char* const pszFilena sys_error (ERR_SEVERE, "Invalid format %s [ImageFile::exportImage]", pszFormat); return false; } - + if (iFormatID == FORMAT_PGM) return writeImagePGM (pszFilename, nxcell, nycell, densmin, densmax); else if (iFormatID == FORMAT_PGMASCII) @@ -1108,7 +1121,7 @@ ImageFile::exportImage (const char* const pszFormat, const char* const pszFilena return writeImagePNG (pszFilename, 8, nxcell, nycell, densmin, densmax); else if (iFormatID == FORMAT_PNG16) return writeImagePNG (pszFilename, 16, nxcell, nycell, densmin, densmax); - + sys_error (ERR_SEVERE, "Invalid format %s [ImageFile::exportImage]", pszFormat); return false; } @@ -1147,7 +1160,7 @@ ImageFile::writeImagePGM (const char* const outfile, int nxcell, int nycell, dou delete rowp; fclose(fp); - + return true; } @@ -1186,7 +1199,7 @@ ImageFile::writeImagePGMASCII (const char* const outfile, int nxcell, int nycell delete rowp; fclose(fp); - + return true; } @@ -1256,7 +1269,7 @@ ImageFile::writeImagePNG (const char* const outfile, int bitdepth, int nxcell, i delete rowp; fclose(fp); - + return true; } #endif @@ -1303,9 +1316,9 @@ ImageFile::writeImageGIF (const char* const outfile, int nxcell, int nycell, dou gdImageGif(gif,out); fclose(out); gdImageDestroy(gif); - + delete rowp; - + return true; } #endif