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