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