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