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.29 2001/01/02 05:34:57 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 (int ix = 0; ix < m_nx; ix++)\r
 104     for (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         std::complex<double> cLHS (vLHS[ix][iy], 0);\r
 470         if (isComplex())\r
 471           cLHS.imag (vLHSImag[ix][iy]);\r
 472         std::complex<double> cRHS (vRHS[ix][iy], 0);\r
 473         if (rRHS.isComplex())\r
 474           cRHS.imag (vRHSImag[ix][iy]);\r
 475         std::complex<double> cResult = cLHS * cRHS;\r
 476         vResult[ix][iy] = cResult.real();\r
 477         vResultImag[ix][iy] = cResult.imag();\r
 478       } else\r
 479         vResult[ix][iy] = vLHS[ix][iy] * vRHS[ix][iy];\r
 480     }\r
 481   }\r
 482   \r
 483   \r
 484   return true;\r
 485 }\r
 486 \r
 487 bool\r
 488 ImageFile::divideImages (const ImageFile& rRHS, ImageFile& result) const\r
 489 {\r
 490   if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) {\r
 491     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");\r
 492     return false;\r
 493   }\r
 494   \r
 495   if (isComplex() || rRHS.isComplex() && ! result.isComplex())\r
 496     result.convertRealToComplex();\r
 497 \r
 498   ImageFileArrayConst vLHS = getArray();\r
 499   ImageFileArrayConst vLHSImag = getImaginaryArray();\r
 500   ImageFileArrayConst vRHS = rRHS.getArray();\r
 501   ImageFileArrayConst vRHSImag = rRHS.getImaginaryArray();\r
 502   ImageFileArray vResult = result.getArray();\r
 503   ImageFileArray vResultImag = result.getImaginaryArray();\r
 504   \r
 505   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
 506     for (unsigned int iy = 0; iy < m_ny; iy++) {\r
 507       if (result.isComplex()) {\r
 508         std::complex<double> cLHS (vLHS[ix][iy], 0);\r
 509         if (isComplex())\r
 510           cLHS.imag (vLHSImag[ix][iy]);\r
 511         std::complex<double> cRHS (vRHS[ix][iy], 0);\r
 512         if (rRHS.isComplex())\r
 513           cRHS.imag (vRHSImag[ix][iy]);\r
 514         std::complex<double> cResult = cLHS / cRHS;\r
 515         vResult[ix][iy] = cResult.real();\r
 516         vResultImag[ix][iy] = cResult.imag();\r
 517       } else {\r
 518         if (vRHS != 0)\r
 519           vResult[ix][iy] = vLHS[ix][iy] / vRHS[ix][iy];\r
 520         else\r
 521           vResult[ix][iy] = 0;\r
 522       }\r
 523     }\r
 524   }\r
 525   \r
 526   \r
 527   return true;\r
 528 }\r
 529 \r
 530 \r
 531 bool\r
 532 ImageFile::invertPixelValues (ImageFile& result) const\r
 533 {\r
 534   if (m_nx != result.nx() || m_ny != result.ny()) {\r
 535     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
 536     return false;\r
 537   }\r
 538   \r
 539   if (isComplex() && ! result.isComplex())\r
 540     result.convertRealToComplex();\r
 541 \r
 542   ImageFileArrayConst vLHS = getArray();\r
 543   ImageFileArray vResult = result.getArray();\r
 544   \r
 545   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
 546     ImageFileColumnConst in = vLHS[ix];\r
 547     ImageFileColumn out = vResult[ix];\r
 548     for (unsigned int iy = 0; iy < m_ny; iy++)\r
 549       *out++ = - *in++;\r
 550   }\r
 551   \r
 552   return true;\r
 553 }\r
 554 \r
 555 bool\r
 556 ImageFile::sqrt (ImageFile& result) const\r
 557 {\r
 558   if (m_nx != result.nx() || m_ny != result.ny()) {\r
 559     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
 560     return false;\r
 561   }\r
 562 \r
 563   if (isComplex() && ! result.isComplex())\r
 564     result.convertRealToComplex();\r
 565   \r
 566   bool bComplexOutput = result.isComplex();\r
 567   ImageFileArrayConst vLHS = getArray();\r
 568   if (! bComplexOutput)   // check if should convert to complex output\r
 569     for (unsigned int ix = 0; ix < m_nx; ix++)\r
 570       for (unsigned int iy = 0; iy < m_ny; iy++)\r
 571         if (! bComplexOutput && vLHS[ix][iy] < 0) {\r
 572           result.convertRealToComplex();\r
 573           bComplexOutput = true;\r
 574           break;\r
 575         }\r
 576 \r
 577   ImageFileArrayConst vLHSImag = getImaginaryArray();\r
 578   ImageFileArray vResult = result.getArray();\r
 579   ImageFileArray vResultImag = result.getImaginaryArray();\r
 580   \r
 581   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
 582     for (unsigned int iy = 0; iy < m_ny; iy++) {\r
 583       if (result.isComplex()) {\r
 584         std::complex<double> cLHS (vLHS[ix][iy], 0);\r
 585         if (isComplex())\r
 586           cLHS.imag (vLHSImag[ix][iy]);\r
 587         std::complex<double> cResult = std::sqrt(cLHS);\r
 588         vResult[ix][iy] = cResult.real();\r
 589         vResultImag[ix][iy] = cResult.imag();\r
 590       } else\r
 591         vResult[ix][iy] = ::sqrt (vLHS[ix][iy]);\r
 592     }\r
 593   }\r
 594   \r
 595   \r
 596   return true;\r
 597 }\r
 598 \r
 599 bool\r
 600 ImageFile::log (ImageFile& result) const\r
 601 {\r
 602   if (m_nx != result.nx() || m_ny != result.ny()) {\r
 603     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
 604     return false;\r
 605   }\r
 606   \r
 607   if (isComplex() && ! result.isComplex())\r
 608     result.convertRealToComplex();\r
 609 \r
 610   ImageFileArrayConst vLHS = getArray();\r
 611   ImageFileArrayConst vLHSImag = getImaginaryArray();\r
 612   ImageFileArray vResult = result.getArray();\r
 613   ImageFileArray vResultImag = result.getImaginaryArray();\r
 614   \r
 615   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
 616     for (unsigned int iy = 0; iy < m_ny; iy++) {\r
 617       if (result.isComplex()) {\r
 618         std::complex<double> cLHS (vLHS[ix][iy], 0);\r
 619         if (isComplex())\r
 620           cLHS.imag (vLHSImag[ix][iy]);\r
 621         std::complex<double> cResult = std::log (cLHS);\r
 622         vResult[ix][iy] = cResult.real();\r
 623         vResultImag[ix][iy] = cResult.imag();\r
 624       } else\r
 625         vResult[ix][iy] = ::log (vLHS[ix][iy]);\r
 626     }\r
 627   }\r
 628   \r
 629   \r
 630   return true;\r
 631 }\r
 632 \r
 633 bool\r
 634 ImageFile::exp (ImageFile& result) const\r
 635 {\r
 636   if (m_nx != result.nx() || m_ny != result.ny()) {\r
 637     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
 638     return false;\r
 639   }\r
 640     \r
 641   if (isComplex() && ! result.isComplex())\r
 642     result.convertRealToComplex();\r
 643 \r
 644   ImageFileArrayConst vLHS = getArray();\r
 645   ImageFileArrayConst vLHSImag = getImaginaryArray();\r
 646   ImageFileArray vResult = result.getArray();\r
 647   ImageFileArray vResultImag = result.getImaginaryArray();\r
 648   \r
 649   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
 650     for (unsigned int iy = 0; iy < m_ny; iy++) {\r
 651       if (result.isComplex()) {\r
 652         std::complex<double> cLHS (vLHS[ix][iy], 0);\r
 653         if (isComplex())\r
 654           cLHS.imag (vLHSImag[ix][iy]);\r
 655         std::complex<double> cResult = std::exp (cLHS);\r
 656         vResult[ix][iy] = cResult.real();\r
 657         vResultImag[ix][iy] = cResult.imag();\r
 658       } else\r
 659         vResult[ix][iy] = ::exp (vLHS[ix][iy]);\r
 660     }\r
 661   }\r
 662   \r
 663   \r
 664   return true;\r
 665 }\r
 666 \r
 667 bool\r
 668 ImageFile::scaleImage (ImageFile& result) const\r
 669 {\r
 670   unsigned int nx = m_nx;\r
 671   unsigned int ny = m_ny;\r
 672   unsigned int newNX = result.nx();\r
 673   unsigned int newNY = result.ny();\r
 674 \r
 675   double dXScale = static_cast<double>(newNX) / static_cast<double>(nx);\r
 676   double dYScale = static_cast<double>(newNY) / static_cast<double>(ny);\r
 677 \r
 678   if (isComplex() && ! result.isComplex())\r
 679     result.convertRealToComplex();\r
 680 \r
 681   ImageFileArrayConst vReal = getArray();\r
 682   ImageFileArrayConst vImag = getImaginaryArray();\r
 683   ImageFileArray vResult = result.getArray();\r
 684   ImageFileArray vResultImag = result.getImaginaryArray();\r
 685   \r
 686   for (unsigned int ix = 0; ix < newNX; ix++) {\r
 687     for (unsigned int iy = 0; iy < newNY; iy++) {\r
 688       unsigned int scaleNX = static_cast<unsigned int> (ix / dXScale);\r
 689       unsigned int scaleNY = static_cast<unsigned int> (iy / dYScale);\r
 690       vResult[ix][iy] = vReal[scaleNX][scaleNY];\r
 691       if (result.isComplex()) {\r
 692         if (isComplex())\r
 693           vResultImag[ix][iy] = vImag[scaleNX][scaleNY];\r
 694         else\r
 695         vResultImag[ix][iy] = 0;\r
 696       }\r
 697     }\r
 698   }\r
 699   \r
 700   return true;\r
 701 }\r
 702 \r
 703 #ifdef HAVE_FFTW\r
 704 bool\r
 705 ImageFile::fft (ImageFile& result) const\r
 706 {\r
 707   if (m_nx != result.nx() || m_ny != result.ny()) {\r
 708     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
 709     return false;\r
 710   }\r
 711 \r
 712   if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {\r
 713     if (! result.convertRealToComplex ())\r
 714       return false;\r
 715   }\r
 716 \r
 717   fftw_complex* in = new fftw_complex [m_nx * m_ny];\r
 718 \r
 719   ImageFileArrayConst vReal = getArray();\r
 720   ImageFileArrayConst vImag = getImaginaryArray();\r
 721 \r
 722   unsigned int ix, iy;\r
 723   unsigned int iArray = 0;\r
 724   for (ix = 0; ix < m_nx; ix++)\r
 725     for (iy = 0; iy < m_ny; iy++) {\r
 726       in[iArray].re = vReal[ix][iy];\r
 727       if (isComplex())\r
 728         in[iArray].im = vImag[ix][iy];\r
 729       else\r
 730         in[iArray].im = 0;\r
 731       iArray++;\r
 732     }\r
 733 \r
 734   fftwnd_plan plan = fftw2d_create_plan (m_nx, m_ny, FFTW_FORWARD, FFTW_IN_PLACE);\r
 735 \r
 736   fftwnd_one (plan, in, NULL);\r
 737 \r
 738   ImageFileArray vRealResult = result.getArray();\r
 739   ImageFileArray vImagResult = result.getImaginaryArray();\r
 740   iArray = 0;\r
 741   unsigned int iScale = m_nx * m_ny;\r
 742   for (ix = 0; ix < m_nx; ix++)\r
 743     for (iy = 0; iy < m_ny; iy++) {\r
 744       vRealResult[ix][iy] = in[iArray].re / iScale;\r
 745       vImagResult[ix][iy] = in[iArray].im / iScale;\r
 746       iArray++;\r
 747     }\r
 748 \r
 749   fftwnd_destroy_plan (plan);\r
 750   delete in;\r
 751 \r
 752 \r
 753   Fourier::shuffleFourierToNaturalOrder (result);\r
 754 \r
 755   return true;\r
 756 }\r
 757 \r
 758 \r
 759 bool\r
 760 ImageFile::ifft (ImageFile& result) const\r
 761 {\r
 762   if (m_nx != result.nx() || m_ny != result.ny()) {\r
 763     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
 764     return false;\r
 765   }\r
 766 \r
 767   if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {\r
 768     if (! result.convertRealToComplex ())\r
 769       return false;\r
 770   }\r
 771 \r
 772   ImageFileArrayConst vReal = getArray();\r
 773   ImageFileArrayConst vImag = getImaginaryArray();\r
 774   ImageFileArray vRealResult = result.getArray();\r
 775   ImageFileArray vImagResult = result.getImaginaryArray();\r
 776   unsigned int ix, iy;\r
 777   for (ix = 0; ix < m_nx; ix++)\r
 778     for (iy = 0; iy < m_ny; iy++) {\r
 779       vRealResult[ix][iy] = vReal[ix][iy];\r
 780       if (isComplex()) \r
 781         vImagResult[ix][iy] = vImag[ix][iy];\r
 782       else\r
 783         vImagResult[ix][iy] = 0;\r
 784     }\r
 785 \r
 786   Fourier::shuffleNaturalToFourierOrder (result);\r
 787 \r
 788   fftw_complex* in = new fftw_complex [m_nx * m_ny];\r
 789 \r
 790   unsigned int iArray = 0;\r
 791   for (ix = 0; ix < m_nx; ix++)\r
 792     for (iy = 0; iy < m_ny; iy++) {\r
 793       in[iArray].re = vRealResult[ix][iy];\r
 794       in[iArray].im = vImagResult[ix][iy];\r
 795       iArray++;\r
 796     }\r
 797 \r
 798   fftwnd_plan plan = fftw2d_create_plan (m_nx, m_ny, FFTW_BACKWARD, FFTW_IN_PLACE);\r
 799 \r
 800   fftwnd_one (plan, in, NULL);\r
 801 \r
 802   iArray = 0;\r
 803   for (ix = 0; ix < m_nx; ix++)\r
 804     for (iy = 0; iy < m_ny; iy++) {\r
 805       vRealResult[ix][iy] = in[iArray].re;\r
 806       vImagResult[ix][iy] = in[iArray].im;\r
 807       iArray++;\r
 808     }\r
 809 \r
 810   fftwnd_destroy_plan (plan);\r
 811 \r
 812   delete in;\r
 813 \r
 814   return true;\r
 815 }\r
 816 #endif // HAVE_FFTW\r
 817 \r
 818 \r
 819   \r
 820 bool\r
 821 ImageFile::fourier (ImageFile& result) const\r
 822 {\r
 823   if (m_nx != result.nx() || m_ny != result.ny()) {\r
 824     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
 825     return false;\r
 826   }\r
 827   \r
 828   if (! result.isComplex())\r
 829     if (! result.convertRealToComplex ())\r
 830       return false;\r
 831   \r
 832   ImageFileArrayConst vLHS = getArray();\r
 833   ImageFileArrayConst vLHSImag = getImaginaryArray();\r
 834   ImageFileArray vRealResult = result.getArray();\r
 835   ImageFileArray vImagResult = result.getImaginaryArray();\r
 836   \r
 837   unsigned int ix, iy;\r
 838 \r
 839   // alloc output matrix\r
 840   CTSimComplex** complexOut = new CTSimComplex* [m_nx];\r
 841   for (ix = 0; ix < m_nx; ix++)\r
 842     complexOut[ix] = new CTSimComplex [m_ny];\r
 843   \r
 844   // fourier each x column\r
 845 #if 1\r
 846   CTSimComplex* pY = new CTSimComplex [m_ny];\r
 847   for (ix = 0; ix < m_nx; ix++) {\r
 848     for (iy = 0; iy < m_ny; iy++) {\r
 849       pY[iy].real (vLHS[ix][iy]);\r
 850       if (isComplex())\r
 851         pY[iy].imag (vLHSImag[ix][iy]);\r
 852       else\r
 853         pY[iy].imag (0);\r
 854     } \r
 855     ProcessSignal::finiteFourierTransform (pY, complexOut[ix], m_ny,  ProcessSignal::FORWARD);\r
 856   }\r
 857 #else\r
 858   double* pY = new double [m_ny];\r
 859   for (ix = 0; ix < m_nx; ix++) {\r
 860     for (iy = 0; iy < m_ny; iy++) {\r
 861       pY[iy] = vLHS[ix][iy];\r
 862     } \r
 863     ProcessSignal::finiteFourierTransform (pY, complexOut[ix], m_ny,  ProcessSignal::FORWARD);\r
 864   }\r
 865 #endif\r
 866   delete [] pY;\r
 867   \r
 868   // fourier each y row\r
 869   CTSimComplex* pX = new CTSimComplex [m_nx];\r
 870   CTSimComplex* complexOutRow = new CTSimComplex [m_nx];\r
 871   for (iy = 0; iy < m_ny; iy++) {\r
 872     for (ix = 0; ix < m_nx; ix++)\r
 873       pX[ix] = complexOut[ix][iy];\r
 874     ProcessSignal::finiteFourierTransform (pX, complexOutRow, m_nx, ProcessSignal::FORWARD);\r
 875     for (ix = 0; ix < m_nx; ix++)\r
 876       complexOut[ix][iy] = complexOutRow[ix];\r
 877   }\r
 878   delete [] pX;\r
 879   delete [] complexOutRow;\r
 880 \r
 881   for (ix = 0; ix < m_nx; ix++)\r
 882     for (iy = 0; iy < m_ny; iy++) {\r
 883       vRealResult[ix][iy] = complexOut[ix][iy].real();\r
 884       vImagResult[ix][iy] = complexOut[ix][iy].imag();\r
 885     }\r
 886     \r
 887   Fourier::shuffleFourierToNaturalOrder (result);\r
 888 \r
 889   // delete complexOut matrix\r
 890   for (ix = 0; ix < m_nx; ix++)\r
 891     delete [] complexOut[ix];\r
 892   delete [] complexOut;\r
 893     \r
 894   return true;\r
 895 }\r
 896 \r
 897 bool\r
 898 ImageFile::inverseFourier (ImageFile& result) const\r
 899 {\r
 900   if (m_nx != result.nx() || m_ny != result.ny()) {\r
 901     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
 902     return false;\r
 903   }\r
 904   \r
 905   if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {\r
 906     if (! result.convertRealToComplex ())\r
 907       return false;\r
 908   }\r
 909   \r
 910   ImageFileArrayConst vLHSReal = getArray();\r
 911   ImageFileArrayConst vLHSImag = getImaginaryArray();\r
 912   ImageFileArray vRealResult = result.getArray();\r
 913   ImageFileArray vImagResult = result.getImaginaryArray();\r
 914   \r
 915   unsigned int ix, iy;\r
 916   // alloc 2d complex output matrix\r
 917   CTSimComplex** complexOut = new CTSimComplex* [m_nx];\r
 918   for (ix = 0; ix < m_nx; ix++)\r
 919     complexOut[ix] = new CTSimComplex [m_ny];\r
 920 \r
 921   // put input image into result\r
 922   for (ix = 0; ix < m_nx; ix++)\r
 923     for (iy = 0; iy < m_ny; iy++) {\r
 924       vRealResult[ix][iy] = vLHSReal[ix][iy];\r
 925       if (isComplex())\r
 926         vImagResult[ix][iy] = vLHSImag[ix][iy];\r
 927       else\r
 928         vImagResult[ix][iy] = 0;\r
 929     }\r
 930 \r
 931   Fourier::shuffleNaturalToFourierOrder (result);\r
 932 \r
 933   // ifourier each x column\r
 934   CTSimComplex* pCol = new CTSimComplex [m_ny];\r
 935   for (ix = 0; ix < m_nx; ix++) {\r
 936     for (iy = 0; iy < m_ny; iy++) {\r
 937       pCol[iy].real (vRealResult[ix][iy]);\r
 938       pCol[iy].imag (vImagResult[ix][iy]);\r
 939     }\r
 940     ProcessSignal::finiteFourierTransform (pCol, complexOut[ix], m_ny,  ProcessSignal::BACKWARD);\r
 941   }\r
 942   delete [] pCol;\r
 943 \r
 944   // ifourier each y row\r
 945   CTSimComplex* complexInRow = new CTSimComplex [m_nx];\r
 946   CTSimComplex* complexOutRow = new CTSimComplex [m_nx];\r
 947   for (iy = 0; iy < m_ny; iy++) {\r
 948     for (ix = 0; ix < m_nx; ix++)\r
 949       complexInRow[ix] = complexOut[ix][iy];\r
 950     ProcessSignal::finiteFourierTransform (complexInRow, complexOutRow, m_nx, ProcessSignal::BACKWARD);\r
 951     for (ix = 0; ix < m_nx; ix++)\r
 952       complexOut[ix][iy] = complexOutRow[ix];\r
 953   }\r
 954   delete [] complexInRow;\r
 955   delete [] complexOutRow;\r
 956 \r
 957   for (ix = 0; ix < m_nx; ix++)\r
 958     for (iy = 0; iy < m_ny; iy++) {\r
 959       vRealResult[ix][iy] = complexOut[ix][iy].real();\r
 960       vImagResult[ix][iy] = complexOut[ix][iy].imag();\r
 961     }\r
 962     \r
 963   // delete complexOut matrix\r
 964   for (ix = 0; ix < m_nx; ix++)\r
 965     delete [] complexOut[ix];\r
 966   delete [] complexOut;\r
 967     \r
 968   return true;\r
 969 }\r
 970 \r
 971 \r
 972 bool\r
 973 ImageFile::magnitude (ImageFile& result) const\r
 974 {\r
 975   if (m_nx != result.nx() || m_ny != result.ny()) {\r
 976     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
 977     return false;\r
 978   }\r
 979   \r
 980   ImageFileArray vReal = getArray();\r
 981   ImageFileArray vImag = getImaginaryArray();\r
 982   ImageFileArray vRealResult = result.getArray();\r
 983 \r
 984   for (unsigned int ix = 0; ix < m_nx; ix++)\r
 985     for (unsigned int iy = 0; iy < m_ny; iy++) {\r
 986       if (isComplex()) \r
 987         vRealResult[ix][iy] = ::sqrt (vReal[ix][iy] * vReal[ix][iy] + vImag[ix][iy] * vImag[ix][iy]);\r
 988       else\r
 989         vRealResult[ix][iy] = vReal[ix][iy];\r
 990     }\r
 991 \r
 992   if (result.isComplex())\r
 993     result.convertComplexToReal();\r
 994   \r
 995     return true;\r
 996 }\r
 997 \r
 998 bool\r
 999 ImageFile::phase (ImageFile& result) const\r
1000 {\r
1001   if (m_nx != result.nx() || m_ny != result.ny()) {\r
1002     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
1003     return false;\r
1004   }\r
1005   \r
1006   ImageFileArray vReal = getArray();\r
1007   ImageFileArray vImag = getImaginaryArray();\r
1008   ImageFileArray vRealResult = result.getArray();\r
1009 \r
1010   for (unsigned int ix = 0; ix < m_nx; ix++)\r
1011     for (unsigned int iy = 0; iy < m_ny; iy++) {\r
1012       if (isComplex())\r
1013         vRealResult[ix][iy] = ::atan2 (vImag[ix][iy], vReal[ix][iy]);\r
1014       else\r
1015         vRealResult[ix][iy] = 0;\r
1016     }\r
1017     \r
1018   if (result.isComplex())\r
1019     result.convertComplexToReal();\r
1020   \r
1021   return true;\r
1022 }\r
1023 \r
1024 bool\r
1025 ImageFile::square (ImageFile& result) const\r
1026 {\r
1027   if (m_nx != result.nx() || m_ny != result.ny()) {\r
1028     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
1029     return false;\r
1030   }\r
1031   \r
1032   if (isComplex() && ! result.isComplex())\r
1033     result.convertRealToComplex();\r
1034 \r
1035   ImageFileArrayConst vLHS = getArray();\r
1036   ImageFileArrayConst vLHSImag = getImaginaryArray();\r
1037   ImageFileArray vResult = result.getArray();\r
1038   ImageFileArray vResultImag = result.getImaginaryArray();\r
1039   \r
1040   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
1041     for (unsigned int iy = 0; iy < m_ny; iy++) {\r
1042       if (result.isComplex()) {\r
1043         std::complex<double> cLHS (vLHS[ix][iy], 0);\r
1044         if (isComplex())\r
1045           cLHS.imag (vLHSImag[ix][iy]);\r
1046         std::complex<double> cResult = cLHS * cLHS;\r
1047         vResult[ix][iy] = cResult.real();\r
1048         vResultImag[ix][iy] = cResult.imag();\r
1049       } else\r
1050         vResult[ix][iy] = vLHS[ix][iy] * vLHS[ix][iy];\r
1051     }\r
1052   }\r
1053   \r
1054   \r
1055   return true;\r
1056 }\r
1057 \r
1058 int\r
1059 ImageFile::convertFormatNameToID (const char* const formatName)\r
1060 {\r
1061   int formatID = FORMAT_INVALID;\r
1062 \r
1063   for (int i = 0; i < s_iFormatCount; i++)\r
1064       if (strcasecmp (formatName, s_aszFormatName[i]) == 0) {\r
1065           formatID = i;\r
1066           break;\r
1067       }\r
1068 \r
1069   return (formatID);\r
1070 }\r
1071 \r
1072 const char*\r
1073 ImageFile::convertFormatIDToName (int formatID)\r
1074 {\r
1075   static const char *formatName = "";\r
1076 \r
1077   if (formatID >= 0 && formatID < s_iFormatCount)\r
1078       return (s_aszFormatName[formatID]);\r
1079 \r
1080   return (formatName);\r
1081 }\r
1082 \r
1083 const char*\r
1084 ImageFile::convertFormatIDToTitle (const int formatID)\r
1085 {\r
1086   static const char *formatTitle = "";\r
1087 \r
1088   if (formatID >= 0 && formatID < s_iFormatCount)\r
1089       return (s_aszFormatTitle[formatID]);\r
1090 \r
1091   return (formatTitle);\r
1092 }\r
1093 \r
1094 bool\r
1095 ImageFile::exportImage (const char* const pszFormat, const char* const pszFilename, int nxcell, int nycell, double densmin, double densmax)\r
1096 {\r
1097   int iFormatID = convertFormatNameToID (pszFormat);\r
1098   if (iFormatID == FORMAT_INVALID) {\r
1099     sys_error (ERR_SEVERE, "Invalid format %s [ImageFile::exportImage]", pszFormat);\r
1100     return false;\r
1101   }\r
1102 \r
1103   if (iFormatID == FORMAT_PGM)\r
1104     return writeImagePGM (pszFilename, nxcell, nycell, densmin, densmax);\r
1105   else if (iFormatID == FORMAT_PGMASCII)\r
1106     return writeImagePGMASCII (pszFilename, nxcell, nycell, densmin, densmax);\r
1107   else if (iFormatID == FORMAT_PNG)\r
1108     return writeImagePNG (pszFilename, 8, nxcell, nycell, densmin, densmax);\r
1109   else if (iFormatID == FORMAT_PNG16)\r
1110     return writeImagePNG (pszFilename, 16, nxcell, nycell, densmin, densmax);\r
1111 \r
1112   sys_error (ERR_SEVERE, "Invalid format %s [ImageFile::exportImage]", pszFormat);\r
1113   return false;\r
1114 }\r
1115 \r
1116 bool
1117 ImageFile::writeImagePGM (const char* const outfile, int nxcell, int nycell, double densmin, double densmax)
1118 {
1119   FILE *fp;
1120   int nx = m_nx;
1121   int ny = m_ny;
1122   ImageFileArray v = getArray();
1123
1124   unsigned char* rowp = new unsigned char [nx * nxcell];
1125
1126   if ((fp = fopen (outfile, "wb")) == NULL)
1127     return false;
1128
1129   fprintf(fp, "P5\n");
1130   fprintf(fp, "%d %d\n", nx, ny);
1131   fprintf(fp, "255\n");
1132
1133   for (int irow = ny - 1; irow >= 0; irow--) {
1134     for (int icol = 0; icol < nx; icol++) {
1135       int pos = icol * nxcell;
1136       double dens = (v[icol][irow] - densmin) / (densmax - densmin);
1137       dens = clamp (dens, 0., 1.);
1138       for (int p = pos; p < pos + nxcell; p++) {
1139         rowp[p] = static_cast<unsigned int> (dens * 255.);
1140       }
1141     }
1142     for (int ir = 0; ir < nycell; ir++) {
1143       for (int ic = 0; ic < nx * nxcell; ic++)
1144         fputc( rowp[ic], fp );
1145     }
1146   }
1147   \r
1148   delete rowp;
1149   fclose(fp);\r
1150 \r
1151   return true;
1152 }
1153
1154 bool
1155 ImageFile::writeImagePGMASCII (const char* const outfile, int nxcell, int nycell, double densmin, double densmax)
1156 {
1157   FILE *fp;
1158   int nx = m_nx;
1159   int ny = m_ny;
1160   ImageFileArray v = getArray();
1161
1162   unsigned char* rowp = new unsigned char [nx * nxcell];
1163
1164   if ((fp = fopen (outfile, "wb")) == NULL)
1165     return false;
1166
1167   fprintf(fp, "P2\n");
1168   fprintf(fp, "%d %d\n", nx, ny);
1169   fprintf(fp, "255\n");
1170
1171   for (int irow = ny - 1; irow >= 0; irow--) {
1172     for (int icol = 0; icol < nx; icol++) {
1173       int pos = icol * nxcell;
1174       double dens = (v[icol][irow] - densmin) / (densmax - densmin);
1175       dens = clamp (dens, 0., 1.);
1176       for (int p = pos; p < pos + nxcell; p++) {
1177         rowp[p] = static_cast<unsigned int> (dens * 255.);
1178       }
1179     }
1180     for (int ir = 0; ir < nycell; ir++) {
1181       for (int ic = 0; ic < nx * nxcell; ic++)
1182         fprintf(fp, "%d ", rowp[ic]);
1183       fprintf(fp, "\n");
1184     }
1185   }
1186   \r
1187   delete rowp;
1188   fclose(fp);\r
1189 \r
1190   return true;
1191 }
1192
1193
1194 #ifdef HAVE_PNG
1195 bool
1196 ImageFile::writeImagePNG (const char* const outfile, int bitdepth, int nxcell, int nycell, double densmin, double densmax)
1197 {
1198   double max_out_level = (1 << bitdepth) - 1;
1199   int nx = m_nx;
1200   int ny = m_ny;
1201   ImageFileArray v = getArray();
1202
1203   unsigned char* rowp = new unsigned char [nx * nxcell * (bitdepth / 8)];
1204
1205   FILE *fp = fopen (outfile, "wb");\r
1206   if (! fp)
1207     return false;
1208
1209   png_structp png_ptr = png_create_write_struct (PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
1210   if (! png_ptr)
1211     return false;
1212
1213   png_infop info_ptr = png_create_info_struct (png_ptr);
1214   if (! info_ptr) {
1215     png_destroy_write_struct (&png_ptr, (png_infopp) NULL);
1216     fclose (fp);
1217     return false;
1218   }
1219
1220   if (setjmp (png_ptr->jmpbuf)) {
1221     png_destroy_write_struct (&png_ptr, &info_ptr);
1222     fclose (fp);
1223     return false;
1224   }
1225
1226   png_init_io(png_ptr, fp);
1227
1228   png_set_IHDR (png_ptr, info_ptr, nx * nxcell, ny * nycell, bitdepth, PNG_COLOR_TYPE_GRAY, PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_DEFAULT);
1229
1230   png_write_info(png_ptr, info_ptr);
1231   for (int irow = ny - 1; irow >= 0; irow--) {
1232     png_bytep row_pointer = rowp;
1233
1234     for (int icol = 0; icol < nx; icol++) {
1235       int pos = icol * nxcell;
1236       double dens = (v[icol][irow] - densmin) / (densmax - densmin);
1237       dens = clamp (dens, 0., 1.);
1238       unsigned int outval = static_cast<unsigned int> (dens * max_out_level);
1239
1240       for (int p = pos; p < pos + nxcell; p++) {
1241         if (bitdepth == 8)
1242           rowp[p] = outval;
1243         else {
1244           int rowpos = p * 2;
1245           rowp[rowpos] = (outval >> 8) & 0xFF;
1246           rowp[rowpos+1] = (outval & 0xFF);
1247         }
1248       }
1249     }
1250     for (int ir = 0; ir < nycell; ir++)
1251       png_write_rows (png_ptr, &row_pointer, 1);
1252   }
1253
1254   png_write_end (png_ptr, info_ptr);
1255   png_destroy_write_struct (&png_ptr, &info_ptr);
1256   delete rowp;\r
1257
1258   fclose(fp);\r
1259 \r
1260   return true;
1261 }
1262 #endif
1263
1264 #ifdef HAVE_GD
1265 #include "gd.h"
1266 static const int N_GRAYSCALE=256;
1267
1268 bool
1269 ImageFile::writeImageGIF (const char* const outfile, int nxcell, int nycell, double densmin, double densmax)
1270 {
1271   int gs_indices[N_GRAYSCALE];
1272   int nx = m_nx;
1273   int ny = m_ny;
1274   ImageFileArray v = getArray();
1275
1276   unsigned char* rowp = new unsigned char [nx * nxcell];
1277
1278   gdImagePtr gif = gdImageCreate(nx * nxcell, ny * nycell);
1279   for (int i = 0; i < N_GRAYSCALE; i++)
1280     gs_indices[i] = gdImageColorAllocate(gif, i, i, i);
1281
1282   int lastrow = ny * nycell - 1;
1283   for (int irow = 0; irow < ny; irow++) {
1284     int rpos = irow * nycell;
1285     for (int ir = rpos; ir < rpos + nycell; ir++) {
1286       for (int icol = 0; icol < nx; icol++) {
1287         int cpos = icol * nxcell;
1288         double dens = (v[icol][irow] - densmin) / (densmax - densmin);
1289         dens = clamp(dens, 0., 1.);
1290         for (int ic = cpos; ic < cpos + nxcell; ic++) {
1291           rowp[ic] = (unsigned int) (dens * (double) (N_GRAYSCALE - 1));
1292           gdImageSetPixel(gif, ic, lastrow - ir, gs_indices[rowp[ic]]);
1293         }
1294       }
1295     }
1296   }
1297
1298   FILE *out;\r
1299   if ((out = fopen (outfile,"w")) == NULL) {
1300     sys_error(ERR_FATAL, "Error opening output file %s for writing", outfile);
1301     return false;
1302   }
1303   gdImageGif(gif,out);
1304   fclose(out);
1305   gdImageDestroy(gif);\r
1306 \r
1307   delete rowp;\r
1308 \r
1309   return true;
1310 }
1311 #endif
1312