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.34 2001/01/04 21:28:41 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] = vReal[scaleNX][scaleNY] +
 705           dXFrac * (vReal[scaleNX+1][scaleNY] - vReal[scaleNX][scaleNY]) +
 706           dYFrac * (vReal[scaleNX][scaleNY+1] - vReal[scaleNX][scaleNY]);
 707         if (result.isComplex()) {
 708           if (isComplex())
 709             vResultImag[ix][iy] = vImag[scaleNX][scaleNY] +
 710             dXFrac * (vImag[scaleNX+1][scaleNY] - vImag[scaleNX][scaleNY]) +
 711             dYFrac * (vImag[scaleNX][scaleNY+1] - vImag[scaleNX][scaleNY]);
 712           else
 713             vResultImag[ix][iy] = 0;
 714         }
 715       }
 716     }
 717   }
 718
 719   return true;
 720 }
 721
 722 #ifdef HAVE_FFTW
 723 bool
 724 ImageFile::fft (ImageFile& result) const
 725 {
 726   if (m_nx != result.nx() || m_ny != result.ny()) {
 727     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
 728     return false;
 729   }
 730
 731   if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
 732     if (! result.convertRealToComplex ())
 733       return false;
 734   }
 735
 736   fftw_complex* in = new fftw_complex [m_nx * m_ny];
 737
 738   ImageFileArrayConst vReal = getArray();
 739   ImageFileArrayConst vImag = getImaginaryArray();
 740
 741   unsigned int ix, iy;
 742   unsigned int iArray = 0;
 743   for (ix = 0; ix < m_nx; ix++)
 744     for (iy = 0; iy < m_ny; iy++) {
 745       in[iArray].re = vReal[ix][iy];
 746       if (isComplex())
 747         in[iArray].im = vImag[ix][iy];
 748       else
 749         in[iArray].im = 0;
 750       iArray++;
 751     }
 752
 753     fftwnd_plan plan = fftw2d_create_plan (m_nx, m_ny, FFTW_FORWARD, FFTW_IN_PLACE);
 754
 755     fftwnd_one (plan, in, NULL);
 756
 757     ImageFileArray vRealResult = result.getArray();
 758     ImageFileArray vImagResult = result.getImaginaryArray();
 759     iArray = 0;
 760     unsigned int iScale = m_nx * m_ny;
 761     for (ix = 0; ix < m_nx; ix++)
 762       for (iy = 0; iy < m_ny; iy++) {
 763         vRealResult[ix][iy] = in[iArray].re / iScale;
 764         vImagResult[ix][iy] = in[iArray].im / iScale;
 765         iArray++;
 766       }
 767
 768       fftwnd_destroy_plan (plan);
 769       delete in;
 770
 771
 772       Fourier::shuffleFourierToNaturalOrder (result);
 773
 774       return true;
 775 }
 776
 777
 778 bool
 779 ImageFile::ifft (ImageFile& result) const
 780 {
 781   if (m_nx != result.nx() || m_ny != result.ny()) {
 782     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
 783     return false;
 784   }
 785
 786   if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
 787     if (! result.convertRealToComplex ())
 788       return false;
 789   }
 790
 791   ImageFileArrayConst vReal = getArray();
 792   ImageFileArrayConst vImag = getImaginaryArray();
 793   ImageFileArray vRealResult = result.getArray();
 794   ImageFileArray vImagResult = result.getImaginaryArray();
 795   unsigned int ix, iy;
 796   for (ix = 0; ix < m_nx; ix++)
 797     for (iy = 0; iy < m_ny; iy++) {
 798       vRealResult[ix][iy] = vReal[ix][iy];
 799       if (isComplex())
 800         vImagResult[ix][iy] = vImag[ix][iy];
 801       else
 802         vImagResult[ix][iy] = 0;
 803     }
 804
 805     Fourier::shuffleNaturalToFourierOrder (result);
 806
 807     fftw_complex* in = new fftw_complex [m_nx * m_ny];
 808
 809     unsigned int iArray = 0;
 810     for (ix = 0; ix < m_nx; ix++)
 811       for (iy = 0; iy < m_ny; iy++) {
 812         in[iArray].re = vRealResult[ix][iy];
 813         in[iArray].im = vImagResult[ix][iy];
 814         iArray++;
 815       }
 816
 817       fftwnd_plan plan = fftw2d_create_plan (m_nx, m_ny, FFTW_BACKWARD, FFTW_IN_PLACE);
 818
 819       fftwnd_one (plan, in, NULL);
 820
 821       iArray = 0;
 822       for (ix = 0; ix < m_nx; ix++)
 823         for (iy = 0; iy < m_ny; iy++) {
 824           vRealResult[ix][iy] = in[iArray].re;
 825           vImagResult[ix][iy] = in[iArray].im;
 826           iArray++;
 827         }
 828
 829         fftwnd_destroy_plan (plan);
 830
 831         delete in;
 832
 833         return true;
 834 }
 835 #endif // HAVE_FFTW
 836
 837
 838
 839 bool
 840 ImageFile::fourier (ImageFile& result) const
 841 {
 842   if (m_nx != result.nx() || m_ny != result.ny()) {
 843     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
 844     return false;
 845   }
 846
 847   if (! result.isComplex())
 848     if (! result.convertRealToComplex ())
 849       return false;
 850
 851     ImageFileArrayConst vLHS = getArray();
 852     ImageFileArrayConst vLHSImag = getImaginaryArray();
 853     ImageFileArray vRealResult = result.getArray();
 854     ImageFileArray vImagResult = result.getImaginaryArray();
 855
 856     unsigned int ix, iy;
 857
 858     // alloc output matrix
 859     CTSimComplex** complexOut = new CTSimComplex* [m_nx];
 860     for (ix = 0; ix < m_nx; ix++)
 861       complexOut[ix] = new CTSimComplex [m_ny];
 862
 863     // fourier each x column
 864     CTSimComplex* pY = new CTSimComplex [m_ny];
 865     for (ix = 0; ix < m_nx; ix++) {
 866       for (iy = 0; iy < m_ny; iy++) {
 867         double dImag = 0;
 868         if (isComplex())
 869           dImag = vLHSImag[ix][iy];
 870         pY[iy] = std::complex<double>(vLHS[ix][iy], dImag);
 871       }
 872       ProcessSignal::finiteFourierTransform (pY, complexOut[ix], m_ny,  ProcessSignal::FORWARD);
 873     }
 874     delete [] pY;
 875
 876     // fourier each y row
 877     CTSimComplex* pX = new CTSimComplex [m_nx];
 878     CTSimComplex* complexOutRow = new CTSimComplex [m_nx];
 879     for (iy = 0; iy < m_ny; iy++) {
 880       for (ix = 0; ix < m_nx; ix++)
 881         pX[ix] = complexOut[ix][iy];
 882       ProcessSignal::finiteFourierTransform (pX, complexOutRow, m_nx, ProcessSignal::FORWARD);
 883       for (ix = 0; ix < m_nx; ix++)
 884         complexOut[ix][iy] = complexOutRow[ix];
 885     }
 886     delete [] pX;
 887     delete [] complexOutRow;
 888
 889     for (ix = 0; ix < m_nx; ix++)
 890       for (iy = 0; iy < m_ny; iy++) {
 891         vRealResult[ix][iy] = complexOut[ix][iy].real();
 892         vImagResult[ix][iy] = complexOut[ix][iy].imag();
 893       }
 894
 895       Fourier::shuffleFourierToNaturalOrder (result);
 896
 897       // delete complexOut matrix
 898       for (ix = 0; ix < m_nx; ix++)
 899         delete [] complexOut[ix];
 900       delete [] complexOut;
 901
 902       return true;
 903 }
 904
 905 bool
 906 ImageFile::inverseFourier (ImageFile& result) const
 907 {
 908   if (m_nx != result.nx() || m_ny != result.ny()) {
 909     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
 910     return false;
 911   }
 912
 913   if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
 914     if (! result.convertRealToComplex ())
 915       return false;
 916   }
 917
 918   ImageFileArrayConst vLHSReal = getArray();
 919   ImageFileArrayConst vLHSImag = getImaginaryArray();
 920   ImageFileArray vRealResult = result.getArray();
 921   ImageFileArray vImagResult = result.getImaginaryArray();
 922
 923   unsigned int ix, iy;
 924   // alloc 2d complex output matrix
 925   CTSimComplex** complexOut = new CTSimComplex* [m_nx];
 926   for (ix = 0; ix < m_nx; ix++)
 927     complexOut[ix] = new CTSimComplex [m_ny];
 928
 929   // put input image into result
 930   for (ix = 0; ix < m_nx; ix++)
 931     for (iy = 0; iy < m_ny; iy++) {
 932       vRealResult[ix][iy] = vLHSReal[ix][iy];
 933       if (isComplex())
 934         vImagResult[ix][iy] = vLHSImag[ix][iy];
 935       else
 936         vImagResult[ix][iy] = 0;
 937     }
 938
 939     Fourier::shuffleNaturalToFourierOrder (result);
 940
 941     // ifourier each x column
 942     CTSimComplex* pCol = new CTSimComplex [m_ny];
 943     for (ix = 0; ix < m_nx; ix++) {
 944       for (iy = 0; iy < m_ny; iy++) {
 945         pCol[iy] = std::complex<double> (vRealResult[ix][iy], vImagResult[ix][iy]);
 946       }
 947       ProcessSignal::finiteFourierTransform (pCol, complexOut[ix], m_ny,  ProcessSignal::BACKWARD);
 948     }
 949     delete [] pCol;
 950
 951     // ifourier each y row
 952     CTSimComplex* complexInRow = new CTSimComplex [m_nx];
 953     CTSimComplex* complexOutRow = new CTSimComplex [m_nx];
 954     for (iy = 0; iy < m_ny; iy++) {
 955       for (ix = 0; ix < m_nx; ix++)
 956         complexInRow[ix] = complexOut[ix][iy];
 957       ProcessSignal::finiteFourierTransform (complexInRow, complexOutRow, m_nx, ProcessSignal::BACKWARD);
 958       for (ix = 0; ix < m_nx; ix++)
 959         complexOut[ix][iy] = complexOutRow[ix];
 960     }
 961     delete [] complexInRow;
 962     delete [] complexOutRow;
 963
 964     for (ix = 0; ix < m_nx; ix++)
 965       for (iy = 0; iy < m_ny; iy++) {
 966         vRealResult[ix][iy] = complexOut[ix][iy].real();
 967         vImagResult[ix][iy] = complexOut[ix][iy].imag();
 968       }
 969
 970       // delete complexOut matrix
 971       for (ix = 0; ix < m_nx; ix++)
 972         delete [] complexOut[ix];
 973       delete [] complexOut;
 974
 975       return true;
 976 }
 977
 978
 979 bool
 980 ImageFile::magnitude (ImageFile& result) const
 981 {
 982   if (m_nx != result.nx() || m_ny != result.ny()) {
 983     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
 984     return false;
 985   }
 986
 987   ImageFileArray vReal = getArray();
 988   ImageFileArray vImag = getImaginaryArray();
 989   ImageFileArray vRealResult = result.getArray();
 990
 991   for (unsigned int ix = 0; ix < m_nx; ix++)
 992     for (unsigned int iy = 0; iy < m_ny; iy++) {
 993       if (isComplex())
 994         vRealResult[ix][iy] = ::sqrt (vReal[ix][iy] * vReal[ix][iy] + vImag[ix][iy] * vImag[ix][iy]);
 995       else
 996         vRealResult[ix][iy] = vReal[ix][iy];
 997     }
 998
 999     if (result.isComplex())
1000       result.convertComplexToReal();
1001
1002     return true;
1003 }
1004
1005 bool
1006 ImageFile::phase (ImageFile& result) const
1007 {
1008   if (m_nx != result.nx() || m_ny != result.ny()) {
1009     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
1010     return false;
1011   }
1012
1013   ImageFileArray vReal = getArray();
1014   ImageFileArray vImag = getImaginaryArray();
1015   ImageFileArray vRealResult = result.getArray();
1016
1017   for (unsigned int ix = 0; ix < m_nx; ix++)
1018     for (unsigned int iy = 0; iy < m_ny; iy++) {
1019       if (isComplex())
1020         vRealResult[ix][iy] = ::atan2 (vImag[ix][iy], vReal[ix][iy]);
1021       else
1022         vRealResult[ix][iy] = 0;
1023     }
1024
1025     if (result.isComplex())
1026       result.convertComplexToReal();
1027
1028     return true;
1029 }
1030
1031 bool
1032 ImageFile::square (ImageFile& result) const
1033 {
1034   if (m_nx != result.nx() || m_ny != result.ny()) {
1035     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
1036     return false;
1037   }
1038
1039   if (isComplex() && ! result.isComplex())
1040     result.convertRealToComplex();
1041
1042   ImageFileArrayConst vLHS = getArray();
1043   ImageFileArrayConst vLHSImag = getImaginaryArray();
1044   ImageFileArray vResult = result.getArray();
1045   ImageFileArray vResultImag = result.getImaginaryArray();
1046
1047   for (unsigned int ix = 0; ix < m_nx; ix++) {
1048     for (unsigned int iy = 0; iy < m_ny; iy++) {
1049       if (result.isComplex()) {
1050         double dImag = 0;
1051         if (isComplex())
1052           dImag = vLHSImag[ix][iy];
1053         std::complex<double> cLHS (vLHS[ix][iy], dImag);
1054         std::complex<double> cResult = cLHS * cLHS;
1055         vResult[ix][iy] = cResult.real();
1056         vResultImag[ix][iy] = cResult.imag();
1057       } else
1058         vResult[ix][iy] = vLHS[ix][iy] * vLHS[ix][iy];
1059     }
1060   }
1061
1062
1063   return true;
1064 }
1065
1066 int
1067 ImageFile::convertFormatNameToID (const char* const formatName)
1068 {
1069   int formatID = FORMAT_INVALID;
1070
1071   for (int i = 0; i < s_iFormatCount; i++)
1072     if (strcasecmp (formatName, s_aszFormatName[i]) == 0) {
1073       formatID = i;
1074       break;
1075     }
1076
1077     return (formatID);
1078 }
1079
1080 const char*
1081 ImageFile::convertFormatIDToName (int formatID)
1082 {
1083   static const char *formatName = "";
1084
1085   if (formatID >= 0 && formatID < s_iFormatCount)
1086     return (s_aszFormatName[formatID]);
1087
1088   return (formatName);
1089 }
1090
1091 const char*
1092 ImageFile::convertFormatIDToTitle (const int formatID)
1093 {
1094   static const char *formatTitle = "";
1095
1096   if (formatID >= 0 && formatID < s_iFormatCount)
1097     return (s_aszFormatTitle[formatID]);
1098
1099   return (formatTitle);
1100 }
1101
1102 bool
1103 ImageFile::exportImage (const char* const pszFormat, const char* const pszFilename, int nxcell, int nycell, double densmin, double densmax)
1104 {
1105   int iFormatID = convertFormatNameToID (pszFormat);
1106   if (iFormatID == FORMAT_INVALID) {
1107     sys_error (ERR_SEVERE, "Invalid format %s [ImageFile::exportImage]", pszFormat);
1108     return false;
1109   }
1110
1111   if (iFormatID == FORMAT_PGM)
1112     return writeImagePGM (pszFilename, nxcell, nycell, densmin, densmax);
1113   else if (iFormatID == FORMAT_PGMASCII)
1114     return writeImagePGMASCII (pszFilename, nxcell, nycell, densmin, densmax);
1115   else if (iFormatID == FORMAT_PNG)
1116     return writeImagePNG (pszFilename, 8, nxcell, nycell, densmin, densmax);
1117   else if (iFormatID == FORMAT_PNG16)
1118     return writeImagePNG (pszFilename, 16, nxcell, nycell, densmin, densmax);
1119
1120   sys_error (ERR_SEVERE, "Invalid format %s [ImageFile::exportImage]", pszFormat);
1121   return false;
1122 }
1123
1124 bool
1125 ImageFile::writeImagePGM (const char* const outfile, int nxcell, int nycell, double densmin, double densmax)
1126 {
1127   FILE *fp;
1128   int nx = m_nx;
1129   int ny = m_ny;
1130   ImageFileArray v = getArray();
1131
1132   unsigned char* rowp = new unsigned char [nx * nxcell];
1133
1134   if ((fp = fopen (outfile, "wb")) == NULL)
1135     return false;
1136
1137   fprintf(fp, "P5\n");
1138   fprintf(fp, "%d %d\n", nx, ny);
1139   fprintf(fp, "255\n");
1140
1141   for (int irow = ny - 1; irow >= 0; irow--) {
1142     for (int icol = 0; icol < nx; icol++) {
1143       int pos = icol * nxcell;
1144       double dens = (v[icol][irow] - densmin) / (densmax - densmin);
1145       dens = clamp (dens, 0., 1.);
1146       for (int p = pos; p < pos + nxcell; p++) {
1147         rowp[p] = static_cast<unsigned int> (dens * 255.);
1148       }
1149     }
1150     for (int ir = 0; ir < nycell; ir++) {
1151       for (int ic = 0; ic < nx * nxcell; ic++)
1152         fputc( rowp[ic], fp );
1153     }
1154   }
1155
1156   delete rowp;
1157   fclose(fp);
1158
1159   return true;
1160 }
1161
1162 bool
1163 ImageFile::writeImagePGMASCII (const char* const outfile, int nxcell, int nycell, double densmin, double densmax)
1164 {
1165   FILE *fp;
1166   int nx = m_nx;
1167   int ny = m_ny;
1168   ImageFileArray v = getArray();
1169
1170   unsigned char* rowp = new unsigned char [nx * nxcell];
1171
1172   if ((fp = fopen (outfile, "wb")) == NULL)
1173     return false;
1174
1175   fprintf(fp, "P2\n");
1176   fprintf(fp, "%d %d\n", nx, ny);
1177   fprintf(fp, "255\n");
1178
1179   for (int irow = ny - 1; irow >= 0; irow--) {
1180     for (int icol = 0; icol < nx; icol++) {
1181       int pos = icol * nxcell;
1182       double dens = (v[icol][irow] - densmin) / (densmax - densmin);
1183       dens = clamp (dens, 0., 1.);
1184       for (int p = pos; p < pos + nxcell; p++) {
1185         rowp[p] = static_cast<unsigned int> (dens * 255.);
1186       }
1187     }
1188     for (int ir = 0; ir < nycell; ir++) {
1189       for (int ic = 0; ic < nx * nxcell; ic++)
1190         fprintf(fp, "%d ", rowp[ic]);
1191       fprintf(fp, "\n");
1192     }
1193   }
1194
1195   delete rowp;
1196   fclose(fp);
1197
1198   return true;
1199 }
1200
1201
1202 #ifdef HAVE_PNG
1203 bool
1204 ImageFile::writeImagePNG (const char* const outfile, int bitdepth, int nxcell, int nycell, double densmin, double densmax)
1205 {
1206   double max_out_level = (1 << bitdepth) - 1;
1207   int nx = m_nx;
1208   int ny = m_ny;
1209   ImageFileArray v = getArray();
1210
1211   unsigned char* rowp = new unsigned char [nx * nxcell * (bitdepth / 8)];
1212
1213   FILE *fp = fopen (outfile, "wb");
1214   if (! fp)
1215     return false;
1216
1217   png_structp png_ptr = png_create_write_struct (PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
1218   if (! png_ptr)
1219     return false;
1220
1221   png_infop info_ptr = png_create_info_struct (png_ptr);
1222   if (! info_ptr) {
1223     png_destroy_write_struct (&png_ptr, (png_infopp) NULL);
1224     fclose (fp);
1225     return false;
1226   }
1227
1228   if (setjmp (png_ptr->jmpbuf)) {
1229     png_destroy_write_struct (&png_ptr, &info_ptr);
1230     fclose (fp);
1231     return false;
1232   }
1233
1234   png_init_io(png_ptr, fp);
1235
1236   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);
1237
1238   png_write_info(png_ptr, info_ptr);
1239   for (int irow = ny - 1; irow >= 0; irow--) {
1240     png_bytep row_pointer = rowp;
1241
1242     for (int icol = 0; icol < nx; icol++) {
1243       int pos = icol * nxcell;
1244       double dens = (v[icol][irow] - densmin) / (densmax - densmin);
1245       dens = clamp (dens, 0., 1.);
1246       unsigned int outval = static_cast<unsigned int> (dens * max_out_level);
1247
1248       for (int p = pos; p < pos + nxcell; p++) {
1249         if (bitdepth == 8)
1250           rowp[p] = outval;
1251         else {
1252           int rowpos = p * 2;
1253           rowp[rowpos+1] = (outval >> 8) & 0xFF;
1254           rowp[rowpos] = (outval & 0xFF);
1255         }
1256       }
1257     }
1258     for (int ir = 0; ir < nycell; ir++)
1259       png_write_rows (png_ptr, &row_pointer, 1);
1260   }
1261
1262   png_write_end (png_ptr, info_ptr);
1263   png_destroy_write_struct (&png_ptr, &info_ptr);
1264   delete rowp;
1265
1266   fclose(fp);
1267
1268   return true;
1269 }
1270 #endif
1271
1272 #ifdef HAVE_GD
1273 #include "gd.h"
1274 static const int N_GRAYSCALE=256;
1275
1276 bool
1277 ImageFile::writeImageGIF (const char* const outfile, int nxcell, int nycell, double densmin, double densmax)
1278 {
1279   int gs_indices[N_GRAYSCALE];
1280   int nx = m_nx;
1281   int ny = m_ny;
1282   ImageFileArray v = getArray();
1283
1284   unsigned char* rowp = new unsigned char [nx * nxcell];
1285
1286   gdImagePtr gif = gdImageCreate(nx * nxcell, ny * nycell);
1287   for (int i = 0; i < N_GRAYSCALE; i++)
1288     gs_indices[i] = gdImageColorAllocate(gif, i, i, i);
1289
1290   int lastrow = ny * nycell - 1;
1291   for (int irow = 0; irow < ny; irow++) {
1292     int rpos = irow * nycell;
1293     for (int ir = rpos; ir < rpos + nycell; ir++) {
1294       for (int icol = 0; icol < nx; icol++) {
1295         int cpos = icol * nxcell;
1296         double dens = (v[icol][irow] - densmin) / (densmax - densmin);
1297         dens = clamp(dens, 0., 1.);
1298         for (int ic = cpos; ic < cpos + nxcell; ic++) {
1299           rowp[ic] = (unsigned int) (dens * (double) (N_GRAYSCALE - 1));
1300           gdImageSetPixel(gif, ic, lastrow - ir, gs_indices[rowp[ic]]);
1301         }
1302       }
1303     }
1304   }
1305
1306   FILE *out;
1307   if ((out = fopen (outfile,"w")) == NULL) {
1308     sys_error(ERR_FATAL, "Error opening output file %s for writing", outfile);
1309     return false;
1310   }
1311   gdImageGif(gif,out);
1312   fclose(out);
1313   gdImageDestroy(gif);
1314
1315   delete rowp;
1316
1317   return true;
1318 }
1319 #endif
1320