libctsim/imagefile.cpp

   1 /*****************************************************************************
   2 ** FILE IDENTIFICATION
   3 **
   4 **      Name:         imagefile.cpp
   5 **      Purpose:      Imagefile classes
   6 **      Programmer:   Kevin Rosenberg
   7 **      Date Started: June 2000
   8 **
   9 **  This is part of the CTSim program
  10 **  Copyright (C) 1983-2000 Kevin Rosenberg
  11 **
  12 **  $Id: imagefile.cpp,v 1.30 2001/01/02 07:18:07 kevin Exp $
  13 **
  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.
  17 **
  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.
  22 **
  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 ******************************************************************************/
  27
  28 #include "ct.h"
  29
  30 const int ImageFile::FORMAT_INVALID = -1;\r
  31 const int ImageFile::FORMAT_PGM = 0;\r
  32 const int ImageFile::FORMAT_PGMASCII = 1;\r
  33 #ifdef HAVE_PNG\r
  34 const int ImageFile::FORMAT_PNG = 2;\r
  35 const int ImageFile::FORMAT_PNG16 = 3;\r
  36 #endif\r
  37 \r
  38 const char* ImageFile::s_aszFormatName[] = \r
  39 {\r
  40   {"pgm"},\r
  41   {"pgmascii"},\r
  42 #ifdef HAVE_PNG\r
  43   {"png"},\r
  44   {"png16"},\r
  45 #endif\r
  46 };\r
  47 \r
  48 const char* ImageFile::s_aszFormatTitle[] = \r
  49 {\r
  50   {"PGM"},\r
  51   {"PGM ASCII"},\r
  52   {"PNG"},\r
  53   {"PNG 16-bit"},\r
  54 };\r
  55 \r
  56 const int ImageFile::s_iFormatCount = sizeof(s_aszFormatName) / sizeof(const char*);\r
  57 \r
  58 \r
  59
  60 F32Image::F32Image (int nx, int ny, int dataType)\r
  61 : Array2dFile (nx, ny, sizeof(kfloat32), Array2dFile::PIXEL_FLOAT32, dataType)\r
  62 {\r
  63 }\r
  64 \r
  65 F32Image::F32Image (void)\r
  66 : Array2dFile()\r
  67 {\r
  68   setPixelFormat (Array2dFile::PIXEL_FLOAT32);\r
  69   setPixelSize (sizeof(kfloat32));\r
  70   setDataType (Array2dFile::DATA_TYPE_REAL);\r
  71 }\r
  72 \r
  73 F64Image::F64Image (int nx, int ny, int dataType)\r
  74 : Array2dFile (nx, ny, sizeof(kfloat64), Array2dFile::PIXEL_FLOAT64, dataType)\r
  75 {\r
  76 }\r
  77 \r
  78 F64Image::F64Image (void)\r
  79 : Array2dFile ()\r
  80 {\r
  81   setPixelFormat (PIXEL_FLOAT64);\r
  82   setPixelSize (sizeof(kfloat64));\r
  83   setDataType (Array2dFile::DATA_TYPE_REAL);\r
  84 }\r
  85
  86 void
  87 ImageFile::filterResponse (const char* const domainName, double bw, const char* const filterName, double filt_param, double dInputScale, double dOutputScale)
  88 {\r
  89   ImageFileArray v = getArray();\r
  90   SignalFilter filter (filterName, domainName, bw, filt_param);\r
  91   \r
  92 #if 1\r
  93   int iXCenter, iYCenter;\r
  94   if (isEven (m_nx))\r
  95     iXCenter = m_nx / 2;\r
  96   else\r
  97     iXCenter = (m_nx - 1) / 2;\r
  98   if (isEven (m_ny))\r
  99     iYCenter = m_ny / 2;\r
 100   else\r
 101     iYCenter = (m_ny - 1) / 2;\r
 102   \r
 103   for (unsigned int ix = 0; ix < m_nx; ix++)\r
 104     for (unsigned int iy = 0; iy < m_ny; iy++) {\r
 105       long lD2 = ((ix - iXCenter) * (ix - iXCenter)) + ((iy - iYCenter) * (iy - iYCenter));\r
 106       double r = ::sqrt (static_cast<double>(lD2)) * dInputScale;\r
 107       v[ix][iy] = filter.response (r) * dOutputScale;\r
 108     }\r
 109 #else
 110     int hx = (m_nx - 1) / 2;
 111     int hy = (m_ny - 1) / 2;
 112
 113     for (int i = -hx; i <= hx; i++) {
 114       for (int j = -hy; j <= hy; j++) {
 115         double r = ::sqrt (i * i + j * j);
 116
 117         v[i+hx][j+hy] = filter.response (r);
 118       }
 119     }\r
 120 #endif
 121 }
 122
 123 int
 124 ImageFile::display (void) const
 125 {
 126   double pmin, pmax;
 127
 128   getMinMax (pmin, pmax);
 129
 130   return (displayScaling (1, pmin, pmax));
 131 }
 132
 133 int
 134 ImageFile::displayScaling (const int scale, const ImageFileValue pmin, const ImageFileValue pmax) const
 135 {
 136   int nx = m_nx;
 137   int ny = m_ny;
 138   ImageFileArrayConst v = getArray();
 139   if (v == NULL || nx == 0 || ny == 0)
 140     return 0;
 141
 142 #if HAVE_G2_H
 143   int* pPens = new int [nx * ny * scale * scale ];
 144
 145   double view_scale = 255 / (pmax - pmin);
 146   int id_X11 = g2_open_X11 (nx * scale, ny * scale);
 147   int grayscale[256];
 148   for (int i = 0; i < 256; i++) {
 149     double cval = i / 255.;
 150     grayscale[i] = g2_ink (id_X11, cval, cval, cval);
 151   }
 152
 153   for (int iy = ny - 1; iy >= 0; iy--) {
 154     int iRowPos = ((ny - 1 - iy) * scale) * (nx * scale);
 155     for (int ix = 0; ix < nx; ix++) {
 156       int cval = static_cast<int>((v[ix][iy] - pmin) * view_scale);
 157       if (cval < 0)
 158         cval = 0;
 159       else if (cval > 255)
 160         cval = 255;
 161       for (int sy = 0; sy < scale; sy++)
 162         for (int sx = 0; sx < scale; sx++)
 163           pPens[iRowPos+(sy * nx * scale)+(sx + (ix * scale))] = grayscale[cval];
 164     }
 165   }
 166
 167   g2_image (id_X11, 0., 0., nx * scale, ny * scale, pPens);
 168
 169   delete pPens;
 170   return (id_X11);
 171 #else
 172   return 0;
 173 #endif
 174 }
 175
 176
 177
 178 // ImageFile::comparativeStatistics    Calculate comparative stats
 179 //
 180 // OUTPUT
 181 //   d   Normalized root mean squared distance measure
 182 //   r   Normalized mean absolute distance measure
 183 //   e   Worst case distance measure
 184 //
 185 // REFERENCES
 186 //  G.T. Herman, Image Reconstruction From Projections, 1980
 187
 188 bool
 189 ImageFile::comparativeStatistics (const ImageFile& imComp, double& d, double& r, double& e) const
 190 {
 191   if (imComp.nx() != m_nx && imComp.ny() != m_ny) {
 192     sys_error (ERR_WARNING, "Image sizes differ [ImageFile::comparativeStatistics]");
 193     return false;
 194   }
 195   ImageFileArrayConst v = getArray();
 196   if (v == NULL || m_nx == 0 || m_ny == 0)
 197     return false;
 198
 199   ImageFileArrayConst vComp = imComp.getArray();
 200
 201   double myMean = 0.;
 202   for (unsigned int ix = 0; ix < m_nx; ix++) {
 203     for (unsigned int iy = 0; iy < m_ny; iy++) {
 204       myMean += v[ix][iy];
 205     }
 206   }
 207   myMean /= (m_nx * m_ny);
 208
 209   double sqErrorSum = 0.;
 210   double absErrorSum = 0.;
 211   double sqDiffFromMeanSum = 0.;
 212   double absValueSum = 0.;
 213   for (unsigned int ix2 = 0; ix2 < m_nx; ix2++) {
 214     for (unsigned int iy = 0; iy < m_ny; iy++) {
 215       double diff = v[ix2][iy] - vComp[ix2][iy];
 216       sqErrorSum += diff * diff;
 217       absErrorSum += fabs(diff);
 218       double diffFromMean = v[ix2][iy] - myMean;
 219       sqDiffFromMeanSum += diffFromMean * diffFromMean;
 220       absValueSum += fabs(v[ix2][iy]);
 221     }
 222   }
 223
 224   d = ::sqrt (sqErrorSum / sqDiffFromMeanSum);
 225   r = absErrorSum / absValueSum;
 226
 227   int hx = m_nx / 2;
 228   int hy = m_ny / 2;
 229   double eMax = -1;
 230   for (int ix3 = 0; ix3 < hx; ix3++) {
 231     for (int iy = 0; iy < hy; iy++) {
 232       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]);
 233       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]);
 234       double error = fabs (avgPixel - avgPixelComp);
 235       if (error > eMax)
 236         eMax = error;
 237     }
 238   }
 239
 240   e = eMax;
 241
 242   return true;
 243 }
 244
 245
 246 bool
 247 ImageFile::printComparativeStatistics (const ImageFile& imComp, std::ostream& os) const
 248 {
 249   double d, r, e;
 250
 251   if (comparativeStatistics (imComp, d, r, e)) {
 252     os << "  Normalized root mean squared distance (d): " << d << std::endl;
 253     os << "      Normalized mean absolute distance (r): " << r << std::endl;
 254     os << "Worst case distance (2x2 pixel average) (e): " << e << std::endl;
 255     return true;
 256   }
 257   return false;
 258 }
 259
 260
 261 void
 262 ImageFile::printStatistics (std::ostream& os) const
 263 {
 264   double min, max, mean, mode, median, stddev;
 265
 266   statistics (min, max, mean, mode, median, stddev);\r
 267   if (isComplex())\r
 268     os << "Real Component Statistics" << std::endl;\r
 269
 270   os << "   min: " << min << std::endl;
 271   os << "   max: " << max << std::endl;
 272   os << "  mean: " << mean << std::endl;
 273   os << "  mode: " << mode << std::endl;
 274   os << "median: " << median << std::endl;
 275   os << "stddev: " << stddev << std::endl;\r
 276   \r
 277   if (isComplex()) {\r
 278     statistics (getImaginaryArray(), min, max, mean, mode, median, stddev);\r
 279     os << std::endl << "Imaginary Component Statistics" << std::endl;
 280     os << "   min: " << min << std::endl;\r
 281     os << "   max: " << max << std::endl;\r
 282     os << "  mean: " << mean << std::endl;\r
 283     os << "  mode: " << mode << std::endl;\r
 284     os << "median: " << median << std::endl;\r
 285     os << "stddev: " << stddev << std::endl;\r
 286   }\r
 287 }
 288
 289
 290 void
 291 ImageFile::statistics (double& min, double& max, double& mean, double& mode, double& median, double& stddev) const
 292 {\r
 293   ImageFileArrayConst v = getArray();\r
 294   statistics (v, min, max, mean, mode, median, stddev);\r
 295 }
 296
 297
 298 void\r
 299 ImageFile::statistics (ImageFileArrayConst v, double& min, double& max, double& mean, double& mode, double& median, double& stddev) const\r
 300 {\r
 301   int nx = m_nx;\r
 302   int ny = m_ny;\r
 303   \r
 304   if (v == NULL || nx == 0 || ny == 0)\r
 305     return;\r
 306   \r
 307   std::vector<double> vecImage;\r
 308   int iVec = 0;\r
 309   vecImage.resize (nx * ny);\r
 310   for (int ix = 0; ix < nx; ix++) {\r
 311     for (int iy = 0; iy < ny; iy++)\r
 312       vecImage[iVec++] = v[ix][iy];\r
 313   }\r
 314   \r
 315   vectorNumericStatistics (vecImage, nx * ny, min, max, mean, mode, median, stddev);\r
 316 }\r
 317 \r
 318 void
 319 ImageFile::getMinMax (double& min, double& max) const
 320 {
 321   int nx = m_nx;
 322   int ny = m_ny;
 323   ImageFileArrayConst v = getArray();
 324
 325   if (v == NULL || nx == 0 || ny == 0)
 326     return;
 327
 328   min = v[0][0];
 329   max = v[0][0];
 330   for (int ix = 0; ix < nx; ix++) {
 331     for (int iy = 0; iy < ny; iy++) {
 332       if (v[ix][iy] > max)
 333         max = v[ix][iy];
 334       if (v[ix][iy] < min)
 335         min = v[ix][iy];
 336     }
 337   }
 338 }
 339 \r
 340 bool\r
 341 ImageFile::convertRealToComplex ()\r
 342 {\r
 343   if (dataType() != Array2dFile::DATA_TYPE_REAL)\r
 344     return false;\r
 345   \r
 346   if (! reallocRealToComplex())\r
 347     return false;\r
 348   \r
 349   ImageFileArray vImag = getImaginaryArray();\r
 350   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
 351     ImageFileColumn vCol = vImag[ix];\r
 352     for (unsigned int iy = 0; iy < m_ny; iy++)\r
 353       *vCol++ = 0;\r
 354   }\r
 355   \r
 356   return true;\r
 357 }\r
 358 \r
 359 bool\r
 360 ImageFile::convertComplexToReal ()\r
 361 {\r
 362   if (dataType() != Array2dFile::DATA_TYPE_COMPLEX)\r
 363     return false;\r
 364   \r
 365   ImageFileArray vReal = getArray();\r
 366   ImageFileArray vImag = getImaginaryArray();\r
 367   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
 368     ImageFileColumn vRealCol = vReal[ix];\r
 369     ImageFileColumn vImagCol = vImag[ix];\r
 370     for (unsigned int iy = 0; iy < m_ny; iy++) {\r
 371       CTSimComplex c (*vRealCol, *vImagCol);\r
 372       *vRealCol++ = std::abs (c);\r
 373       vImagCol++;\r
 374     }\r
 375   }\r
 376   \r
 377   return reallocComplexToReal();\r
 378 }\r
 379 \r
 380 bool\r
 381 ImageFile::subtractImages (const ImageFile& rRHS, ImageFile& result) const\r
 382 {\r
 383   if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) {\r
 384     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");\r
 385     return false;\r
 386   }\r
 387   \r
 388   if (isComplex() || rRHS.isComplex() && ! result.isComplex())\r
 389     result.convertRealToComplex();\r
 390   \r
 391   ImageFileArrayConst vLHS = getArray();\r
 392   ImageFileArrayConst vLHSImag = getImaginaryArray();\r
 393   ImageFileArrayConst vRHS = rRHS.getArray();\r
 394   ImageFileArrayConst vRHSImag = rRHS.getImaginaryArray();\r
 395   ImageFileArray vResult = result.getArray();\r
 396   ImageFileArray vResultImag = result.getImaginaryArray();\r
 397   \r
 398   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
 399     for (unsigned int iy = 0; iy < m_ny; iy++) {\r
 400       vResult[ix][iy] = vLHS[ix][iy] - vRHS[ix][iy];\r
 401       if (result.isComplex()) {\r
 402         vResultImag[ix][iy] = 0;\r
 403         if (isComplex())\r
 404           vResultImag[ix][iy] += vLHSImag[ix][iy];\r
 405         if (rRHS.isComplex())\r
 406           vResultImag[ix][iy] -= vRHSImag[ix][iy];\r
 407       }\r
 408     }\r
 409   }\r
 410   \r
 411   return true;\r
 412 }\r
 413
 414 bool\r
 415 ImageFile::addImages (const ImageFile& rRHS, ImageFile& result) const\r
 416 {\r
 417   if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) {\r
 418     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");\r
 419     return false;\r
 420   }\r
 421   \r
 422   if (isComplex() || rRHS.isComplex() && ! result.isComplex())\r
 423     result.convertRealToComplex();\r
 424   \r
 425   ImageFileArrayConst vLHS = getArray();\r
 426   ImageFileArrayConst vLHSImag = getImaginaryArray();\r
 427   ImageFileArrayConst vRHS = rRHS.getArray();\r
 428   ImageFileArrayConst vRHSImag = rRHS.getImaginaryArray();\r
 429   ImageFileArray vResult = result.getArray();\r
 430   ImageFileArray vResultImag = result.getImaginaryArray();\r
 431   \r
 432   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
 433     for (unsigned int iy = 0; iy < m_ny; iy++) {\r
 434       vResult[ix][iy] = vLHS[ix][iy] + vRHS[ix][iy];\r
 435       if (result.isComplex()) {\r
 436         vResultImag[ix][iy] = 0;\r
 437         if (isComplex())\r
 438           vResultImag[ix][iy] += vLHSImag[ix][iy];\r
 439         if (rRHS.isComplex())\r
 440           vResultImag[ix][iy] += vRHSImag[ix][iy];\r
 441       }\r
 442     }\r
 443   }\r
 444   \r
 445   return true;\r
 446 }\r
 447 \r
 448 bool\r
 449 ImageFile::multiplyImages (const ImageFile& rRHS, ImageFile& result) const\r
 450 {\r
 451   if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) {\r
 452     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");\r
 453     return false;\r
 454   }\r
 455   \r
 456   if (isComplex() || rRHS.isComplex() && ! result.isComplex())\r
 457     result.convertRealToComplex();\r
 458   \r
 459   ImageFileArrayConst vLHS = getArray();\r
 460   ImageFileArrayConst vLHSImag = getImaginaryArray();\r
 461   ImageFileArrayConst vRHS = rRHS.getArray();\r
 462   ImageFileArrayConst vRHSImag = rRHS.getImaginaryArray();\r
 463   ImageFileArray vResult = result.getArray();\r
 464   ImageFileArray vResultImag = result.getImaginaryArray();\r
 465   \r
 466   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
 467     for (unsigned int iy = 0; iy < m_ny; iy++) {\r
 468       if (result.isComplex()) {\r
 469         double dImag = 0;\r
 470         if (isComplex())\r
 471           dImag = vLHSImag[ix][iy];\r
 472         std::complex<double> cLHS (vLHS[ix][iy], dImag);\r
 473         dImag = 0;\r
 474         if (rRHS.isComplex())\r
 475           dImag = vRHSImag[ix][iy];\r
 476         std::complex<double> cRHS (vRHS[ix][iy], dImag);\r
 477         std::complex<double> cResult = cLHS * cRHS;\r
 478         vResult[ix][iy] = cResult.real();\r
 479         vResultImag[ix][iy] = cResult.imag();\r
 480       } else\r
 481         vResult[ix][iy] = vLHS[ix][iy] * vRHS[ix][iy];\r
 482     }\r
 483   }\r
 484   \r
 485   \r
 486   return true;\r
 487 }\r
 488 \r
 489 bool\r
 490 ImageFile::divideImages (const ImageFile& rRHS, ImageFile& result) const\r
 491 {\r
 492   if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) {\r
 493     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");\r
 494     return false;\r
 495   }\r
 496   \r
 497   if (isComplex() || rRHS.isComplex() && ! result.isComplex())\r
 498     result.convertRealToComplex();\r
 499   \r
 500   ImageFileArrayConst vLHS = getArray();\r
 501   ImageFileArrayConst vLHSImag = getImaginaryArray();\r
 502   ImageFileArrayConst vRHS = rRHS.getArray();\r
 503   ImageFileArrayConst vRHSImag = rRHS.getImaginaryArray();\r
 504   ImageFileArray vResult = result.getArray();\r
 505   ImageFileArray vResultImag = result.getImaginaryArray();\r
 506   \r
 507   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
 508     for (unsigned int iy = 0; iy < m_ny; iy++) {\r
 509       if (result.isComplex()) {\r
 510         double dImag = 0;\r
 511         if (isComplex())\r
 512           dImag = vLHSImag[ix][iy];\r
 513         std::complex<double> cLHS (vLHS[ix][iy], dImag);\r
 514         dImag = 0;\r
 515         if (rRHS.isComplex())\r
 516           dImag = vRHSImag[ix][iy];\r
 517         std::complex<double> cRHS (vRHS[ix][iy], dImag);\r
 518         std::complex<double> cResult = cLHS / cRHS;\r
 519         vResult[ix][iy] = cResult.real();\r
 520         vResultImag[ix][iy] = cResult.imag();\r
 521       } else {\r
 522         if (vRHS != 0)\r
 523           vResult[ix][iy] = vLHS[ix][iy] / vRHS[ix][iy];\r
 524         else\r
 525           vResult[ix][iy] = 0;\r
 526       }\r
 527     }\r
 528   }\r
 529   \r
 530   return true;\r
 531 }\r
 532 \r
 533 \r
 534 bool\r
 535 ImageFile::invertPixelValues (ImageFile& result) const\r
 536 {\r
 537   if (m_nx != result.nx() || m_ny != result.ny()) {\r
 538     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
 539     return false;\r
 540   }\r
 541   \r
 542   if (isComplex() && ! result.isComplex())\r
 543     result.convertRealToComplex();\r
 544   \r
 545   ImageFileArrayConst vLHS = getArray();\r
 546   ImageFileArray vResult = result.getArray();\r
 547   \r
 548   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
 549     ImageFileColumnConst in = vLHS[ix];\r
 550     ImageFileColumn out = vResult[ix];\r
 551     for (unsigned int iy = 0; iy < m_ny; iy++)\r
 552       *out++ = - *in++;\r
 553   }\r
 554   \r
 555   return true;\r
 556 }\r
 557 \r
 558 bool\r
 559 ImageFile::sqrt (ImageFile& result) const\r
 560 {\r
 561   if (m_nx != result.nx() || m_ny != result.ny()) {\r
 562     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
 563     return false;\r
 564   }\r
 565   \r
 566   if (isComplex() && ! result.isComplex())\r
 567     result.convertRealToComplex();\r
 568   \r
 569   bool bComplexOutput = result.isComplex();\r
 570   ImageFileArrayConst vLHS = getArray();\r
 571   if (! bComplexOutput)   // check if should convert to complex output\r
 572     for (unsigned int ix = 0; ix < m_nx; ix++)\r
 573       for (unsigned int iy = 0; iy < m_ny; iy++)\r
 574         if (! bComplexOutput && vLHS[ix][iy] < 0) {\r
 575           result.convertRealToComplex();\r
 576           bComplexOutput = true;\r
 577           break;\r
 578         }\r
 579         \r
 580         ImageFileArrayConst vLHSImag = getImaginaryArray();\r
 581         ImageFileArray vResult = result.getArray();\r
 582         ImageFileArray vResultImag = result.getImaginaryArray();\r
 583         \r
 584         for (unsigned int ix = 0; ix < m_nx; ix++) {\r
 585           for (unsigned int iy = 0; iy < m_ny; iy++) {\r
 586             if (result.isComplex()) {\r
 587               double dImag = 0;\r
 588               if (isComplex())\r
 589                 dImag = vLHSImag[ix][iy];\r
 590               std::complex<double> cLHS (vLHS[ix][iy], dImag);\r
 591               std::complex<double> cResult = std::sqrt(cLHS);\r
 592               vResult[ix][iy] = cResult.real();\r
 593               vResultImag[ix][iy] = cResult.imag();\r
 594             } else\r
 595               vResult[ix][iy] = ::sqrt (vLHS[ix][iy]);\r
 596           }\r
 597         }\r
 598         \r
 599         \r
 600         return true;\r
 601 }\r
 602 \r
 603 bool\r
 604 ImageFile::log (ImageFile& result) const\r
 605 {\r
 606   if (m_nx != result.nx() || m_ny != result.ny()) {\r
 607     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
 608     return false;\r
 609   }\r
 610   \r
 611   if (isComplex() && ! result.isComplex())\r
 612     result.convertRealToComplex();\r
 613   \r
 614   ImageFileArrayConst vLHS = getArray();\r
 615   ImageFileArrayConst vLHSImag = getImaginaryArray();\r
 616   ImageFileArray vResult = result.getArray();\r
 617   ImageFileArray vResultImag = result.getImaginaryArray();\r
 618   \r
 619   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
 620     for (unsigned int iy = 0; iy < m_ny; iy++) {\r
 621       if (result.isComplex()) {\r
 622         double dImag = 0;\r
 623         if (isComplex())\r
 624           dImag = vLHSImag[ix][iy];\r
 625         std::complex<double> cLHS (vLHS[ix][iy], dImag);\r
 626         std::complex<double> cResult = std::log (cLHS);\r
 627         vResult[ix][iy] = cResult.real();\r
 628         vResultImag[ix][iy] = cResult.imag();\r
 629       } else\r
 630         vResult[ix][iy] = ::log (vLHS[ix][iy]);\r
 631     }\r
 632   }\r
 633   \r
 634   \r
 635   return true;\r
 636 }\r
 637 \r
 638 bool\r
 639 ImageFile::exp (ImageFile& result) const\r
 640 {\r
 641   if (m_nx != result.nx() || m_ny != result.ny()) {\r
 642     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
 643     return false;\r
 644   }\r
 645   \r
 646   if (isComplex() && ! result.isComplex())\r
 647     result.convertRealToComplex();\r
 648   \r
 649   ImageFileArrayConst vLHS = getArray();\r
 650   ImageFileArrayConst vLHSImag = getImaginaryArray();\r
 651   ImageFileArray vResult = result.getArray();\r
 652   ImageFileArray vResultImag = result.getImaginaryArray();\r
 653   \r
 654   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
 655     for (unsigned int iy = 0; iy < m_ny; iy++) {\r
 656       if (result.isComplex()) {\r
 657         double dImag = 0;\r
 658         if (isComplex())\r
 659           dImag = vLHSImag[ix][iy];\r
 660         std::complex<double> cLHS (vLHS[ix][iy], dImag);\r
 661         std::complex<double> cResult = std::exp (cLHS);\r
 662         vResult[ix][iy] = cResult.real();\r
 663         vResultImag[ix][iy] = cResult.imag();\r
 664       } else\r
 665         vResult[ix][iy] = ::exp (vLHS[ix][iy]);\r
 666     }\r
 667   }\r
 668   \r
 669   \r
 670   return true;\r
 671 }\r
 672 \r
 673 bool\r
 674 ImageFile::scaleImage (ImageFile& result) const\r
 675 {\r
 676   unsigned int nx = m_nx;\r
 677   unsigned int ny = m_ny;\r
 678   unsigned int newNX = result.nx();\r
 679   unsigned int newNY = result.ny();\r
 680   \r
 681   double dXScale = static_cast<double>(newNX) / static_cast<double>(nx);\r
 682   double dYScale = static_cast<double>(newNY) / static_cast<double>(ny);\r
 683   \r
 684   if (isComplex() && ! result.isComplex())\r
 685     result.convertRealToComplex();\r
 686   \r
 687   ImageFileArrayConst vReal = getArray();\r
 688   ImageFileArrayConst vImag = getImaginaryArray();\r
 689   ImageFileArray vResult = result.getArray();\r
 690   ImageFileArray vResultImag = result.getImaginaryArray();\r
 691   \r
 692   for (unsigned int ix = 0; ix < newNX; ix++) {\r
 693     for (unsigned int iy = 0; iy < newNY; iy++) {\r
 694       double dXPos = ix / dXScale;\r
 695       double dYPos = iy / dYScale;\r
 696       unsigned int scaleNX = static_cast<unsigned int> (dXPos);\r
 697       unsigned int scaleNY = static_cast<unsigned int> (dYPos);\r
 698       double dXFrac = dXPos - scaleNX;\r
 699       double dYFrac = dYPos - scaleNY;\r
 700       if (scaleNX >= nx - 1 || scaleNY >= ny - 1) {\r
 701         vResult[ix][iy] = vReal[scaleNX][scaleNY];\r
 702         if (result.isComplex()) {\r
 703           if (isComplex())\r
 704             vResultImag[ix][iy] = vImag[scaleNX][scaleNY];\r
 705           else\r
 706             vResultImag[ix][iy] = 0;\r
 707         }\r
 708       } else {\r
 709         vResult[ix][iy] = vReal[scaleNX][scaleNY] + \r
 710           dXFrac * (vReal[scaleNX+1][scaleNY] - vReal[scaleNX][scaleNY]) + \r
 711           dYFrac * (vReal[scaleNX][scaleNY+1] - vReal[scaleNX][scaleNY]);\r
 712         if (result.isComplex()) {\r
 713           if (isComplex())\r
 714             vResultImag[ix][iy] = vImag[scaleNX][scaleNY] + \r
 715             dXFrac * (vImag[scaleNX+1][scaleNY] - vImag[scaleNX][scaleNY]) + \r
 716             dYFrac * (vImag[scaleNX][scaleNY+1] - vImag[scaleNX][scaleNY]);\r
 717           else\r
 718             vResultImag[ix][iy] = 0;\r
 719         }\r
 720       }\r
 721     }\r
 722   }\r
 723   \r
 724   return true;\r
 725 }\r
 726 \r
 727 #ifdef HAVE_FFTW\r
 728 bool\r
 729 ImageFile::fft (ImageFile& result) const\r
 730 {\r
 731   if (m_nx != result.nx() || m_ny != result.ny()) {\r
 732     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
 733     return false;\r
 734   }\r
 735   \r
 736   if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {\r
 737     if (! result.convertRealToComplex ())\r
 738       return false;\r
 739   }\r
 740   \r
 741   fftw_complex* in = new fftw_complex [m_nx * m_ny];\r
 742   \r
 743   ImageFileArrayConst vReal = getArray();\r
 744   ImageFileArrayConst vImag = getImaginaryArray();\r
 745   \r
 746   unsigned int ix, iy;\r
 747   unsigned int iArray = 0;\r
 748   for (ix = 0; ix < m_nx; ix++)\r
 749     for (iy = 0; iy < m_ny; iy++) {\r
 750       in[iArray].re = vReal[ix][iy];\r
 751       if (isComplex())\r
 752         in[iArray].im = vImag[ix][iy];\r
 753       else\r
 754         in[iArray].im = 0;\r
 755       iArray++;\r
 756     }\r
 757     \r
 758     fftwnd_plan plan = fftw2d_create_plan (m_nx, m_ny, FFTW_FORWARD, FFTW_IN_PLACE);\r
 759     \r
 760     fftwnd_one (plan, in, NULL);\r
 761     \r
 762     ImageFileArray vRealResult = result.getArray();\r
 763     ImageFileArray vImagResult = result.getImaginaryArray();\r
 764     iArray = 0;\r
 765     unsigned int iScale = m_nx * m_ny;\r
 766     for (ix = 0; ix < m_nx; ix++)\r
 767       for (iy = 0; iy < m_ny; iy++) {\r
 768         vRealResult[ix][iy] = in[iArray].re / iScale;\r
 769         vImagResult[ix][iy] = in[iArray].im / iScale;\r
 770         iArray++;\r
 771       }\r
 772       \r
 773       fftwnd_destroy_plan (plan);\r
 774       delete in;\r
 775       \r
 776       \r
 777       Fourier::shuffleFourierToNaturalOrder (result);\r
 778       \r
 779       return true;\r
 780 }\r
 781 \r
 782 \r
 783 bool\r
 784 ImageFile::ifft (ImageFile& result) const\r
 785 {\r
 786   if (m_nx != result.nx() || m_ny != result.ny()) {\r
 787     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
 788     return false;\r
 789   }\r
 790   \r
 791   if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {\r
 792     if (! result.convertRealToComplex ())\r
 793       return false;\r
 794   }\r
 795   \r
 796   ImageFileArrayConst vReal = getArray();\r
 797   ImageFileArrayConst vImag = getImaginaryArray();\r
 798   ImageFileArray vRealResult = result.getArray();\r
 799   ImageFileArray vImagResult = result.getImaginaryArray();\r
 800   unsigned int ix, iy;\r
 801   for (ix = 0; ix < m_nx; ix++)\r
 802     for (iy = 0; iy < m_ny; iy++) {\r
 803       vRealResult[ix][iy] = vReal[ix][iy];\r
 804       if (isComplex()) \r
 805         vImagResult[ix][iy] = vImag[ix][iy];\r
 806       else\r
 807         vImagResult[ix][iy] = 0;\r
 808     }\r
 809     \r
 810     Fourier::shuffleNaturalToFourierOrder (result);\r
 811     \r
 812     fftw_complex* in = new fftw_complex [m_nx * m_ny];\r
 813     \r
 814     unsigned int iArray = 0;\r
 815     for (ix = 0; ix < m_nx; ix++)\r
 816       for (iy = 0; iy < m_ny; iy++) {\r
 817         in[iArray].re = vRealResult[ix][iy];\r
 818         in[iArray].im = vImagResult[ix][iy];\r
 819         iArray++;\r
 820       }\r
 821       \r
 822       fftwnd_plan plan = fftw2d_create_plan (m_nx, m_ny, FFTW_BACKWARD, FFTW_IN_PLACE);\r
 823       \r
 824       fftwnd_one (plan, in, NULL);\r
 825       \r
 826       iArray = 0;\r
 827       for (ix = 0; ix < m_nx; ix++)\r
 828         for (iy = 0; iy < m_ny; iy++) {\r
 829           vRealResult[ix][iy] = in[iArray].re;\r
 830           vImagResult[ix][iy] = in[iArray].im;\r
 831           iArray++;\r
 832         }\r
 833         \r
 834         fftwnd_destroy_plan (plan);\r
 835         \r
 836         delete in;\r
 837         \r
 838         return true;\r
 839 }\r
 840 #endif // HAVE_FFTW\r
 841 \r
 842 \r
 843 \r
 844 bool\r
 845 ImageFile::fourier (ImageFile& result) const\r
 846 {\r
 847   if (m_nx != result.nx() || m_ny != result.ny()) {\r
 848     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
 849     return false;\r
 850   }\r
 851   \r
 852   if (! result.isComplex())\r
 853     if (! result.convertRealToComplex ())\r
 854       return false;\r
 855     \r
 856     ImageFileArrayConst vLHS = getArray();\r
 857     ImageFileArrayConst vLHSImag = getImaginaryArray();\r
 858     ImageFileArray vRealResult = result.getArray();\r
 859     ImageFileArray vImagResult = result.getImaginaryArray();\r
 860     \r
 861     unsigned int ix, iy;\r
 862     \r
 863     // alloc output matrix\r
 864     CTSimComplex** complexOut = new CTSimComplex* [m_nx];\r
 865     for (ix = 0; ix < m_nx; ix++)\r
 866       complexOut[ix] = new CTSimComplex [m_ny];\r
 867     \r
 868     // fourier each x column\r
 869 #if 1\r
 870     CTSimComplex* pY = new CTSimComplex [m_ny];\r
 871     for (ix = 0; ix < m_nx; ix++) {\r
 872       for (iy = 0; iy < m_ny; iy++) {\r
 873         double dImag = 0;\r
 874         if (isComplex())\r
 875           dImag = vLHSImag[ix][iy];\r
 876         pY[iy] = complex<double>(vLHS[ix][iy], dImag);\r
 877       } \r
 878       ProcessSignal::finiteFourierTransform (pY, complexOut[ix], m_ny,  ProcessSignal::FORWARD);\r
 879     }\r
 880 #else\r
 881     double* pY = new double [m_ny];\r
 882     for (ix = 0; ix < m_nx; ix++) {\r
 883       for (iy = 0; iy < m_ny; iy++) {\r
 884         pY[iy] = vLHS[ix][iy];\r
 885       } \r
 886       ProcessSignal::finiteFourierTransform (pY, complexOut[ix], m_ny,  ProcessSignal::FORWARD);\r
 887     }\r
 888 #endif\r
 889     delete [] pY;\r
 890     \r
 891     // fourier each y row\r
 892     CTSimComplex* pX = new CTSimComplex [m_nx];\r
 893     CTSimComplex* complexOutRow = new CTSimComplex [m_nx];\r
 894     for (iy = 0; iy < m_ny; iy++) {\r
 895       for (ix = 0; ix < m_nx; ix++)\r
 896         pX[ix] = complexOut[ix][iy];\r
 897       ProcessSignal::finiteFourierTransform (pX, complexOutRow, m_nx, ProcessSignal::FORWARD);\r
 898       for (ix = 0; ix < m_nx; ix++)\r
 899         complexOut[ix][iy] = complexOutRow[ix];\r
 900     }\r
 901     delete [] pX;\r
 902     delete [] complexOutRow;\r
 903     \r
 904     for (ix = 0; ix < m_nx; ix++)\r
 905       for (iy = 0; iy < m_ny; iy++) {\r
 906         vRealResult[ix][iy] = complexOut[ix][iy].real();\r
 907         vImagResult[ix][iy] = complexOut[ix][iy].imag();\r
 908       }\r
 909       \r
 910       Fourier::shuffleFourierToNaturalOrder (result);\r
 911       \r
 912       // delete complexOut matrix\r
 913       for (ix = 0; ix < m_nx; ix++)\r
 914         delete [] complexOut[ix];\r
 915       delete [] complexOut;\r
 916       \r
 917       return true;\r
 918 }\r
 919 \r
 920 bool\r
 921 ImageFile::inverseFourier (ImageFile& result) const\r
 922 {\r
 923   if (m_nx != result.nx() || m_ny != result.ny()) {\r
 924     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
 925     return false;\r
 926   }\r
 927   \r
 928   if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {\r
 929     if (! result.convertRealToComplex ())\r
 930       return false;\r
 931   }\r
 932   \r
 933   ImageFileArrayConst vLHSReal = getArray();\r
 934   ImageFileArrayConst vLHSImag = getImaginaryArray();\r
 935   ImageFileArray vRealResult = result.getArray();\r
 936   ImageFileArray vImagResult = result.getImaginaryArray();\r
 937   \r
 938   unsigned int ix, iy;\r
 939   // alloc 2d complex output matrix\r
 940   CTSimComplex** complexOut = new CTSimComplex* [m_nx];\r
 941   for (ix = 0; ix < m_nx; ix++)\r
 942     complexOut[ix] = new CTSimComplex [m_ny];\r
 943   \r
 944   // put input image into result\r
 945   for (ix = 0; ix < m_nx; ix++)\r
 946     for (iy = 0; iy < m_ny; iy++) {\r
 947       vRealResult[ix][iy] = vLHSReal[ix][iy];\r
 948       if (isComplex())\r
 949         vImagResult[ix][iy] = vLHSImag[ix][iy];\r
 950       else\r
 951         vImagResult[ix][iy] = 0;\r
 952     }\r
 953     \r
 954     Fourier::shuffleNaturalToFourierOrder (result);\r
 955     \r
 956     // ifourier each x column\r
 957     CTSimComplex* pCol = new CTSimComplex [m_ny];\r
 958     for (ix = 0; ix < m_nx; ix++) {\r
 959       for (iy = 0; iy < m_ny; iy++) {\r
 960         pCol[iy] = std::complex<double> (vRealResult[ix][iy], vImagResult[ix][iy]);\r
 961       }\r
 962       ProcessSignal::finiteFourierTransform (pCol, complexOut[ix], m_ny,  ProcessSignal::BACKWARD);\r
 963     }\r
 964     delete [] pCol;\r
 965     \r
 966     // ifourier each y row\r
 967     CTSimComplex* complexInRow = new CTSimComplex [m_nx];\r
 968     CTSimComplex* complexOutRow = new CTSimComplex [m_nx];\r
 969     for (iy = 0; iy < m_ny; iy++) {\r
 970       for (ix = 0; ix < m_nx; ix++)\r
 971         complexInRow[ix] = complexOut[ix][iy];\r
 972       ProcessSignal::finiteFourierTransform (complexInRow, complexOutRow, m_nx, ProcessSignal::BACKWARD);\r
 973       for (ix = 0; ix < m_nx; ix++)\r
 974         complexOut[ix][iy] = complexOutRow[ix];\r
 975     }\r
 976     delete [] complexInRow;\r
 977     delete [] complexOutRow;\r
 978     \r
 979     for (ix = 0; ix < m_nx; ix++)\r
 980       for (iy = 0; iy < m_ny; iy++) {\r
 981         vRealResult[ix][iy] = complexOut[ix][iy].real();\r
 982         vImagResult[ix][iy] = complexOut[ix][iy].imag();\r
 983       }\r
 984       \r
 985       // delete complexOut matrix\r
 986       for (ix = 0; ix < m_nx; ix++)\r
 987         delete [] complexOut[ix];\r
 988       delete [] complexOut;\r
 989       \r
 990       return true;\r
 991 }\r
 992 \r
 993 \r
 994 bool\r
 995 ImageFile::magnitude (ImageFile& result) const\r
 996 {\r
 997   if (m_nx != result.nx() || m_ny != result.ny()) {\r
 998     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
 999     return false;\r
1000   }\r
1001   \r
1002   ImageFileArray vReal = getArray();\r
1003   ImageFileArray vImag = getImaginaryArray();\r
1004   ImageFileArray vRealResult = result.getArray();\r
1005   \r
1006   for (unsigned int ix = 0; ix < m_nx; ix++)\r
1007     for (unsigned int iy = 0; iy < m_ny; iy++) {\r
1008       if (isComplex()) \r
1009         vRealResult[ix][iy] = ::sqrt (vReal[ix][iy] * vReal[ix][iy] + vImag[ix][iy] * vImag[ix][iy]);\r
1010       else\r
1011         vRealResult[ix][iy] = vReal[ix][iy];\r
1012     }\r
1013     \r
1014     if (result.isComplex())\r
1015       result.convertComplexToReal();\r
1016     \r
1017     return true;\r
1018 }\r
1019 \r
1020 bool\r
1021 ImageFile::phase (ImageFile& result) const\r
1022 {\r
1023   if (m_nx != result.nx() || m_ny != result.ny()) {\r
1024     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
1025     return false;\r
1026   }\r
1027   \r
1028   ImageFileArray vReal = getArray();\r
1029   ImageFileArray vImag = getImaginaryArray();\r
1030   ImageFileArray vRealResult = result.getArray();\r
1031   \r
1032   for (unsigned int ix = 0; ix < m_nx; ix++)\r
1033     for (unsigned int iy = 0; iy < m_ny; iy++) {\r
1034       if (isComplex())\r
1035         vRealResult[ix][iy] = ::atan2 (vImag[ix][iy], vReal[ix][iy]);\r
1036       else\r
1037         vRealResult[ix][iy] = 0;\r
1038     }\r
1039     \r
1040     if (result.isComplex())\r
1041       result.convertComplexToReal();\r
1042     \r
1043     return true;\r
1044 }\r
1045 \r
1046 bool\r
1047 ImageFile::square (ImageFile& result) const\r
1048 {\r
1049   if (m_nx != result.nx() || m_ny != result.ny()) {\r
1050     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
1051     return false;\r
1052   }\r
1053   \r
1054   if (isComplex() && ! result.isComplex())\r
1055     result.convertRealToComplex();\r
1056   \r
1057   ImageFileArrayConst vLHS = getArray();\r
1058   ImageFileArrayConst vLHSImag = getImaginaryArray();\r
1059   ImageFileArray vResult = result.getArray();\r
1060   ImageFileArray vResultImag = result.getImaginaryArray();\r
1061   \r
1062   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
1063     for (unsigned int iy = 0; iy < m_ny; iy++) {\r
1064       if (result.isComplex()) {\r
1065         double dImag = 0;\r
1066         if (isComplex())\r
1067           dImag = vLHSImag[ix][iy];\r
1068         std::complex<double> cLHS (vLHS[ix][iy], dImag);\r
1069         std::complex<double> cResult = cLHS * cLHS;\r
1070         vResult[ix][iy] = cResult.real();\r
1071         vResultImag[ix][iy] = cResult.imag();\r
1072       } else\r
1073         vResult[ix][iy] = vLHS[ix][iy] * vLHS[ix][iy];\r
1074     }\r
1075   }\r
1076   \r
1077   \r
1078   return true;\r
1079 }\r
1080 \r
1081 int\r
1082 ImageFile::convertFormatNameToID (const char* const formatName)\r
1083 {\r
1084   int formatID = FORMAT_INVALID;\r
1085   \r
1086   for (int i = 0; i < s_iFormatCount; i++)\r
1087     if (strcasecmp (formatName, s_aszFormatName[i]) == 0) {\r
1088       formatID = i;\r
1089       break;\r
1090     }\r
1091     \r
1092     return (formatID);\r
1093 }\r
1094 \r
1095 const char*\r
1096 ImageFile::convertFormatIDToName (int formatID)\r
1097 {\r
1098   static const char *formatName = "";\r
1099   \r
1100   if (formatID >= 0 && formatID < s_iFormatCount)\r
1101     return (s_aszFormatName[formatID]);\r
1102   \r
1103   return (formatName);\r
1104 }\r
1105 \r
1106 const char*\r
1107 ImageFile::convertFormatIDToTitle (const int formatID)\r
1108 {\r
1109   static const char *formatTitle = "";\r
1110   \r
1111   if (formatID >= 0 && formatID < s_iFormatCount)\r
1112     return (s_aszFormatTitle[formatID]);\r
1113   \r
1114   return (formatTitle);\r
1115 }\r
1116 \r
1117 bool\r
1118 ImageFile::exportImage (const char* const pszFormat, const char* const pszFilename, int nxcell, int nycell, double densmin, double densmax)\r
1119 {\r
1120   int iFormatID = convertFormatNameToID (pszFormat);\r
1121   if (iFormatID == FORMAT_INVALID) {\r
1122     sys_error (ERR_SEVERE, "Invalid format %s [ImageFile::exportImage]", pszFormat);\r
1123     return false;\r
1124   }\r
1125   \r
1126   if (iFormatID == FORMAT_PGM)\r
1127     return writeImagePGM (pszFilename, nxcell, nycell, densmin, densmax);\r
1128   else if (iFormatID == FORMAT_PGMASCII)\r
1129     return writeImagePGMASCII (pszFilename, nxcell, nycell, densmin, densmax);\r
1130   else if (iFormatID == FORMAT_PNG)\r
1131     return writeImagePNG (pszFilename, 8, nxcell, nycell, densmin, densmax);\r
1132   else if (iFormatID == FORMAT_PNG16)\r
1133     return writeImagePNG (pszFilename, 16, nxcell, nycell, densmin, densmax);\r
1134   \r
1135   sys_error (ERR_SEVERE, "Invalid format %s [ImageFile::exportImage]", pszFormat);\r
1136   return false;\r
1137 }\r
1138 \r
1139 bool
1140 ImageFile::writeImagePGM (const char* const outfile, int nxcell, int nycell, double densmin, double densmax)
1141 {
1142   FILE *fp;
1143   int nx = m_nx;
1144   int ny = m_ny;
1145   ImageFileArray v = getArray();
1146
1147   unsigned char* rowp = new unsigned char [nx * nxcell];
1148
1149   if ((fp = fopen (outfile, "wb")) == NULL)
1150     return false;
1151
1152   fprintf(fp, "P5\n");
1153   fprintf(fp, "%d %d\n", nx, ny);
1154   fprintf(fp, "255\n");
1155
1156   for (int irow = ny - 1; irow >= 0; irow--) {
1157     for (int icol = 0; icol < nx; icol++) {
1158       int pos = icol * nxcell;
1159       double dens = (v[icol][irow] - densmin) / (densmax - densmin);
1160       dens = clamp (dens, 0., 1.);
1161       for (int p = pos; p < pos + nxcell; p++) {
1162         rowp[p] = static_cast<unsigned int> (dens * 255.);
1163       }
1164     }
1165     for (int ir = 0; ir < nycell; ir++) {
1166       for (int ic = 0; ic < nx * nxcell; ic++)
1167         fputc( rowp[ic], fp );
1168     }
1169   }
1170   \r
1171   delete rowp;
1172   fclose(fp);\r
1173   \r
1174   return true;
1175 }
1176
1177 bool
1178 ImageFile::writeImagePGMASCII (const char* const outfile, int nxcell, int nycell, double densmin, double densmax)
1179 {
1180   FILE *fp;
1181   int nx = m_nx;
1182   int ny = m_ny;
1183   ImageFileArray v = getArray();
1184
1185   unsigned char* rowp = new unsigned char [nx * nxcell];
1186
1187   if ((fp = fopen (outfile, "wb")) == NULL)
1188     return false;
1189
1190   fprintf(fp, "P2\n");
1191   fprintf(fp, "%d %d\n", nx, ny);
1192   fprintf(fp, "255\n");
1193
1194   for (int irow = ny - 1; irow >= 0; irow--) {
1195     for (int icol = 0; icol < nx; icol++) {
1196       int pos = icol * nxcell;
1197       double dens = (v[icol][irow] - densmin) / (densmax - densmin);
1198       dens = clamp (dens, 0., 1.);
1199       for (int p = pos; p < pos + nxcell; p++) {
1200         rowp[p] = static_cast<unsigned int> (dens * 255.);
1201       }
1202     }
1203     for (int ir = 0; ir < nycell; ir++) {
1204       for (int ic = 0; ic < nx * nxcell; ic++)
1205         fprintf(fp, "%d ", rowp[ic]);
1206       fprintf(fp, "\n");
1207     }
1208   }
1209   \r
1210   delete rowp;
1211   fclose(fp);\r
1212   \r
1213   return true;
1214 }
1215
1216
1217 #ifdef HAVE_PNG
1218 bool
1219 ImageFile::writeImagePNG (const char* const outfile, int bitdepth, int nxcell, int nycell, double densmin, double densmax)
1220 {
1221   double max_out_level = (1 << bitdepth) - 1;
1222   int nx = m_nx;
1223   int ny = m_ny;
1224   ImageFileArray v = getArray();
1225
1226   unsigned char* rowp = new unsigned char [nx * nxcell * (bitdepth / 8)];
1227
1228   FILE *fp = fopen (outfile, "wb");\r
1229   if (! fp)
1230     return false;
1231
1232   png_structp png_ptr = png_create_write_struct (PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
1233   if (! png_ptr)
1234     return false;
1235
1236   png_infop info_ptr = png_create_info_struct (png_ptr);
1237   if (! info_ptr) {
1238     png_destroy_write_struct (&png_ptr, (png_infopp) NULL);
1239     fclose (fp);
1240     return false;
1241   }
1242
1243   if (setjmp (png_ptr->jmpbuf)) {
1244     png_destroy_write_struct (&png_ptr, &info_ptr);
1245     fclose (fp);
1246     return false;
1247   }
1248
1249   png_init_io(png_ptr, fp);
1250
1251   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);
1252
1253   png_write_info(png_ptr, info_ptr);
1254   for (int irow = ny - 1; irow >= 0; irow--) {
1255     png_bytep row_pointer = rowp;
1256
1257     for (int icol = 0; icol < nx; icol++) {
1258       int pos = icol * nxcell;
1259       double dens = (v[icol][irow] - densmin) / (densmax - densmin);
1260       dens = clamp (dens, 0., 1.);
1261       unsigned int outval = static_cast<unsigned int> (dens * max_out_level);
1262
1263       for (int p = pos; p < pos + nxcell; p++) {
1264         if (bitdepth == 8)
1265           rowp[p] = outval;
1266         else {
1267           int rowpos = p * 2;
1268           rowp[rowpos] = (outval >> 8) & 0xFF;
1269           rowp[rowpos+1] = (outval & 0xFF);
1270         }
1271       }
1272     }
1273     for (int ir = 0; ir < nycell; ir++)
1274       png_write_rows (png_ptr, &row_pointer, 1);
1275   }
1276
1277   png_write_end (png_ptr, info_ptr);
1278   png_destroy_write_struct (&png_ptr, &info_ptr);
1279   delete rowp;\r
1280
1281   fclose(fp);\r
1282   \r
1283   return true;
1284 }
1285 #endif
1286
1287 #ifdef HAVE_GD
1288 #include "gd.h"
1289 static const int N_GRAYSCALE=256;
1290
1291 bool
1292 ImageFile::writeImageGIF (const char* const outfile, int nxcell, int nycell, double densmin, double densmax)
1293 {
1294   int gs_indices[N_GRAYSCALE];
1295   int nx = m_nx;
1296   int ny = m_ny;
1297   ImageFileArray v = getArray();
1298
1299   unsigned char* rowp = new unsigned char [nx * nxcell];
1300
1301   gdImagePtr gif = gdImageCreate(nx * nxcell, ny * nycell);
1302   for (int i = 0; i < N_GRAYSCALE; i++)
1303     gs_indices[i] = gdImageColorAllocate(gif, i, i, i);
1304
1305   int lastrow = ny * nycell - 1;
1306   for (int irow = 0; irow < ny; irow++) {
1307     int rpos = irow * nycell;
1308     for (int ir = rpos; ir < rpos + nycell; ir++) {
1309       for (int icol = 0; icol < nx; icol++) {
1310         int cpos = icol * nxcell;
1311         double dens = (v[icol][irow] - densmin) / (densmax - densmin);
1312         dens = clamp(dens, 0., 1.);
1313         for (int ic = cpos; ic < cpos + nxcell; ic++) {
1314           rowp[ic] = (unsigned int) (dens * (double) (N_GRAYSCALE - 1));
1315           gdImageSetPixel(gif, ic, lastrow - ir, gs_indices[rowp[ic]]);
1316         }
1317       }
1318     }
1319   }
1320
1321   FILE *out;\r
1322   if ((out = fopen (outfile,"w")) == NULL) {
1323     sys_error(ERR_FATAL, "Error opening output file %s for writing", outfile);
1324     return false;
1325   }
1326   gdImageGif(gif,out);
1327   fclose(out);
1328   gdImageDestroy(gif);\r
1329   \r
1330   delete rowp;\r
1331   \r
1332   return true;
1333 }
1334 #endif
1335