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