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.27 2000/12/29 20:15:37 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 {
  60         int hx = (m_nx - 1) / 2;
  61         int hy = (m_ny - 1) / 2;
  62         ImageFileArray v = getArray();
  63         SignalFilter filter (filterName, domainName, bw, filt_param);
  64
  65         for (int i = -hx; i <= hx; i++) {
  66                 for (int j = -hy; j <= hy; j++) {
  67       double r = ::sqrt (i * i + j * j);
  68
  69                         v[i+hx][j+hy] = filter.response (r);
  70                 }
  71         }
  72 }
  73
  74 int
  75 ImageFile::display (void) const
  76 {
  77     double pmin, pmax;
  78
  79     getMinMax (pmin, pmax);
  80
  81     return (displayScaling (1, pmin, pmax));
  82 }
  83
  84 int
  85 ImageFile::displayScaling (const int scale, const ImageFileValue pmin, const ImageFileValue pmax) const
  86 {
  87     int nx = m_nx;
  88     int ny = m_ny;
  89     ImageFileArrayConst v = getArray();
  90     if (v == NULL || nx == 0 || ny == 0)
  91                 return 0;
  92
  93 #if HAVE_G2_H
  94     int* pPens = new int [nx * ny * scale * scale ];
  95
  96     double view_scale = 255 / (pmax - pmin);
  97     int id_X11 = g2_open_X11 (nx * scale, ny * scale);
  98     int grayscale[256];
  99     for (int i = 0; i < 256; i++) {
 100                 double cval = i / 255.;
 101                 grayscale[i] = g2_ink (id_X11, cval, cval, cval);
 102     }
 103
 104     for (int iy = ny - 1; iy >= 0; iy--) {
 105                 int iRowPos = ((ny - 1 - iy) * scale) * (nx * scale);
 106                 for (int ix = 0; ix < nx; ix++) {
 107                         int cval = static_cast<int>((v[ix][iy] - pmin) * view_scale);
 108                         if (cval < 0)
 109                                 cval = 0;
 110                         else if (cval > 255)
 111                                 cval = 255;
 112                         for (int sy = 0; sy < scale; sy++)
 113                                 for (int sx = 0; sx < scale; sx++)
 114                                         pPens[iRowPos+(sy * nx * scale)+(sx + (ix * scale))] = grayscale[cval];
 115                 }
 116     }
 117
 118     g2_image (id_X11, 0., 0., nx * scale, ny * scale, pPens);
 119
 120     delete pPens;
 121     return (id_X11);
 122 #else
 123     return 0;
 124 #endif
 125 }
 126
 127
 128
 129 // ImageFile::comparativeStatistics    Calculate comparative stats
 130 //
 131 // OUTPUT
 132 //   d   Normalized root mean squared distance measure
 133 //   r   Normalized mean absolute distance measure
 134 //   e   Worst case distance measure
 135 //
 136 // REFERENCES
 137 //  G.T. Herman, Image Reconstruction From Projections, 1980
 138
 139 bool
 140 ImageFile::comparativeStatistics (const ImageFile& imComp, double& d, double& r, double& e) const
 141 {
 142     if (imComp.nx() != m_nx && imComp.ny() != m_ny) {
 143                 sys_error (ERR_WARNING, "Image sizes differ [ImageFile::comparativeStatistics]");
 144                 return false;
 145     }
 146     ImageFileArrayConst v = getArray();
 147     if (v == NULL || m_nx == 0 || m_ny == 0)
 148                 return false;
 149
 150     ImageFileArrayConst vComp = imComp.getArray();
 151
 152     double myMean = 0.;
 153     for (unsigned int ix = 0; ix < m_nx; ix++) {
 154                 for (unsigned int iy = 0; iy < m_ny; iy++) {
 155                         myMean += v[ix][iy];
 156                 }
 157     }
 158     myMean /= (m_nx * m_ny);
 159
 160     double sqErrorSum = 0.;
 161     double absErrorSum = 0.;
 162     double sqDiffFromMeanSum = 0.;
 163     double absValueSum = 0.;
 164     for (unsigned int ix2 = 0; ix2 < m_nx; ix2++) {
 165                 for (unsigned int iy = 0; iy < m_ny; iy++) {
 166                         double diff = v[ix2][iy] - vComp[ix2][iy];
 167                         sqErrorSum += diff * diff;
 168                         absErrorSum += fabs(diff);
 169                         double diffFromMean = v[ix2][iy] - myMean;
 170                         sqDiffFromMeanSum += diffFromMean * diffFromMean;
 171                         absValueSum += fabs(v[ix2][iy]);
 172                 }
 173     }
 174
 175     d = ::sqrt (sqErrorSum / sqDiffFromMeanSum);
 176     r = absErrorSum / absValueSum;
 177
 178     int hx = m_nx / 2;
 179     int hy = m_ny / 2;
 180     double eMax = -1;
 181     for (int ix3 = 0; ix3 < hx; ix3++) {
 182                 for (int iy = 0; iy < hy; iy++) {
 183                         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]);
 184                         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]);
 185                         double error = fabs (avgPixel - avgPixelComp);
 186                         if (error > eMax)
 187                                 eMax = error;
 188                 }
 189     }
 190
 191     e = eMax;
 192
 193     return true;
 194 }
 195
 196
 197 bool
 198 ImageFile::printComparativeStatistics (const ImageFile& imComp, std::ostream& os) const
 199 {
 200         double d, r, e;
 201
 202         if (comparativeStatistics (imComp, d, r, e)) {
 203                 os << "  Normalized root mean squared distance (d): " << d << std::endl;
 204                 os << "      Normalized mean absolute distance (r): " << r << std::endl;
 205                 os << "Worst case distance (2x2 pixel average) (e): " << e << std::endl;
 206                 return true;
 207         }
 208         return false;
 209 }
 210
 211
 212 void
 213 ImageFile::printStatistics (std::ostream& os) const
 214 {
 215     double min, max, mean, mode, median, stddev;
 216
 217     statistics (min, max, mean, mode, median, stddev);\r
 218     if (dataType() == Array2dFile::DATA_TYPE_COMPLEX)\r
 219       os << "Real Component Statistics" << std::endl;\r
 220
 221     os << "   min: " << min << std::endl;
 222     os << "   max: " << max << std::endl;
 223     os << "  mean: " << mean << std::endl;
 224     os << "  mode: " << mode << std::endl;
 225     os << "median: " << median << std::endl;
 226     os << "stddev: " << stddev << std::endl;\r
 227 \r
 228     if (dataType() == Array2dFile::DATA_TYPE_COMPLEX) {\r
 229       statistics (getImaginaryArray(), min, max, mean, mode, median, stddev);\r
 230       os << std::endl << "Imaginary Component Statistics" << std::endl;
 231     os << "   min: " << min << std::endl;\r
 232     os << "   max: " << max << std::endl;\r
 233     os << "  mean: " << mean << std::endl;\r
 234     os << "  mode: " << mode << std::endl;\r
 235     os << "median: " << median << std::endl;\r
 236     os << "stddev: " << stddev << std::endl;\r
 237     }\r
 238 }
 239
 240
 241 void
 242 ImageFile::statistics (double& min, double& max, double& mean, double& mode, double& median, double& stddev) const
 243 {\r
 244   ImageFileArrayConst v = getArray();\r
 245   statistics (v, min, max, mean, mode, median, stddev);\r
 246 }
 247
 248
 249 void\r
 250 ImageFile::statistics (ImageFileArrayConst& v, double& min, double& max, double& mean, double& mode, double& median, double& stddev) const\r
 251 {\r
 252     int nx = m_nx;\r
 253     int ny = m_ny;\r
 254     \r
 255     if (v == NULL || nx == 0 || ny == 0)\r
 256                 return;\r
 257 \r
 258         std::vector<double> vecImage;\r
 259         int iVec = 0;\r
 260         vecImage.resize (nx * ny);\r
 261     for (int ix = 0; ix < nx; ix++) {\r
 262                 for (int iy = 0; iy < ny; iy++)\r
 263                         vecImage[iVec++] = v[ix][iy];\r
 264         }\r
 265 \r
 266         vectorNumericStatistics (vecImage, nx * ny, min, max, mean, mode, median, stddev);\r
 267 }\r
 268 \r
 269 void
 270 ImageFile::getMinMax (double& min, double& max) const
 271 {
 272     int nx = m_nx;
 273     int ny = m_ny;
 274     ImageFileArrayConst v = getArray();
 275
 276     if (v == NULL || nx == 0 || ny == 0)
 277                 return;
 278
 279     min = v[0][0];
 280     max = v[0][0];
 281     for (int ix = 0; ix < nx; ix++) {
 282                 for (int iy = 0; iy < ny; iy++) {
 283                         if (v[ix][iy] > max)
 284                                 max = v[ix][iy];
 285                         if (v[ix][iy] < min)
 286                                 min = v[ix][iy];
 287                 }
 288     }
 289 }
 290 \r
 291 bool\r
 292 ImageFile::convertRealToComplex ()\r
 293 {\r
 294   if (dataType() != Array2dFile::DATA_TYPE_REAL)\r
 295     return false;\r
 296 \r
 297   if (! reallocRealToComplex())\r
 298     return false;\r
 299 \r
 300   ImageFileArray vImag = getImaginaryArray();\r
 301   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
 302     ImageFileColumn vCol = vImag[ix];\r
 303     for (unsigned int iy = 0; iy < m_ny; iy++)\r
 304       *vCol++ = 0;\r
 305   }\r
 306 \r
 307   return true;\r
 308 }\r
 309 \r
 310 bool\r
 311 ImageFile::convertComplexToReal ()\r
 312 {\r
 313   if (dataType() != Array2dFile::DATA_TYPE_COMPLEX)\r
 314     return false;\r
 315 \r
 316   ImageFileArray vReal = getArray();\r
 317   ImageFileArray vImag = getImaginaryArray();\r
 318   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
 319     ImageFileColumn vRealCol = vReal[ix];\r
 320     ImageFileColumn vImagCol = vImag[ix];\r
 321     for (unsigned int iy = 0; iy < m_ny; iy++) {\r
 322       std::complex<double> c (*vRealCol, *vImagCol);\r
 323       *vRealCol++ = std::abs (c);\r
 324       vImagCol++;\r
 325     }\r
 326   }\r
 327 \r
 328   return reallocComplexToReal();\r
 329 }\r
 330 \r
 331 bool\r
 332 ImageFile::subtractImages (const ImageFile& rRHS, ImageFile& result) const\r
 333 {\r
 334   if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) {\r
 335     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");\r
 336     return false;\r
 337   }\r
 338 \r
 339   ImageFileArrayConst vLHS = getArray();\r
 340   ImageFileArrayConst vRHS = rRHS.getArray();\r
 341   ImageFileArray vResult = result.getArray();\r
 342 \r
 343   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
 344     ImageFileColumnConst in1 = vLHS[ix];\r
 345     ImageFileColumnConst in2 = vRHS[ix];\r
 346     ImageFileColumn out = vResult[ix];\r
 347     for (unsigned int iy = 0; iy < m_ny; iy++)\r
 348         *out++ = *in1++ - *in2++;\r
 349   }\r
 350 \r
 351     return true;\r
 352 }\r
 353
 354 bool\r
 355 ImageFile::addImages (const ImageFile& rRHS, ImageFile& result) const\r
 356 {\r
 357   if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) {\r
 358     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");\r
 359     return false;\r
 360   }\r
 361 \r
 362   ImageFileArrayConst vLHS = getArray();\r
 363   ImageFileArrayConst vRHS = rRHS.getArray();\r
 364   ImageFileArray vResult = result.getArray();\r
 365 \r
 366   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
 367     ImageFileColumnConst in1 = vLHS[ix];\r
 368     ImageFileColumnConst in2 = vRHS[ix];\r
 369     ImageFileColumn out = vResult[ix];\r
 370     for (unsigned int iy = 0; iy < m_ny; iy++)\r
 371         *out++ = *in1++ + *in2++;\r
 372   }\r
 373 \r
 374     return true;\r
 375 }\r
 376 \r
 377 bool\r
 378 ImageFile::multiplyImages (const ImageFile& rRHS, ImageFile& result) const\r
 379 {\r
 380   if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) {\r
 381     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");\r
 382     return false;\r
 383   }\r
 384 \r
 385   ImageFileArrayConst vLHS = getArray();\r
 386   ImageFileArrayConst vRHS = rRHS.getArray();\r
 387   ImageFileArray vResult = result.getArray();\r
 388 \r
 389   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
 390     ImageFileColumnConst in1 = vLHS[ix];\r
 391     ImageFileColumnConst in2 = vRHS[ix];\r
 392     ImageFileColumn out = vResult[ix];\r
 393     for (unsigned int iy = 0; iy < m_ny; iy++)\r
 394         *out++ = *in1++ * *in2++;\r
 395   }\r
 396 \r
 397     return true;\r
 398 }\r
 399 \r
 400 bool\r
 401 ImageFile::divideImages (const ImageFile& rRHS, ImageFile& result) const\r
 402 {\r
 403   if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) {\r
 404     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");\r
 405     return false;\r
 406   }\r
 407 \r
 408   ImageFileArrayConst vLHS = getArray();\r
 409   ImageFileArrayConst vRHS = rRHS.getArray();\r
 410   ImageFileArray vResult = result.getArray();\r
 411 \r
 412   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
 413     ImageFileColumnConst in1 = vLHS[ix];\r
 414     ImageFileColumnConst in2 = vRHS[ix];\r
 415     ImageFileColumn out = vResult[ix];\r
 416     for (unsigned int iy = 0; iy < m_ny; iy++) {\r
 417       if (*in2 != 0.)\r
 418         *out++ = *in1++ / *in2++;\r
 419       else\r
 420         *out++ = 0;\r
 421     }\r
 422   }\r
 423 \r
 424     return true;\r
 425 }\r
 426 \r
 427 \r
 428 bool\r
 429 ImageFile::invertPixelValues (ImageFile& result) const\r
 430 {\r
 431   if (m_nx != result.nx() || m_ny != result.ny()) {\r
 432     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
 433     return false;\r
 434   }\r
 435 \r
 436   ImageFileArrayConst vLHS = getArray();\r
 437   ImageFileArray vResult = result.getArray();\r
 438 \r
 439   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
 440     ImageFileColumnConst in = vLHS[ix];\r
 441     ImageFileColumn out = vResult[ix];\r
 442     for (unsigned int iy = 0; iy < m_ny; iy++)\r
 443         *out++ = - *in++;\r
 444   }\r
 445 \r
 446   return true;\r
 447 }\r
 448 \r
 449 bool\r
 450 ImageFile::sqrt (ImageFile& result) const\r
 451 {\r
 452   if (m_nx != result.nx() || m_ny != result.ny()) {\r
 453     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
 454     return false;\r
 455   }\r
 456 \r
 457   ImageFileArrayConst vLHS = getArray();\r
 458   ImageFileArray vResult = result.getArray();\r
 459 \r
 460   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
 461     ImageFileColumnConst in = vLHS[ix];\r
 462     ImageFileColumn out = vResult[ix];\r
 463     for (unsigned int iy = 0; iy < m_ny; iy++)\r
 464       if (*in < 0)\r
 465         *out++ = -::sqrt(-*in++);\r
 466       else\r
 467         *out++ = ::sqrt(*in++);\r
 468   }\r
 469 \r
 470   return true;\r
 471 }\r
 472 \r
 473 bool\r
 474 ImageFile::log (ImageFile& result) const\r
 475 {\r
 476   if (m_nx != result.nx() || m_ny != result.ny()) {\r
 477     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
 478     return false;\r
 479   }\r
 480 \r
 481   ImageFileArrayConst vLHS = getArray();\r
 482   ImageFileArray vResult = result.getArray();\r
 483 \r
 484   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
 485     ImageFileColumnConst in = vLHS[ix];\r
 486     ImageFileColumn out = vResult[ix];\r
 487     for (unsigned int iy = 0; iy < m_ny; iy++)\r
 488       if (*in <= 0)\r
 489         *out++ = 0;\r
 490       else\r
 491         *out++ = ::log(*in++);\r
 492   }\r
 493 \r
 494   return true;\r
 495 }\r
 496 \r
 497 bool\r
 498 ImageFile::exp (ImageFile& result) const\r
 499 {\r
 500   if (m_nx != result.nx() || m_ny != result.ny()) {\r
 501     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
 502     return false;\r
 503   }\r
 504 \r
 505   ImageFileArrayConst vLHS = getArray();\r
 506   ImageFileArray vResult = result.getArray();\r
 507 \r
 508   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
 509     ImageFileColumnConst in = vLHS[ix];\r
 510     ImageFileColumn out = vResult[ix];\r
 511     for (unsigned int iy = 0; iy < m_ny; iy++)\r
 512       *out++ = ::exp (*in++);\r
 513   }\r
 514 \r
 515   return true;\r
 516 }\r
 517 \r
 518 bool\r
 519 ImageFile::inverseFourier (ImageFile& result) const\r
 520 {\r
 521   if (m_nx != result.nx() || m_ny != result.ny()) {\r
 522     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
 523     return false;\r
 524   }\r
 525 \r
 526   return true;\r
 527 }\r
 528 \r
 529 bool\r
 530 ImageFile::fourier (ImageFile& result) const\r
 531 {\r
 532   if (m_nx != result.nx() || m_ny != result.ny()) {\r
 533     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
 534     return false;\r
 535   }\r
 536 \r
 537   if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {\r
 538     if (! result.reallocRealToComplex ())\r
 539       return false;\r
 540   }\r
 541 \r
 542   ImageFileArrayConst vLHS = getArray();\r
 543   ImageFileArray vRealResult = result.getArray();\r
 544   ImageFileArray vImagResult = result.getImaginaryArray();\r
 545 \r
 546   unsigned int ix, iy;\r
 547   double* pY = new double [m_ny];\r
 548   std::complex<double>** complexOut = new std::complex<double>* [m_nx];\r
 549   for (ix = 0; ix < m_nx; ix++)\r
 550     complexOut[ix] = new std::complex<double> [m_ny];\r
 551 \r
 552   for (ix = 0; ix < m_nx; ix++) {\r
 553     for (iy = 0; iy < m_ny; iy++)\r
 554       pY[iy] = vLHS[ix][iy];\r
 555     ProcessSignal::finiteFourierTransform (pY, complexOut[ix], m_ny,  ProcessSignal::FORWARD);\r
 556   }\r
 557   delete pY;\r
 558 \r
 559   std::complex<double>* pX = new std::complex<double> [m_nx];\r
 560   std::complex<double>* complexOutCol = new std::complex<double> [m_nx];\r
 561   for (iy = 0; iy < m_ny; iy++) {\r
 562     for (ix = 0; ix < m_nx; ix++)\r
 563       pX[ix] = complexOut[ix][iy];\r
 564     ProcessSignal::finiteFourierTransform (pX, complexOutCol, m_nx, ProcessSignal::FORWARD);\r
 565     for (ix = 0; ix < m_nx; ix++)\r
 566       complexOut[ix][iy] = complexOutCol[ix];\r
 567   }\r
 568   delete [] pX;\r
 569 \r
 570   for (ix = 0; ix < m_nx; ix++)\r
 571     ProcessSignal::shuffleFourierToNaturalOrder (complexOut[ix], m_ny);\r
 572 \r
 573   \r
 574   for (iy = 0; iy < m_ny; iy++) {\r
 575     for (ix = 0; ix < m_nx; ix++)\r
 576       complexOutCol[ix] = complexOut[ix][iy];\r
 577     ProcessSignal::shuffleFourierToNaturalOrder (complexOutCol, m_nx);;\r
 578     for (ix = 0; ix < m_nx; ix++)\r
 579       complexOut[ix][iy] = complexOutCol[ix];\r
 580 \r
 581   }\r
 582   delete [] complexOutCol;\r
 583 \r
 584   for (ix = 0; ix < m_nx; ix++)\r
 585     for (iy = 0; iy < m_ny; iy++) {\r
 586       vRealResult[ix][iy] = complexOut[ix][iy].real();\r
 587       vImagResult[ix][iy] = complexOut[ix][iy].imag();\r
 588     }\r
 589 \r
 590   for (ix = 0; ix < m_nx; ix++)\r
 591     delete [] complexOut[ix];\r
 592 \r
 593   delete [] complexOut;\r
 594 \r
 595   return true;\r
 596 }\r
 597 \r
 598 bool\r
 599 ImageFile::magnitude (ImageFile& result) const\r
 600 {\r
 601   if (m_nx != result.nx() || m_ny != result.ny()) {\r
 602     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
 603     return false;\r
 604   }\r
 605 \r
 606   ImageFileArray vReal = getArray();\r
 607   ImageFileArray vImag = NULL;\r
 608   if (dataType() == Array2dFile::DATA_TYPE_COMPLEX)\r
 609     vImag = getImaginaryArray();\r
 610 \r
 611   ImageFileArray vRealResult = result.getArray();\r
 612   if (result.dataType() == Array2dFile::DATA_TYPE_COMPLEX) {\r
 613     ImageFileArray vImagResult = result.getImaginaryArray();\r
 614     for (unsigned int ix = 0; ix < m_nx; ix++)\r
 615       for (unsigned int iy = 0; iy < m_ny; iy++)\r
 616         vImagResult[ix][iy] = 0;\r
 617   }\r
 618 \r
 619   for (unsigned int ix = 0; ix < m_nx; ix++)\r
 620     for (unsigned int iy = 0; iy < m_ny; iy++) {\r
 621       if (dataType() == Array2dFile::DATA_TYPE_COMPLEX) \r
 622         vRealResult[ix][iy] = ::sqrt (vReal[ix][iy] * vReal[ix][iy] + vImag[ix][iy] * vImag[ix][iy]);\r
 623       else\r
 624         vRealResult[ix][iy] = vReal[ix][iy];\r
 625     }\r
 626 \r
 627   return true;\r
 628 }\r
 629 \r
 630 bool\r
 631 ImageFile::phase (ImageFile& result) const\r
 632 {\r
 633   if (m_nx != result.nx() || m_ny != result.ny()) {\r
 634     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
 635     return false;\r
 636   }\r
 637 \r
 638   ImageFileArray vReal = getArray();\r
 639   ImageFileArray vImag = NULL;\r
 640   if (dataType() == Array2dFile::DATA_TYPE_COMPLEX)\r
 641     vImag = getImaginaryArray();\r
 642 \r
 643   ImageFileArray vRealResult = result.getArray();\r
 644   if (result.dataType() == Array2dFile::DATA_TYPE_COMPLEX) {\r
 645     ImageFileArray vImagResult = result.getImaginaryArray();\r
 646     for (unsigned int ix = 0; ix < m_nx; ix++)\r
 647       for (unsigned int iy = 0; iy < m_ny; iy++)\r
 648         vImagResult[ix][iy] = 0;\r
 649   }\r
 650 \r
 651   for (unsigned int ix = 0; ix < m_nx; ix++)\r
 652     for (unsigned int iy = 0; iy < m_ny; iy++) {\r
 653       if (dataType() == Array2dFile::DATA_TYPE_COMPLEX) \r
 654         vRealResult[ix][iy] = ::atan (vImag[ix][iy] / vReal[ix][iy]);\r
 655       else\r
 656         vRealResult[ix][iy] = vReal[ix][iy];\r
 657     }\r
 658 \r
 659   return true;\r
 660 }\r
 661 \r
 662 bool\r
 663 ImageFile::square (ImageFile& result) const\r
 664 {\r
 665   if (m_nx != result.nx() || m_ny != result.ny()) {\r
 666     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");\r
 667     return false;\r
 668   }\r
 669 \r
 670   ImageFileArrayConst vLHS = getArray();\r
 671   ImageFileArray vResult = result.getArray();\r
 672 \r
 673   for (unsigned int ix = 0; ix < m_nx; ix++) {\r
 674     ImageFileColumnConst in = vLHS[ix];\r
 675     ImageFileColumn out = vResult[ix];\r
 676     for (unsigned int iy = 0; iy < m_ny; iy++) {\r
 677         *out++ = *in * *in;\r
 678         in++;\r
 679     }\r
 680   }\r
 681 \r
 682   return true;\r
 683 }\r
 684 \r
 685 \r
 686 void
 687 ImageFile::writeImagePGM (const char *outfile, int nxcell, int nycell, double densmin, double densmax)
 688 {
 689         FILE *fp;
 690         int nx = m_nx;
 691         int ny = m_ny;
 692         ImageFileArray v = getArray();
 693
 694         unsigned char* rowp = new unsigned char [nx * nxcell];
 695
 696         if ((fp = fopen (outfile, "wb")) == NULL)
 697                 return;
 698
 699         fprintf(fp, "P5\n");
 700         fprintf(fp, "%d %d\n", nx, ny);
 701         fprintf(fp, "255\n");
 702
 703         for (int irow = ny - 1; irow >= 0; irow--) {
 704                 for (int icol = 0; icol < nx; icol++) {
 705                         int pos = icol * nxcell;
 706                         double dens = (v[icol][irow] - densmin) / (densmax - densmin);
 707                         dens = clamp (dens, 0., 1.);
 708                         for (int p = pos; p < pos + nxcell; p++) {
 709                                 rowp[p] = static_cast<unsigned int> (dens * 255.);
 710                         }
 711                 }
 712                 for (int ir = 0; ir < nycell; ir++) {
 713                         for (int ic = 0; ic < nx * nxcell; ic++)
 714                                 fputc( rowp[ic], fp );
 715                 }
 716         }
 717         \r
 718         delete rowp;
 719         fclose(fp);
 720 }
 721
 722 void
 723 ImageFile::writeImagePGMASCII (const char *outfile, int nxcell, int nycell, double densmin, double densmax)
 724 {
 725         FILE *fp;
 726         int nx = m_nx;
 727         int ny = m_ny;
 728         ImageFileArray v = getArray();
 729
 730         unsigned char* rowp = new unsigned char [nx * nxcell];
 731
 732         if ((fp = fopen (outfile, "wb")) == NULL)
 733                 return;
 734
 735         fprintf(fp, "P2\n");
 736         fprintf(fp, "%d %d\n", nx, ny);
 737         fprintf(fp, "255\n");
 738
 739         for (int irow = ny - 1; irow >= 0; irow--) {
 740                 for (int icol = 0; icol < nx; icol++) {
 741                         int pos = icol * nxcell;
 742                         double dens = (v[icol][irow] - densmin) / (densmax - densmin);
 743                         dens = clamp (dens, 0., 1.);
 744                         for (int p = pos; p < pos + nxcell; p++) {
 745                                 rowp[p] = static_cast<unsigned int> (dens * 255.);
 746                         }
 747                 }
 748                 for (int ir = 0; ir < nycell; ir++) {
 749                         for (int ic = 0; ic < nx * nxcell; ic++)
 750                                 fprintf(fp, "%d ", rowp[ic]);
 751                         fprintf(fp, "\n");
 752                 }
 753         }
 754         \r
 755         delete rowp;
 756         fclose(fp);
 757 }
 758
 759
 760 #ifdef HAVE_PNG
 761 void
 762 ImageFile::writeImagePNG (const char *outfile, int bitdepth, int nxcell, int nycell, double densmin, double densmax)
 763 {
 764         FILE *fp;
 765         png_structp png_ptr;
 766         png_infop info_ptr;
 767         double max_out_level = (1 << bitdepth) - 1;
 768         int nx = m_nx;
 769         int ny = m_ny;
 770         ImageFileArray v = getArray();
 771
 772         unsigned char* rowp = new unsigned char [nx * nxcell * (bitdepth / 8)];
 773
 774         if ((fp = fopen (outfile, "wb")) == NULL)
 775                 return;
 776
 777         png_ptr = png_create_write_struct (PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
 778         if (! png_ptr)
 779                 return;
 780
 781         info_ptr = png_create_info_struct(png_ptr);
 782         if (! info_ptr) {
 783                 png_destroy_write_struct(&png_ptr, (png_infopp) NULL);
 784                 fclose(fp);
 785                 return;
 786         }
 787
 788         if (setjmp(png_ptr->jmpbuf)) {
 789                 png_destroy_write_struct(&png_ptr, &info_ptr);
 790                 fclose(fp);
 791                 return;
 792         }
 793
 794         png_init_io(png_ptr, fp);
 795
 796         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);
 797
 798         png_write_info(png_ptr, info_ptr);
 799         for (int irow = ny - 1; irow >= 0; irow--) {
 800                 png_bytep row_pointer = rowp;
 801
 802                 for (int icol = 0; icol < nx; icol++) {
 803                         int pos = icol * nxcell;
 804                         double dens = (v[icol][irow] - densmin) / (densmax - densmin);
 805                         dens = clamp (dens, 0., 1.);
 806                         unsigned int outval = static_cast<unsigned int> (dens * max_out_level);
 807
 808                         for (int p = pos; p < pos + nxcell; p++) {
 809                                 if (bitdepth == 8)
 810                                         rowp[p] = outval;
 811                                 else {
 812                                         int rowpos = p * 2;
 813                                         rowp[rowpos] = (outval >> 8) & 0xFF;
 814                                         rowp[rowpos+1] = (outval & 0xFF);
 815                                 }
 816                         }
 817                 }
 818                 for (int ir = 0; ir < nycell; ir++)
 819                         png_write_rows (png_ptr, &row_pointer, 1);
 820         }
 821
 822         png_write_end(png_ptr, info_ptr);
 823         png_destroy_write_struct(&png_ptr, &info_ptr);
 824         delete rowp;\r
 825
 826         fclose(fp);
 827 }
 828 #endif
 829
 830 #ifdef HAVE_GD
 831 #include "gd.h"
 832 static const int N_GRAYSCALE=256;
 833
 834 void
 835 ImageFile::writeImageGIF (const char *outfile, int nxcell, int nycell, double densmin, double densmax)
 836 {
 837         gdImagePtr gif;
 838         FILE *out;
 839         int gs_indices[N_GRAYSCALE];
 840         int nx = m_nx;
 841         int ny = m_ny;
 842         ImageFileArray v = getArray();
 843
 844         unsigned char rowp [nx * nxcell];
 845         if (rowp == NULL)
 846                 return;
 847
 848         gif = gdImageCreate(nx * nxcell, ny * nycell);
 849         for (int i = 0; i < N_GRAYSCALE; i++)
 850                 gs_indices[i] = gdImageColorAllocate(gif, i, i, i);
 851
 852         int lastrow = ny * nycell - 1;
 853         for (int irow = 0; irow < ny; irow++) {
 854                 int rpos = irow * nycell;
 855                 for (int ir = rpos; ir < rpos + nycell; ir++) {
 856                         for (int icol = 0; icol < nx; icol++) {
 857                                 int cpos = icol * nxcell;
 858                                 double dens = (v[icol][irow] - densmin) / (densmax - densmin);
 859                                 dens = clamp(dens, 0., 1.);
 860                                 for (int ic = cpos; ic < cpos + nxcell; ic++) {
 861                                         rowp[ic] = (unsigned int) (dens * (double) (N_GRAYSCALE - 1));
 862                                         gdImageSetPixel(gif, ic, lastrow - ir, gs_indices[rowp[ic]]);
 863                                 }
 864                         }
 865                 }
 866         }
 867
 868         if ((out = fopen(outfile,"w")) == NULL) {
 869                 sys_error(ERR_FATAL, "Error opening output file %s for writing", outfile);
 870                 return (1);
 871         }
 872         gdImageGif(gif,out);
 873         fclose(out);
 874         gdImageDestroy(gif);
 875 }
 876 #endif
 877