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