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-2001 Kevin Rosenberg
12 ** $Id: imagefile.cpp,v 1.45 2001/09/24 09:40:42 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 ******************************************************************************/
33 const double ImageFile::s_dRedGrayscaleFactor = 0.299;
34 const double ImageFile::s_dGreenGrayscaleFactor = 0.587;
35 const double ImageFile::s_dBlueGrayscaleFactor = 0.114;
38 const int ImageFile::EXPORT_FORMAT_INVALID = -1;
39 const int ImageFile::EXPORT_FORMAT_TEXT = 0;
40 const int ImageFile::EXPORT_FORMAT_PGM = 1;
41 const int ImageFile::EXPORT_FORMAT_PGMASCII = 2;
43 const int ImageFile::EXPORT_FORMAT_PNG = 3;
44 const int ImageFile::EXPORT_FORMAT_PNG16 = 4;
47 const int ImageFile::EXPORT_FORMAT_DICOM = 5;
49 const int ImageFile::EXPORT_FORMAT_RAW = 6;
51 const char* ImageFile::s_aszExportFormatName[] =
65 const char* ImageFile::s_aszExportFormatTitle[] =
76 const int ImageFile::s_iExportFormatCount = sizeof(s_aszExportFormatName) / sizeof(const char*);
79 const int ImageFile::IMPORT_FORMAT_INVALID = -1;
80 const int ImageFile::IMPORT_FORMAT_PPM = 0;
82 const int ImageFile::IMPORT_FORMAT_PNG = 1;
85 const int ImageFile::IMPORT_FORMAT_DICOM = 2;
89 const char* ImageFile::s_aszImportFormatName[] =
100 const char* ImageFile::s_aszImportFormatTitle[] =
104 #ifdef HAVE_CTN_DICOM
108 const int ImageFile::s_iImportFormatCount = sizeof(s_aszImportFormatName) / sizeof(const char*);
112 F32Image::F32Image (int nx, int ny, int dataType)
113 : Array2dFile (nx, ny, sizeof(kfloat32), Array2dFile::PIXEL_FLOAT32, dataType)
117 F32Image::F32Image (void)
120 setPixelFormat (Array2dFile::PIXEL_FLOAT32);
121 setPixelSize (sizeof(kfloat32));
122 setDataType (Array2dFile::DATA_TYPE_REAL);
125 F64Image::F64Image (int nx, int ny, int dataType)
126 : Array2dFile (nx, ny, sizeof(kfloat64), Array2dFile::PIXEL_FLOAT64, dataType)
130 F64Image::F64Image (void)
133 setPixelFormat (PIXEL_FLOAT64);
134 setPixelSize (sizeof(kfloat64));
135 setDataType (Array2dFile::DATA_TYPE_REAL);
139 ImageFile::getCenterCoordinates (unsigned int& iXCenter, unsigned int& iYCenter)
144 iXCenter = (m_nx - 1) / 2;
149 iYCenter = (m_ny - 1) / 2;
154 ImageFile::filterResponse (const char* const domainName, double bw, const char* const filterName,
155 double filt_param, double dInputScale, double dOutputScale)
157 ImageFileArray v = getArray();
158 SignalFilter filter (filterName, domainName, bw, filt_param);
160 unsigned int iXCenter, iYCenter;
161 getCenterCoordinates (iXCenter, iYCenter);
163 for (unsigned int ix = 0; ix < m_nx; ix++)
164 for (unsigned int iy = 0; iy < m_ny; iy++) {
165 long lD2 = ((ix - iXCenter) * (ix - iXCenter)) + ((iy - iYCenter) * (iy - iYCenter));
166 double r = ::sqrt (static_cast<double>(lD2)) * dInputScale;
167 v[ix][iy] = filter.response (r) * dOutputScale;
172 // ImageFile::comparativeStatistics Calculate comparative stats
175 // d Normalized root mean squared distance measure
176 // r Normalized mean absolute distance measure
177 // e Worst case distance measure
180 // G.T. Herman, Image Reconstruction From Projections, 1980
183 ImageFile::comparativeStatistics (const ImageFile& imComp, double& d, double& r, double& e) const
185 if (imComp.nx() != m_nx && imComp.ny() != m_ny) {
186 sys_error (ERR_WARNING, "Image sizes differ [ImageFile::comparativeStatistics]");
189 ImageFileArrayConst v = getArray();
190 if (v == NULL || m_nx == 0 || m_ny == 0)
193 ImageFileArrayConst vComp = imComp.getArray();
196 for (unsigned int ix = 0; ix < m_nx; ix++) {
197 for (unsigned int iy = 0; iy < m_ny; iy++) {
201 myMean /= (m_nx * m_ny);
203 double sqErrorSum = 0.;
204 double absErrorSum = 0.;
205 double sqDiffFromMeanSum = 0.;
206 double absValueSum = 0.;
207 for (unsigned int ix2 = 0; ix2 < m_nx; ix2++) {
208 for (unsigned int iy = 0; iy < m_ny; iy++) {
209 double diff = v[ix2][iy] - vComp[ix2][iy];
210 sqErrorSum += diff * diff;
211 absErrorSum += fabs(diff);
212 double diffFromMean = v[ix2][iy] - myMean;
213 sqDiffFromMeanSum += diffFromMean * diffFromMean;
214 absValueSum += fabs(v[ix2][iy]);
218 d = ::sqrt (sqErrorSum / sqDiffFromMeanSum);
219 r = absErrorSum / absValueSum;
224 for (int ix3 = 0; ix3 < hx; ix3++) {
225 for (int iy = 0; iy < hy; iy++) {
226 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]);
227 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]);
228 double error = fabs (avgPixel - avgPixelComp);
241 ImageFile::printComparativeStatistics (const ImageFile& imComp, std::ostream& os) const
245 if (comparativeStatistics (imComp, d, r, e)) {
246 os << " Normalized root mean squared distance (d): " << d << std::endl;
247 os << " Normalized mean absolute distance (r): " << r << std::endl;
248 os << "Worst case distance (2x2 pixel average) (e): " << e << std::endl;
256 ImageFile::printStatistics (std::ostream& os) const
258 double min, max, mean, mode, median, stddev;
260 statistics (min, max, mean, mode, median, stddev);
262 os << "Real Component Statistics" << std::endl;
264 os << " min: " << min << std::endl;
265 os << " max: " << max << std::endl;
266 os << " mean: " << mean << std::endl;
267 os << " mode: " << mode << std::endl;
268 os << "median: " << median << std::endl;
269 os << "stddev: " << stddev << std::endl;
272 statistics (getImaginaryArray(), min, max, mean, mode, median, stddev);
273 os << std::endl << "Imaginary Component Statistics" << std::endl;
274 os << " min: " << min << std::endl;
275 os << " max: " << max << std::endl;
276 os << " mean: " << mean << std::endl;
277 os << " mode: " << mode << std::endl;
278 os << "median: " << median << std::endl;
279 os << "stddev: " << stddev << std::endl;
285 ImageFile::statistics (double& min, double& max, double& mean, double& mode, double& median, double& stddev) const
287 ImageFileArrayConst v = getArray();
288 statistics (v, min, max, mean, mode, median, stddev);
293 ImageFile::statistics (ImageFileArrayConst v, double& min, double& max, double& mean, double& mode, double& median, double& stddev) const
298 if (v == NULL || nx == 0 || ny == 0)
301 std::vector<double> vecImage;
303 vecImage.resize (nx * ny);
304 for (int ix = 0; ix < nx; ix++) {
305 for (int iy = 0; iy < ny; iy++)
306 vecImage[iVec++] = v[ix][iy];
309 vectorNumericStatistics (vecImage, nx * ny, min, max, mean, mode, median, stddev);
313 ImageFile::getMinMax (double& min, double& max) const
317 ImageFileArrayConst v = getArray();
319 if (v == NULL || nx == 0 || ny == 0)
324 for (int ix = 0; ix < nx; ix++) {
325 for (int iy = 0; iy < ny; iy++) {
335 ImageFile::convertRealToComplex ()
337 if (dataType() != Array2dFile::DATA_TYPE_REAL)
340 if (! reallocRealToComplex())
343 ImageFileArray vImag = getImaginaryArray();
344 for (unsigned int ix = 0; ix < m_nx; ix++) {
345 ImageFileColumn vCol = vImag[ix];
346 for (unsigned int iy = 0; iy < m_ny; iy++)
354 ImageFile::convertComplexToReal ()
356 if (dataType() != Array2dFile::DATA_TYPE_COMPLEX)
359 ImageFileArray vReal = getArray();
360 ImageFileArray vImag = getImaginaryArray();
361 for (unsigned int ix = 0; ix < m_nx; ix++) {
362 ImageFileColumn vRealCol = vReal[ix];
363 ImageFileColumn vImagCol = vImag[ix];
364 for (unsigned int iy = 0; iy < m_ny; iy++) {
365 CTSimComplex c (*vRealCol, *vImagCol);
366 *vRealCol++ = std::abs (c);
371 return reallocComplexToReal();
375 ImageFile::subtractImages (const ImageFile& rRHS, ImageFile& result) const
377 if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) {
378 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");
382 if (isComplex() || rRHS.isComplex() && ! result.isComplex())
383 result.convertRealToComplex();
385 ImageFileArrayConst vLHS = getArray();
386 ImageFileArrayConst vLHSImag = getImaginaryArray();
387 ImageFileArrayConst vRHS = rRHS.getArray();
388 ImageFileArrayConst vRHSImag = rRHS.getImaginaryArray();
389 ImageFileArray vResult = result.getArray();
390 ImageFileArray vResultImag = result.getImaginaryArray();
392 for (unsigned int ix = 0; ix < m_nx; ix++) {
393 for (unsigned int iy = 0; iy < m_ny; iy++) {
394 vResult[ix][iy] = vLHS[ix][iy] - vRHS[ix][iy];
395 if (result.isComplex()) {
396 vResultImag[ix][iy] = 0;
398 vResultImag[ix][iy] += vLHSImag[ix][iy];
399 if (rRHS.isComplex())
400 vResultImag[ix][iy] -= vRHSImag[ix][iy];
409 ImageFile::addImages (const ImageFile& rRHS, ImageFile& result) const
411 if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) {
412 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");
416 if (isComplex() || rRHS.isComplex() && ! result.isComplex())
417 result.convertRealToComplex();
419 ImageFileArrayConst vLHS = getArray();
420 ImageFileArrayConst vLHSImag = getImaginaryArray();
421 ImageFileArrayConst vRHS = rRHS.getArray();
422 ImageFileArrayConst vRHSImag = rRHS.getImaginaryArray();
423 ImageFileArray vResult = result.getArray();
424 ImageFileArray vResultImag = result.getImaginaryArray();
426 for (unsigned int ix = 0; ix < m_nx; ix++) {
427 for (unsigned int iy = 0; iy < m_ny; iy++) {
428 vResult[ix][iy] = vLHS[ix][iy] + vRHS[ix][iy];
429 if (result.isComplex()) {
430 vResultImag[ix][iy] = 0;
432 vResultImag[ix][iy] += vLHSImag[ix][iy];
433 if (rRHS.isComplex())
434 vResultImag[ix][iy] += vRHSImag[ix][iy];
443 ImageFile::multiplyImages (const ImageFile& rRHS, ImageFile& result) const
445 if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) {
446 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");
450 if (isComplex() || rRHS.isComplex() && ! result.isComplex())
451 result.convertRealToComplex();
453 ImageFileArrayConst vLHS = getArray();
454 ImageFileArrayConst vLHSImag = getImaginaryArray();
455 ImageFileArrayConst vRHS = rRHS.getArray();
456 ImageFileArrayConst vRHSImag = rRHS.getImaginaryArray();
457 ImageFileArray vResult = result.getArray();
458 ImageFileArray vResultImag = result.getImaginaryArray();
460 for (unsigned int ix = 0; ix < m_nx; ix++) {
461 for (unsigned int iy = 0; iy < m_ny; iy++) {
462 if (result.isComplex()) {
465 dImag = vLHSImag[ix][iy];
466 std::complex<double> cLHS (vLHS[ix][iy], dImag);
468 if (rRHS.isComplex())
469 dImag = vRHSImag[ix][iy];
470 std::complex<double> cRHS (vRHS[ix][iy], dImag);
471 std::complex<double> cResult = cLHS * cRHS;
472 vResult[ix][iy] = cResult.real();
473 vResultImag[ix][iy] = cResult.imag();
475 vResult[ix][iy] = vLHS[ix][iy] * vRHS[ix][iy];
484 ImageFile::divideImages (const ImageFile& rRHS, ImageFile& result) const
486 if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) {
487 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");
491 if (isComplex() || rRHS.isComplex() && ! result.isComplex())
492 result.convertRealToComplex();
494 ImageFileArrayConst vLHS = getArray();
495 ImageFileArrayConst vLHSImag = getImaginaryArray();
496 ImageFileArrayConst vRHS = rRHS.getArray();
497 ImageFileArrayConst vRHSImag = rRHS.getImaginaryArray();
498 ImageFileArray vResult = result.getArray();
499 ImageFileArray vResultImag = result.getImaginaryArray();
501 for (unsigned int ix = 0; ix < m_nx; ix++) {
502 for (unsigned int iy = 0; iy < m_ny; iy++) {
503 if (result.isComplex()) {
506 dImag = vLHSImag[ix][iy];
507 std::complex<double> cLHS (vLHS[ix][iy], dImag);
509 if (rRHS.isComplex())
510 dImag = vRHSImag[ix][iy];
511 std::complex<double> cRHS (vRHS[ix][iy], dImag);
512 std::complex<double> cResult = cLHS / cRHS;
513 vResult[ix][iy] = cResult.real();
514 vResultImag[ix][iy] = cResult.imag();
517 vResult[ix][iy] = vLHS[ix][iy] / vRHS[ix][iy];
529 ImageFile::invertPixelValues (ImageFile& result) const
531 if (m_nx != result.nx() || m_ny != result.ny()) {
532 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
536 if (isComplex() && ! result.isComplex())
537 result.convertRealToComplex();
539 ImageFileArrayConst vLHS = getArray();
540 ImageFileArray vResult = result.getArray();
542 for (unsigned int ix = 0; ix < m_nx; ix++) {
543 ImageFileColumnConst in = vLHS[ix];
544 ImageFileColumn out = vResult[ix];
545 for (unsigned int iy = 0; iy < m_ny; iy++)
553 ImageFile::sqrt (ImageFile& result) const
555 if (m_nx != result.nx() || m_ny != result.ny()) {
556 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
560 if (isComplex() && ! result.isComplex())
561 result.convertRealToComplex();
563 bool bComplexOutput = result.isComplex();
564 ImageFileArrayConst vLHS = getArray();
565 if (! bComplexOutput) // check if should convert to complex output
566 for (unsigned int ix = 0; ix < m_nx; ix++)
567 for (unsigned int iy = 0; iy < m_ny; iy++)
568 if (! bComplexOutput && vLHS[ix][iy] < 0) {
569 result.convertRealToComplex();
570 bComplexOutput = true;
574 ImageFileArrayConst vLHSImag = getImaginaryArray();
575 ImageFileArray vResult = result.getArray();
576 ImageFileArray vResultImag = result.getImaginaryArray();
578 for (unsigned int ix = 0; ix < m_nx; ix++) {
579 for (unsigned int iy = 0; iy < m_ny; iy++) {
580 if (result.isComplex()) {
583 dImag = vLHSImag[ix][iy];
584 std::complex<double> cLHS (vLHS[ix][iy], dImag);
585 std::complex<double> cResult = std::sqrt(cLHS);
586 vResult[ix][iy] = cResult.real();
587 vResultImag[ix][iy] = cResult.imag();
589 vResult[ix][iy] = ::sqrt (vLHS[ix][iy]);
598 ImageFile::log (ImageFile& result) const
600 if (m_nx != result.nx() || m_ny != result.ny()) {
601 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::log]");
605 if (isComplex() && ! result.isComplex())
606 result.convertRealToComplex();
608 ImageFileArrayConst vLHS = getArray();
609 ImageFileArrayConst vLHSImag = getImaginaryArray();
610 ImageFileArray vResult = result.getArray();
611 ImageFileArray vResultImag = result.getImaginaryArray();
613 for (unsigned int ix = 0; ix < m_nx; ix++) {
614 for (unsigned int iy = 0; iy < m_ny; iy++) {
615 if (result.isComplex()) {
618 dImag = vLHSImag[ix][iy];
619 std::complex<double> cLHS (vLHS[ix][iy], dImag);
620 std::complex<double> cResult = std::log (cLHS);
621 vResult[ix][iy] = cResult.real();
622 vResultImag[ix][iy] = cResult.imag();
624 if (vLHS[ix][iy] > 0)
625 vResult[ix][iy] = ::log (vLHS[ix][iy]);
637 ImageFile::exp (ImageFile& result) const
639 if (m_nx != result.nx() || m_ny != result.ny()) {
640 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
644 if (isComplex() && ! result.isComplex())
645 result.convertRealToComplex();
647 ImageFileArrayConst vLHS = getArray();
648 ImageFileArrayConst vLHSImag = getImaginaryArray();
649 ImageFileArray vResult = result.getArray();
650 ImageFileArray vResultImag = result.getImaginaryArray();
652 for (unsigned int ix = 0; ix < m_nx; ix++) {
653 for (unsigned int iy = 0; iy < m_ny; iy++) {
654 if (result.isComplex()) {
657 dImag = vLHSImag[ix][iy];
658 std::complex<double> cLHS (vLHS[ix][iy], dImag);
659 std::complex<double> cResult = std::exp (cLHS);
660 vResult[ix][iy] = cResult.real();
661 vResultImag[ix][iy] = cResult.imag();
663 vResult[ix][iy] = ::exp (vLHS[ix][iy]);
672 ImageFile::scaleImage (ImageFile& result) const
674 unsigned int nx = m_nx;
675 unsigned int ny = m_ny;
676 unsigned int newNX = result.nx();
677 unsigned int newNY = result.ny();
679 double dXScale = static_cast<double>(newNX) / static_cast<double>(nx);
680 double dYScale = static_cast<double>(newNY) / static_cast<double>(ny);
682 if (isComplex() && ! result.isComplex())
683 result.convertRealToComplex();
685 ImageFileArrayConst vReal = getArray();
686 ImageFileArrayConst vImag = getImaginaryArray();
687 ImageFileArray vResult = result.getArray();
688 ImageFileArray vResultImag = result.getImaginaryArray();
690 BilinearInterpolator<ImageFileValue> realInterp (vReal, nx, ny);
691 BilinearInterpolator<ImageFileValue> imagInterp (vImag, nx, ny);
693 for (unsigned int ix = 0; ix < newNX; ix++) {
694 for (unsigned int iy = 0; iy < newNY; iy++) {
695 double dXPos = ix / dXScale;
696 double dYPos = iy / dYScale;
697 vResult[ix][iy] = realInterp.interpolate (dXPos, dYPos);
698 if (result.isComplex())
700 vResultImag[ix][iy] = imagInterp.interpolate (dXPos, dYPos);
702 vResultImag[ix][iy] = 0;
711 ImageFile::fft (ImageFile& result) const
713 if (m_nx != result.nx() || m_ny != result.ny()) {
714 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
718 if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
719 if (! result.convertRealToComplex ())
723 fftw_complex* in = new fftw_complex [m_nx * m_ny];
725 ImageFileArrayConst vReal = getArray();
726 ImageFileArrayConst vImag = getImaginaryArray();
729 unsigned int iArray = 0;
730 for (ix = 0; ix < m_nx; ix++) {
731 for (iy = 0; iy < m_ny; iy++) {
732 in[iArray].re = vReal[ix][iy];
734 in[iArray].im = vImag[ix][iy];
741 fftwnd_plan plan = fftw2d_create_plan (m_nx, m_ny, FFTW_FORWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
742 fftwnd_one (plan, in, NULL);
744 ImageFileArray vRealResult = result.getArray();
745 ImageFileArray vImagResult = result.getImaginaryArray();
747 unsigned int iScale = m_nx * m_ny;
748 for (ix = 0; ix < m_nx; ix++) {
749 for (iy = 0; iy < m_ny; iy++) {
750 vRealResult[ix][iy] = in[iArray].re / iScale;
751 vImagResult[ix][iy] = in[iArray].im / iScale;
756 fftwnd_destroy_plan (plan);
758 Fourier::shuffleFourierToNaturalOrder (result);
765 ImageFile::ifft (ImageFile& result) const
767 if (m_nx != result.nx() || m_ny != result.ny()) {
768 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
772 if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
773 if (! result.convertRealToComplex ())
777 ImageFileArrayConst vReal = getArray();
778 ImageFileArrayConst vImag = getImaginaryArray();
779 ImageFileArray vRealResult = result.getArray();
780 ImageFileArray vImagResult = result.getImaginaryArray();
782 for (ix = 0; ix < m_nx; ix++) {
783 for (iy = 0; iy < m_ny; iy++) {
784 vRealResult[ix][iy] = vReal[ix][iy];
786 vImagResult[ix][iy] = vImag[ix][iy];
788 vImagResult[ix][iy] = 0;
792 Fourier::shuffleNaturalToFourierOrder (result);
794 fftw_complex* in = new fftw_complex [m_nx * m_ny];
796 unsigned int iArray = 0;
797 for (ix = 0; ix < m_nx; ix++) {
798 for (iy = 0; iy < m_ny; iy++) {
799 in[iArray].re = vRealResult[ix][iy];
800 in[iArray].im = vImagResult[ix][iy];
805 fftwnd_plan plan = fftw2d_create_plan (m_nx, m_ny, FFTW_BACKWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
807 fftwnd_one (plan, in, NULL);
810 for (ix = 0; ix < m_nx; ix++) {
811 for (iy = 0; iy < m_ny; iy++) {
812 vRealResult[ix][iy] = in[iArray].re;
813 vImagResult[ix][iy] = in[iArray].im;
817 fftwnd_destroy_plan (plan);
825 ImageFile::fftRows (ImageFile& result) const
827 if (m_nx != result.nx() || m_ny != result.ny()) {
828 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::fftRows]");
832 if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
833 if (! result.convertRealToComplex ())
837 ImageFileArrayConst vReal = getArray();
838 ImageFileArrayConst vImag = getImaginaryArray();
840 fftw_plan plan = fftw_create_plan (m_nx, FFTW_FORWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
842 fftw_complex* in = new fftw_complex [m_nx];
843 std::complex<double>* pcRow = new std::complex<double> [m_nx];
844 for (int iy = 0; iy < m_ny; iy++) {
846 for (ix = 0; ix < m_nx; ix++) {
847 in[ix].re = vReal[ix][iy];
849 in[ix].im = vImag[ix][iy];
854 fftw_one (plan, in, NULL);
856 for (ix = 0; ix < m_nx; ix++)
857 pcRow[ix] = std::complex<double>(in[ix].re, in[ix].im);
859 Fourier::shuffleFourierToNaturalOrder (pcRow, m_nx);
860 for (ix = 0; ix < m_nx; ix++) {
861 vReal[ix][iy] = pcRow[ix].real() / m_nx;
862 vImag[ix][iy] = pcRow[ix].imag() / m_nx;
867 fftw_destroy_plan (plan);
874 ImageFile::ifftRows (ImageFile& result) const
876 if (m_nx != result.nx() || m_ny != result.ny()) {
877 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::fftRows]");
881 if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
882 if (! result.convertRealToComplex ())
886 fftw_complex* in = new fftw_complex [m_nx];
888 ImageFileArrayConst vReal = getArray();
889 ImageFileArrayConst vImag = getImaginaryArray();
891 fftw_plan plan = fftw_create_plan (m_nx, FFTW_BACKWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
892 std::complex<double>* pcRow = new std::complex<double> [m_nx];
895 unsigned int iArray = 0;
896 for (iy = 0; iy < m_ny; iy++) {
897 for (ix = 0; ix < m_nx; ix++) {
900 dImag = vImag[ix][iy];
901 pcRow[ix] = std::complex<double> (vReal[ix][iy], dImag);
904 Fourier::shuffleNaturalToFourierOrder (pcRow, m_nx);
906 for (ix = 0; ix < m_nx; ix++) {
907 in[ix].re = pcRow[ix].real();
908 in[ix].im = pcRow[ix].imag();
911 fftw_one (plan, in, NULL);
913 for (ix = 0; ix < m_nx; ix++) {
914 vReal[ix][iy] = in[ix].re;
915 vImag[ix][iy] = in[ix].im;
920 fftw_destroy_plan (plan);
927 ImageFile::fftCols (ImageFile& result) const
929 if (m_nx != result.nx() || m_ny != result.ny()) {
930 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::fftRows]");
934 if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
935 if (! result.convertRealToComplex ())
939 ImageFileArrayConst vReal = getArray();
940 ImageFileArrayConst vImag = getImaginaryArray();
942 fftw_plan plan = fftw_create_plan (m_ny, FFTW_FORWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
944 std::complex<double>* pcCol = new std::complex<double> [m_ny];
945 fftw_complex* in = new fftw_complex [m_ny];
946 for (int ix = 0; ix < m_nx; ix++) {
948 for (iy = 0; iy < m_ny; iy++) {
949 in[iy].re = vReal[ix][iy];
951 in[iy].im = vImag[ix][iy];
956 fftw_one (plan, in, NULL);
958 for (iy = 0; iy < m_ny; iy++)
959 pcCol[iy] = std::complex<double>(in[iy].re, in[iy].im);
961 Fourier::shuffleFourierToNaturalOrder (pcCol, m_ny);
962 for (iy = 0; iy < m_ny; iy++) {
963 vReal[ix][iy] = pcCol[iy].real() / m_ny;
964 vImag[ix][iy] = pcCol[iy].imag() / m_ny;
969 fftw_destroy_plan (plan);
976 ImageFile::ifftCols (ImageFile& result) const
978 if (m_nx != result.nx() || m_ny != result.ny()) {
979 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::fftRows]");
983 if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
984 if (! result.convertRealToComplex ())
988 fftw_complex* in = new fftw_complex [m_ny];
990 ImageFileArrayConst vReal = getArray();
991 ImageFileArrayConst vImag = getImaginaryArray();
993 fftw_plan plan = fftw_create_plan (m_ny, FFTW_BACKWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
994 std::complex<double>* pcCol = new std::complex<double> [m_ny];
997 unsigned int iArray = 0;
998 for (ix = 0; ix < m_nx; ix++) {
999 for (iy = 0; iy < m_ny; iy++) {
1002 dImag = vImag[ix][iy];
1003 pcCol[iy] = std::complex<double> (vReal[ix][iy], dImag);
1006 Fourier::shuffleNaturalToFourierOrder (pcCol, m_ny);
1008 for (iy = 0; iy < m_ny; iy++) {
1009 in[iy].re = pcCol[iy].real();
1010 in[iy].im = pcCol[iy].imag();
1013 fftw_one (plan, in, NULL);
1015 for (iy = 0; iy < m_ny; iy++) {
1016 vReal[ix][iy] = in[iy].re;
1017 vImag[ix][iy] = in[iy].im;
1022 fftw_destroy_plan (plan);
1033 ImageFile::fourier (ImageFile& result) const
1035 if (m_nx != result.nx() || m_ny != result.ny()) {
1036 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
1040 if (! result.isComplex())
1041 if (! result.convertRealToComplex ())
1044 ImageFileArrayConst vLHS = getArray();
1045 ImageFileArrayConst vLHSImag = getImaginaryArray();
1046 ImageFileArray vRealResult = result.getArray();
1047 ImageFileArray vImagResult = result.getImaginaryArray();
1049 unsigned int ix, iy;
1051 // alloc output matrix
1052 CTSimComplex** complexOut = new CTSimComplex* [m_nx];
1053 for (ix = 0; ix < m_nx; ix++)
1054 complexOut[ix] = new CTSimComplex [m_ny];
1056 // fourier each x column
1057 CTSimComplex* pY = new CTSimComplex [m_ny];
1058 for (ix = 0; ix < m_nx; ix++) {
1059 for (iy = 0; iy < m_ny; iy++) {
1062 dImag = vLHSImag[ix][iy];
1063 pY[iy] = std::complex<double>(vLHS[ix][iy], dImag);
1065 ProcessSignal::finiteFourierTransform (pY, complexOut[ix], m_ny, ProcessSignal::FORWARD);
1069 // fourier each y row
1070 CTSimComplex* pX = new CTSimComplex [m_nx];
1071 CTSimComplex* complexOutRow = new CTSimComplex [m_nx];
1072 for (iy = 0; iy < m_ny; iy++) {
1073 for (ix = 0; ix < m_nx; ix++)
1074 pX[ix] = complexOut[ix][iy];
1075 ProcessSignal::finiteFourierTransform (pX, complexOutRow, m_nx, ProcessSignal::FORWARD);
1076 for (ix = 0; ix < m_nx; ix++)
1077 complexOut[ix][iy] = complexOutRow[ix];
1080 delete [] complexOutRow;
1082 for (ix = 0; ix < m_nx; ix++)
1083 for (iy = 0; iy < m_ny; iy++) {
1084 vRealResult[ix][iy] = complexOut[ix][iy].real();
1085 vImagResult[ix][iy] = complexOut[ix][iy].imag();
1088 Fourier::shuffleFourierToNaturalOrder (result);
1090 // delete complexOut matrix
1091 for (ix = 0; ix < m_nx; ix++)
1092 delete [] complexOut[ix];
1093 delete [] complexOut;
1099 ImageFile::inverseFourier (ImageFile& result) const
1101 if (m_nx != result.nx() || m_ny != result.ny()) {
1102 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
1106 if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
1107 if (! result.convertRealToComplex ())
1111 ImageFileArrayConst vLHSReal = getArray();
1112 ImageFileArrayConst vLHSImag = getImaginaryArray();
1113 ImageFileArray vRealResult = result.getArray();
1114 ImageFileArray vImagResult = result.getImaginaryArray();
1116 unsigned int ix, iy;
1117 // alloc 2d complex output matrix
1118 CTSimComplex** complexOut = new CTSimComplex* [m_nx];
1119 for (ix = 0; ix < m_nx; ix++)
1120 complexOut[ix] = new CTSimComplex [m_ny];
1122 // put input image into result
1123 for (ix = 0; ix < m_nx; ix++)
1124 for (iy = 0; iy < m_ny; iy++) {
1125 vRealResult[ix][iy] = vLHSReal[ix][iy];
1127 vImagResult[ix][iy] = vLHSImag[ix][iy];
1129 vImagResult[ix][iy] = 0;
1132 Fourier::shuffleNaturalToFourierOrder (result);
1134 // ifourier each x column
1135 CTSimComplex* pCol = new CTSimComplex [m_ny];
1136 for (ix = 0; ix < m_nx; ix++) {
1137 for (iy = 0; iy < m_ny; iy++) {
1138 pCol[iy] = std::complex<double> (vRealResult[ix][iy], vImagResult[ix][iy]);
1140 ProcessSignal::finiteFourierTransform (pCol, complexOut[ix], m_ny, ProcessSignal::BACKWARD);
1144 // ifourier each y row
1145 CTSimComplex* complexInRow = new CTSimComplex [m_nx];
1146 CTSimComplex* complexOutRow = new CTSimComplex [m_nx];
1147 for (iy = 0; iy < m_ny; iy++) {
1148 for (ix = 0; ix < m_nx; ix++)
1149 complexInRow[ix] = complexOut[ix][iy];
1150 ProcessSignal::finiteFourierTransform (complexInRow, complexOutRow, m_nx, ProcessSignal::BACKWARD);
1151 for (ix = 0; ix < m_nx; ix++)
1152 complexOut[ix][iy] = complexOutRow[ix];
1154 delete [] complexInRow;
1155 delete [] complexOutRow;
1157 for (ix = 0; ix < m_nx; ix++)
1158 for (iy = 0; iy < m_ny; iy++) {
1159 vRealResult[ix][iy] = complexOut[ix][iy].real();
1160 vImagResult[ix][iy] = complexOut[ix][iy].imag();
1163 // delete complexOut matrix
1164 for (ix = 0; ix < m_nx; ix++)
1165 delete [] complexOut[ix];
1166 delete [] complexOut;
1173 ImageFile::magnitude (ImageFile& result) const
1175 if (m_nx != result.nx() || m_ny != result.ny()) {
1176 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
1180 ImageFileArray vReal = getArray();
1181 ImageFileArray vImag = getImaginaryArray();
1182 ImageFileArray vRealResult = result.getArray();
1184 for (unsigned int ix = 0; ix < m_nx; ix++)
1185 for (unsigned int iy = 0; iy < m_ny; iy++) {
1187 vRealResult[ix][iy] = ::sqrt (vReal[ix][iy] * vReal[ix][iy] + vImag[ix][iy] * vImag[ix][iy]);
1189 vRealResult[ix][iy] = ::fabs(vReal[ix][iy]);
1192 if (result.isComplex())
1193 result.reallocComplexToReal();
1199 ImageFile::phase (ImageFile& result) const
1201 if (m_nx != result.nx() || m_ny != result.ny()) {
1202 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
1206 ImageFileArray vReal = getArray();
1207 ImageFileArray vImag = getImaginaryArray();
1208 ImageFileArray vRealResult = result.getArray();
1210 for (unsigned int ix = 0; ix < m_nx; ix++) {
1211 for (unsigned int iy = 0; iy < m_ny; iy++) {
1213 vRealResult[ix][iy] = ::atan2 (vImag[ix][iy], vReal[ix][iy]);
1215 vRealResult[ix][iy] = 0;
1218 if (result.isComplex())
1219 result.reallocComplexToReal();
1225 ImageFile::real (ImageFile& result) const
1227 if (m_nx != result.nx() || m_ny != result.ny()) {
1228 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
1232 ImageFileArray vReal = getArray();
1233 ImageFileArray vRealResult = result.getArray();
1235 for (unsigned int ix = 0; ix < m_nx; ix++) {
1236 for (unsigned int iy = 0; iy < m_ny; iy++) {
1237 vRealResult[ix][iy] = vReal[ix][iy];
1241 if (result.isComplex())
1242 result.reallocComplexToReal();
1248 ImageFile::imaginary (ImageFile& result) const
1250 if (m_nx != result.nx() || m_ny != result.ny()) {
1251 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
1255 ImageFileArray vImag = getArray();
1256 ImageFileArray vRealResult = result.getArray();
1258 for (unsigned int ix = 0; ix < m_nx; ix++) {
1259 for (unsigned int iy = 0; iy < m_ny; iy++) {
1261 vRealResult[ix][iy] = vImag[ix][iy];
1263 vRealResult[ix][iy] = 0;
1267 if (result.isComplex())
1268 result.reallocComplexToReal();
1274 ImageFile::square (ImageFile& result) const
1276 if (m_nx != result.nx() || m_ny != result.ny()) {
1277 sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
1281 if (isComplex() && ! result.isComplex())
1282 result.convertRealToComplex();
1284 ImageFileArrayConst vLHS = getArray();
1285 ImageFileArrayConst vLHSImag = getImaginaryArray();
1286 ImageFileArray vResult = result.getArray();
1287 ImageFileArray vResultImag = result.getImaginaryArray();
1289 for (unsigned int ix = 0; ix < m_nx; ix++) {
1290 for (unsigned int iy = 0; iy < m_ny; iy++) {
1291 if (result.isComplex()) {
1294 dImag = vLHSImag[ix][iy];
1295 std::complex<double> cLHS (vLHS[ix][iy], dImag);
1296 std::complex<double> cResult = cLHS * cLHS;
1297 vResult[ix][iy] = cResult.real();
1298 vResultImag[ix][iy] = cResult.imag();
1300 vResult[ix][iy] = vLHS[ix][iy] * vLHS[ix][iy];
1308 ImageFile::convertExportFormatNameToID (const char* const formatName)
1310 int formatID = EXPORT_FORMAT_INVALID;
1312 for (int i = 0; i < s_iExportFormatCount; i++)
1313 if (strcasecmp (formatName, s_aszExportFormatName[i]) == 0) {
1322 ImageFile::convertExportFormatIDToName (int formatID)
1324 static const char *formatName = "";
1326 if (formatID >= 0 && formatID < s_iExportFormatCount)
1327 return (s_aszExportFormatName[formatID]);
1329 return (formatName);
1333 ImageFile::convertExportFormatIDToTitle (const int formatID)
1335 static const char *formatTitle = "";
1337 if (formatID >= 0 && formatID < s_iExportFormatCount)
1338 return (s_aszExportFormatTitle[formatID]);
1340 return (formatTitle);
1344 ImageFile::convertImportFormatNameToID (const char* const formatName)
1346 int formatID = IMPORT_FORMAT_INVALID;
1348 for (int i = 0; i < s_iImportFormatCount; i++)
1349 if (strcasecmp (formatName, s_aszImportFormatName[i]) == 0) {
1358 ImageFile::convertImportFormatIDToName (int formatID)
1360 static const char *formatName = "";
1362 if (formatID >= 0 && formatID < s_iImportFormatCount)
1363 return (s_aszImportFormatName[formatID]);
1365 return (formatName);
1369 ImageFile::convertImportFormatIDToTitle (const int formatID)
1371 static const char *formatTitle = "";
1373 if (formatID >= 0 && formatID < s_iImportFormatCount)
1374 return (s_aszImportFormatTitle[formatID]);
1376 return (formatTitle);
1380 ImageFile::importImage (const char* const pszFormat, const char* const pszFilename)
1382 int iFormatID = convertImportFormatNameToID (pszFormat);
1384 if (iFormatID == IMPORT_FORMAT_PPM)
1385 return readImagePPM (pszFilename);
1387 else if (iFormatID == IMPORT_FORMAT_PNG)
1388 return readImagePNG (pszFilename);
1391 sys_error (ERR_SEVERE, "Invalid format %s [ImageFile::importImage]", pszFormat);
1396 ImageFile::skipSpacePPM (FILE* fp)
1399 while (isspace (c) || c == '#') {
1400 if (c == '#') { // comment until end of line
1402 while (c != 13 && c != 10)
1413 ImageFile::readImagePPM (const char* const pszFile)
1415 FILE* fp = fopen (pszFile, "r");
1416 if ((fp = fopen (pszFile, "r")) == NULL)
1418 char cSignature = toupper(fgetc(fp));
1419 if (cSignature != 'P') {
1423 cSignature = fgetc(fp);
1424 if (cSignature == '5' || cSignature == '6') { // binary modes
1426 fp = fopen(pszFile, "rb"); // reopen in binary mode
1429 } else if (cSignature != '2' && cSignature != '3') {
1434 int nRows, nCols, iMaxValue;
1436 if (fscanf (fp, "%d", &nCols) != 1) {
1441 if (fscanf (fp, "%d", &nRows) != 1) {
1446 if (fscanf (fp, "%d", &iMaxValue) != 1) {
1450 setArraySize (nRows, nCols);
1452 if (cSignature == '5' || cSignature == '6') { // binary modes
1456 if (c != 10) // read msdos 13-10 newline
1460 skipSpacePPM (fp); // ascii may have comments
1462 bool bMonochromeImage = false;
1463 double dMaxValue = iMaxValue;
1464 double dMaxValue3 = iMaxValue * 3;
1466 ImageFileArray v = getArray();
1467 for (int iy = nRows - 1; iy >= 0; iy--) {
1468 for (int ix = 0; ix < nCols; ix++) {
1469 int iGS, iR, iG, iB;
1471 switch (cSignature) {
1473 if (fscanf(fp, "%d ", &iGS) != 1) {
1477 v[ix][iy] = iGS / dMaxValue;
1485 v[ix][iy] = iGS / dMaxValue;
1488 if (fscanf (fp, "%d %d %d ", &iR, &iG, &iB) != 3) {
1492 if (ix == 0 && iy == 0 && (iR == iG && iG == iB))
1493 bMonochromeImage = true;
1494 if (bMonochromeImage)
1495 v[ix][iy] = (iR + iG + iB) / dMaxValue3;
1497 dR = iR / dMaxValue;
1498 dG = iG / dMaxValue;
1499 dB = iB / dMaxValue;
1500 v[ix][iy] = colorToGrayscale (dR, dG, dB);
1512 if (ix == 0 && iy == 0 && (iR == iG && iG == iB))
1513 bMonochromeImage = true;
1515 if (bMonochromeImage)
1516 v[ix][iy] = (iR + iG + iB) / dMaxValue3;
1518 dR = iR / dMaxValue;
1519 dG = iG / dMaxValue;
1520 dB = iB / dMaxValue;
1521 v[ix][iy] = colorToGrayscale (dR, dG, dB);
1534 ImageFile::readImagePNG (const char* const pszFile)
1536 FILE* fp = fopen(pszFile, "rb");
1539 unsigned char header[8];
1540 fread (header, 1, 8, fp);
1541 if (png_sig_cmp (header, 0, 8)) {
1546 png_structp png_ptr = png_create_read_struct (PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
1552 png_infop info_ptr = png_create_info_struct(png_ptr);
1554 png_destroy_read_struct(&png_ptr, (png_infopp)NULL, (png_infopp)NULL);
1559 png_infop end_info = png_create_info_struct(png_ptr);
1561 png_destroy_read_struct(&png_ptr, &info_ptr, (png_infopp)NULL);
1566 if (setjmp(png_ptr->jmpbuf)) {
1567 png_destroy_read_struct(&png_ptr, &info_ptr, &end_info);
1572 png_init_io(png_ptr, fp);
1573 png_set_sig_bytes(png_ptr, 8);
1574 png_read_info(png_ptr, info_ptr);
1576 int width = png_get_image_width (png_ptr, info_ptr);
1577 int height = png_get_image_height (png_ptr, info_ptr);
1578 int bit_depth = png_get_bit_depth (png_ptr, info_ptr);
1579 int color_type = png_get_color_type (png_ptr, info_ptr);
1581 if (color_type == PNG_COLOR_TYPE_PALETTE && bit_depth <= 8)
1582 png_set_expand(png_ptr);
1584 if (color_type == PNG_COLOR_TYPE_GRAY && bit_depth < 8)
1585 png_set_expand(png_ptr);
1588 png_set_packing(png_ptr);
1590 if (color_type & PNG_COLOR_MASK_ALPHA)
1591 png_set_strip_alpha(png_ptr);
1593 if (bit_depth == 16)
1594 png_set_swap(png_ptr); // convert to little-endian format
1596 png_read_update_info(png_ptr, info_ptr); // update with transformations
1597 int rowbytes = png_get_rowbytes (png_ptr, info_ptr);
1598 bit_depth = png_get_bit_depth (png_ptr, info_ptr);
1599 color_type = png_get_color_type (png_ptr, info_ptr);
1601 png_bytep* row_pointers = new png_bytep [height];
1603 for (i = 0; i < height; i++)
1604 row_pointers[i] = new unsigned char [rowbytes];
1606 png_read_image(png_ptr, row_pointers);
1608 setArraySize (width, height);
1609 ImageFileArray v = getArray();
1610 for (int iy = 0; iy < height; iy++) {
1611 for (int ix = 0; ix < width; ix++) {
1613 if (color_type == PNG_COLOR_TYPE_GRAY) {
1615 dV = row_pointers[iy][ix] / 255.;
1616 else if (bit_depth == 16) {
1618 dV = (row_pointers[iy][iBase] + (row_pointers[iy][iBase+1] << 8)) / 65536.;
1620 } else if (color_type == PNG_COLOR_TYPE_RGB) {
1621 if (bit_depth == 8) {
1623 double dR = row_pointers[iy][iBase] / 255.;
1624 double dG = row_pointers[iy][iBase+1] / 255.;
1625 double dB = row_pointers[iy][iBase+2] / 255.;
1626 dV = colorToGrayscale (dR, dG, dR);
1629 v[ix][height-iy-1] = dV;
1633 png_read_end(png_ptr, end_info);
1634 png_destroy_read_struct(&png_ptr, &info_ptr, &end_info);
1636 for (i = 0; i < height; i++)
1637 delete row_pointers[i];
1638 delete row_pointers;
1646 ImageFile::exportImage (const char* const pszFormat, const char* const pszFilename, int nxcell, int nycell, double densmin, double densmax)
1648 int iFormatID = convertExportFormatNameToID (pszFormat);
1650 if (iFormatID == EXPORT_FORMAT_PGM)
1651 return writeImagePGM (pszFilename, nxcell, nycell, densmin, densmax);
1652 else if (iFormatID == EXPORT_FORMAT_PGMASCII)
1653 return writeImagePGMASCII (pszFilename, nxcell, nycell, densmin, densmax);
1654 else if (iFormatID == EXPORT_FORMAT_TEXT)
1655 return writeImageText (pszFilename);
1657 else if (iFormatID == EXPORT_FORMAT_PNG)
1658 return writeImagePNG (pszFilename, 8, nxcell, nycell, densmin, densmax);
1659 else if (iFormatID == EXPORT_FORMAT_PNG16)
1660 return writeImagePNG (pszFilename, 16, nxcell, nycell, densmin, densmax);
1662 #ifdef HAVE_CTN_DICOM
1663 else if (iFormatID == EXPORT_FORMAT_DICOM) {
1664 DicomExporter dicomExport (this);
1665 bool bSuccess = dicomExport.writeFile (pszFilename);
1667 sys_error (ERR_SEVERE, dicomExport.failMessage().c_str());
1671 else if (iFormatID == EXPORT_FORMAT_RAW)
1672 return writeImageRaw(pszFilename, nxcell, nycell);
1675 sys_error (ERR_SEVERE, "Invalid format %s [ImageFile::exportImage]", pszFormat);
1681 ImageFile::writeImagePGM (const char* const outfile, int nxcell, int nycell, double densmin, double densmax)
1686 ImageFileArray v = getArray();
1688 unsigned char* rowp = new unsigned char [nx * nxcell];
1690 if ((fp = fopen (outfile, "wb")) == NULL)
1693 fprintf(fp, "P5\n");
1694 fprintf(fp, "%d %d\n", nx, ny);
1695 fprintf(fp, "255\n");
1697 for (int irow = ny - 1; irow >= 0; irow--) {
1698 for (int icol = 0; icol < nx; icol++) {
1699 int pos = icol * nxcell;
1700 double dens = (v[icol][irow] - densmin) / (densmax - densmin);
1701 dens = clamp (dens, 0., 1.);
1702 for (int p = pos; p < pos + nxcell; p++) {
1703 rowp[p] = static_cast<unsigned int> (dens * 255.);
1706 for (int ir = 0; ir < nycell; ir++) {
1707 for (int ic = 0; ic < nx * nxcell; ic++)
1708 fputc( rowp[ic], fp );
1719 ImageFile::writeImagePGMASCII (const char* const outfile, int nxcell, int nycell, double densmin, double densmax)
1724 ImageFileArray v = getArray();
1726 unsigned char* rowp = new unsigned char [nx * nxcell];
1728 if ((fp = fopen (outfile, "wb")) == NULL)
1731 fprintf(fp, "P2\n");
1732 fprintf(fp, "%d %d\n", nx, ny);
1733 fprintf(fp, "255\n");
1735 for (int irow = ny - 1; irow >= 0; irow--) {
1736 for (int icol = 0; icol < nx; icol++) {
1737 int pos = icol * nxcell;
1738 double dens = (v[icol][irow] - densmin) / (densmax - densmin);
1739 dens = clamp (dens, 0., 1.);
1740 for (int p = pos; p < pos + nxcell; p++) {
1741 rowp[p] = static_cast<unsigned int> (dens * 255.);
1744 for (int ir = 0; ir < nycell; ir++) {
1745 for (int ic = 0; ic < nx * nxcell; ic++)
1746 fprintf(fp, "%d ", rowp[ic]);
1758 ImageFile::writeImageText (const char* const outfile)
1763 ImageFileArray v = getArray();
1764 ImageFileArray vImag = getImaginaryArray();
1766 if ((fp = fopen (outfile, "w")) == NULL)
1769 for (int irow = ny - 1; irow >= 0; irow--) {
1770 for (int icol = 0; icol < nx; icol++) {
1772 if (vImag[icol][irow] >= 0)
1773 fprintf (fp, "%.9g+%.9gi ", v[icol][irow], vImag[icol][irow]);
1775 fprintf (fp, "%.9g-%.9gi ", v[icol][irow], -vImag[icol][irow]);
1777 fprintf (fp, "%12.8g ", v[icol][irow]);
1790 ImageFile::writeImagePNG (const char* const outfile, int bitdepth, int nxcell, int nycell, double densmin, double densmax)
1792 double max_out_level = (1 << bitdepth) - 1;
1795 ImageFileArray v = getArray();
1797 unsigned char* rowp = new unsigned char [nx * nxcell * (bitdepth / 8)];
1799 FILE *fp = fopen (outfile, "wb");
1803 png_structp png_ptr = png_create_write_struct (PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
1807 png_infop info_ptr = png_create_info_struct (png_ptr);
1809 png_destroy_write_struct (&png_ptr, (png_infopp) NULL);
1814 if (setjmp (png_ptr->jmpbuf)) {
1815 png_destroy_write_struct (&png_ptr, &info_ptr);
1820 png_init_io(png_ptr, fp);
1822 png_set_IHDR (png_ptr, info_ptr, nx * nxcell, ny * nycell, bitdepth, PNG_COLOR_TYPE_GRAY, PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT);
1824 png_write_info(png_ptr, info_ptr);
1825 for (int irow = ny - 1; irow >= 0; irow--) {
1826 png_bytep row_pointer = rowp;
1828 for (int icol = 0; icol < nx; icol++) {
1829 int pos = icol * nxcell;
1830 double dens = (v[icol][irow] - densmin) / (densmax - densmin);
1831 dens = clamp (dens, 0., 1.);
1832 unsigned int outval = static_cast<unsigned int> (dens * max_out_level);
1834 for (int p = pos; p < pos + nxcell; p++) {
1839 rowp[rowpos+1] = (outval >> 8) & 0xFF;
1840 rowp[rowpos] = (outval & 0xFF);
1844 for (int ir = 0; ir < nycell; ir++)
1845 png_write_rows (png_ptr, &row_pointer, 1);
1848 png_write_end (png_ptr, info_ptr);
1849 png_destroy_write_struct (&png_ptr, &info_ptr);
1860 static const int N_GRAYSCALE=256;
1863 ImageFile::writeImageGIF (const char* const outfile, int nxcell, int nycell, double densmin, double densmax)
1865 int gs_indices[N_GRAYSCALE];
1868 ImageFileArray v = getArray();
1870 unsigned char* rowp = new unsigned char [nx * nxcell];
1872 gdImagePtr gif = gdImageCreate(nx * nxcell, ny * nycell);
1873 for (int i = 0; i < N_GRAYSCALE; i++)
1874 gs_indices[i] = gdImageColorAllocate(gif, i, i, i);
1876 int lastrow = ny * nycell - 1;
1877 for (int irow = 0; irow < ny; irow++) {
1878 int rpos = irow * nycell;
1879 for (int ir = rpos; ir < rpos + nycell; ir++) {
1880 for (int icol = 0; icol < nx; icol++) {
1881 int cpos = icol * nxcell;
1882 double dens = (v[icol][irow] - densmin) / (densmax - densmin);
1883 dens = clamp(dens, 0., 1.);
1884 for (int ic = cpos; ic < cpos + nxcell; ic++) {
1885 rowp[ic] = (unsigned int) (dens * (double) (N_GRAYSCALE - 1));
1886 gdImageSetPixel(gif, ic, lastrow - ir, gs_indices[rowp[ic]]);
1893 if ((out = fopen (outfile,"w")) == NULL) {
1894 sys_error(ERR_SEVERE, "Error opening output file %s for writing", outfile);
1897 gdImageGif(gif,out);
1899 gdImageDestroy(gif);
1908 ImageFile::writeImageRaw (const char* const outfile, int nxcell, int nycell)
1913 ImageFileArray v = getArray();
1915 if ((fp = fopen (outfile, "wb")) == NULL)
1918 for (int irow = ny - 1; irow >= 0; irow--) {
1919 for (int icol = 0; icol < nx; icol++) {
1920 float dens = v[icol][irow];
1921 fwrite(&dens, sizeof(float), 1, fp);