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