** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: imagefile.cpp,v 1.33 2001/01/02 16:02:13 kevin Exp $
+** $Id: imagefile.cpp,v 1.35 2001/01/06 15:33:15 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
}
void
-ImageFile::filterResponse (const char* const domainName, double bw, const char* const filterName, double filt_param, double dInputScale, double dOutputScale)
+ImageFile::getCenterCoordinates (unsigned int& iXCenter, unsigned int& iYCenter)
{
- ImageFileArray v = getArray();
- SignalFilter filter (filterName, domainName, bw, filt_param);
-
- int iXCenter, 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));
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]);
+ 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] = vImag[scaleNX][scaleNY] +
- dXFrac * (vImag[scaleNX+1][scaleNY] - vImag[scaleNX][scaleNY]) +
- dYFrac * (vImag[scaleNX][scaleNY+1] - vImag[scaleNX][scaleNY]);
+ 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;
}
+
+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;
+ }
+
+ fftw_complex* in = new fftw_complex [m_nx];
+
+ ImageFileArrayConst vReal = getArray();
+ ImageFileArrayConst vImag = getImaginaryArray();
+
+ fftw_plan plan = fftw_create_plan (m_nx, FFTW_FORWARD, FFTW_IN_PLACE);
+ 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++) {
+ in[ix].re = vReal[ix][iy];
+ if (isComplex())
+ in[ix].im = vImag[ix][iy];
+ else
+ in[ix].im = 0;
+ }
+
+ fftw_one (plan, in, NULL);
+
+ for (ix = 0; ix < m_nx; ix++)
+ pcRow[ix] = std::complex<double>(in[ix].re, in[ix].im);
+
+ Fourier::shuffleFourierToNaturalOrder (pcRow, m_nx);
+ for (ix = 0; ix < m_nx; ix++) {
+ vReal[ix][iy] = pcRow[ix].real();
+ vImag[ix][iy] = pcRow[ix].imag();
+ }
+ }
+ delete [] pcRow;
+
+ fftw_destroy_plan (plan);
+ delete 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;
+ }
+
+ fftw_complex* in = new fftw_complex [m_nx];
+
+ ImageFileArrayConst vReal = getArray();
+ ImageFileArrayConst vImag = getImaginaryArray();
+
+ fftw_plan plan = fftw_create_plan (m_nx, FFTW_BACKWARD, FFTW_IN_PLACE);
+ 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].re = pcRow[ix].real();
+ in[ix].im = pcRow[ix].imag();
+ }
+
+ fftw_one (plan, in, NULL);
+
+ for (ix = 0; ix < m_nx; ix++) {
+ vReal[ix][iy] = in[ix].re;
+ vImag[ix][iy] = in[ix].im;
+ }
+ }
+ delete [] pcRow;
+
+ fftw_destroy_plan (plan);
+ delete in;
+
+ return true;
+}
+
+bool
+ImageFile::fftCols (ImageFile& result) const
+{
+ return false;
+}
+
+bool
+ImageFile::ifftCols (ImageFile& result) const
+{
+ return false;
+}
+
#endif // HAVE_FFTW