libctsim/projections.cpp

   1 /*****************************************************************************
   2 ** FILE IDENTIFICATION
   3 **
   4 **   Name:         projections.cpp         Projection data classes
   5 **   Programmer:   Kevin Rosenberg
   6 **   Date Started: Aug 84
   7 **
   8 **  This is part of the CTSim program
   9 **  Copyright (C) 1983-2000 Kevin Rosenberg
  10 **
  11 **  $Id: projections.cpp,v 1.23 2000/08/27 20:32:55 kevin Exp $
  12 **
  13 **  This program is free software; you can redistribute it and/or modify
  14 **  it under the terms of the GNU General Public License (version 2) as
  15 **  published by the Free Software Foundation.
  16 **
  17 **  This program is distributed in the hope that it will be useful,
  18 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
  19 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  20 **  GNU General Public License for more details.
  21 **
  22 **  You should have received a copy of the GNU General Public License
  23 **  along with this program; if not, write to the Free Software
  24 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  25 ******************************************************************************/
  26
  27 #include "ct.h"
  28
  29
  30 /* NAME
  31  *    Projections               Constructor for projections matrix storage
  32  *
  33  * SYNOPSIS
  34  *    proj = projections_create (filename, nView, nDet)
  35  *    Projections& proj         Allocated projections structure & matrix
  36  *    int nView                 Number of rotated view
  37  *    int nDet                  Number of detectors
  38  *
  39  */
  40
  41 Projections::Projections (const Scanner& scanner)
  42   : m_projData(0)
  43 {
  44   initFromScanner (scanner);
  45 }
  46
  47
  48 Projections::Projections (const int nView, const int nDet)
  49   : m_projData(0)
  50 {
  51   init (nView, nDet);
  52 }
  53
  54 Projections::Projections (void)
  55   : m_projData(0)
  56 {
  57   init (0, 0);
  58 }
  59
  60 Projections::~Projections (void)
  61 {
  62   deleteProjData();
  63 }
  64
  65
  66 void
  67 Projections::init (const int nView, const int nDet)
  68 {
  69   m_nView = nView;
  70   m_nDet = nDet;
  71   newProjData ();
  72 }
  73
  74 void
  75 Projections::initFromScanner (const Scanner& scanner)
  76 {
  77   deleteProjData();
  78   init (scanner.nView(), scanner.nDet());
  79
  80   m_phmLen = scanner.phmLen();
  81   m_rotInc = scanner.rotInc();
  82   m_detInc = scanner.detInc();
  83   m_focalLength = scanner.focalLength();
  84   m_fieldOfView = scanner.fieldOfView();
  85   m_rotStart = 0;
  86   m_detStart =  -(scanner.detLen() / 2) + (scanner.detInc() / 2);
  87 }
  88
  89 void
  90 Projections::setNView (int nView)  // used by MPI to reduce # of views
  91 {
  92   deleteProjData();
  93   init (nView, m_nDet);
  94 }
  95
  96 // NAME
  97 // newProjData
  98
  99 void
 100 Projections::newProjData (void)
 101 {
 102   if (m_projData)
 103     sys_error(ERR_WARNING, "m_projData != NULL [newProjData]");
 104
 105   if (m_nView > 0 && m_nDet) {
 106     m_projData = new DetectorArray* [m_nView];
 107
 108     for (int i = 0; i < m_nView; i++)
 109       m_projData[i] = new DetectorArray (m_nDet);
 110   }
 111 }
 112
 113
 114 /* NAME
 115  *   projections_free                   Free memory allocated to projections
 116  *
 117  * SYNOPSIS
 118  *   projections_free(proj)
 119  *   Projections& proj                          Projectionss to be deallocated
 120  */
 121
 122 void
 123 Projections::deleteProjData (void)
 124 {
 125   if (m_projData != NULL) {
 126     for (int i = 0; i < m_nView; i++)
 127       delete m_projData[i];
 128
 129     delete m_projData;
 130     m_projData = NULL;
 131   }
 132 }
 133
 134
 135 /* NAME
 136  *    Projections::headerWwrite         Write data header for projections file
 137  *
 138  */
 139
 140 bool
 141 Projections::headerWrite (fnetorderstream& fs)
 142 {
 143   kuint16 _hsize = m_headerSize;
 144   kuint16 _signature = m_signature;
 145   kuint32 _nView = m_nView;
 146   kuint32 _nDet = m_nDet;
 147   kuint32 _geom = m_geometry;
 148   kuint16 _remarksize = m_remark.length();
 149   kuint16 _year = m_year;
 150   kuint16 _month = m_month;
 151   kuint16 _day = m_day;
 152   kuint16 _hour = m_hour;
 153   kuint16 _minute = m_minute;
 154   kuint16 _second = m_second;
 155
 156   kfloat64 _calcTime = m_calcTime;
 157   kfloat64 _rotStart = m_rotStart;
 158   kfloat64 _rotInc = m_rotInc;
 159   kfloat64 _detStart = m_detStart;
 160   kfloat64 _detInc = m_detInc;
 161   kfloat64 _phmLen = m_phmLen;
 162   kfloat64 _fieldOfView = m_fieldOfView;
 163   kfloat64 _focalLength = m_focalLength;
 164
 165   fs.seekp(0);
 166   if (! fs)
 167     return false;
 168
 169   fs.writeInt16 (_hsize);
 170   fs.writeInt16 (_signature);
 171   fs.writeInt32 (_nView);
 172   fs.writeInt32 (_nDet);
 173   fs.writeInt32 (_geom);
 174   fs.writeFloat64 (_calcTime);
 175   fs.writeFloat64 (_rotStart);
 176   fs.writeFloat64 (_rotInc);
 177   fs.writeFloat64 (_detStart);
 178   fs.writeFloat64 (_detInc);
 179   fs.writeFloat64 (_phmLen);
 180   fs.writeFloat64 (_focalLength);
 181   fs.writeFloat64 (_fieldOfView);
 182   fs.writeInt16 (_year);
 183   fs.writeInt16 (_month);
 184   fs.writeInt16 (_day);
 185   fs.writeInt16 (_hour);
 186   fs.writeInt16 (_minute);
 187   fs.writeInt16 (_second);
 188   fs.writeInt16 (_remarksize);
 189   fs.write (m_remark.c_str(), _remarksize);
 190
 191   m_headerSize = fs.tellp();
 192   _hsize = m_headerSize;
 193   fs.seekp(0);
 194   fs.writeInt16 (_hsize);
 195   if (! fs)
 196       return false;
 197
 198   return true;
 199 }
 200
 201 /* NAME
 202  *    projections_read_header         Read data header for projections file
 203  *
 204  */
 205 bool
 206 Projections::headerRead (fnetorderstream& fs)
 207 {
 208   kuint16 _hsize, _signature, _year, _month, _day, _hour, _minute, _second, _remarksize = 0;
 209   kuint32 _nView, _nDet, _geom;
 210   kfloat64 _calcTime, _rotStart, _rotInc, _detStart, _detInc, _phmLen, _focalLength, _fieldOfView;
 211
 212   fs.seekg(0);
 213   if (! fs)
 214       return false;
 215
 216   fs.readInt16 (_hsize);
 217   fs.readInt16 (_signature);
 218   fs.readInt32 (_nView);
 219   fs.readInt32 (_nDet);
 220   fs.readInt32 (_geom);
 221   fs.readFloat64 (_calcTime);
 222   fs.readFloat64 (_rotStart);
 223   fs.readFloat64 (_rotInc);
 224   fs.readFloat64 (_detStart);
 225   fs.readFloat64 (_detInc);
 226   fs.readFloat64 (_phmLen);
 227   fs.readFloat64 (_focalLength);
 228   fs.readFloat64 (_fieldOfView);
 229   fs.readInt16 (_year);
 230   fs.readInt16 (_month);
 231   fs.readInt16 (_day);
 232   fs.readInt16 (_hour);
 233   fs.readInt16 (_minute);
 234   fs.readInt16 (_second);
 235   fs.readInt16 (_remarksize);
 236
 237   if (! fs) {
 238       sys_error (ERR_SEVERE, "Error reading header information , _remarksize=%d [projections_read_header]", _remarksize);
 239       return false;
 240   }
 241
 242   if (_signature != m_signature) {
 243     sys_error (ERR_SEVERE, "File %s does not have a valid projection file signature", m_filename.c_str());
 244     return false;
 245   }
 246
 247   char remarkStorage[_remarksize+1];
 248   fs.read (remarkStorage, _remarksize);
 249   if (! fs) {
 250     sys_error (ERR_SEVERE, "Error reading remark, _remarksize = %d", _remarksize);
 251     return false;
 252   }
 253   remarkStorage[_remarksize] = 0;
 254   m_remark = remarkStorage;
 255
 256   off_t _hsizeread = fs.tellg();
 257   if (!fs || _hsizeread != _hsize) {
 258       sys_error (ERR_WARNING, "File header size read %ld != file header size stored %ld [read_projections_header]\n_remarksize=%ld", (long int) _hsizeread, _hsize, _remarksize);
 259       return false;
 260   }
 261
 262   m_headerSize = _hsize;
 263   m_nView = _nView;
 264   m_nDet = _nDet;
 265   m_geometry = _geom;
 266   m_calcTime = _calcTime;
 267   m_rotStart = _rotStart;
 268   m_rotInc = _rotInc;
 269   m_detStart = _detStart;
 270   m_detInc = _detInc;
 271   m_phmLen = _phmLen;
 272   m_focalLength = _focalLength;
 273   m_fieldOfView = _fieldOfView;
 274   m_year = _year;
 275   m_month = _month;
 276   m_day = _day;
 277   m_hour = _hour;
 278   m_minute = _minute;
 279   m_second = _second;
 280
 281   m_label.setLabelType (Array2dFileLabel::L_HISTORY);
 282   m_label.setLabelString (m_remark);
 283   m_label.setCalcTime (m_calcTime);
 284   m_label.setDateTime (m_year, m_month, m_day, m_hour, m_minute, m_second);
 285
 286   return true;
 287 }
 288
 289 bool
 290 Projections::read (const char* filename)
 291 {
 292   frnetorderstream fileRead (filename, ios::in | ios::binary);
 293   m_filename = filename;
 294
 295   if (! fileRead)
 296     return false;
 297
 298   if (! headerRead (fileRead))
 299     return false;
 300
 301   deleteProjData ();
 302   newProjData();
 303
 304   for (int i = 0; i < m_nView; i++) {
 305     if (! detarrayRead (fileRead, *m_projData[i], i))
 306       break;
 307   }
 308
 309   fileRead.close();
 310   return true;
 311 }
 312
 313
 314 bool
 315 Projections::write (const char* filename)
 316 {
 317   frnetorderstream fs (filename, ios::out | ios::binary | ios::trunc | ios::ate);
 318   m_filename = filename;
 319   if (! fs) {
 320     sys_error (ERR_SEVERE, "Error opening file %s for output [projections_create]", filename);
 321     return false;
 322   }
 323
 324   time_t t = time(NULL);
 325   tm* lt = localtime(&t);
 326   m_year = lt->tm_year;
 327   m_month = lt->tm_mon;
 328   m_day = lt->tm_mday;
 329   m_hour = lt->tm_hour;
 330   m_minute = lt->tm_min;
 331   m_second = lt->tm_sec;
 332
 333   if (! headerWrite (fs))
 334       return false;
 335
 336   if (m_projData != NULL) {
 337     for (int i = 0; i < m_nView; i++) {
 338       if (! detarrayWrite (fs, *m_projData[i], i))
 339         break;
 340     }
 341   }
 342   if (! fs)
 343     return false;
 344
 345  fs.close();
 346
 347   return true;
 348 }
 349
 350 /* NAME
 351  *   detarrayRead               Read a Detector Array structure from the disk
 352  *
 353  * SYNOPSIS
 354  *   detarrayRead (proj, darray, view_num)
 355  *   DETARRAY *darray           Detector array storage location to be filled
 356  *   int      view_num          View number to read
 357  */
 358
 359 bool
 360 Projections::detarrayRead (fnetorderstream& fs, DetectorArray& darray, const int iview)
 361 {
 362   const int detval_bytes = darray.nDet() * sizeof(kfloat32);
 363   const int detheader_bytes = sizeof(kfloat64) /* view_angle */ + sizeof(kint32) /* nDet */;
 364   const int view_bytes = detheader_bytes + detval_bytes;
 365   const off_t start_data = m_headerSize + (iview * view_bytes);
 366   DetectorValue* detval_ptr = darray.detValues();
 367   kfloat64 view_angle;
 368   kuint32 nDet;
 369
 370   fs.seekg (start_data);
 371
 372   fs.readFloat64 (view_angle);
 373   fs.readInt32 (nDet);
 374   darray.setViewAngle (view_angle);
 375   //  darray.setNDet ( nDet);
 376
 377   for (unsigned int i = 0; i < nDet; i++) {
 378       kfloat32 detval;
 379       fs.readFloat32 (detval);
 380       detval_ptr[i] = detval;
 381   }
 382   if (! fs)
 383     return false;
 384
 385   return true;
 386 }
 387
 388
 389 /* NAME
 390  *   detarrayWrite                      Write detector array data to the disk
 391  *
 392  * SYNOPSIS
 393  *   detarrayWrite (darray, view_num)
 394  *   DETARRAY *darray                   Detector array data to be written
 395  *   int      view_num                  View number to write
 396  *
 397  * DESCRIPTION
 398  *       This routine writes the detarray data from the disk sequentially to
 399  *    the file that was opened with open_projections().  Data is written in
 400  *    binary format.
 401  */
 402
 403 bool
 404 Projections::detarrayWrite (fnetorderstream& fs, const DetectorArray& darray, const int iview)
 405 {
 406   const int detval_bytes = darray.nDet() * sizeof(float);
 407   const int detheader_bytes = sizeof(kfloat64) /* view_angle */ + sizeof(kint32) /* nDet */;
 408   const int view_bytes = detheader_bytes + detval_bytes;
 409   const off_t start_data = m_headerSize + (iview * view_bytes);
 410   const DetectorValue* const detval_ptr = darray.detValues();
 411   kfloat64 view_angle = darray.viewAngle();
 412   kuint32 nDet = darray.nDet();
 413
 414   fs.seekp (start_data);
 415   if (! fs) {
 416     sys_error (ERR_SEVERE, "Error seeking detectory array [detarrayWrite]");
 417     return false;
 418   }
 419
 420   fs.writeFloat64 (view_angle);
 421   fs.writeInt32 (nDet);
 422
 423   for (unsigned int i = 0; i < nDet; i++) {
 424     kfloat32 detval = detval_ptr[i];
 425     fs.writeFloat32 (detval);
 426   }
 427
 428   if (! fs)
 429     return (false);
 430
 431   return true;
 432 }
 433
 434 /* NAME
 435  *   prt_projections                    Print projections data
 436  *
 437  * SYNOPSIS
 438  *   prt_projections (proj)
 439  *   Projections& proj                  Projection data to be printed
 440  */
 441
 442 void
 443 Projections::printProjectionData (void)
 444 {
 445   printf("Projections Print\n\n");
 446   printf("Description: %s\n", m_remark.c_str());
 447   printf("nView = %8d          nDet = %8d\n", m_nView, m_nDet);
 448   printf("rotStart = %8.4f     rotInc = %8.4f\n", m_rotStart, m_rotInc);
 449   printf("detStart = %8.4f     detInc = %8.4f\n", m_detStart, m_detInc);
 450   printf("focalLength = %8.4f  fieldOfView = %8.4f\n", m_focalLength, m_fieldOfView);
 451   if (m_projData != NULL) {
 452       for (int ir = 0; ir < m_nView; ir++) {
 453           DetectorValue* detval = m_projData[ir]->detValues();
 454           for (int id = 0; id < m_projData[ir]->nDet(); id++)
 455               printf("%8.4f  ", detval[id]);
 456           printf("\n");
 457       }
 458   }
 459 }
 460
 461 void
 462 Projections::printScanInfo (void) const
 463 {
 464   printf ("Number of detectors: %d\n", m_nDet);
 465   printf ("    Number of views: %d\n", m_nView);
 466   printf ("             Remark: %s\n", m_remark.c_str());
 467   printf ("             phmLen: %f\n", m_phmLen);
 468   printf ("           detStart: %f\n", m_detStart);
 469   printf ("             detInc: %f\n", m_detInc);
 470   printf ("           rotStart: %f\n", m_rotStart);
 471   printf ("             rotInc: %f\n", m_rotInc);
 472 }
 473
 474
 475
 476 /* NAME
 477  *   Projections::reconstruct      Reconstruct Image from Projections
 478  *
 479  * SYNOPSIS
 480  *   im = proj.reconstruct (im, filt_type, filt_param, interp_type)
 481  *   IMAGE *im                  Output image
 482  *   int filt_type              Type of convolution filter to use
 483  *   double filt_param          Filter specific parameter
 484  *                              Currently, used only with Hamming filters
 485  *   int interp_type            Type of interpolation method to use
 486  *
 487  * ALGORITHM
 488  *
 489  *      Calculate one-dimensional filter in spatial domain
 490  *      Allocate & clear (zero) the 2d output image array
 491  *      For each projection view
 492  *          Convolve raysum array with filter
 493  *          Backproject raysums and add (summate) to image array
 494  *      end
 495  */
 496
 497 bool
 498 Projections::reconstruct (ImageFile& im, const char* const filterName, double filt_param, const char* const filterMethodName, const int zeropad, const char* filterGenerationName, const char* const interpName, int interpFactor, const char* const backprojectName, const int trace) const
 499 {
 500   double detInc = m_detInc;
 501   int n_filteredProj = m_nDet * interpFactor;
 502   double filteredProj [n_filteredProj];   // filtered projections
 503
 504 #ifdef HAVE_BSPLINE_INTERP
 505   int spline_order = 0, zoom_factor = 0;
 506   if (interp_type == I_BSPLINE) {
 507     zoom_factor = interpFactor;
 508     spline_order = 3;
 509     zoom_factor = 3;
 510     n_filteredProj = (m_nDet - 1) * (zoom_factor + 1) + 1;
 511   }
 512 #endif
 513
 514   double filterBW = 1. / detInc;
 515   ProcessSignal processSignal (filterName, filterMethodName, filterBW, m_detInc, m_nDet, filt_param, "spatial", filterGenerationName, zeropad, interpFactor, trace, m_geometry);
 516
 517   if (processSignal.fail()) {
 518       sys_error (ERR_SEVERE, "%s [Projections::reconstruct]", processSignal.failMessage().c_str());
 519       return false;
 520   }
 521
 522   if (trace)
 523     cout << "Reconstruct: filter="<<filterName<< ", interp="<<interpName<<", backproject="<<backprojectName<<endl;
 524
 525 #if HAVE_SGP
 526   int nVecFilter = processSignal.getNFilterPoints();
 527   double plot_xaxis [nVecFilter];                       // array for plotting
 528
 529   if (trace > Trace::TRACE_CONSOLE && nVecFilter > 0)  {
 530     int i;
 531     double f;
 532     double filterInc = processSignal.getFilterIncrement();
 533     for (i = 0, f = processSignal.getFilterMin(); i < nVecFilter; i++, f += filterInc)
 534       plot_xaxis[i] = f;
 535
 536     if (processSignal.getFilter()) {
 537       SGPDriver sgpDriver ("Filter Function");
 538       SGP sgp (sgpDriver);
 539       EZPlot ezplot (sgp);
 540
 541       ezplot.ezset ("title Filter Response");
 542       ezplot.addCurve (plot_xaxis, processSignal.getFilter(), nVecFilter);
 543       ezplot.plot();
 544       cio_put_str ("Press any key to continue");
 545       cio_kb_getc ();
 546     }
 547   }
 548   if (trace >= Trace::TRACE_CONSOLE) {
 549     printf ("nview=%d, ndet=%d, det_start=%.4f, detInc=%.4f\n", m_nView, m_nDet, m_detStart, m_detInc);
 550   }
 551 #endif  //HAVE_SGP
 552
 553   Backprojector bj (*this, im, backprojectName, interpName, interpFactor);
 554   if (bj.fail()) {
 555     sys_error (ERR_SEVERE, "%s [Projections::reconstruct]", bj.failMessage().c_str());
 556     return false;
 557   }
 558
 559   for (int iview = 0; iview < m_nView; iview++)  {
 560     if (trace >= Trace::TRACE_CONSOLE)
 561       printf ("Reconstructing view %d (last = %d)\n",  iview, m_nView - 1);
 562
 563     const DetectorArray& darray = getDetectorArray (iview);
 564     const DetectorValue* detval = darray.detValues();
 565
 566     processSignal.filterSignal (detval, filteredProj, iview);
 567
 568 #ifdef HAVE_BSPLINE_INTERP
 569     if (interp_type == I_BSPLINE)
 570         bspline (m_nDet, zoom_factor, spline_order, filteredProj, filteredProj);
 571
 572 #ifdef HAVE_SGP
 573     if (trace >= Trace::TRACE_PLOT && interp_type == I_BSPLINE) {
 574         bspline (m_nDet, zoom_factor, spline_order, filteredProj, filteredProj);
 575       ezplot_1d (filteredProj, n_filteredProj);
 576     }
 577 #endif
 578 #endif
 579
 580     bj.BackprojectView (filteredProj, darray.viewAngle());
 581
 582 #ifdef HAVE_SGP
 583     if (trace >= Trace::TRACE_PLOT) {
 584       SGPDriver sgpDriverProj ("Projection");
 585       SGP sgpProj (sgpDriverProj);
 586       EZPlot ezplotProj (sgpProj);
 587
 588       ezplotProj.ezset  ("clear");
 589       ezplotProj.ezset ("title Filtered Projection");
 590       ezplotProj.ezset  ("xticks major 5.");
 591       ezplotProj.ezset  ("xlabel ");
 592       ezplotProj.ezset  ("ylabel ");
 593       ezplotProj.ezset  ("yporigin .5.");
 594       ezplotProj.ezset  ("ylength .5.");
 595       ezplotProj.ezset  ("box.");
 596       ezplotProj.ezset  ("grid.");
 597       ezplotProj.addCurve (plot_xaxis, detval, m_nDet);
 598       ezplotProj.plot();
 599       ezplotProj.ezset  ("clear");
 600       ezplotProj.ezset  ("xticks major 5.");
 601       ezplotProj.ezset  ("xlabel ");
 602       ezplotProj.ezset  ("ylabel ");
 603       ezplotProj.ezset  ("ylength .5.");
 604       ezplotProj.ezset ("box");
 605       ezplotProj.ezset ("grid");
 606       ezplotProj.addCurve (plot_xaxis, filteredProj,  n_filteredProj);
 607       ezplotProj.plot();
 608
 609       cout << "Do you want to exit with current pic (y/n)? " << flush;
 610       char str[256];
 611       fgets(str, sizeof(str), stdin);
 612       if (tolower(str[0]) == 'y') {
 613         break;
 614       }
 615     }
 616 #endif  //HAVE_SGP
 617   }
 618
 619   return true;
 620 }
 621