1 /*****************************************************************************
5 ** Purpose: Imagefile classes
6 ** Programmer: Kevin Rosenberg
7 ** Date Started: June 2000
9 ** This is part of the CTSim program
10 ** Copyright (C) 1983-2000 Kevin Rosenberg
12 ** $Id: imagefile.cpp,v 1.28 2001/01/01 10:14:34 kevin Exp $
14 ** This program is free software; you can redistribute it and/or modify
15 ** it under the terms of the GNU General Public License (version 2) as
16 ** published by the Free Software Foundation.
18 ** This program is distributed in the hope that it will be useful,
19 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
20 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 ** GNU General Public License for more details.
23 ** You should have received a copy of the GNU General Public License
24 ** along with this program; if not, write to the Free Software
25 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26 ******************************************************************************/
31 F32Image::F32Image (int nx, int ny, int dataType)
\r
32 : Array2dFile (nx, ny, sizeof(kfloat32), Array2dFile::PIXEL_FLOAT32, dataType)
\r
36 F32Image::F32Image (void)
\r
39 setPixelFormat (Array2dFile::PIXEL_FLOAT32);
\r
40 setPixelSize (sizeof(kfloat32));
\r
41 setDataType (Array2dFile::DATA_TYPE_REAL);
\r
44 F64Image::F64Image (int nx, int ny, int dataType)
\r
45 : Array2dFile (nx, ny, sizeof(kfloat64), Array2dFile::PIXEL_FLOAT64, dataType)
\r
49 F64Image::F64Image (void)
\r
52 setPixelFormat (PIXEL_FLOAT64);
\r
53 setPixelSize (sizeof(kfloat64));
\r
54 setDataType (Array2dFile::DATA_TYPE_REAL);
\r
58 ImageFile::filterResponse (const char* const domainName, double bw, const char* const filterName, double filt_param)
60 ImageFileArray v = getArray();
\r
61 SignalFilter filter (filterName, domainName, bw, filt_param);
\r
64 int iXCenter, iYCenter;
\r
66 iXCenter = m_nx / 2;
\r
68 iXCenter = (m_nx - 1) / 2;
\r
70 iYCenter = m_ny / 2;
\r
72 iYCenter = (m_ny - 1) / 2;
\r
74 for (int ix = 0; ix < m_nx; ix++)
\r
75 for (int iy = 0; iy < m_ny; iy++) {
\r
76 long lD2 = (ix - iXCenter) * (ix - iXCenter) + (iy - iYCenter) * (iy - iYCenter);
\r
77 double r = ::sqrt (static_cast<double>(lD2));
\r
78 v[ix][iy] = filter.response (r);
\r
81 int hx = (m_nx - 1) / 2;
82 int hy = (m_ny - 1) / 2;
84 for (int i = -hx; i <= hx; i++) {
85 for (int j = -hy; j <= hy; j++) {
86 double r = ::sqrt (i * i + j * j);
88 v[i+hx][j+hy] = filter.response (r);
95 ImageFile::display (void) const
99 getMinMax (pmin, pmax);
101 return (displayScaling (1, pmin, pmax));
105 ImageFile::displayScaling (const int scale, const ImageFileValue pmin, const ImageFileValue pmax) const
109 ImageFileArrayConst v = getArray();
110 if (v == NULL || nx == 0 || ny == 0)
114 int* pPens = new int [nx * ny * scale * scale ];
116 double view_scale = 255 / (pmax - pmin);
117 int id_X11 = g2_open_X11 (nx * scale, ny * scale);
119 for (int i = 0; i < 256; i++) {
120 double cval = i / 255.;
121 grayscale[i] = g2_ink (id_X11, cval, cval, cval);
124 for (int iy = ny - 1; iy >= 0; iy--) {
125 int iRowPos = ((ny - 1 - iy) * scale) * (nx * scale);
126 for (int ix = 0; ix < nx; ix++) {
127 int cval = static_cast<int>((v[ix][iy] - pmin) * view_scale);
132 for (int sy = 0; sy < scale; sy++)
133 for (int sx = 0; sx < scale; sx++)
134 pPens[iRowPos+(sy * nx * scale)+(sx + (ix * scale))] = grayscale[cval];
138 g2_image (id_X11, 0., 0., nx * scale, ny * scale, pPens);
149 // ImageFile::comparativeStatistics Calculate comparative stats
152 // d Normalized root mean squared distance measure
153 // r Normalized mean absolute distance measure
154 // e Worst case distance measure
157 // G.T. Herman, Image Reconstruction From Projections, 1980
160 ImageFile::comparativeStatistics (const ImageFile& imComp, double& d, double& r, double& e) const
162 if (imComp.nx() != m_nx && imComp.ny() != m_ny) {
163 sys_error (ERR_WARNING, "Image sizes differ [ImageFile::comparativeStatistics]");
166 ImageFileArrayConst v = getArray();
167 if (v == NULL || m_nx == 0 || m_ny == 0)
170 ImageFileArrayConst vComp = imComp.getArray();
173 for (unsigned int ix = 0; ix < m_nx; ix++) {
174 for (unsigned int iy = 0; iy < m_ny; iy++) {
178 myMean /= (m_nx * m_ny);
180 double sqErrorSum = 0.;
181 double absErrorSum = 0.;
182 double sqDiffFromMeanSum = 0.;
183 double absValueSum = 0.;
184 for (unsigned int ix2 = 0; ix2 < m_nx; ix2++) {
185 for (unsigned int iy = 0; iy < m_ny; iy++) {
186 double diff = v[ix2][iy] - vComp[ix2][iy];
187 sqErrorSum += diff * diff;
188 absErrorSum += fabs(diff);
189 double diffFromMean = v[ix2][iy] - myMean;
190 sqDiffFromMeanSum += diffFromMean * diffFromMean;
191 absValueSum += fabs(v[ix2][iy]);
195 d = ::sqrt (sqErrorSum / sqDiffFromMeanSum);
196 r = absErrorSum / absValueSum;
201 for (int ix3 = 0; ix3 < hx; ix3++) {
202 for (int iy = 0; iy < hy; iy++) {
203 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]);
204 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]);
205 double error = fabs (avgPixel - avgPixelComp);
218 ImageFile::printComparativeStatistics (const ImageFile& imComp, std::ostream& os) const
222 if (comparativeStatistics (imComp, d, r, e)) {
223 os << " Normalized root mean squared distance (d): " << d << std::endl;
224 os << " Normalized mean absolute distance (r): " << r << std::endl;
225 os << "Worst case distance (2x2 pixel average) (e): " << e << std::endl;
233 ImageFile::printStatistics (std::ostream& os) const
235 double min, max, mean, mode, median, stddev;
237 statistics (min, max, mean, mode, median, stddev);
\r
239 os << "Real Component Statistics" << std::endl;
\r
241 os << " min: " << min << std::endl;
242 os << " max: " << max << std::endl;
243 os << " mean: " << mean << std::endl;
244 os << " mode: " << mode << std::endl;
245 os << "median: " << median << std::endl;
246 os << "stddev: " << stddev << std::endl;
\r
249 statistics (getImaginaryArray(), min, max, mean, mode, median, stddev);
\r
250 os << std::endl << "Imaginary Component Statistics" << std::endl;
251 os << " min: " << min << std::endl;
\r
252 os << " max: " << max << std::endl;
\r
253 os << " mean: " << mean << std::endl;
\r
254 os << " mode: " << mode << std::endl;
\r
255 os << "median: " << median << std::endl;
\r
256 os << "stddev: " << stddev << std::endl;
\r
262 ImageFile::statistics (double& min, double& max, double& mean, double& mode, double& median, double& stddev) const
264 ImageFileArrayConst v = getArray();
\r
265 statistics (v, min, max, mean, mode, median, stddev);
\r
270 ImageFile::statistics (ImageFileArrayConst v, double& min, double& max, double& mean, double& mode, double& median, double& stddev) const
\r
275 if (v == NULL || nx == 0 || ny == 0)
\r
278 std::vector<double> vecImage;
\r
280 vecImage.resize (nx * ny);
\r
281 for (int ix = 0; ix < nx; ix++) {
\r
282 for (int iy = 0; iy < ny; iy++)
\r
283 vecImage[iVec++] = v[ix][iy];
\r
286 vectorNumericStatistics (vecImage, nx * ny, min, max, mean, mode, median, stddev);
\r
290 ImageFile::getMinMax (double& min, double& max) const
294 ImageFileArrayConst v = getArray();
296 if (v == NULL || nx == 0 || ny == 0)
301 for (int ix = 0; ix < nx; ix++) {
302 for (int iy = 0; iy < ny; iy++) {
312 ImageFile::convertRealToComplex ()
\r
314 if (dataType() != Array2dFile::DATA_TYPE_REAL)
\r
317 if (! reallocRealToComplex())
\r
320 ImageFileArray vImag = getImaginaryArray();
\r
321 for (unsigned int ix = 0; ix < m_nx; ix++) {
\r
322 ImageFileColumn vCol = vImag[ix];
\r
323 for (unsigned int iy = 0; iy < m_ny; iy++)
\r
331 ImageFile::convertComplexToReal ()
\r
333 if (dataType() != Array2dFile::DATA_TYPE_COMPLEX)
\r
336 ImageFileArray vReal = getArray();
\r
337 ImageFileArray vImag = getImaginaryArray();
\r
338 for (unsigned int ix = 0; ix < m_nx; ix++) {
\r
339 ImageFileColumn vRealCol = vReal[ix];
\r
340 ImageFileColumn vImagCol = vImag[ix];
\r
341 for (unsigned int iy = 0; iy < m_ny; iy++) {
\r
342 CTSimComplex c (*vRealCol, *vImagCol);
\r
343 *vRealCol++ = std::abs (c);
\r
348 return reallocComplexToReal();
\r
352 ImageFile::subtractImages (const ImageFile& rRHS, ImageFile& result) const
\r
354 if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) {
\r
355 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");
\r
359 ImageFileArrayConst vLHS = getArray();
\r
360 ImageFileArrayConst vRHS = rRHS.getArray();
\r
361 ImageFileArray vResult = result.getArray();
\r
363 for (unsigned int ix = 0; ix < m_nx; ix++) {
\r
364 ImageFileColumnConst in1 = vLHS[ix];
\r
365 ImageFileColumnConst in2 = vRHS[ix];
\r
366 ImageFileColumn out = vResult[ix];
\r
367 for (unsigned int iy = 0; iy < m_ny; iy++)
\r
368 *out++ = *in1++ - *in2++;
\r
375 ImageFile::addImages (const ImageFile& rRHS, ImageFile& result) const
\r
377 if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) {
\r
378 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");
\r
382 ImageFileArrayConst vLHS = getArray();
\r
383 ImageFileArrayConst vRHS = rRHS.getArray();
\r
384 ImageFileArray vResult = result.getArray();
\r
386 for (unsigned int ix = 0; ix < m_nx; ix++) {
\r
387 ImageFileColumnConst in1 = vLHS[ix];
\r
388 ImageFileColumnConst in2 = vRHS[ix];
\r
389 ImageFileColumn out = vResult[ix];
\r
390 for (unsigned int iy = 0; iy < m_ny; iy++)
\r
391 *out++ = *in1++ + *in2++;
\r
398 ImageFile::multiplyImages (const ImageFile& rRHS, ImageFile& result) const
\r
400 if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) {
\r
401 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");
\r
405 ImageFileArrayConst vLHS = getArray();
\r
406 ImageFileArrayConst vRHS = rRHS.getArray();
\r
407 ImageFileArray vResult = result.getArray();
\r
409 for (unsigned int ix = 0; ix < m_nx; ix++) {
\r
410 ImageFileColumnConst in1 = vLHS[ix];
\r
411 ImageFileColumnConst in2 = vRHS[ix];
\r
412 ImageFileColumn out = vResult[ix];
\r
413 for (unsigned int iy = 0; iy < m_ny; iy++)
\r
414 *out++ = *in1++ * *in2++;
\r
421 ImageFile::divideImages (const ImageFile& rRHS, ImageFile& result) const
\r
423 if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) {
\r
424 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");
\r
428 ImageFileArrayConst vLHS = getArray();
\r
429 ImageFileArrayConst vRHS = rRHS.getArray();
\r
430 ImageFileArray vResult = result.getArray();
\r
432 for (unsigned int ix = 0; ix < m_nx; ix++) {
\r
433 ImageFileColumnConst in1 = vLHS[ix];
\r
434 ImageFileColumnConst in2 = vRHS[ix];
\r
435 ImageFileColumn out = vResult[ix];
\r
436 for (unsigned int iy = 0; iy < m_ny; iy++) {
\r
438 *out++ = *in1++ / *in2++;
\r
449 ImageFile::invertPixelValues (ImageFile& result) const
\r
451 if (m_nx != result.nx() || m_ny != result.ny()) {
\r
452 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
\r
456 ImageFileArrayConst vLHS = getArray();
\r
457 ImageFileArray vResult = result.getArray();
\r
459 for (unsigned int ix = 0; ix < m_nx; ix++) {
\r
460 ImageFileColumnConst in = vLHS[ix];
\r
461 ImageFileColumn out = vResult[ix];
\r
462 for (unsigned int iy = 0; iy < m_ny; iy++)
\r
470 ImageFile::sqrt (ImageFile& result) const
\r
472 if (m_nx != result.nx() || m_ny != result.ny()) {
\r
473 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
\r
477 ImageFileArrayConst vLHS = getArray();
\r
478 ImageFileArray vResult = result.getArray();
\r
480 for (unsigned int ix = 0; ix < m_nx; ix++) {
\r
481 ImageFileColumnConst in = vLHS[ix];
\r
482 ImageFileColumn out = vResult[ix];
\r
483 for (unsigned int iy = 0; iy < m_ny; iy++)
\r
485 *out++ = -::sqrt(-*in++);
\r
487 *out++ = ::sqrt(*in++);
\r
494 ImageFile::log (ImageFile& result) const
\r
496 if (m_nx != result.nx() || m_ny != result.ny()) {
\r
497 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
\r
501 ImageFileArrayConst vLHS = getArray();
\r
502 ImageFileArray vResult = result.getArray();
\r
504 for (unsigned int ix = 0; ix < m_nx; ix++) {
\r
505 ImageFileColumnConst in = vLHS[ix];
\r
506 ImageFileColumn out = vResult[ix];
\r
507 for (unsigned int iy = 0; iy < m_ny; iy++)
\r
511 *out++ = ::log(*in++);
\r
518 ImageFile::exp (ImageFile& result) const
\r
520 if (m_nx != result.nx() || m_ny != result.ny()) {
\r
521 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
\r
525 ImageFileArrayConst vLHS = getArray();
\r
526 ImageFileArray vResult = result.getArray();
\r
528 for (unsigned int ix = 0; ix < m_nx; ix++) {
\r
529 ImageFileColumnConst in = vLHS[ix];
\r
530 ImageFileColumn out = vResult[ix];
\r
531 for (unsigned int iy = 0; iy < m_ny; iy++)
\r
532 *out++ = ::exp (*in++);
\r
538 namespace ProcessImage {
\r
541 shuffleFourierToNaturalOrder (ImageFile& im)
\r
543 ImageFileArray vReal = im.getArray();
\r
544 ImageFileArray vImag = im.getImaginaryArray();
\r
545 unsigned int ix, iy;
\r
546 unsigned int nx = im.nx();
\r
547 unsigned int ny = im.ny();
\r
549 // shuffle each column
\r
550 for (ix = 0; ix < nx; ix++) {
\r
551 ProcessSignal::shuffleFourierToNaturalOrder (vReal[ix], ny);
\r
552 if (im.isComplex())
\r
553 ProcessSignal::shuffleFourierToNaturalOrder (vImag[ix], ny);
\r
556 // shuffle each row
\r
557 float* pRow = new float [nx];
\r
558 for (iy = 0; iy < ny; iy++) {
\r
559 for (ix = 0; ix < nx; ix++)
\r
560 pRow[ix] = vReal[ix][iy];
\r
561 ProcessSignal::shuffleFourierToNaturalOrder (pRow, nx);
\r
562 for (ix = 0; ix < nx; ix++)
\r
563 vReal[ix][iy] = pRow[ix];
\r
564 if (im.isComplex()) {
\r
565 for (ix = 0; ix < nx; ix++)
\r
566 pRow[ix] = vImag[ix][iy];
\r
567 ProcessSignal::shuffleFourierToNaturalOrder (pRow, nx);;
\r
568 for (ix = 0; ix < nx; ix++)
\r
569 vImag[ix][iy] = pRow[ix];
\r
576 shuffleNaturalToFourierOrder (ImageFile& im)
\r
578 ImageFileArray vReal = im.getArray();
\r
579 ImageFileArray vImag = im.getImaginaryArray();
\r
580 unsigned int ix, iy;
\r
581 unsigned int nx = im.nx();
\r
582 unsigned int ny = im.ny();
\r
584 // shuffle each x column
\r
585 for (ix = 0; ix < nx; ix++) {
\r
586 ProcessSignal::shuffleNaturalToFourierOrder (vReal[ix], ny);
\r
587 if (im.isComplex())
\r
588 ProcessSignal::shuffleNaturalToFourierOrder (vImag[ix], ny);
\r
591 // shuffle each y row
\r
592 float* pRow = new float [nx];
\r
593 for (iy = 0; iy < ny; iy++) {
\r
594 for (ix = 0; ix < nx; ix++)
\r
595 pRow[ix] = vReal[ix][iy];
\r
596 ProcessSignal::shuffleNaturalToFourierOrder (pRow, nx);
\r
597 for (ix = 0; ix < nx; ix++)
\r
598 vReal[ix][iy] = pRow[ix];
\r
599 if (im.isComplex()) {
\r
600 for (ix = 0; ix < nx; ix++)
\r
601 pRow[ix] = vImag[ix][iy];
\r
602 ProcessSignal::shuffleNaturalToFourierOrder (pRow, nx);
\r
603 for (ix = 0; ix < nx; ix++)
\r
604 vImag[ix][iy] = pRow[ix];
\r
611 }; // namespace ProcessIamge
\r
615 ImageFile::fft (ImageFile& result) const
\r
617 if (m_nx != result.nx() || m_ny != result.ny()) {
\r
618 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
\r
622 if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
\r
623 if (! result.reallocRealToComplex ())
\r
627 fftw_complex* in = new fftw_complex [m_nx * m_ny];
\r
629 ImageFileArrayConst vReal = getArray();
\r
630 ImageFileArrayConst vImag = getImaginaryArray();
\r
632 unsigned int ix, iy;
\r
633 unsigned int iArray = 0;
\r
634 for (ix = 0; ix < m_nx; ix++)
\r
635 for (iy = 0; iy < m_ny; iy++) {
\r
636 in[iArray].re = vReal[ix][iy];
\r
638 in[iArray].im = vImag[ix][iy];
\r
644 fftwnd_plan plan = fftw2d_create_plan (m_nx, m_ny, FFTW_FORWARD, FFTW_IN_PLACE);
\r
646 fftwnd_one (plan, in, NULL);
\r
648 ImageFileArray vRealResult = result.getArray();
\r
649 ImageFileArray vImagResult = result.getImaginaryArray();
\r
651 unsigned int iScale = m_nx * m_ny;
\r
652 for (ix = 0; ix < m_nx; ix++)
\r
653 for (iy = 0; iy < m_ny; iy++) {
\r
654 vRealResult[ix][iy] = in[iArray].re / iScale;
\r
655 vImagResult[ix][iy] = in[iArray].im / iScale;
\r
659 fftwnd_destroy_plan (plan);
\r
663 ProcessImage::shuffleFourierToNaturalOrder (result);
\r
670 ImageFile::ifft (ImageFile& result) const
\r
672 if (m_nx != result.nx() || m_ny != result.ny()) {
\r
673 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
\r
677 if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
\r
678 if (! result.reallocRealToComplex ())
\r
682 ImageFileArrayConst vReal = getArray();
\r
683 ImageFileArrayConst vImag = getImaginaryArray();
\r
684 ImageFileArray vRealResult = result.getArray();
\r
685 ImageFileArray vImagResult = result.getImaginaryArray();
\r
686 unsigned int ix, iy;
\r
687 for (ix = 0; ix < m_nx; ix++)
\r
688 for (iy = 0; iy < m_ny; iy++) {
\r
689 vRealResult[ix][iy] = vReal[ix][iy];
\r
691 vImagResult[ix][iy] = vImag[ix][iy];
\r
693 vImagResult[ix][iy] = 0;
\r
696 ProcessImage::shuffleNaturalToFourierOrder (result);
\r
698 fftw_complex* in = new fftw_complex [m_nx * m_ny];
\r
700 unsigned int iArray = 0;
\r
701 for (ix = 0; ix < m_nx; ix++)
\r
702 for (iy = 0; iy < m_ny; iy++) {
\r
703 in[iArray].re = vRealResult[ix][iy];
\r
704 in[iArray].im = vImagResult[ix][iy];
\r
708 fftwnd_plan plan = fftw2d_create_plan (m_nx, m_ny, FFTW_BACKWARD, FFTW_IN_PLACE);
\r
710 fftwnd_one (plan, in, NULL);
\r
713 for (ix = 0; ix < m_nx; ix++)
\r
714 for (iy = 0; iy < m_ny; iy++) {
\r
715 vRealResult[ix][iy] = in[iArray].re;
\r
716 vImagResult[ix][iy] = in[iArray].im;
\r
720 fftwnd_destroy_plan (plan);
\r
726 #endif // HAVE_FFTW
\r
731 ImageFile::fourier (ImageFile& result) const
\r
733 if (m_nx != result.nx() || m_ny != result.ny()) {
\r
734 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
\r
738 if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
\r
739 if (! result.reallocRealToComplex ())
\r
743 ImageFileArrayConst vLHS = getArray();
\r
744 ImageFileArrayConst vLHSImag = getImaginaryArray();
\r
745 ImageFileArray vRealResult = result.getArray();
\r
746 ImageFileArray vImagResult = result.getImaginaryArray();
\r
748 unsigned int ix, iy;
\r
750 // alloc output matrix
\r
751 CTSimComplex** complexOut = new CTSimComplex* [m_nx];
\r
752 for (ix = 0; ix < m_nx; ix++)
\r
753 complexOut[ix] = new CTSimComplex [m_ny];
\r
755 // fourier each x column
\r
757 CTSimComplex* pY = new CTSimComplex [m_ny];
\r
758 for (ix = 0; ix < m_nx; ix++) {
\r
759 for (iy = 0; iy < m_ny; iy++) {
\r
760 pY[iy].real (vLHS[ix][iy]);
\r
762 pY[iy].imag (vLHSImag[ix][iy]);
\r
766 ProcessSignal::finiteFourierTransform (pY, complexOut[ix], m_ny, ProcessSignal::FORWARD);
\r
769 double* pY = new double [m_ny];
\r
770 for (ix = 0; ix < m_nx; ix++) {
\r
771 for (iy = 0; iy < m_ny; iy++) {
\r
772 pY[iy] = vLHS[ix][iy];
\r
774 ProcessSignal::finiteFourierTransform (pY, complexOut[ix], m_ny, ProcessSignal::FORWARD);
\r
779 // fourier each y row
\r
780 CTSimComplex* pX = new CTSimComplex [m_nx];
\r
781 CTSimComplex* complexOutRow = new CTSimComplex [m_nx];
\r
782 for (iy = 0; iy < m_ny; iy++) {
\r
783 for (ix = 0; ix < m_nx; ix++)
\r
784 pX[ix] = complexOut[ix][iy];
\r
785 ProcessSignal::finiteFourierTransform (pX, complexOutRow, m_nx, ProcessSignal::FORWARD);
\r
786 for (ix = 0; ix < m_nx; ix++)
\r
787 complexOut[ix][iy] = complexOutRow[ix];
\r
791 // shuffle each column
\r
792 for (ix = 0; ix < m_nx; ix++)
\r
793 ProcessSignal::shuffleFourierToNaturalOrder (complexOut[ix], m_ny);
\r
795 // shuffle each row
\r
796 for (iy = 0; iy < m_ny; iy++) {
\r
797 for (ix = 0; ix < m_nx; ix++)
\r
798 complexOutRow[ix] = complexOut[ix][iy];
\r
799 ProcessSignal::shuffleFourierToNaturalOrder (complexOutRow, m_nx);;
\r
800 for (ix = 0; ix < m_nx; ix++)
\r
801 complexOut[ix][iy] = complexOutRow[ix];
\r
804 delete [] complexOutRow;
\r
806 for (ix = 0; ix < m_nx; ix++)
\r
807 for (iy = 0; iy < m_ny; iy++) {
\r
808 vRealResult[ix][iy] = complexOut[ix][iy].real();
\r
809 vImagResult[ix][iy] = complexOut[ix][iy].imag();
\r
812 for (ix = 0; ix < m_nx; ix++)
\r
813 delete [] complexOut[ix];
\r
814 delete [] complexOut;
\r
820 ImageFile::inverseFourier (ImageFile& result) const
\r
822 if (m_nx != result.nx() || m_ny != result.ny()) {
\r
823 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
\r
827 if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
\r
828 if (! result.reallocRealToComplex ())
\r
832 ImageFileArrayConst vLHSReal = getArray();
\r
833 ImageFileArrayConst vLHSImag = getImaginaryArray();
\r
835 unsigned int ix, iy;
\r
837 // alloc 2d complex output matrix
\r
838 CTSimComplex** complexOut = new CTSimComplex* [m_nx];
\r
839 for (ix = 0; ix < m_nx; ix++)
\r
840 complexOut[ix] = new CTSimComplex [m_ny];
\r
842 // put input image into complexOut
\r
843 for (ix = 0; ix < m_nx; ix++)
\r
844 for (iy = 0; iy < m_ny; iy++) {
\r
845 complexOut[ix][iy].real (vLHSReal[ix][iy]);
\r
847 complexOut[ix][iy].imag (vLHSImag[ix][iy]);
\r
849 complexOut[ix][iy].imag (0);
\r
852 // shuffle each x column
\r
853 for (ix = 0; ix < m_nx; ix++)
\r
854 ProcessSignal::shuffleNaturalToFourierOrder (complexOut[ix], m_ny);
\r
856 // shuffle each y row
\r
857 CTSimComplex* pComplexRow = new CTSimComplex [m_nx];
\r
858 for (iy = 0; iy < m_ny; iy++) {
\r
859 for (ix = 0; ix < m_nx; ix++)
\r
860 pComplexRow[ix] = complexOut[ix][iy];
\r
861 ProcessSignal::shuffleNaturalToFourierOrder (pComplexRow, m_nx);
\r
862 for (ix = 0; ix < m_nx; ix++)
\r
863 complexOut[ix][iy] = pComplexRow[ix];
\r
865 delete [] pComplexRow;
\r
867 // ifourier each x column
\r
868 CTSimComplex* pCol = new CTSimComplex [m_ny];
\r
869 for (ix = 0; ix < m_nx; ix++) {
\r
870 for (iy = 0; iy < m_ny; iy++)
\r
871 pCol[iy] = complexOut[ix][iy];
\r
872 ProcessSignal::finiteFourierTransform (pCol, complexOut[ix], m_ny, ProcessSignal::BACKWARD);
\r
876 // ifourier each y row
\r
877 CTSimComplex* complexInRow = new CTSimComplex [m_nx];
\r
878 CTSimComplex* complexOutRow = new CTSimComplex [m_nx];
\r
879 for (iy = 0; iy < m_ny; iy++) {
\r
880 for (ix = 0; ix < m_nx; ix++)
\r
881 complexInRow[ix] = complexOut[ix][iy];
\r
882 ProcessSignal::finiteFourierTransform (complexInRow, complexOutRow, m_nx, ProcessSignal::BACKWARD);
\r
883 for (ix = 0; ix < m_nx; ix++)
\r
884 complexOut[ix][iy] = complexOutRow[ix];
\r
886 delete [] complexInRow;
\r
887 delete [] complexOutRow;
\r
889 ImageFileArray vRealResult = result.getArray();
\r
890 ImageFileArray vImagResult = result.getImaginaryArray();
\r
891 for (ix = 0; ix < m_nx; ix++)
\r
892 for (iy = 0; iy < m_ny; iy++) {
\r
893 vRealResult[ix][iy] = complexOut[ix][iy].real();
\r
894 vImagResult[ix][iy] = complexOut[ix][iy].imag();
\r
897 // delete complexOut matrix
\r
898 for (ix = 0; ix < m_nx; ix++)
\r
899 delete [] complexOut[ix];
\r
900 delete [] complexOut;
\r
907 ImageFile::magnitude (ImageFile& result) const
\r
909 if (m_nx != result.nx() || m_ny != result.ny()) {
\r
910 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
\r
914 ImageFileArray vReal = getArray();
\r
915 ImageFileArray vImag = getImaginaryArray();
\r
917 ImageFileArray vRealResult = result.getArray();
\r
918 for (unsigned int ix = 0; ix < m_nx; ix++)
\r
919 for (unsigned int iy = 0; iy < m_ny; iy++) {
\r
921 vRealResult[ix][iy] = ::sqrt (vReal[ix][iy] * vReal[ix][iy] + vImag[ix][iy] * vImag[ix][iy]);
\r
923 vRealResult[ix][iy] = vReal[ix][iy];
\r
926 if (result.isComplex())
\r
927 result.reallocComplexToReal();
\r
933 ImageFile::phase (ImageFile& result) const
\r
935 if (m_nx != result.nx() || m_ny != result.ny()) {
\r
936 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
\r
940 ImageFileArray vReal = getArray();
\r
941 ImageFileArray vImag = getImaginaryArray();
\r
943 ImageFileArray vRealResult = result.getArray();
\r
944 for (unsigned int ix = 0; ix < m_nx; ix++)
\r
945 for (unsigned int iy = 0; iy < m_ny; iy++) {
\r
947 vRealResult[ix][iy] = ::atan2 (vImag[ix][iy], vReal[ix][iy]);
\r
949 vRealResult[ix][iy] = 0;
\r
952 if (result.isComplex())
\r
953 result.reallocComplexToReal();
\r
959 ImageFile::square (ImageFile& result) const
\r
961 if (m_nx != result.nx() || m_ny != result.ny()) {
\r
962 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
\r
966 ImageFileArrayConst vLHS = getArray();
\r
967 ImageFileArray vResult = result.getArray();
\r
969 for (unsigned int ix = 0; ix < m_nx; ix++) {
\r
970 ImageFileColumnConst in = vLHS[ix];
\r
971 ImageFileColumn out = vResult[ix];
\r
972 for (unsigned int iy = 0; iy < m_ny; iy++) {
\r
973 *out++ = *in * *in;
\r
983 ImageFile::writeImagePGM (const char *outfile, int nxcell, int nycell, double densmin, double densmax)
988 ImageFileArray v = getArray();
990 unsigned char* rowp = new unsigned char [nx * nxcell];
992 if ((fp = fopen (outfile, "wb")) == NULL)
996 fprintf(fp, "%d %d\n", nx, ny);
997 fprintf(fp, "255\n");
999 for (int irow = ny - 1; irow >= 0; irow--) {
1000 for (int icol = 0; icol < nx; icol++) {
1001 int pos = icol * nxcell;
1002 double dens = (v[icol][irow] - densmin) / (densmax - densmin);
1003 dens = clamp (dens, 0., 1.);
1004 for (int p = pos; p < pos + nxcell; p++) {
1005 rowp[p] = static_cast<unsigned int> (dens * 255.);
1008 for (int ir = 0; ir < nycell; ir++) {
1009 for (int ic = 0; ic < nx * nxcell; ic++)
1010 fputc( rowp[ic], fp );
1019 ImageFile::writeImagePGMASCII (const char *outfile, int nxcell, int nycell, double densmin, double densmax)
1024 ImageFileArray v = getArray();
1026 unsigned char* rowp = new unsigned char [nx * nxcell];
1028 if ((fp = fopen (outfile, "wb")) == NULL)
1031 fprintf(fp, "P2\n");
1032 fprintf(fp, "%d %d\n", nx, ny);
1033 fprintf(fp, "255\n");
1035 for (int irow = ny - 1; irow >= 0; irow--) {
1036 for (int icol = 0; icol < nx; icol++) {
1037 int pos = icol * nxcell;
1038 double dens = (v[icol][irow] - densmin) / (densmax - densmin);
1039 dens = clamp (dens, 0., 1.);
1040 for (int p = pos; p < pos + nxcell; p++) {
1041 rowp[p] = static_cast<unsigned int> (dens * 255.);
1044 for (int ir = 0; ir < nycell; ir++) {
1045 for (int ic = 0; ic < nx * nxcell; ic++)
1046 fprintf(fp, "%d ", rowp[ic]);
1058 ImageFile::writeImagePNG (const char *outfile, int bitdepth, int nxcell, int nycell, double densmin, double densmax)
1061 png_structp png_ptr;
1063 double max_out_level = (1 << bitdepth) - 1;
1066 ImageFileArray v = getArray();
1068 unsigned char* rowp = new unsigned char [nx * nxcell * (bitdepth / 8)];
1070 if ((fp = fopen (outfile, "wb")) == NULL)
1073 png_ptr = png_create_write_struct (PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
1077 info_ptr = png_create_info_struct(png_ptr);
1079 png_destroy_write_struct(&png_ptr, (png_infopp) NULL);
1084 if (setjmp(png_ptr->jmpbuf)) {
1085 png_destroy_write_struct(&png_ptr, &info_ptr);
1090 png_init_io(png_ptr, fp);
1092 png_set_IHDR(png_ptr, info_ptr, nx * nxcell, ny * nycell, bitdepth, PNG_COLOR_TYPE_GRAY, PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_DEFAULT);
1094 png_write_info(png_ptr, info_ptr);
1095 for (int irow = ny - 1; irow >= 0; irow--) {
1096 png_bytep row_pointer = rowp;
1098 for (int icol = 0; icol < nx; icol++) {
1099 int pos = icol * nxcell;
1100 double dens = (v[icol][irow] - densmin) / (densmax - densmin);
1101 dens = clamp (dens, 0., 1.);
1102 unsigned int outval = static_cast<unsigned int> (dens * max_out_level);
1104 for (int p = pos; p < pos + nxcell; p++) {
1109 rowp[rowpos] = (outval >> 8) & 0xFF;
1110 rowp[rowpos+1] = (outval & 0xFF);
1114 for (int ir = 0; ir < nycell; ir++)
1115 png_write_rows (png_ptr, &row_pointer, 1);
1118 png_write_end(png_ptr, info_ptr);
1119 png_destroy_write_struct(&png_ptr, &info_ptr);
1128 static const int N_GRAYSCALE=256;
1131 ImageFile::writeImageGIF (const char *outfile, int nxcell, int nycell, double densmin, double densmax)
1135 int gs_indices[N_GRAYSCALE];
1138 ImageFileArray v = getArray();
1140 unsigned char rowp [nx * nxcell];
1144 gif = gdImageCreate(nx * nxcell, ny * nycell);
1145 for (int i = 0; i < N_GRAYSCALE; i++)
1146 gs_indices[i] = gdImageColorAllocate(gif, i, i, i);
1148 int lastrow = ny * nycell - 1;
1149 for (int irow = 0; irow < ny; irow++) {
1150 int rpos = irow * nycell;
1151 for (int ir = rpos; ir < rpos + nycell; ir++) {
1152 for (int icol = 0; icol < nx; icol++) {
1153 int cpos = icol * nxcell;
1154 double dens = (v[icol][irow] - densmin) / (densmax - densmin);
1155 dens = clamp(dens, 0., 1.);
1156 for (int ic = cpos; ic < cpos + nxcell; ic++) {
1157 rowp[ic] = (unsigned int) (dens * (double) (N_GRAYSCALE - 1));
1158 gdImageSetPixel(gif, ic, lastrow - ir, gs_indices[rowp[ic]]);
1164 if ((out = fopen(outfile,"w")) == NULL) {
1165 sys_error(ERR_FATAL, "Error opening output file %s for writing", outfile);
1168 gdImageGif(gif,out);
1170 gdImageDestroy(gif);