libctsim/imagefile.cpp

   1 /*****************************************************************************
   2 ** FILE IDENTIFICATION
   3 **
   4 **      Name:         imagefile.cpp
   5 **      Purpose:      Imagefile classes
   6 **      Programmer:   Kevin Rosenberg
   7 **      Date Started: June 2000
   8 **
   9 **  This is part of the CTSim program
  10 **  Copyright (C) 1983-2000 Kevin Rosenberg
  11 **
  12 **  $Id: imagefile.cpp,v 1.4 2000/06/25 17:32:24 kevin Exp $
  13 **
  14 **  This program is free software; you can redistribute it and/or modify
  15 **  it under the terms of the GNU General Public License (version 2) as
  16 **  published by the Free Software Foundation.
  17 **
  18 **  This program is distributed in the hope that it will be useful,
  19 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
  20 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  21 **  GNU General Public License for more details.
  22 **
  23 **  You should have received a copy of the GNU General Public License
  24 **  along with this program; if not, write to the Free Software
  25 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  26 ******************************************************************************/
  27
  28 #include "ct.h"
  29
  30
  31 ///////////////////////////////////////////////////////////////////////////
  32 // CLASS IMPLEMENTATION
  33 //
  34 //     Name: Array2dFileLabel
  35 //     Purpose: Labels for Array2dFiles
  36 ///////////////////////////////////////////////////////////////////////////
  37
  38 void
  39 Array2dFileLabel::init (void)
  40 {
  41     m_calcTime = 0;
  42     m_labelType = L_EMPTY;
  43     TIMEDATE td;
  44     td_get_tmdt (&td);
  45     m_year = td.d.year;
  46     m_month = td.d.month;
  47     m_day = td.d.date;
  48     m_hour = td.t.hour;
  49     m_minute = td.t.minute;
  50     m_second = td.t.second;
  51 }
  52
  53 Array2dFileLabel::Array2dFileLabel()
  54 {
  55     init();
  56 }
  57
  58 Array2dFileLabel::Array2dFileLabel(const char* const str, double ctime = 0.)
  59   : m_strLabel (str)
  60 {
  61     init();
  62
  63     m_labelType = L_USER;
  64     m_calcTime = ctime;
  65 }
  66
  67 Array2dFileLabel::Array2dFileLabel(const int type, const char* const str, double ctime = 0.)
  68   :  m_strLabel (str)
  69 {
  70     init();
  71
  72     m_labelType = type;
  73     m_calcTime = ctime;
  74 }
  75
  76 Array2dFileLabel::~Array2dFileLabel()
  77 {
  78 }
  79
  80 const string&
  81 Array2dFileLabel::getDateString (void) const
  82 {
  83   ostringstream oss;
  84   oss <<  m_month <<"/"<< m_day <<"/"<< m_year << " " << m_hour <<":"<<  m_minute <<":"<< m_second;
  85   m_strDate = oss.str();
  86   return m_strDate;
  87 }
  88
  89
  90 /* FILE
  91  *   image.c                            Routines for managing images
  92  */
  93
  94
  95 void
  96 ImageFile::filterResponse (const char* const domainName, double bw, const char* const filterName, double filt_param)
  97 {
  98   int hx = (mNX - 1) / 2;
  99   int hy = (mNY - 1) / 2;
 100   ImageFileArray v = getArray();
 101   SignalFilter filter (filterName, domainName, bw, filt_param);
 102
 103   for (int i = -hx; i <= hx; i++) {
 104     for (int j = -hy; j <= hy; j++) {
 105       double r = sqrt (i * i + j * j);
 106
 107       v[i+hx][j+hy] = filter.response (r);
 108     }
 109   }
 110 }
 111
 112 int
 113 ImageFile::display (void)
 114 {
 115     ImageFileValue pmin, pmax;
 116
 117     getPixelValueRange (pmin, pmax);
 118
 119     return (displayScaling (1, pmin, pmax));
 120 }
 121
 122 int
 123 ImageFile::displayScaling (const int scale, const ImageFileValue pmin, const ImageFileValue pmax)
 124 {
 125     int grayscale[256];
 126     int nx = mNX;
 127     int ny = mNY;
 128     ImageFileArray v = getArray();
 129
 130 #if HAVE_G2_H
 131     int pens [nx * ny * scale * scale ];
 132
 133     double view_scale = 255 / (pmax - pmin);
 134     int id_X11 = g2_open_X11 (nx * scale, ny * scale);
 135
 136     for (int i = 0; i < 256; i++) {
 137         double cval = i / 255.;
 138         grayscale[i] = g2_ink (id_X11, cval, cval, cval);
 139     }
 140
 141     for (int i= 0, iy = ny - 1; iy >= 0; iy--) {
 142         for (int ix = 0; ix < nx; ix++) {
 143             int cval = static_cast<int>((v[ix][iy] - pmin) * view_scale);
 144             if (cval < 0)
 145               cval = 0;
 146             else if (cval > 255)
 147               cval = 255;
 148             pens[i++] = grayscale[cval];
 149         }
 150     }
 151
 152     g2_image (id_X11, 0., 0., nx, ny, pens);
 153
 154     return (id_X11);
 155
 156 #endif
 157 }
 158
 159
 160
 161 // ImageFile::comparativeStatistics    Calculate comparative stats
 162 //
 163 // OUTPUT
 164 //   d   Normalized root mean squared distance measure
 165 //   r   Normalized mean absolute distance measure
 166 //   e   Worst case distance measure
 167 //
 168 // REFERENCES
 169 //  G.T. Herman, Image Reconstruction From Projections, 1980
 170
 171 bool
 172 ImageFile::comparativeStatistics (const ImageFile& imComp, double& d, double& r, double& e) const
 173 {
 174     if (imComp.nx() != mNX && imComp.ny() != mNY) {
 175         sys_error (ERR_WARNING, "Image sizes differ [ImageFile::comparativeStatistics]");
 176         return false;
 177     }
 178     ImageFileArray v = getArray();
 179     ImageFileArray vComp = imComp.getArray();
 180
 181     double myMean = 0.;
 182     for (int ix = 0; ix < mNX; ix++) {
 183         for (int iy = 0; iy < mNY; iy++) {
 184             myMean += v[ix][iy];
 185         }
 186     }
 187     myMean /= (mNX * mNY);
 188
 189     double sqErrorSum = 0.;
 190     double absErrorSum = 0.;
 191     double sqDiffFromMeanSum = 0.;
 192     double absValueSum = 0.;
 193     for (int ix = 0; ix < mNX; ix++) {
 194         for (int iy = 0; iy < mNY; iy++) {
 195             double diff = v[ix][iy] - vComp[ix][iy];
 196             sqErrorSum += diff * diff;
 197             absErrorSum += fabs(diff);
 198             double diffFromMean = v[ix][iy] - myMean;
 199             sqDiffFromMeanSum += diffFromMean * diffFromMean;
 200             absValueSum += fabs(v[ix][iy]);
 201         }
 202     }
 203
 204     d = sqrt (sqErrorSum / sqDiffFromMeanSum);
 205     r = absErrorSum / absValueSum;
 206
 207     int hx = mNX / 2;
 208     int hy = mNY / 2;
 209     double eMax = -1;
 210     for (int ix = 0; ix < hx; ix++) {
 211       for (int iy = 0; iy < hy; iy++) {
 212         double avgPixel = 0.25 * (v[2*ix][2*iy] + v[2*ix+1][2*iy] + v[2*ix][2*iy+1] + v[2*ix+1][2*iy+1]);
 213         double avgPixelComp = 0.25 * (vComp[2*ix][2*iy] + vComp[2*ix+1][2*iy] + vComp[2*ix][2*iy+1] + vComp[2*ix+1][2*iy+1]);
 214         double error = fabs (avgPixel - avgPixelComp);
 215         if (error > eMax)
 216           eMax = error;
 217       }
 218     }
 219
 220     e = eMax;
 221
 222     return true;
 223 }
 224
 225
 226 bool
 227 ImageFile::printComparativeStatistics (const ImageFile& imComp, ostream& os) const
 228 {
 229   double d, r, e;
 230
 231   if (comparativeStatistics (imComp, d, r, e)) {
 232     os << "  Normalized root mean squared distance (d): " << d << endl;
 233     os << "      Normalized mean absolute distance (r): " << r << endl;
 234     os << "Worst case distance (2x2 pixel average) (e): " << e << endl;
 235   }
 236 }
 237
 238
 239 void
 240 ImageFile::printStatistics (ostream& os) const
 241 {
 242     double min, max, mean, mode, median, stddev;
 243
 244     statistics (min, max, mean, mode, median, stddev);
 245
 246     os << "   min: " << min << endl;
 247     os << "   max: " << max << endl;
 248     os << "  mean: " << mean << endl;
 249     os << "  mode: " << mode << endl;
 250     os << "median: " << median << endl;
 251     os << "stddev: " << stddev << endl;
 252 }
 253
 254
 255 void
 256 ImageFile::statistics (double& min, double& max, double& mean, double& mode, double& median, double& stddev) const
 257 {
 258     int nx = mNX;
 259     int ny = mNY;
 260     ImageFileArray v = getArray();
 261
 262     mean = 0;
 263     min = v[0][0];
 264     max = v[0][0];
 265     for (int ix = 0; ix < nx; ix++) {
 266       for (int iy = 0; iy < ny; iy++) {
 267         if (v[ix][iy] > max)
 268           max = v[ix][iy];
 269         if (v[ix][iy] < min)
 270           min = v[ix][iy];
 271         mean += v[ix][iy];
 272       }
 273     }
 274     mean /= (nx * ny);
 275
 276     static const int nbin = 256;
 277     int hist[ nbin ] = {0};
 278     double spread = max - min;
 279     mode = 0;
 280     stddev = 0;
 281     for (int ix = 0; ix < nx; ix++) {
 282       for (int iy = 0; iy < ny; iy++) {
 283         int b = static_cast<int>((((v[ix][iy] - min) / spread) * (nbin - 1)) + 0.5);
 284         hist[b]++;
 285         double diff = (v[ix][iy] - mean);
 286         stddev += diff * diff;
 287         }
 288     }
 289     stddev = sqrt(stddev / (nx * ny));
 290
 291     int max_binindex = 0;
 292     int max_bin = -1;
 293     for (int ibin = 0; ibin < nbin; ibin++) {
 294       if (hist[ibin] > max_bin) {
 295         max_bin = hist[ibin];
 296         max_binindex = ibin;
 297       }
 298     }
 299
 300     mode = (max_binindex * spread / (nbin - 1)) + min;
 301
 302     median = 0.;
 303 }
 304
 305