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-2001 Kevin Rosenberg
  11 **
  12 **  $Id: imagefile.cpp,v 1.41 2001/03/18 18:08:25 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 #ifdef HAVE_CTN_DICOM
  30 #include "ctndicom.h"
  31 #endif
  32
  33 const double ImageFile::s_dRedGrayscaleFactor = 0.299;
  34 const double ImageFile::s_dGreenGrayscaleFactor = 0.587;
  35 const double ImageFile::s_dBlueGrayscaleFactor = 0.114;
  36
  37
  38 const int ImageFile::EXPORT_FORMAT_INVALID = -1;
  39 const int ImageFile::EXPORT_FORMAT_TEXT = 0;
  40 const int ImageFile::EXPORT_FORMAT_PGM = 1;
  41 const int ImageFile::EXPORT_FORMAT_PGMASCII = 2;
  42 #ifdef HAVE_PNG
  43 const int ImageFile::EXPORT_FORMAT_PNG = 3;
  44 const int ImageFile::EXPORT_FORMAT_PNG16 = 4;
  45 #endif
  46 #ifdef HAVE_CTN_DICOM
  47 const int ImageFile::EXPORT_FORMAT_DICOM = 5;
  48 #endif
  49
  50 const char* ImageFile::s_aszExportFormatName[] =
  51 {
  52   {"text"},
  53   {"pgm"},
  54   {"pgmascii"},
  55 #ifdef HAVE_PNG
  56   {"png"},
  57   {"png16"},
  58 #endif
  59 #ifdef HAVE_CTN_DICOM
  60   {"dicom"},
  61 #endif
  62 };
  63
  64 const char* ImageFile::s_aszExportFormatTitle[] =
  65 {
  66   {"Text"},
  67   {"PGM"},
  68   {"PGM ASCII"},
  69   {"PNG"},
  70   {"PNG 16-bit"},
  71 #ifdef HAVE_CTN_DICOM
  72   {"Dicom"},
  73 #endif
  74 };
  75 const int ImageFile::s_iExportFormatCount = sizeof(s_aszExportFormatName) / sizeof(const char*);
  76
  77
  78 const int ImageFile::IMPORT_FORMAT_INVALID = -1;
  79 const int ImageFile::IMPORT_FORMAT_PPM = 0;
  80 #ifdef HAVE_PNG
  81 const int ImageFile::IMPORT_FORMAT_PNG = 1;
  82 #endif
  83 #ifdef HAVE_CTN_DICOM
  84 const int ImageFile::IMPORT_FORMAT_DICOM = 2;
  85 #endif
  86
  87
  88 const char* ImageFile::s_aszImportFormatName[] =
  89 {
  90   {"ppm"},
  91 #ifdef HAVE_PNG
  92   {"png"},
  93 #endif
  94 #ifdef HAVE_CTN_DICOM
  95   {"dicom"},
  96 #endif
  97 };
  98
  99 const char* ImageFile::s_aszImportFormatTitle[] =
 100 {
 101   {"PPM"},
 102   {"PNG"},
 103 #ifdef HAVE_CTN_DICOM
 104   {"Dicom"},
 105 #endif
 106 };
 107 const int ImageFile::s_iImportFormatCount = sizeof(s_aszImportFormatName) / sizeof(const char*);
 108
 109
 110
 111 F32Image::F32Image (int nx, int ny, int dataType)
 112 : Array2dFile (nx, ny, sizeof(kfloat32), Array2dFile::PIXEL_FLOAT32, dataType)
 113 {
 114 }
 115
 116 F32Image::F32Image (void)
 117 : Array2dFile()
 118 {
 119   setPixelFormat (Array2dFile::PIXEL_FLOAT32);
 120   setPixelSize (sizeof(kfloat32));
 121   setDataType (Array2dFile::DATA_TYPE_REAL);
 122 }
 123
 124 F64Image::F64Image (int nx, int ny, int dataType)
 125 : Array2dFile (nx, ny, sizeof(kfloat64), Array2dFile::PIXEL_FLOAT64, dataType)
 126 {
 127 }
 128
 129 F64Image::F64Image (void)
 130 : Array2dFile ()
 131 {
 132   setPixelFormat (PIXEL_FLOAT64);
 133   setPixelSize (sizeof(kfloat64));
 134   setDataType (Array2dFile::DATA_TYPE_REAL);
 135 }
 136
 137 void
 138 ImageFile::getCenterCoordinates (unsigned int& iXCenter, unsigned int& iYCenter)
 139 {
 140   if (isEven (m_nx))
 141     iXCenter = m_nx / 2;
 142   else
 143     iXCenter = (m_nx - 1) / 2;
 144
 145   if (isEven (m_ny))
 146     iYCenter = m_ny / 2;
 147   else
 148     iYCenter = (m_ny - 1) / 2;
 149 }
 150
 151
 152 void
 153 ImageFile::filterResponse (const char* const domainName, double bw, const char* const filterName, double filt_param, double dInputScale, double dOutputScale)
 154 {
 155   ImageFileArray v = getArray();
 156   SignalFilter filter (filterName, domainName, bw, filt_param);
 157
 158   unsigned int iXCenter, iYCenter;
 159   getCenterCoordinates (iXCenter, iYCenter);
 160
 161   for (unsigned int ix = 0; ix < m_nx; ix++)
 162     for (unsigned int iy = 0; iy < m_ny; iy++) {
 163       long lD2 = ((ix - iXCenter) * (ix - iXCenter)) + ((iy - iYCenter) * (iy - iYCenter));
 164       double r = ::sqrt (static_cast<double>(lD2)) * dInputScale;
 165       v[ix][iy] = filter.response (r) * dOutputScale;
 166     }
 167 }
 168
 169 int
 170 ImageFile::display (void) const
 171 {
 172   double pmin, pmax;
 173
 174   getMinMax (pmin, pmax);
 175
 176   return (displayScaling (1, pmin, pmax));
 177 }
 178
 179 int
 180 ImageFile::displayScaling (const int scale, const ImageFileValue pmin, const ImageFileValue pmax) const
 181 {
 182   int nx = m_nx;
 183   int ny = m_ny;
 184   ImageFileArrayConst v = getArray();
 185   if (v == NULL || nx == 0 || ny == 0)
 186     return 0;
 187
 188 #if HAVE_G2_H
 189   int* pPens = new int [nx * ny * scale * scale ];
 190
 191   double view_scale = 255 / (pmax - pmin);
 192   int id_X11 = g2_open_X11 (nx * scale, ny * scale);
 193   int grayscale[256];
 194   for (int i = 0; i < 256; i++) {
 195     double cval = i / 255.;
 196     grayscale[i] = g2_ink (id_X11, cval, cval, cval);
 197   }
 198
 199   for (int iy = ny - 1; iy >= 0; iy--) {
 200     int iRowPos = ((ny - 1 - iy) * scale) * (nx * scale);
 201     for (int ix = 0; ix < nx; ix++) {
 202       int cval = static_cast<int>((v[ix][iy] - pmin) * view_scale);
 203       if (cval < 0)
 204         cval = 0;
 205       else if (cval > 255)
 206         cval = 255;
 207       for (int sy = 0; sy < scale; sy++)
 208         for (int sx = 0; sx < scale; sx++)
 209           pPens[iRowPos+(sy * nx * scale)+(sx + (ix * scale))] = grayscale[cval];
 210     }
 211   }
 212
 213   g2_image (id_X11, 0., 0., nx * scale, ny * scale, pPens);
 214
 215   delete pPens;
 216   return (id_X11);
 217 #else
 218   return 0;
 219 #endif
 220 }
 221
 222
 223
 224 // ImageFile::comparativeStatistics    Calculate comparative stats
 225 //
 226 // OUTPUT
 227 //   d   Normalized root mean squared distance measure
 228 //   r   Normalized mean absolute distance measure
 229 //   e   Worst case distance measure
 230 //
 231 // REFERENCES
 232 //  G.T. Herman, Image Reconstruction From Projections, 1980
 233
 234 bool
 235 ImageFile::comparativeStatistics (const ImageFile& imComp, double& d, double& r, double& e) const
 236 {
 237   if (imComp.nx() != m_nx && imComp.ny() != m_ny) {
 238     sys_error (ERR_WARNING, "Image sizes differ [ImageFile::comparativeStatistics]");
 239     return false;
 240   }
 241   ImageFileArrayConst v = getArray();
 242   if (v == NULL || m_nx == 0 || m_ny == 0)
 243     return false;
 244
 245   ImageFileArrayConst vComp = imComp.getArray();
 246
 247   double myMean = 0.;
 248   for (unsigned int ix = 0; ix < m_nx; ix++) {
 249     for (unsigned int iy = 0; iy < m_ny; iy++) {
 250       myMean += v[ix][iy];
 251     }
 252   }
 253   myMean /= (m_nx * m_ny);
 254
 255   double sqErrorSum = 0.;
 256   double absErrorSum = 0.;
 257   double sqDiffFromMeanSum = 0.;
 258   double absValueSum = 0.;
 259   for (unsigned int ix2 = 0; ix2 < m_nx; ix2++) {
 260     for (unsigned int iy = 0; iy < m_ny; iy++) {
 261       double diff = v[ix2][iy] - vComp[ix2][iy];
 262       sqErrorSum += diff * diff;
 263       absErrorSum += fabs(diff);
 264       double diffFromMean = v[ix2][iy] - myMean;
 265       sqDiffFromMeanSum += diffFromMean * diffFromMean;
 266       absValueSum += fabs(v[ix2][iy]);
 267     }
 268   }
 269
 270   d = ::sqrt (sqErrorSum / sqDiffFromMeanSum);
 271   r = absErrorSum / absValueSum;
 272
 273   int hx = m_nx / 2;
 274   int hy = m_ny / 2;
 275   double eMax = -1;
 276   for (int ix3 = 0; ix3 < hx; ix3++) {
 277     for (int iy = 0; iy < hy; iy++) {
 278       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]);
 279       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]);
 280       double error = fabs (avgPixel - avgPixelComp);
 281       if (error > eMax)
 282         eMax = error;
 283     }
 284   }
 285
 286   e = eMax;
 287
 288   return true;
 289 }
 290
 291
 292 bool
 293 ImageFile::printComparativeStatistics (const ImageFile& imComp, std::ostream& os) const
 294 {
 295   double d, r, e;
 296
 297   if (comparativeStatistics (imComp, d, r, e)) {
 298     os << "  Normalized root mean squared distance (d): " << d << std::endl;
 299     os << "      Normalized mean absolute distance (r): " << r << std::endl;
 300     os << "Worst case distance (2x2 pixel average) (e): " << e << std::endl;
 301     return true;
 302   }
 303   return false;
 304 }
 305
 306
 307 void
 308 ImageFile::printStatistics (std::ostream& os) const
 309 {
 310   double min, max, mean, mode, median, stddev;
 311
 312   statistics (min, max, mean, mode, median, stddev);
 313   if (isComplex())
 314     os << "Real Component Statistics" << std::endl;
 315
 316   os << "   min: " << min << std::endl;
 317   os << "   max: " << max << std::endl;
 318   os << "  mean: " << mean << std::endl;
 319   os << "  mode: " << mode << std::endl;
 320   os << "median: " << median << std::endl;
 321   os << "stddev: " << stddev << std::endl;
 322
 323   if (isComplex()) {
 324     statistics (getImaginaryArray(), min, max, mean, mode, median, stddev);
 325     os << std::endl << "Imaginary Component Statistics" << std::endl;
 326     os << "   min: " << min << std::endl;
 327     os << "   max: " << max << std::endl;
 328     os << "  mean: " << mean << std::endl;
 329     os << "  mode: " << mode << std::endl;
 330     os << "median: " << median << std::endl;
 331     os << "stddev: " << stddev << std::endl;
 332   }
 333 }
 334
 335
 336 void
 337 ImageFile::statistics (double& min, double& max, double& mean, double& mode, double& median, double& stddev) const
 338 {
 339   ImageFileArrayConst v = getArray();
 340   statistics (v, min, max, mean, mode, median, stddev);
 341 }
 342
 343
 344 void
 345 ImageFile::statistics (ImageFileArrayConst v, double& min, double& max, double& mean, double& mode, double& median, double& stddev) const
 346 {
 347   int nx = m_nx;
 348   int ny = m_ny;
 349
 350   if (v == NULL || nx == 0 || ny == 0)
 351     return;
 352
 353   std::vector<double> vecImage;
 354   int iVec = 0;
 355   vecImage.resize (nx * ny);
 356   for (int ix = 0; ix < nx; ix++) {
 357     for (int iy = 0; iy < ny; iy++)
 358       vecImage[iVec++] = v[ix][iy];
 359   }
 360
 361   vectorNumericStatistics (vecImage, nx * ny, min, max, mean, mode, median, stddev);
 362 }
 363
 364 void
 365 ImageFile::getMinMax (double& min, double& max) const
 366 {
 367   int nx = m_nx;
 368   int ny = m_ny;
 369   ImageFileArrayConst v = getArray();
 370
 371   if (v == NULL || nx == 0 || ny == 0)
 372     return;
 373
 374   min = v[0][0];
 375   max = v[0][0];
 376   for (int ix = 0; ix < nx; ix++) {
 377     for (int iy = 0; iy < ny; iy++) {
 378       if (v[ix][iy] > max)
 379         max = v[ix][iy];
 380       if (v[ix][iy] < min)
 381         min = v[ix][iy];
 382     }
 383   }
 384 }
 385
 386 bool
 387 ImageFile::convertRealToComplex ()
 388 {
 389   if (dataType() != Array2dFile::DATA_TYPE_REAL)
 390     return false;
 391
 392   if (! reallocRealToComplex())
 393     return false;
 394
 395   ImageFileArray vImag = getImaginaryArray();
 396   for (unsigned int ix = 0; ix < m_nx; ix++) {
 397     ImageFileColumn vCol = vImag[ix];
 398     for (unsigned int iy = 0; iy < m_ny; iy++)
 399       *vCol++ = 0;
 400   }
 401
 402   return true;
 403 }
 404
 405 bool
 406 ImageFile::convertComplexToReal ()
 407 {
 408   if (dataType() != Array2dFile::DATA_TYPE_COMPLEX)
 409     return false;
 410
 411   ImageFileArray vReal = getArray();
 412   ImageFileArray vImag = getImaginaryArray();
 413   for (unsigned int ix = 0; ix < m_nx; ix++) {
 414     ImageFileColumn vRealCol = vReal[ix];
 415     ImageFileColumn vImagCol = vImag[ix];
 416     for (unsigned int iy = 0; iy < m_ny; iy++) {
 417       CTSimComplex c (*vRealCol, *vImagCol);
 418       *vRealCol++ = std::abs (c);
 419       vImagCol++;
 420     }
 421   }
 422
 423   return reallocComplexToReal();
 424 }
 425
 426 bool
 427 ImageFile::subtractImages (const ImageFile& rRHS, ImageFile& result) const
 428 {
 429   if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) {
 430     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");
 431     return false;
 432   }
 433
 434   if (isComplex() || rRHS.isComplex() && ! result.isComplex())
 435     result.convertRealToComplex();
 436
 437   ImageFileArrayConst vLHS = getArray();
 438   ImageFileArrayConst vLHSImag = getImaginaryArray();
 439   ImageFileArrayConst vRHS = rRHS.getArray();
 440   ImageFileArrayConst vRHSImag = rRHS.getImaginaryArray();
 441   ImageFileArray vResult = result.getArray();
 442   ImageFileArray vResultImag = result.getImaginaryArray();
 443
 444   for (unsigned int ix = 0; ix < m_nx; ix++) {
 445     for (unsigned int iy = 0; iy < m_ny; iy++) {
 446       vResult[ix][iy] = vLHS[ix][iy] - vRHS[ix][iy];
 447       if (result.isComplex()) {
 448         vResultImag[ix][iy] = 0;
 449         if (isComplex())
 450           vResultImag[ix][iy] += vLHSImag[ix][iy];
 451         if (rRHS.isComplex())
 452           vResultImag[ix][iy] -= vRHSImag[ix][iy];
 453       }
 454     }
 455   }
 456
 457   return true;
 458 }
 459
 460 bool
 461 ImageFile::addImages (const ImageFile& rRHS, ImageFile& result) const
 462 {
 463   if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) {
 464     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");
 465     return false;
 466   }
 467
 468   if (isComplex() || rRHS.isComplex() && ! result.isComplex())
 469     result.convertRealToComplex();
 470
 471   ImageFileArrayConst vLHS = getArray();
 472   ImageFileArrayConst vLHSImag = getImaginaryArray();
 473   ImageFileArrayConst vRHS = rRHS.getArray();
 474   ImageFileArrayConst vRHSImag = rRHS.getImaginaryArray();
 475   ImageFileArray vResult = result.getArray();
 476   ImageFileArray vResultImag = result.getImaginaryArray();
 477
 478   for (unsigned int ix = 0; ix < m_nx; ix++) {
 479     for (unsigned int iy = 0; iy < m_ny; iy++) {
 480       vResult[ix][iy] = vLHS[ix][iy] + vRHS[ix][iy];
 481       if (result.isComplex()) {
 482         vResultImag[ix][iy] = 0;
 483         if (isComplex())
 484           vResultImag[ix][iy] += vLHSImag[ix][iy];
 485         if (rRHS.isComplex())
 486           vResultImag[ix][iy] += vRHSImag[ix][iy];
 487       }
 488     }
 489   }
 490
 491   return true;
 492 }
 493
 494 bool
 495 ImageFile::multiplyImages (const ImageFile& rRHS, ImageFile& result) const
 496 {
 497   if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) {
 498     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");
 499     return false;
 500   }
 501
 502   if (isComplex() || rRHS.isComplex() && ! result.isComplex())
 503     result.convertRealToComplex();
 504
 505   ImageFileArrayConst vLHS = getArray();
 506   ImageFileArrayConst vLHSImag = getImaginaryArray();
 507   ImageFileArrayConst vRHS = rRHS.getArray();
 508   ImageFileArrayConst vRHSImag = rRHS.getImaginaryArray();
 509   ImageFileArray vResult = result.getArray();
 510   ImageFileArray vResultImag = result.getImaginaryArray();
 511
 512   for (unsigned int ix = 0; ix < m_nx; ix++) {
 513     for (unsigned int iy = 0; iy < m_ny; iy++) {
 514       if (result.isComplex()) {
 515         double dImag = 0;
 516         if (isComplex())
 517           dImag = vLHSImag[ix][iy];
 518         std::complex<double> cLHS (vLHS[ix][iy], dImag);
 519         dImag = 0;
 520         if (rRHS.isComplex())
 521           dImag = vRHSImag[ix][iy];
 522         std::complex<double> cRHS (vRHS[ix][iy], dImag);
 523         std::complex<double> cResult = cLHS * cRHS;
 524         vResult[ix][iy] = cResult.real();
 525         vResultImag[ix][iy] = cResult.imag();
 526       } else
 527         vResult[ix][iy] = vLHS[ix][iy] * vRHS[ix][iy];
 528     }
 529   }
 530
 531
 532   return true;
 533 }
 534
 535 bool
 536 ImageFile::divideImages (const ImageFile& rRHS, ImageFile& result) const
 537 {
 538   if (m_nx != rRHS.nx() || m_ny != rRHS.ny() || m_nx != result.nx() || m_ny != result.ny()) {
 539     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::subtractImage]");
 540     return false;
 541   }
 542
 543   if (isComplex() || rRHS.isComplex() && ! result.isComplex())
 544     result.convertRealToComplex();
 545
 546   ImageFileArrayConst vLHS = getArray();
 547   ImageFileArrayConst vLHSImag = getImaginaryArray();
 548   ImageFileArrayConst vRHS = rRHS.getArray();
 549   ImageFileArrayConst vRHSImag = rRHS.getImaginaryArray();
 550   ImageFileArray vResult = result.getArray();
 551   ImageFileArray vResultImag = result.getImaginaryArray();
 552
 553   for (unsigned int ix = 0; ix < m_nx; ix++) {
 554     for (unsigned int iy = 0; iy < m_ny; iy++) {
 555       if (result.isComplex()) {
 556         double dImag = 0;
 557         if (isComplex())
 558           dImag = vLHSImag[ix][iy];
 559         std::complex<double> cLHS (vLHS[ix][iy], dImag);
 560         dImag = 0;
 561         if (rRHS.isComplex())
 562           dImag = vRHSImag[ix][iy];
 563         std::complex<double> cRHS (vRHS[ix][iy], dImag);
 564         std::complex<double> cResult = cLHS / cRHS;
 565         vResult[ix][iy] = cResult.real();
 566         vResultImag[ix][iy] = cResult.imag();
 567       } else {
 568         if (vRHS != 0)
 569           vResult[ix][iy] = vLHS[ix][iy] / vRHS[ix][iy];
 570         else
 571           vResult[ix][iy] = 0;
 572       }
 573     }
 574   }
 575
 576   return true;
 577 }
 578
 579
 580 bool
 581 ImageFile::invertPixelValues (ImageFile& result) const
 582 {
 583   if (m_nx != result.nx() || m_ny != result.ny()) {
 584     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
 585     return false;
 586   }
 587
 588   if (isComplex() && ! result.isComplex())
 589     result.convertRealToComplex();
 590
 591   ImageFileArrayConst vLHS = getArray();
 592   ImageFileArray vResult = result.getArray();
 593
 594   for (unsigned int ix = 0; ix < m_nx; ix++) {
 595     ImageFileColumnConst in = vLHS[ix];
 596     ImageFileColumn out = vResult[ix];
 597     for (unsigned int iy = 0; iy < m_ny; iy++)
 598       *out++ = - *in++;
 599   }
 600
 601   return true;
 602 }
 603
 604 bool
 605 ImageFile::sqrt (ImageFile& result) const
 606 {
 607   if (m_nx != result.nx() || m_ny != result.ny()) {
 608     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
 609     return false;
 610   }
 611
 612   if (isComplex() && ! result.isComplex())
 613     result.convertRealToComplex();
 614
 615   bool bComplexOutput = result.isComplex();
 616   ImageFileArrayConst vLHS = getArray();
 617   if (! bComplexOutput)   // check if should convert to complex output
 618     for (unsigned int ix = 0; ix < m_nx; ix++)
 619       for (unsigned int iy = 0; iy < m_ny; iy++)
 620         if (! bComplexOutput && vLHS[ix][iy] < 0) {
 621           result.convertRealToComplex();
 622           bComplexOutput = true;
 623           break;
 624         }
 625
 626         ImageFileArrayConst vLHSImag = getImaginaryArray();
 627         ImageFileArray vResult = result.getArray();
 628         ImageFileArray vResultImag = result.getImaginaryArray();
 629
 630         for (unsigned int ix = 0; ix < m_nx; ix++) {
 631           for (unsigned int iy = 0; iy < m_ny; iy++) {
 632             if (result.isComplex()) {
 633               double dImag = 0;
 634               if (isComplex())
 635                 dImag = vLHSImag[ix][iy];
 636               std::complex<double> cLHS (vLHS[ix][iy], dImag);
 637               std::complex<double> cResult = std::sqrt(cLHS);
 638               vResult[ix][iy] = cResult.real();
 639               vResultImag[ix][iy] = cResult.imag();
 640             } else
 641               vResult[ix][iy] = ::sqrt (vLHS[ix][iy]);
 642           }
 643         }
 644
 645
 646         return true;
 647 }
 648
 649 bool
 650 ImageFile::log (ImageFile& result) const
 651 {
 652   if (m_nx != result.nx() || m_ny != result.ny()) {
 653     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::log]");
 654     return false;
 655   }
 656
 657   if (isComplex() && ! result.isComplex())
 658     result.convertRealToComplex();
 659
 660   ImageFileArrayConst vLHS = getArray();
 661   ImageFileArrayConst vLHSImag = getImaginaryArray();
 662   ImageFileArray vResult = result.getArray();
 663   ImageFileArray vResultImag = result.getImaginaryArray();
 664
 665   for (unsigned int ix = 0; ix < m_nx; ix++) {
 666     for (unsigned int iy = 0; iy < m_ny; iy++) {
 667       if (result.isComplex()) {
 668         double dImag = 0;
 669         if (isComplex())
 670           dImag = vLHSImag[ix][iy];
 671         std::complex<double> cLHS (vLHS[ix][iy], dImag);
 672         std::complex<double> cResult = std::log (cLHS);
 673         vResult[ix][iy] = cResult.real();
 674         vResultImag[ix][iy] = cResult.imag();
 675       } else {
 676         if (vLHS[ix][iy] > 0)
 677           vResult[ix][iy] = ::log (vLHS[ix][iy]);
 678         else
 679           vResult[ix][iy] = 0;
 680       }
 681     }
 682   }
 683
 684
 685   return true;
 686 }
 687
 688 bool
 689 ImageFile::exp (ImageFile& result) const
 690 {
 691   if (m_nx != result.nx() || m_ny != result.ny()) {
 692     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
 693     return false;
 694   }
 695
 696   if (isComplex() && ! result.isComplex())
 697     result.convertRealToComplex();
 698
 699   ImageFileArrayConst vLHS = getArray();
 700   ImageFileArrayConst vLHSImag = getImaginaryArray();
 701   ImageFileArray vResult = result.getArray();
 702   ImageFileArray vResultImag = result.getImaginaryArray();
 703
 704   for (unsigned int ix = 0; ix < m_nx; ix++) {
 705     for (unsigned int iy = 0; iy < m_ny; iy++) {
 706       if (result.isComplex()) {
 707         double dImag = 0;
 708         if (isComplex())
 709           dImag = vLHSImag[ix][iy];
 710         std::complex<double> cLHS (vLHS[ix][iy], dImag);
 711         std::complex<double> cResult = std::exp (cLHS);
 712         vResult[ix][iy] = cResult.real();
 713         vResultImag[ix][iy] = cResult.imag();
 714       } else
 715         vResult[ix][iy] = ::exp (vLHS[ix][iy]);
 716     }
 717   }
 718
 719
 720   return true;
 721 }
 722
 723 bool
 724 ImageFile::scaleImage (ImageFile& result) const
 725 {
 726   unsigned int nx = m_nx;
 727   unsigned int ny = m_ny;
 728   unsigned int newNX = result.nx();
 729   unsigned int newNY = result.ny();
 730
 731   double dXScale = static_cast<double>(newNX) / static_cast<double>(nx);
 732   double dYScale = static_cast<double>(newNY) / static_cast<double>(ny);
 733
 734   if (isComplex() && ! result.isComplex())
 735     result.convertRealToComplex();
 736
 737   ImageFileArrayConst vReal = getArray();
 738   ImageFileArrayConst vImag = getImaginaryArray();
 739   ImageFileArray vResult = result.getArray();
 740   ImageFileArray vResultImag = result.getImaginaryArray();
 741
 742   for (unsigned int ix = 0; ix < newNX; ix++) {
 743     for (unsigned int iy = 0; iy < newNY; iy++) {
 744       double dXPos = ix / dXScale;
 745       double dYPos = iy / dYScale;
 746       unsigned int scaleNX = static_cast<unsigned int> (dXPos);
 747       unsigned int scaleNY = static_cast<unsigned int> (dYPos);
 748       double dXFrac = dXPos - scaleNX;
 749       double dYFrac = dYPos - scaleNY;
 750       if (scaleNX >= nx - 1 || scaleNY >= ny - 1) {
 751         vResult[ix][iy] = vReal[scaleNX][scaleNY];
 752         if (result.isComplex()) {
 753           if (isComplex())
 754             vResultImag[ix][iy] = vImag[scaleNX][scaleNY];
 755           else
 756             vResultImag[ix][iy] = 0;
 757         }
 758       } else {
 759         vResult[ix][iy] = (1 - dXFrac) * (1 - dYFrac) * vReal[scaleNX][scaleNY] +
 760           dXFrac * (1 - dYFrac) * vReal[scaleNX+1][scaleNY] +
 761           dYFrac * (1 - dXFrac) * vReal[scaleNX][scaleNY+1] +
 762           dXFrac * dYFrac * vReal[scaleNX+1][scaleNY+1];
 763         if (result.isComplex()) {
 764           if (isComplex())
 765             vResultImag[ix][iy] = (1 - dXFrac) * (1 - dYFrac) * vImag[scaleNX][scaleNY] +
 766             dXFrac * (1 - dYFrac) * vImag[scaleNX+1][scaleNY] +
 767             dYFrac * (1 - dXFrac) * vImag[scaleNX][scaleNY+1] +
 768             dXFrac * dYFrac * vImag[scaleNX+1][scaleNY+1];
 769           else
 770             vResultImag[ix][iy] = 0;
 771         }
 772       }
 773     }
 774   }
 775
 776   return true;
 777 }
 778
 779 #ifdef HAVE_FFTW
 780 bool
 781 ImageFile::fft (ImageFile& result) const
 782 {
 783   if (m_nx != result.nx() || m_ny != result.ny()) {
 784     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
 785     return false;
 786   }
 787
 788   if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
 789     if (! result.convertRealToComplex ())
 790       return false;
 791   }
 792
 793   fftw_complex* in = new fftw_complex [m_nx * m_ny];
 794
 795   ImageFileArrayConst vReal = getArray();
 796   ImageFileArrayConst vImag = getImaginaryArray();
 797
 798   unsigned int ix, iy;
 799   unsigned int iArray = 0;
 800   for (ix = 0; ix < m_nx; ix++)
 801     for (iy = 0; iy < m_ny; iy++) {
 802       in[iArray].re = vReal[ix][iy];
 803       if (isComplex())
 804         in[iArray].im = vImag[ix][iy];
 805       else
 806         in[iArray].im = 0;
 807       iArray++;
 808     }
 809
 810     fftwnd_plan plan = fftw2d_create_plan (m_nx, m_ny, FFTW_FORWARD, FFTW_IN_PLACE);
 811
 812     fftwnd_one (plan, in, NULL);
 813
 814     ImageFileArray vRealResult = result.getArray();
 815     ImageFileArray vImagResult = result.getImaginaryArray();
 816     iArray = 0;
 817     unsigned int iScale = m_nx * m_ny;
 818     for (ix = 0; ix < m_nx; ix++)
 819       for (iy = 0; iy < m_ny; iy++) {
 820         vRealResult[ix][iy] = in[iArray].re / iScale;
 821         vImagResult[ix][iy] = in[iArray].im / iScale;
 822         iArray++;
 823       }
 824
 825       fftwnd_destroy_plan (plan);
 826       delete in;
 827
 828
 829       Fourier::shuffleFourierToNaturalOrder (result);
 830
 831       return true;
 832 }
 833
 834
 835 bool
 836 ImageFile::ifft (ImageFile& result) const
 837 {
 838   if (m_nx != result.nx() || m_ny != result.ny()) {
 839     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
 840     return false;
 841   }
 842
 843   if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
 844     if (! result.convertRealToComplex ())
 845       return false;
 846   }
 847
 848   ImageFileArrayConst vReal = getArray();
 849   ImageFileArrayConst vImag = getImaginaryArray();
 850   ImageFileArray vRealResult = result.getArray();
 851   ImageFileArray vImagResult = result.getImaginaryArray();
 852   unsigned int ix, iy;
 853   for (ix = 0; ix < m_nx; ix++)
 854     for (iy = 0; iy < m_ny; iy++) {
 855       vRealResult[ix][iy] = vReal[ix][iy];
 856       if (isComplex())
 857         vImagResult[ix][iy] = vImag[ix][iy];
 858       else
 859         vImagResult[ix][iy] = 0;
 860     }
 861
 862     Fourier::shuffleNaturalToFourierOrder (result);
 863
 864     fftw_complex* in = new fftw_complex [m_nx * m_ny];
 865
 866     unsigned int iArray = 0;
 867     for (ix = 0; ix < m_nx; ix++)
 868       for (iy = 0; iy < m_ny; iy++) {
 869         in[iArray].re = vRealResult[ix][iy];
 870         in[iArray].im = vImagResult[ix][iy];
 871         iArray++;
 872       }
 873
 874       fftwnd_plan plan = fftw2d_create_plan (m_nx, m_ny, FFTW_BACKWARD, FFTW_IN_PLACE);
 875
 876       fftwnd_one (plan, in, NULL);
 877
 878       iArray = 0;
 879       for (ix = 0; ix < m_nx; ix++)
 880         for (iy = 0; iy < m_ny; iy++) {
 881           vRealResult[ix][iy] = in[iArray].re;
 882           vImagResult[ix][iy] = in[iArray].im;
 883           iArray++;
 884         }
 885
 886         fftwnd_destroy_plan (plan);
 887
 888         delete in;
 889
 890         return true;
 891 }
 892
 893 bool
 894 ImageFile::fftRows (ImageFile& result) const
 895 {
 896   if (m_nx != result.nx() || m_ny != result.ny()) {
 897     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::fftRows]");
 898     return false;
 899   }
 900
 901   if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
 902     if (! result.convertRealToComplex ())
 903       return false;
 904   }
 905
 906   fftw_complex* in = new fftw_complex [m_nx];
 907
 908   ImageFileArrayConst vReal = getArray();
 909   ImageFileArrayConst vImag = getImaginaryArray();
 910
 911   fftw_plan plan = fftw_create_plan (m_nx, FFTW_FORWARD, FFTW_IN_PLACE);
 912   std::complex<double>* pcRow = new std::complex<double> [m_nx];
 913
 914   unsigned int ix, iy;
 915   unsigned int iArray = 0;
 916   for (iy = 0; iy < m_ny; iy++) {
 917     for (ix = 0; ix < m_nx; ix++) {
 918       in[ix].re = vReal[ix][iy];
 919       if (isComplex())
 920         in[ix].im = vImag[ix][iy];
 921       else
 922         in[ix].im = 0;
 923     }
 924
 925     fftw_one (plan, in, NULL);
 926
 927     for (ix = 0; ix < m_nx; ix++)
 928       pcRow[ix] = std::complex<double>(in[ix].re, in[ix].im);
 929
 930     Fourier::shuffleFourierToNaturalOrder (pcRow, m_nx);
 931     for (ix = 0; ix < m_nx; ix++) {
 932       vReal[ix][iy] = pcRow[ix].real() / m_nx;
 933       vImag[ix][iy] = pcRow[ix].imag() / m_nx;
 934     }
 935   }
 936   delete [] pcRow;
 937
 938   fftw_destroy_plan (plan);
 939   delete in;
 940
 941   return true;
 942 }
 943
 944 bool
 945 ImageFile::ifftRows (ImageFile& result) const
 946 {
 947   if (m_nx != result.nx() || m_ny != result.ny()) {
 948     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::fftRows]");
 949     return false;
 950   }
 951
 952   if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
 953     if (! result.convertRealToComplex ())
 954       return false;
 955   }
 956
 957   fftw_complex* in = new fftw_complex [m_nx];
 958
 959   ImageFileArrayConst vReal = getArray();
 960   ImageFileArrayConst vImag = getImaginaryArray();
 961
 962   fftw_plan plan = fftw_create_plan (m_nx, FFTW_BACKWARD, FFTW_IN_PLACE);
 963   std::complex<double>* pcRow = new std::complex<double> [m_nx];
 964
 965   unsigned int ix, iy;
 966   unsigned int iArray = 0;
 967   for (iy = 0; iy < m_ny; iy++) {
 968     for (ix = 0; ix < m_nx; ix++) {
 969       double dImag = 0;
 970       if (isComplex())
 971         dImag = vImag[ix][iy];
 972       pcRow[ix] = std::complex<double> (vReal[ix][iy], dImag);
 973     }
 974
 975     Fourier::shuffleNaturalToFourierOrder (pcRow, m_nx);
 976
 977     for (ix = 0; ix < m_nx; ix++) {
 978       in[ix].re = pcRow[ix].real();
 979       in[ix].im = pcRow[ix].imag();
 980     }
 981
 982     fftw_one (plan, in, NULL);
 983
 984     for (ix = 0; ix < m_nx; ix++) {
 985       vReal[ix][iy] = in[ix].re;
 986       vImag[ix][iy] = in[ix].im;
 987     }
 988   }
 989   delete [] pcRow;
 990
 991   fftw_destroy_plan (plan);
 992   delete in;
 993
 994   return true;
 995 }
 996
 997 bool
 998 ImageFile::fftCols (ImageFile& result) const
 999 {
1000   if (m_nx != result.nx() || m_ny != result.ny()) {
1001     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::fftRows]");
1002     return false;
1003   }
1004
1005   if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
1006     if (! result.convertRealToComplex ())
1007       return false;
1008   }
1009
1010   fftw_complex* in = new fftw_complex [m_ny];
1011
1012   ImageFileArrayConst vReal = getArray();
1013   ImageFileArrayConst vImag = getImaginaryArray();
1014
1015   fftw_plan plan = fftw_create_plan (m_ny, FFTW_FORWARD, FFTW_IN_PLACE);
1016   std::complex<double>* pcCol = new std::complex<double> [m_ny];
1017
1018   unsigned int ix, iy;
1019   unsigned int iArray = 0;
1020   for (ix = 0; ix < m_nx; ix++) {
1021     for (iy = 0; iy < m_ny; iy++) {
1022       in[iy].re = vReal[ix][iy];
1023       if (isComplex())
1024         in[iy].im = vImag[ix][iy];
1025       else
1026         in[iy].im = 0;
1027     }
1028
1029     fftw_one (plan, in, NULL);
1030
1031     for (iy = 0; iy < m_ny; iy++)
1032       pcCol[iy] = std::complex<double>(in[iy].re, in[iy].im);
1033
1034     Fourier::shuffleFourierToNaturalOrder (pcCol, m_ny);
1035     for (iy = 0; iy < m_ny; iy++) {
1036       vReal[ix][iy] = pcCol[iy].real() / m_ny;
1037       vImag[ix][iy] = pcCol[iy].imag() / m_ny;
1038     }
1039   }
1040   delete [] pcCol;
1041
1042   fftw_destroy_plan (plan);
1043   delete in;
1044
1045   return true;
1046 }
1047
1048 bool
1049 ImageFile::ifftCols (ImageFile& result) const
1050 {
1051   if (m_nx != result.nx() || m_ny != result.ny()) {
1052     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::fftRows]");
1053     return false;
1054   }
1055
1056   if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
1057     if (! result.convertRealToComplex ())
1058       return false;
1059   }
1060
1061   fftw_complex* in = new fftw_complex [m_ny];
1062
1063   ImageFileArrayConst vReal = getArray();
1064   ImageFileArrayConst vImag = getImaginaryArray();
1065
1066   fftw_plan plan = fftw_create_plan (m_ny, FFTW_BACKWARD, FFTW_IN_PLACE);
1067   std::complex<double>* pcCol = new std::complex<double> [m_ny];
1068
1069   unsigned int ix, iy;
1070   unsigned int iArray = 0;
1071   for (ix = 0; ix < m_nx; ix++) {
1072     for (iy = 0; iy < m_ny; iy++) {
1073       double dImag = 0;
1074       if (isComplex())
1075         dImag = vImag[ix][iy];
1076       pcCol[iy] = std::complex<double> (vReal[ix][iy], dImag);
1077     }
1078
1079     Fourier::shuffleNaturalToFourierOrder (pcCol, m_ny);
1080
1081     for (iy = 0; iy < m_ny; iy++) {
1082       in[iy].re = pcCol[iy].real();
1083       in[iy].im = pcCol[iy].imag();
1084     }
1085
1086     fftw_one (plan, in, NULL);
1087
1088     for (iy = 0; iy < m_ny; iy++) {
1089       vReal[ix][iy] = in[iy].re;
1090       vImag[ix][iy] = in[iy].im;
1091     }
1092   }
1093   delete [] pcCol;
1094
1095   fftw_destroy_plan (plan);
1096   delete in;
1097
1098   return true;
1099 }
1100
1101 #endif // HAVE_FFTW
1102
1103
1104
1105 bool
1106 ImageFile::fourier (ImageFile& result) const
1107 {
1108   if (m_nx != result.nx() || m_ny != result.ny()) {
1109     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
1110     return false;
1111   }
1112
1113   if (! result.isComplex())
1114     if (! result.convertRealToComplex ())
1115       return false;
1116
1117     ImageFileArrayConst vLHS = getArray();
1118     ImageFileArrayConst vLHSImag = getImaginaryArray();
1119     ImageFileArray vRealResult = result.getArray();
1120     ImageFileArray vImagResult = result.getImaginaryArray();
1121
1122     unsigned int ix, iy;
1123
1124     // alloc output matrix
1125     CTSimComplex** complexOut = new CTSimComplex* [m_nx];
1126     for (ix = 0; ix < m_nx; ix++)
1127       complexOut[ix] = new CTSimComplex [m_ny];
1128
1129     // fourier each x column
1130     CTSimComplex* pY = new CTSimComplex [m_ny];
1131     for (ix = 0; ix < m_nx; ix++) {
1132       for (iy = 0; iy < m_ny; iy++) {
1133         double dImag = 0;
1134         if (isComplex())
1135           dImag = vLHSImag[ix][iy];
1136         pY[iy] = std::complex<double>(vLHS[ix][iy], dImag);
1137       }
1138       ProcessSignal::finiteFourierTransform (pY, complexOut[ix], m_ny,  ProcessSignal::FORWARD);
1139     }
1140     delete [] pY;
1141
1142     // fourier each y row
1143     CTSimComplex* pX = new CTSimComplex [m_nx];
1144     CTSimComplex* complexOutRow = new CTSimComplex [m_nx];
1145     for (iy = 0; iy < m_ny; iy++) {
1146       for (ix = 0; ix < m_nx; ix++)
1147         pX[ix] = complexOut[ix][iy];
1148       ProcessSignal::finiteFourierTransform (pX, complexOutRow, m_nx, ProcessSignal::FORWARD);
1149       for (ix = 0; ix < m_nx; ix++)
1150         complexOut[ix][iy] = complexOutRow[ix];
1151     }
1152     delete [] pX;
1153     delete [] complexOutRow;
1154
1155     for (ix = 0; ix < m_nx; ix++)
1156       for (iy = 0; iy < m_ny; iy++) {
1157         vRealResult[ix][iy] = complexOut[ix][iy].real();
1158         vImagResult[ix][iy] = complexOut[ix][iy].imag();
1159       }
1160
1161       Fourier::shuffleFourierToNaturalOrder (result);
1162
1163       // delete complexOut matrix
1164       for (ix = 0; ix < m_nx; ix++)
1165         delete [] complexOut[ix];
1166       delete [] complexOut;
1167
1168       return true;
1169 }
1170
1171 bool
1172 ImageFile::inverseFourier (ImageFile& result) const
1173 {
1174   if (m_nx != result.nx() || m_ny != result.ny()) {
1175     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
1176     return false;
1177   }
1178
1179   if (result.dataType() == Array2dFile::DATA_TYPE_REAL) {
1180     if (! result.convertRealToComplex ())
1181       return false;
1182   }
1183
1184   ImageFileArrayConst vLHSReal = getArray();
1185   ImageFileArrayConst vLHSImag = getImaginaryArray();
1186   ImageFileArray vRealResult = result.getArray();
1187   ImageFileArray vImagResult = result.getImaginaryArray();
1188
1189   unsigned int ix, iy;
1190   // alloc 2d complex output matrix
1191   CTSimComplex** complexOut = new CTSimComplex* [m_nx];
1192   for (ix = 0; ix < m_nx; ix++)
1193     complexOut[ix] = new CTSimComplex [m_ny];
1194
1195   // put input image into result
1196   for (ix = 0; ix < m_nx; ix++)
1197     for (iy = 0; iy < m_ny; iy++) {
1198       vRealResult[ix][iy] = vLHSReal[ix][iy];
1199       if (isComplex())
1200         vImagResult[ix][iy] = vLHSImag[ix][iy];
1201       else
1202         vImagResult[ix][iy] = 0;
1203     }
1204
1205     Fourier::shuffleNaturalToFourierOrder (result);
1206
1207     // ifourier each x column
1208     CTSimComplex* pCol = new CTSimComplex [m_ny];
1209     for (ix = 0; ix < m_nx; ix++) {
1210       for (iy = 0; iy < m_ny; iy++) {
1211         pCol[iy] = std::complex<double> (vRealResult[ix][iy], vImagResult[ix][iy]);
1212       }
1213       ProcessSignal::finiteFourierTransform (pCol, complexOut[ix], m_ny,  ProcessSignal::BACKWARD);
1214     }
1215     delete [] pCol;
1216
1217     // ifourier each y row
1218     CTSimComplex* complexInRow = new CTSimComplex [m_nx];
1219     CTSimComplex* complexOutRow = new CTSimComplex [m_nx];
1220     for (iy = 0; iy < m_ny; iy++) {
1221       for (ix = 0; ix < m_nx; ix++)
1222         complexInRow[ix] = complexOut[ix][iy];
1223       ProcessSignal::finiteFourierTransform (complexInRow, complexOutRow, m_nx, ProcessSignal::BACKWARD);
1224       for (ix = 0; ix < m_nx; ix++)
1225         complexOut[ix][iy] = complexOutRow[ix];
1226     }
1227     delete [] complexInRow;
1228     delete [] complexOutRow;
1229
1230     for (ix = 0; ix < m_nx; ix++)
1231       for (iy = 0; iy < m_ny; iy++) {
1232         vRealResult[ix][iy] = complexOut[ix][iy].real();
1233         vImagResult[ix][iy] = complexOut[ix][iy].imag();
1234       }
1235
1236       // delete complexOut matrix
1237       for (ix = 0; ix < m_nx; ix++)
1238         delete [] complexOut[ix];
1239       delete [] complexOut;
1240
1241       return true;
1242 }
1243
1244
1245 bool
1246 ImageFile::magnitude (ImageFile& result) const
1247 {
1248   if (m_nx != result.nx() || m_ny != result.ny()) {
1249     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
1250     return false;
1251   }
1252
1253   ImageFileArray vReal = getArray();
1254   ImageFileArray vImag = getImaginaryArray();
1255   ImageFileArray vRealResult = result.getArray();
1256
1257   for (unsigned int ix = 0; ix < m_nx; ix++)
1258     for (unsigned int iy = 0; iy < m_ny; iy++) {
1259       if (isComplex())
1260         vRealResult[ix][iy] = ::sqrt (vReal[ix][iy] * vReal[ix][iy] + vImag[ix][iy] * vImag[ix][iy]);
1261       else
1262         vRealResult[ix][iy] = vReal[ix][iy];
1263     }
1264
1265     if (result.isComplex())
1266       result.convertComplexToReal();
1267
1268     return true;
1269 }
1270
1271 bool
1272 ImageFile::phase (ImageFile& result) const
1273 {
1274   if (m_nx != result.nx() || m_ny != result.ny()) {
1275     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
1276     return false;
1277   }
1278
1279   ImageFileArray vReal = getArray();
1280   ImageFileArray vImag = getImaginaryArray();
1281   ImageFileArray vRealResult = result.getArray();
1282
1283   for (unsigned int ix = 0; ix < m_nx; ix++)
1284     for (unsigned int iy = 0; iy < m_ny; iy++) {
1285       if (isComplex())
1286         vRealResult[ix][iy] = ::atan2 (vImag[ix][iy], vReal[ix][iy]);
1287       else
1288         vRealResult[ix][iy] = 0;
1289     }
1290
1291     if (result.isComplex())
1292       result.convertComplexToReal();
1293
1294     return true;
1295 }
1296
1297 bool
1298 ImageFile::square (ImageFile& result) const
1299 {
1300   if (m_nx != result.nx() || m_ny != result.ny()) {
1301     sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
1302     return false;
1303   }
1304
1305   if (isComplex() && ! result.isComplex())
1306     result.convertRealToComplex();
1307
1308   ImageFileArrayConst vLHS = getArray();
1309   ImageFileArrayConst vLHSImag = getImaginaryArray();
1310   ImageFileArray vResult = result.getArray();
1311   ImageFileArray vResultImag = result.getImaginaryArray();
1312
1313   for (unsigned int ix = 0; ix < m_nx; ix++) {
1314     for (unsigned int iy = 0; iy < m_ny; iy++) {
1315       if (result.isComplex()) {
1316         double dImag = 0;
1317         if (isComplex())
1318           dImag = vLHSImag[ix][iy];
1319         std::complex<double> cLHS (vLHS[ix][iy], dImag);
1320         std::complex<double> cResult = cLHS * cLHS;
1321         vResult[ix][iy] = cResult.real();
1322         vResultImag[ix][iy] = cResult.imag();
1323       } else
1324         vResult[ix][iy] = vLHS[ix][iy] * vLHS[ix][iy];
1325     }
1326   }
1327
1328
1329   return true;
1330 }
1331
1332 int
1333 ImageFile::convertExportFormatNameToID (const char* const formatName)
1334 {
1335   int formatID = EXPORT_FORMAT_INVALID;
1336
1337   for (int i = 0; i < s_iExportFormatCount; i++)
1338     if (strcasecmp (formatName, s_aszExportFormatName[i]) == 0) {
1339       formatID = i;
1340       break;
1341     }
1342
1343     return (formatID);
1344 }
1345
1346 const char*
1347 ImageFile::convertExportFormatIDToName (int formatID)
1348 {
1349   static const char *formatName = "";
1350
1351   if (formatID >= 0 && formatID < s_iExportFormatCount)
1352     return (s_aszExportFormatName[formatID]);
1353
1354   return (formatName);
1355 }
1356
1357 const char*
1358 ImageFile::convertExportFormatIDToTitle (const int formatID)
1359 {
1360   static const char *formatTitle = "";
1361
1362   if (formatID >= 0 && formatID < s_iExportFormatCount)
1363     return (s_aszExportFormatTitle[formatID]);
1364
1365   return (formatTitle);
1366 }
1367
1368 int
1369 ImageFile::convertImportFormatNameToID (const char* const formatName)
1370 {
1371   int formatID = IMPORT_FORMAT_INVALID;
1372
1373   for (int i = 0; i < s_iImportFormatCount; i++)
1374     if (strcasecmp (formatName, s_aszImportFormatName[i]) == 0) {
1375       formatID = i;
1376       break;
1377     }
1378
1379     return (formatID);
1380 }
1381
1382 const char*
1383 ImageFile::convertImportFormatIDToName (int formatID)
1384 {
1385   static const char *formatName = "";
1386
1387   if (formatID >= 0 && formatID < s_iImportFormatCount)
1388     return (s_aszImportFormatName[formatID]);
1389
1390   return (formatName);
1391 }
1392
1393 const char*
1394 ImageFile::convertImportFormatIDToTitle (const int formatID)
1395 {
1396   static const char *formatTitle = "";
1397
1398   if (formatID >= 0 && formatID < s_iImportFormatCount)
1399     return (s_aszImportFormatTitle[formatID]);
1400
1401   return (formatTitle);
1402 }
1403
1404 bool
1405 ImageFile::importImage (const char* const pszFormat, const char* const pszFilename)
1406 {
1407   int iFormatID = convertImportFormatNameToID (pszFormat);
1408
1409   if (iFormatID == IMPORT_FORMAT_PPM)
1410     return readImagePPM (pszFilename);
1411 #ifdef HAVE_PNG
1412   else if (iFormatID == IMPORT_FORMAT_PNG)
1413     return readImagePNG (pszFilename);
1414 #endif
1415
1416   sys_error (ERR_SEVERE, "Invalid format %s [ImageFile::importImage]", pszFormat);
1417   return false;
1418 }
1419
1420 void
1421 ImageFile::skipSpacePPM (FILE* fp)
1422 {
1423   int c = fgetc (fp);
1424   while (isspace (c) || c == '#') {
1425     if (c == '#') {   // comment until end of line
1426       c = fgetc(fp);
1427       while (c != 13 && c != 10)
1428         c = fgetc(fp);
1429     }
1430     else
1431       c = fgetc(fp);
1432   }
1433
1434   ungetc (c, fp);
1435 }
1436
1437 bool
1438 ImageFile::readImagePPM (const char* const pszFile)
1439 {
1440   FILE* fp = fopen (pszFile, "r");
1441   if ((fp = fopen (pszFile, "r")) == NULL)
1442     return false;
1443   char cSignature = toupper(fgetc(fp));
1444   if (cSignature != 'P') {
1445     fclose(fp);
1446     return false;
1447   }
1448   cSignature = fgetc(fp);
1449   if (cSignature == '5' || cSignature == '6') { // binary modes
1450     fclose(fp);
1451     fp = fopen(pszFile, "rb"); // reopen in binary mode
1452     fgetc(fp);
1453     fgetc(fp);
1454   } else if (cSignature != '2' && cSignature != '3') {
1455     fclose(fp);
1456     return false;
1457   }
1458
1459   int nRows, nCols, iMaxValue;
1460   skipSpacePPM (fp);
1461   if (fscanf (fp, "%d", &nCols) != 1) {
1462     fclose(fp);
1463     return false;
1464   }
1465   skipSpacePPM (fp);
1466   if (fscanf (fp, "%d", &nRows) != 1) {
1467     fclose(fp);
1468     return false;
1469   }
1470   skipSpacePPM (fp);
1471   if (fscanf (fp, "%d", &iMaxValue) != 1) {
1472     fclose(fp);
1473     return false;
1474   }
1475   setArraySize (nRows, nCols);
1476
1477   if (cSignature == '5' || cSignature == '6') { // binary modes
1478     int c = fgetc(fp);
1479     if (c == 13) {
1480       c = fgetc(fp);
1481       if (c != 10)  // read msdos 13-10 newline
1482         ungetc(c, fp);
1483     }
1484   } else
1485     skipSpacePPM (fp); // ascii may have comments
1486
1487   double dMaxValue = iMaxValue;
1488   ImageFileArray v = getArray();
1489   for (int iy = nRows - 1; iy >= 0; iy--) {
1490     for (int ix = 0; ix < nCols; ix++) {
1491       int iGS, iR, iG, iB;
1492       double dR, dG, dB;
1493       switch (cSignature) {
1494       case '2':
1495         if (fscanf(fp, "%d ", &iGS) != 1) {
1496           fclose(fp);
1497           return false;
1498         }
1499         v[ix][iy] = iGS / dMaxValue;
1500         break;
1501       case '5':
1502         iGS = fgetc(fp);
1503         if (iGS == EOF) {
1504           fclose(fp);
1505           return false;
1506         }
1507         v[ix][iy] = iGS / dMaxValue;
1508         break;
1509       case '3':
1510         if (fscanf (fp, "%d %d %d ", &iR, &iG, &iB) != 3) {
1511           fclose(fp);
1512           return false;
1513         }
1514         dR = iR / dMaxValue;
1515         dG = iG / dMaxValue;
1516         dB = iB / dMaxValue;
1517         v[ix][iy] = colorToGrayscale (dR, dG, dB);
1518         break;
1519       case '6':
1520         iR = fgetc(fp);
1521         iG = fgetc(fp);
1522         iB = fgetc(fp);
1523         if (iB == EOF) {
1524           fclose(fp);
1525           return false;
1526         }
1527         dR = iR / dMaxValue;
1528         dG = iG / dMaxValue;
1529         dB = iB / dMaxValue;
1530         v[ix][iy] = colorToGrayscale (dR, dG, dB);
1531         break;
1532       }
1533     }
1534   }
1535
1536   fclose(fp);
1537   return true;
1538 }
1539
1540 #ifdef HAVE_PNG
1541 bool
1542 ImageFile::readImagePNG (const char* const pszFile)
1543 {
1544   FILE* fp = fopen(pszFile, "rb");
1545   if (!fp)
1546     return false;
1547   unsigned char header[8];
1548   fread (header, 1, 8, fp);
1549   if (png_sig_cmp (header, 0, 8)) {
1550     fclose (fp);
1551     return false;
1552   }
1553
1554   png_structp png_ptr = png_create_read_struct (PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
1555   if (!png_ptr) {
1556     fclose(fp);
1557     return false;
1558   }
1559
1560   png_infop info_ptr = png_create_info_struct(png_ptr);
1561   if (!info_ptr) {
1562     png_destroy_read_struct(&png_ptr, (png_infopp)NULL, (png_infopp)NULL);
1563     fclose(fp);
1564     return false;
1565   }
1566
1567   png_infop end_info = png_create_info_struct(png_ptr);
1568   if (!end_info) {
1569     png_destroy_read_struct(&png_ptr, &info_ptr, (png_infopp)NULL);
1570     fclose(fp);
1571     return false;
1572   }
1573
1574   if (setjmp(png_ptr->jmpbuf)) {
1575     png_destroy_read_struct(&png_ptr, &info_ptr, &end_info);
1576     fclose(fp);
1577     return false;
1578   }
1579
1580   png_init_io(png_ptr, fp);
1581   png_set_sig_bytes(png_ptr, 8);
1582   png_read_info(png_ptr, info_ptr);
1583
1584   int width = png_get_image_width (png_ptr, info_ptr);
1585   int height = png_get_image_height (png_ptr, info_ptr);
1586   int bit_depth = png_get_bit_depth (png_ptr, info_ptr);
1587   int color_type = png_get_color_type (png_ptr, info_ptr);
1588
1589   if (color_type == PNG_COLOR_TYPE_PALETTE && bit_depth <= 8)
1590     png_set_expand(png_ptr);
1591
1592   if (color_type == PNG_COLOR_TYPE_GRAY && bit_depth < 8)
1593     png_set_expand(png_ptr);
1594
1595   if (bit_depth < 8)
1596     png_set_packing(png_ptr);
1597
1598   if (color_type & PNG_COLOR_MASK_ALPHA)
1599     png_set_strip_alpha(png_ptr);
1600
1601   if (bit_depth == 16)
1602     png_set_swap(png_ptr); // convert to little-endian format
1603
1604   png_read_update_info(png_ptr, info_ptr); // update with transformations
1605   int rowbytes = png_get_rowbytes (png_ptr, info_ptr);
1606   bit_depth = png_get_bit_depth (png_ptr, info_ptr);
1607   color_type = png_get_color_type (png_ptr, info_ptr);
1608
1609   png_bytep* row_pointers = new png_bytep [height];
1610   int i;
1611   for (i = 0; i < height; i++)
1612     row_pointers[i] = new unsigned char [rowbytes];
1613
1614   png_read_image(png_ptr, row_pointers);
1615
1616   setArraySize (width, height);
1617   ImageFileArray v = getArray();
1618   for (int iy = 0; iy < height; iy++) {
1619     for (int ix = 0; ix < width; ix++) {
1620       double dV;
1621       if (color_type == PNG_COLOR_TYPE_GRAY) {
1622         if (bit_depth == 8)
1623           dV = row_pointers[iy][ix] / 255.;
1624         else if (bit_depth == 16) {
1625           int iBase = ix * 2;
1626           dV = (row_pointers[iy][iBase] + (row_pointers[iy][iBase+1] << 8)) / 65536.;
1627         }
1628       } else if (color_type == PNG_COLOR_TYPE_RGB) {
1629         if (bit_depth == 8) {
1630           int iBase = ix * 3;
1631           double dR = row_pointers[iy][iBase] / 255.;
1632           double dG = row_pointers[iy][iBase+1] / 255.;
1633           double dB = row_pointers[iy][iBase+2] / 255.;
1634           dV = colorToGrayscale (dR, dG, dR);
1635         }
1636       }
1637       v[ix][height-iy-1] = dV;
1638     }
1639   }
1640
1641   png_read_end(png_ptr, end_info);
1642   png_destroy_read_struct(&png_ptr, &info_ptr, &end_info);
1643
1644   for (i = 0; i < height; i++)
1645     delete row_pointers[i];
1646   delete row_pointers;
1647
1648   fclose (fp);
1649   return true;
1650 }
1651 #endif
1652
1653 bool
1654 ImageFile::exportImage (const char* const pszFormat, const char* const pszFilename, int nxcell, int nycell, double densmin, double densmax)
1655 {
1656   int iFormatID = convertExportFormatNameToID (pszFormat);
1657
1658   if (iFormatID == EXPORT_FORMAT_PGM)
1659     return writeImagePGM (pszFilename, nxcell, nycell, densmin, densmax);
1660   else if (iFormatID == EXPORT_FORMAT_PGMASCII)
1661     return writeImagePGMASCII (pszFilename, nxcell, nycell, densmin, densmax);
1662   else if (iFormatID == EXPORT_FORMAT_TEXT)
1663     return writeImageText (pszFilename);
1664 #ifdef HAVE_PNG
1665   else if (iFormatID == EXPORT_FORMAT_PNG)
1666     return writeImagePNG (pszFilename, 8, nxcell, nycell, densmin, densmax);
1667   else if (iFormatID == EXPORT_FORMAT_PNG16)
1668     return writeImagePNG (pszFilename, 16, nxcell, nycell, densmin, densmax);
1669 #endif
1670 #ifdef HAVE_CTN_DICOM
1671   else if (iFormatID == EXPORT_FORMAT_DICOM) {
1672     DicomExporter dicomExport (this);
1673     bool bSuccess = dicomExport.writeFile (pszFilename);
1674     if (! bSuccess)
1675       sys_error (ERR_SEVERE, dicomExport.failMessage().c_str());
1676     return bSuccess;
1677   }
1678 #endif
1679
1680   sys_error (ERR_SEVERE, "Invalid format %s [ImageFile::exportImage]", pszFormat);
1681   return false;
1682 }
1683
1684
1685 bool
1686 ImageFile::writeImagePGM (const char* const outfile, int nxcell, int nycell, double densmin, double densmax)
1687 {
1688   FILE *fp;
1689   int nx = m_nx;
1690   int ny = m_ny;
1691   ImageFileArray v = getArray();
1692
1693   unsigned char* rowp = new unsigned char [nx * nxcell];
1694
1695   if ((fp = fopen (outfile, "wb")) == NULL)
1696     return false;
1697
1698   fprintf(fp, "P5\n");
1699   fprintf(fp, "%d %d\n", nx, ny);
1700   fprintf(fp, "255\n");
1701
1702   for (int irow = ny - 1; irow >= 0; irow--) {
1703     for (int icol = 0; icol < nx; icol++) {
1704       int pos = icol * nxcell;
1705       double dens = (v[icol][irow] - densmin) / (densmax - densmin);
1706       dens = clamp (dens, 0., 1.);
1707       for (int p = pos; p < pos + nxcell; p++) {
1708         rowp[p] = static_cast<unsigned int> (dens * 255.);
1709       }
1710     }
1711     for (int ir = 0; ir < nycell; ir++) {
1712       for (int ic = 0; ic < nx * nxcell; ic++)
1713         fputc( rowp[ic], fp );
1714     }
1715   }
1716
1717   delete rowp;
1718   fclose(fp);
1719
1720   return true;
1721 }
1722
1723 bool
1724 ImageFile::writeImagePGMASCII (const char* const outfile, int nxcell, int nycell, double densmin, double densmax)
1725 {
1726   FILE *fp;
1727   int nx = m_nx;
1728   int ny = m_ny;
1729   ImageFileArray v = getArray();
1730
1731   unsigned char* rowp = new unsigned char [nx * nxcell];
1732
1733   if ((fp = fopen (outfile, "wb")) == NULL)
1734     return false;
1735
1736   fprintf(fp, "P2\n");
1737   fprintf(fp, "%d %d\n", nx, ny);
1738   fprintf(fp, "255\n");
1739
1740   for (int irow = ny - 1; irow >= 0; irow--) {
1741     for (int icol = 0; icol < nx; icol++) {
1742       int pos = icol * nxcell;
1743       double dens = (v[icol][irow] - densmin) / (densmax - densmin);
1744       dens = clamp (dens, 0., 1.);
1745       for (int p = pos; p < pos + nxcell; p++) {
1746         rowp[p] = static_cast<unsigned int> (dens * 255.);
1747       }
1748     }
1749     for (int ir = 0; ir < nycell; ir++) {
1750       for (int ic = 0; ic < nx * nxcell; ic++)
1751         fprintf(fp, "%d ", rowp[ic]);
1752       fprintf(fp, "\n");
1753     }
1754   }
1755
1756   delete rowp;
1757   fclose(fp);
1758
1759   return true;
1760 }
1761
1762 bool
1763 ImageFile::writeImageText (const char* const outfile)
1764 {
1765   FILE *fp;
1766   int nx = m_nx;
1767   int ny = m_ny;
1768   ImageFileArray v = getArray();
1769   ImageFileArray vImag = getImaginaryArray();
1770
1771   if ((fp = fopen (outfile, "w")) == NULL)
1772     return false;
1773
1774   for (int irow = ny - 1; irow >= 0; irow--) {
1775     for (int icol = 0; icol < nx; icol++) {
1776       if (isComplex()) {
1777         if (vImag[icol][irow] >= 0)
1778           fprintf (fp, "%.9g+%.9gi ", v[icol][irow], vImag[icol][irow]);
1779         else
1780           fprintf (fp, "%.9g-%.9gi ", v[icol][irow], -vImag[icol][irow]);
1781       } else
1782         fprintf (fp, "%12.8g ", v[icol][irow]);
1783     }
1784     fprintf(fp, "\n");
1785   }
1786
1787   fclose(fp);
1788
1789   return true;
1790 }
1791
1792
1793 #ifdef HAVE_PNG
1794 bool
1795 ImageFile::writeImagePNG (const char* const outfile, int bitdepth, int nxcell, int nycell, double densmin, double densmax)
1796 {
1797   double max_out_level = (1 << bitdepth) - 1;
1798   int nx = m_nx;
1799   int ny = m_ny;
1800   ImageFileArray v = getArray();
1801
1802   unsigned char* rowp = new unsigned char [nx * nxcell * (bitdepth / 8)];
1803
1804   FILE *fp = fopen (outfile, "wb");
1805   if (! fp)
1806     return false;
1807
1808   png_structp png_ptr = png_create_write_struct (PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
1809   if (! png_ptr)
1810     return false;
1811
1812   png_infop info_ptr = png_create_info_struct (png_ptr);
1813   if (! info_ptr) {
1814     png_destroy_write_struct (&png_ptr, (png_infopp) NULL);
1815     fclose (fp);
1816     return false;
1817   }
1818
1819   if (setjmp (png_ptr->jmpbuf)) {
1820     png_destroy_write_struct (&png_ptr, &info_ptr);
1821     fclose (fp);
1822     return false;
1823   }
1824
1825   png_init_io(png_ptr, fp);
1826
1827   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);
1828
1829   png_write_info(png_ptr, info_ptr);
1830   for (int irow = ny - 1; irow >= 0; irow--) {
1831     png_bytep row_pointer = rowp;
1832
1833     for (int icol = 0; icol < nx; icol++) {
1834       int pos = icol * nxcell;
1835       double dens = (v[icol][irow] - densmin) / (densmax - densmin);
1836       dens = clamp (dens, 0., 1.);
1837       unsigned int outval = static_cast<unsigned int> (dens * max_out_level);
1838
1839       for (int p = pos; p < pos + nxcell; p++) {
1840         if (bitdepth == 8)
1841           rowp[p] = outval;
1842         else {
1843           int rowpos = p * 2;
1844           rowp[rowpos+1] = (outval >> 8) & 0xFF;
1845           rowp[rowpos] = (outval & 0xFF);
1846         }
1847       }
1848     }
1849     for (int ir = 0; ir < nycell; ir++)
1850       png_write_rows (png_ptr, &row_pointer, 1);
1851   }
1852
1853   png_write_end (png_ptr, info_ptr);
1854   png_destroy_write_struct (&png_ptr, &info_ptr);
1855   delete rowp;
1856
1857   fclose(fp);
1858
1859   return true;
1860 }
1861 #endif
1862
1863 #ifdef HAVE_GD
1864 #include "gd.h"
1865 static const int N_GRAYSCALE=256;
1866
1867 bool
1868 ImageFile::writeImageGIF (const char* const outfile, int nxcell, int nycell, double densmin, double densmax)
1869 {
1870   int gs_indices[N_GRAYSCALE];
1871   int nx = m_nx;
1872   int ny = m_ny;
1873   ImageFileArray v = getArray();
1874
1875   unsigned char* rowp = new unsigned char [nx * nxcell];
1876
1877   gdImagePtr gif = gdImageCreate(nx * nxcell, ny * nycell);
1878   for (int i = 0; i < N_GRAYSCALE; i++)
1879     gs_indices[i] = gdImageColorAllocate(gif, i, i, i);
1880
1881   int lastrow = ny * nycell - 1;
1882   for (int irow = 0; irow < ny; irow++) {
1883     int rpos = irow * nycell;
1884     for (int ir = rpos; ir < rpos + nycell; ir++) {
1885       for (int icol = 0; icol < nx; icol++) {
1886         int cpos = icol * nxcell;
1887         double dens = (v[icol][irow] - densmin) / (densmax - densmin);
1888         dens = clamp(dens, 0., 1.);
1889         for (int ic = cpos; ic < cpos + nxcell; ic++) {
1890           rowp[ic] = (unsigned int) (dens * (double) (N_GRAYSCALE - 1));
1891           gdImageSetPixel(gif, ic, lastrow - ir, gs_indices[rowp[ic]]);
1892         }
1893       }
1894     }
1895   }
1896
1897   FILE *out;
1898   if ((out = fopen (outfile,"w")) == NULL) {
1899     sys_error(ERR_SEVERE, "Error opening output file %s for writing", outfile);
1900     return false;
1901   }
1902   gdImageGif(gif,out);
1903   fclose(out);
1904   gdImageDestroy(gif);
1905
1906   delete rowp;
1907
1908   return true;
1909 }
1910 #endif
1911