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