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