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.31 2001/01/02 08:15:02 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     CTSimComplex* pY = new CTSimComplex [m_ny];\r
 870     for (ix = 0; ix < m_nx; ix++) {\r
 871       for (iy = 0; iy < m_ny; iy++) {\r
 872         double dImag = 0;\r
 873         if (isComplex())\r
 874           dImag = vLHSImag[ix][iy];\r
 875         pY[iy] = std::complex<double>(vLHS[ix][iy], dImag);\r
 876       } \r
 877       ProcessSignal::finiteFourierTransform (pY, complexOut[ix], m_ny,  ProcessSignal::FORWARD);\r
 878     }\r
 879     delete [] pY;\r
 880     \r
 881     // fourier each y row\r
 882     CTSimComplex* pX = new CTSimComplex [m_nx];\r
 883     CTSimComplex* complexOutRow = new CTSimComplex [m_nx];\r
 884     for (iy = 0; iy < m_ny; iy++) {\r
 885       for (ix = 0; ix < m_nx; ix++)\r
 886         pX[ix] = complexOut[ix][iy];\r
 887       ProcessSignal::finiteFourierTransform (pX, complexOutRow, m_nx, ProcessSignal::FORWARD);\r
 888       for (ix = 0; ix < m_nx; ix++)\r
 889         complexOut[ix][iy] = complexOutRow[ix];\r
 890     }\r
 891     delete [] pX;\r
 892     delete [] complexOutRow;\r
 893     \r
 894     for (ix = 0; ix < m_nx; ix++)\r
 895       for (iy = 0; iy < m_ny; iy++) {\r
 896         vRealResult[ix][iy] = complexOut[ix][iy].real();\r
 897         vImagResult[ix][iy] = complexOut[ix][iy].imag();\r
 898       }\r
 899       \r
 900       Fourier::shuffleFourierToNaturalOrder (result);\r
 901       \r
 902       // delete complexOut matrix\r
 903       for (ix = 0; ix < m_nx; ix++)\r
 904         delete [] complexOut[ix];\r
 905       delete [] complexOut;\r
 906       \r
 907       return true;\r
 908 }\r
 909 \r
 910 bool\r
 911 ImageFile::inverseFourier (ImageFile& result) const\r
 912 {\r
 913   if (m_nx != result.nx() || m_ny != result.ny()) {\r
 914     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
 915     return false;\r
 916   }\r
 917   \r
 918   if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {\r
 919     if (! result.convertRealToComplex ())\r
 920       return false;\r
 921   }\r
 922   \r
 923   ImageFileArrayConst vLHSReal = getArray();\r
 924   ImageFileArrayConst vLHSImag = getImaginaryArray();\r
 925   ImageFileArray vRealResult = result.getArray();\r
 926   ImageFileArray vImagResult = result.getImaginaryArray();\r
 927   \r
 928   unsigned int ix, iy;\r
 929   // alloc 2d complex output matrix\r
 930   CTSimComplex** complexOut = new CTSimComplex* [m_nx];\r
 931   for (ix = 0; ix < m_nx; ix++)\r
 932     complexOut[ix] = new CTSimComplex [m_ny];\r
 933   \r
 934   // put input image into result\r
 935   for (ix = 0; ix < m_nx; ix++)\r
 936     for (iy = 0; iy < m_ny; iy++) {\r
 937       vRealResult[ix][iy] = vLHSReal[ix][iy];\r
 938       if (isComplex())\r
 939         vImagResult[ix][iy] = vLHSImag[ix][iy];\r
 940       else\r
 941         vImagResult[ix][iy] = 0;\r
 942     }\r
 943     \r
 944     Fourier::shuffleNaturalToFourierOrder (result);\r
 945     \r
 946     // ifourier each x column\r
 947     CTSimComplex* pCol = new CTSimComplex [m_ny];\r
 948     for (ix = 0; ix < m_nx; ix++) {\r
 949       for (iy = 0; iy < m_ny; iy++) {\r
 950         pCol[iy] = std::complex<double> (vRealResult[ix][iy], vImagResult[ix][iy]);\r
 951       }\r
 952       ProcessSignal::finiteFourierTransform (pCol, complexOut[ix], m_ny,  ProcessSignal::BACKWARD);\r
 953     }\r
 954     delete [] pCol;\r
 955     \r
 956     // ifourier each y row\r
 957     CTSimComplex* complexInRow = new CTSimComplex [m_nx];\r
 958     CTSimComplex* complexOutRow = new CTSimComplex [m_nx];\r
 959     for (iy = 0; iy < m_ny; iy++) {\r
 960       for (ix = 0; ix < m_nx; ix++)\r
 961         complexInRow[ix] = complexOut[ix][iy];\r
 962       ProcessSignal::finiteFourierTransform (complexInRow, complexOutRow, m_nx, ProcessSignal::BACKWARD);\r
 963       for (ix = 0; ix < m_nx; ix++)\r
 964         complexOut[ix][iy] = complexOutRow[ix];\r
 965     }\r
 966     delete [] complexInRow;\r
 967     delete [] complexOutRow;\r
 968     \r
 969     for (ix = 0; ix < m_nx; ix++)\r
 970       for (iy = 0; iy < m_ny; iy++) {\r
 971         vRealResult[ix][iy] = complexOut[ix][iy].real();\r
 972         vImagResult[ix][iy] = complexOut[ix][iy].imag();\r
 973       }\r
 974       \r
 975       // delete complexOut matrix\r
 976       for (ix = 0; ix < m_nx; ix++)\r
 977         delete [] complexOut[ix];\r
 978       delete [] complexOut;\r
 979       \r
 980       return true;\r
 981 }\r
 982 \r
 983 \r
 984 bool\r
 985 ImageFile::magnitude (ImageFile& result) const\r
 986 {\r
 987   if (m_nx != result.nx() || m_ny != result.ny()) {\r
 988     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
 989     return false;\r
 990   }\r
 991   \r
 992   ImageFileArray vReal = getArray();\r
 993   ImageFileArray vImag = getImaginaryArray();\r
 994   ImageFileArray vRealResult = result.getArray();\r
 995   \r
 996   for (unsigned int ix = 0; ix < m_nx; ix++)\r
 997     for (unsigned int iy = 0; iy < m_ny; iy++) {\r
 998       if (isComplex()) \r
 999         vRealResult[ix][iy] = ::sqrt (vReal[ix][iy] * vReal[ix][iy] + vImag[ix][iy] * vImag[ix][iy]);\r
1000       else\r
1001         vRealResult[ix][iy] = vReal[ix][iy];\r
1002     }\r
1003     \r
1004     if (result.isComplex())\r
1005       result.convertComplexToReal();\r
1006     \r
1007     return true;\r
1008 }\r
1009 \r
1010 bool\r
1011 ImageFile::phase (ImageFile& result) const\r
1012 {\r
1013   if (m_nx != result.nx() || m_ny != result.ny()) {\r
1014     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
1015     return false;\r
1016   }\r
1017   \r
1018   ImageFileArray vReal = getArray();\r
1019   ImageFileArray vImag = getImaginaryArray();\r
1020   ImageFileArray vRealResult = result.getArray();\r
1021   \r
1022   for (unsigned int ix = 0; ix < m_nx; ix++)\r
1023     for (unsigned int iy = 0; iy < m_ny; iy++) {\r
1024       if (isComplex())\r
1025         vRealResult[ix][iy] = ::atan2 (vImag[ix][iy], vReal[ix][iy]);\r
1026       else\r
1027         vRealResult[ix][iy] = 0;\r
1028     }\r
1029     \r
1030     if (result.isComplex())\r
1031       result.convertComplexToReal();\r
1032     \r
1033     return true;\r
1034 }\r
1035 \r
1036 bool\r
1037 ImageFile::square (ImageFile& result) const\r
1038 {\r
1039   if (m_nx != result.nx() || m_ny != result.ny()) {\r
1040     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
1041     return false;\r
1042   }\r
1043   \r
1044   if (isComplex() && ! result.isComplex())\r
1045     result.convertRealToComplex();\r
1046   \r
1047   ImageFileArrayConst vLHS = getArray();\r
1048   ImageFileArrayConst vLHSImag = getImaginaryArray();\r
1049   ImageFileArray vResult = result.getArray();\r
1050   ImageFileArray vResultImag = result.getImaginaryArray();\r
1051   \r
1052   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
1053     for (unsigned int iy = 0; iy < m_ny; iy++) {\r
1054       if (result.isComplex()) {\r
1055         double dImag = 0;\r
1056         if (isComplex())\r
1057           dImag = vLHSImag[ix][iy];\r
1058         std::complex<double> cLHS (vLHS[ix][iy], dImag);\r
1059         std::complex<double> cResult = cLHS * cLHS;\r
1060         vResult[ix][iy] = cResult.real();\r
1061         vResultImag[ix][iy] = cResult.imag();\r
1062       } else\r
1063         vResult[ix][iy] = vLHS[ix][iy] * vLHS[ix][iy];\r
1064     }\r
1065   }\r
1066   \r
1067   \r
1068   return true;\r
1069 }\r
1070 \r
1071 int\r
1072 ImageFile::convertFormatNameToID (const char* const formatName)\r
1073 {\r
1074   int formatID = FORMAT_INVALID;\r
1075   \r
1076   for (int i = 0; i < s_iFormatCount; i++)\r
1077     if (strcasecmp (formatName, s_aszFormatName[i]) == 0) {\r
1078       formatID = i;\r
1079       break;\r
1080     }\r
1081     \r
1082     return (formatID);\r
1083 }\r
1084 \r
1085 const char*\r
1086 ImageFile::convertFormatIDToName (int formatID)\r
1087 {\r
1088   static const char *formatName = "";\r
1089   \r
1090   if (formatID >= 0 && formatID < s_iFormatCount)\r
1091     return (s_aszFormatName[formatID]);\r
1092   \r
1093   return (formatName);\r
1094 }\r
1095 \r
1096 const char*\r
1097 ImageFile::convertFormatIDToTitle (const int formatID)\r
1098 {\r
1099   static const char *formatTitle = "";\r
1100   \r
1101   if (formatID >= 0 && formatID < s_iFormatCount)\r
1102     return (s_aszFormatTitle[formatID]);\r
1103   \r
1104   return (formatTitle);\r
1105 }\r
1106 \r
1107 bool\r
1108 ImageFile::exportImage (const char* const pszFormat, const char* const pszFilename, int nxcell, int nycell, double densmin, double densmax)\r
1109 {\r
1110   int iFormatID = convertFormatNameToID (pszFormat);\r
1111   if (iFormatID == FORMAT_INVALID) {\r
1112     sys_error (ERR_SEVERE, "Invalid format %s [ImageFile::exportImage]", pszFormat);\r
1113     return false;\r
1114   }\r
1115   \r
1116   if (iFormatID == FORMAT_PGM)\r
1117     return writeImagePGM (pszFilename, nxcell, nycell, densmin, densmax);\r
1118   else if (iFormatID == FORMAT_PGMASCII)\r
1119     return writeImagePGMASCII (pszFilename, nxcell, nycell, densmin, densmax);\r
1120   else if (iFormatID == FORMAT_PNG)\r
1121     return writeImagePNG (pszFilename, 8, nxcell, nycell, densmin, densmax);\r
1122   else if (iFormatID == FORMAT_PNG16)\r
1123     return writeImagePNG (pszFilename, 16, nxcell, nycell, densmin, densmax);\r
1124   \r
1125   sys_error (ERR_SEVERE, "Invalid format %s [ImageFile::exportImage]", pszFormat);\r
1126   return false;\r
1127 }\r
1128 \r
1129 bool
1130 ImageFile::writeImagePGM (const char* const outfile, int nxcell, int nycell, double densmin, double densmax)
1131 {
1132   FILE *fp;
1133   int nx = m_nx;
1134   int ny = m_ny;
1135   ImageFileArray v = getArray();
1136
1137   unsigned char* rowp = new unsigned char [nx * nxcell];
1138
1139   if ((fp = fopen (outfile, "wb")) == NULL)
1140     return false;
1141
1142   fprintf(fp, "P5\n");
1143   fprintf(fp, "%d %d\n", nx, ny);
1144   fprintf(fp, "255\n");
1145
1146   for (int irow = ny - 1; irow >= 0; irow--) {
1147     for (int icol = 0; icol < nx; icol++) {
1148       int pos = icol * nxcell;
1149       double dens = (v[icol][irow] - densmin) / (densmax - densmin);
1150       dens = clamp (dens, 0., 1.);
1151       for (int p = pos; p < pos + nxcell; p++) {
1152         rowp[p] = static_cast<unsigned int> (dens * 255.);
1153       }
1154     }
1155     for (int ir = 0; ir < nycell; ir++) {
1156       for (int ic = 0; ic < nx * nxcell; ic++)
1157         fputc( rowp[ic], fp );
1158     }
1159   }
1160   \r
1161   delete rowp;
1162   fclose(fp);\r
1163   \r
1164   return true;
1165 }
1166
1167 bool
1168 ImageFile::writeImagePGMASCII (const char* const outfile, int nxcell, int nycell, double densmin, double densmax)
1169 {
1170   FILE *fp;
1171   int nx = m_nx;
1172   int ny = m_ny;
1173   ImageFileArray v = getArray();
1174
1175   unsigned char* rowp = new unsigned char [nx * nxcell];
1176
1177   if ((fp = fopen (outfile, "wb")) == NULL)
1178     return false;
1179
1180   fprintf(fp, "P2\n");
1181   fprintf(fp, "%d %d\n", nx, ny);
1182   fprintf(fp, "255\n");
1183
1184   for (int irow = ny - 1; irow >= 0; irow--) {
1185     for (int icol = 0; icol < nx; icol++) {
1186       int pos = icol * nxcell;
1187       double dens = (v[icol][irow] - densmin) / (densmax - densmin);
1188       dens = clamp (dens, 0., 1.);
1189       for (int p = pos; p < pos + nxcell; p++) {
1190         rowp[p] = static_cast<unsigned int> (dens * 255.);
1191       }
1192     }
1193     for (int ir = 0; ir < nycell; ir++) {
1194       for (int ic = 0; ic < nx * nxcell; ic++)
1195         fprintf(fp, "%d ", rowp[ic]);
1196       fprintf(fp, "\n");
1197     }
1198   }
1199   \r
1200   delete rowp;
1201   fclose(fp);\r
1202   \r
1203   return true;
1204 }
1205
1206
1207 #ifdef HAVE_PNG
1208 bool
1209 ImageFile::writeImagePNG (const char* const outfile, int bitdepth, int nxcell, int nycell, double densmin, double densmax)
1210 {
1211   double max_out_level = (1 << bitdepth) - 1;
1212   int nx = m_nx;
1213   int ny = m_ny;
1214   ImageFileArray v = getArray();
1215
1216   unsigned char* rowp = new unsigned char [nx * nxcell * (bitdepth / 8)];
1217
1218   FILE *fp = fopen (outfile, "wb");\r
1219   if (! fp)
1220     return false;
1221
1222   png_structp png_ptr = png_create_write_struct (PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
1223   if (! png_ptr)
1224     return false;
1225
1226   png_infop info_ptr = png_create_info_struct (png_ptr);
1227   if (! info_ptr) {
1228     png_destroy_write_struct (&png_ptr, (png_infopp) NULL);
1229     fclose (fp);
1230     return false;
1231   }
1232
1233   if (setjmp (png_ptr->jmpbuf)) {
1234     png_destroy_write_struct (&png_ptr, &info_ptr);
1235     fclose (fp);
1236     return false;
1237   }
1238
1239   png_init_io(png_ptr, fp);
1240
1241   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);
1242
1243   png_write_info(png_ptr, info_ptr);
1244   for (int irow = ny - 1; irow >= 0; irow--) {
1245     png_bytep row_pointer = rowp;
1246
1247     for (int icol = 0; icol < nx; icol++) {
1248       int pos = icol * nxcell;
1249       double dens = (v[icol][irow] - densmin) / (densmax - densmin);
1250       dens = clamp (dens, 0., 1.);
1251       unsigned int outval = static_cast<unsigned int> (dens * max_out_level);
1252
1253       for (int p = pos; p < pos + nxcell; p++) {
1254         if (bitdepth == 8)
1255           rowp[p] = outval;
1256         else {
1257           int rowpos = p * 2;
1258           rowp[rowpos] = (outval >> 8) & 0xFF;
1259           rowp[rowpos+1] = (outval & 0xFF);
1260         }
1261       }
1262     }
1263     for (int ir = 0; ir < nycell; ir++)
1264       png_write_rows (png_ptr, &row_pointer, 1);
1265   }
1266
1267   png_write_end (png_ptr, info_ptr);
1268   png_destroy_write_struct (&png_ptr, &info_ptr);
1269   delete rowp;\r
1270
1271   fclose(fp);\r
1272   \r
1273   return true;
1274 }
1275 #endif
1276
1277 #ifdef HAVE_GD
1278 #include "gd.h"
1279 static const int N_GRAYSCALE=256;
1280
1281 bool
1282 ImageFile::writeImageGIF (const char* const outfile, int nxcell, int nycell, double densmin, double densmax)
1283 {
1284   int gs_indices[N_GRAYSCALE];
1285   int nx = m_nx;
1286   int ny = m_ny;
1287   ImageFileArray v = getArray();
1288
1289   unsigned char* rowp = new unsigned char [nx * nxcell];
1290
1291   gdImagePtr gif = gdImageCreate(nx * nxcell, ny * nycell);
1292   for (int i = 0; i < N_GRAYSCALE; i++)
1293     gs_indices[i] = gdImageColorAllocate(gif, i, i, i);
1294
1295   int lastrow = ny * nycell - 1;
1296   for (int irow = 0; irow < ny; irow++) {
1297     int rpos = irow * nycell;
1298     for (int ir = rpos; ir < rpos + nycell; ir++) {
1299       for (int icol = 0; icol < nx; icol++) {
1300         int cpos = icol * nxcell;
1301         double dens = (v[icol][irow] - densmin) / (densmax - densmin);
1302         dens = clamp(dens, 0., 1.);
1303         for (int ic = cpos; ic < cpos + nxcell; ic++) {
1304           rowp[ic] = (unsigned int) (dens * (double) (N_GRAYSCALE - 1));
1305           gdImageSetPixel(gif, ic, lastrow - ir, gs_indices[rowp[ic]]);
1306         }
1307       }
1308     }
1309   }
1310
1311   FILE *out;\r
1312   if ((out = fopen (outfile,"w")) == NULL) {
1313     sys_error(ERR_FATAL, "Error opening output file %s for writing", outfile);
1314     return false;
1315   }
1316   gdImageGif(gif,out);
1317   fclose(out);
1318   gdImageDestroy(gif);\r
1319   \r
1320   delete rowp;\r
1321   \r
1322   return true;
1323 }
1324 #endif
1325