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.25 2000/09/02 05:10:39 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_label.setLabelType (Array2dFileLabel::L_HISTORY);
  70   m_nView = nView;
  71   m_nDet = nDet;
  72   newProjData ();
  73 }
  74
  75 void
  76 Projections::initFromScanner (const Scanner& scanner)
  77 {
  78   m_label.setLabelType (Array2dFileLabel::L_HISTORY);
  79   deleteProjData();
  80   init (scanner.nView(), scanner.nDet());
  81
  82   m_phmLen = scanner.phmLen();
  83   m_rotInc = scanner.rotInc();
  84   m_detInc = scanner.detInc();
  85   m_geometry = scanner.geometry();
  86   m_focalLength = scanner.focalLength();
  87   m_fieldOfView = scanner.fieldOfView();
  88   m_rotStart = 0;
  89   m_detStart =  -(scanner.detLen() / 2);
  90 #if 0
  91   if (m_geometry == Scanner::GEOMETRY_EQUILINEAR) {
  92     m_detInc /= 2;
  93     cout << "Kludge: detInc /= 2 in Projections::initFromScanner" << endl;
  94   }
  95 #endif
  96 }
  97
  98 void
  99 Projections::setNView (int nView)  // used by MPI to reduce # of views
 100 {
 101   deleteProjData();
 102   init (nView, m_nDet);
 103 }
 104
 105 // NAME
 106 // newProjData
 107
 108 void
 109 Projections::newProjData (void)
 110 {
 111   if (m_projData)
 112     sys_error(ERR_WARNING, "m_projData != NULL [newProjData]");
 113
 114   if (m_nView > 0 && m_nDet) {
 115     m_projData = new DetectorArray* [m_nView];
 116
 117     for (int i = 0; i < m_nView; i++)
 118       m_projData[i] = new DetectorArray (m_nDet);
 119   }
 120 }
 121
 122
 123 /* NAME
 124  *   projections_free                   Free memory allocated to projections
 125  *
 126  * SYNOPSIS
 127  *   projections_free(proj)
 128  *   Projections& proj                          Projectionss to be deallocated
 129  */
 130
 131 void
 132 Projections::deleteProjData (void)
 133 {
 134   if (m_projData != NULL) {
 135     for (int i = 0; i < m_nView; i++)
 136       delete m_projData[i];
 137
 138     delete m_projData;
 139     m_projData = NULL;
 140   }
 141 }
 142
 143
 144 /* NAME
 145  *    Projections::headerWwrite         Write data header for projections file
 146  *
 147  */
 148
 149 bool
 150 Projections::headerWrite (fnetorderstream& fs)
 151 {
 152   kuint16 _hsize = m_headerSize;
 153   kuint16 _signature = m_signature;
 154   kuint32 _nView = m_nView;
 155   kuint32 _nDet = m_nDet;
 156   kuint32 _geom = m_geometry;
 157   kuint16 _remarksize = m_remark.length();
 158   kuint16 _year = m_year;
 159   kuint16 _month = m_month;
 160   kuint16 _day = m_day;
 161   kuint16 _hour = m_hour;
 162   kuint16 _minute = m_minute;
 163   kuint16 _second = m_second;
 164
 165   kfloat64 _calcTime = m_calcTime;
 166   kfloat64 _rotStart = m_rotStart;
 167   kfloat64 _rotInc = m_rotInc;
 168   kfloat64 _detStart = m_detStart;
 169   kfloat64 _detInc = m_detInc;
 170   kfloat64 _phmLen = m_phmLen;
 171   kfloat64 _fieldOfView = m_fieldOfView;
 172   kfloat64 _focalLength = m_focalLength;
 173
 174   fs.seekp(0);
 175   if (! fs)
 176     return false;
 177
 178   fs.writeInt16 (_hsize);
 179   fs.writeInt16 (_signature);
 180   fs.writeInt32 (_nView);
 181   fs.writeInt32 (_nDet);
 182   fs.writeInt32 (_geom);
 183   fs.writeFloat64 (_calcTime);
 184   fs.writeFloat64 (_rotStart);
 185   fs.writeFloat64 (_rotInc);
 186   fs.writeFloat64 (_detStart);
 187   fs.writeFloat64 (_detInc);
 188   fs.writeFloat64 (_phmLen);
 189   fs.writeFloat64 (_focalLength);
 190   fs.writeFloat64 (_fieldOfView);
 191   fs.writeInt16 (_year);
 192   fs.writeInt16 (_month);
 193   fs.writeInt16 (_day);
 194   fs.writeInt16 (_hour);
 195   fs.writeInt16 (_minute);
 196   fs.writeInt16 (_second);
 197   fs.writeInt16 (_remarksize);
 198   fs.write (m_remark.c_str(), _remarksize);
 199
 200   m_headerSize = fs.tellp();
 201   _hsize = m_headerSize;
 202   fs.seekp(0);
 203   fs.writeInt16 (_hsize);
 204   if (! fs)
 205       return false;
 206
 207   return true;
 208 }
 209
 210 /* NAME
 211  *    projections_read_header         Read data header for projections file
 212  *
 213  */
 214 bool
 215 Projections::headerRead (fnetorderstream& fs)
 216 {
 217   kuint16 _hsize, _signature, _year, _month, _day, _hour, _minute, _second, _remarksize = 0;
 218   kuint32 _nView, _nDet, _geom;
 219   kfloat64 _calcTime, _rotStart, _rotInc, _detStart, _detInc, _phmLen, _focalLength, _fieldOfView;
 220
 221   fs.seekg(0);
 222   if (! fs)
 223       return false;
 224
 225   fs.readInt16 (_hsize);
 226   fs.readInt16 (_signature);
 227   fs.readInt32 (_nView);
 228   fs.readInt32 (_nDet);
 229   fs.readInt32 (_geom);
 230   fs.readFloat64 (_calcTime);
 231   fs.readFloat64 (_rotStart);
 232   fs.readFloat64 (_rotInc);
 233   fs.readFloat64 (_detStart);
 234   fs.readFloat64 (_detInc);
 235   fs.readFloat64 (_phmLen);
 236   fs.readFloat64 (_focalLength);
 237   fs.readFloat64 (_fieldOfView);
 238   fs.readInt16 (_year);
 239   fs.readInt16 (_month);
 240   fs.readInt16 (_day);
 241   fs.readInt16 (_hour);
 242   fs.readInt16 (_minute);
 243   fs.readInt16 (_second);
 244   fs.readInt16 (_remarksize);
 245
 246   if (! fs) {
 247       sys_error (ERR_SEVERE, "Error reading header information , _remarksize=%d [projections_read_header]", _remarksize);
 248       return false;
 249   }
 250
 251   if (_signature != m_signature) {
 252     sys_error (ERR_SEVERE, "File %s does not have a valid projection file signature", m_filename.c_str());
 253     return false;
 254   }
 255
 256   char remarkStorage[_remarksize+1];
 257   fs.read (remarkStorage, _remarksize);
 258   if (! fs) {
 259     sys_error (ERR_SEVERE, "Error reading remark, _remarksize = %d", _remarksize);
 260     return false;
 261   }
 262   remarkStorage[_remarksize] = 0;
 263   m_remark = remarkStorage;
 264
 265   off_t _hsizeread = fs.tellg();
 266   if (!fs || _hsizeread != _hsize) {
 267       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);
 268       return false;
 269   }
 270
 271   m_headerSize = _hsize;
 272   m_nView = _nView;
 273   m_nDet = _nDet;
 274   m_geometry = _geom;
 275   m_calcTime = _calcTime;
 276   m_rotStart = _rotStart;
 277   m_rotInc = _rotInc;
 278   m_detStart = _detStart;
 279   m_detInc = _detInc;
 280   m_phmLen = _phmLen;
 281   m_focalLength = _focalLength;
 282   m_fieldOfView = _fieldOfView;
 283   m_year = _year;
 284   m_month = _month;
 285   m_day = _day;
 286   m_hour = _hour;
 287   m_minute = _minute;
 288   m_second = _second;
 289
 290   m_label.setLabelType (Array2dFileLabel::L_HISTORY);
 291   m_label.setLabelString (m_remark);
 292   m_label.setCalcTime (m_calcTime);
 293   m_label.setDateTime (m_year, m_month, m_day, m_hour, m_minute, m_second);
 294
 295   return true;
 296 }
 297
 298 bool
 299 Projections::read (const string& filename)
 300 {
 301   return read (filename.c_str());
 302 }
 303
 304
 305 bool
 306 Projections::read (const char* filename)
 307 {
 308   frnetorderstream fileRead (filename, ios::in | ios::binary);
 309   m_filename = filename;
 310
 311   if (! fileRead)
 312     return false;
 313
 314   if (! headerRead (fileRead))
 315     return false;
 316
 317   deleteProjData ();
 318   newProjData();
 319
 320   for (int i = 0; i < m_nView; i++) {
 321     if (! detarrayRead (fileRead, *m_projData[i], i))
 322       break;
 323   }
 324
 325   fileRead.close();
 326   return true;
 327 }
 328
 329
 330 bool
 331 Projections::write (const string& filename)
 332 {
 333   return write (filename.c_str());
 334 }
 335
 336 bool
 337 Projections::write (const char* filename)
 338 {
 339   frnetorderstream fs (filename, ios::out | ios::binary | ios::trunc | ios::ate);
 340   m_filename = filename;
 341   if (! fs) {
 342     sys_error (ERR_SEVERE, "Error opening file %s for output [projections_create]", filename);
 343     return false;
 344   }
 345
 346   time_t t = time(NULL);
 347   tm* lt = localtime(&t);
 348   m_year = lt->tm_year;
 349   m_month = lt->tm_mon;
 350   m_day = lt->tm_mday;
 351   m_hour = lt->tm_hour;
 352   m_minute = lt->tm_min;
 353   m_second = lt->tm_sec;
 354
 355   if (! headerWrite (fs))
 356       return false;
 357
 358   if (m_projData != NULL) {
 359     for (int i = 0; i < m_nView; i++) {
 360       if (! detarrayWrite (fs, *m_projData[i], i))
 361         break;
 362     }
 363   }
 364   if (! fs)
 365     return false;
 366
 367  fs.close();
 368
 369   return true;
 370 }
 371
 372 /* NAME
 373  *   detarrayRead               Read a Detector Array structure from the disk
 374  *
 375  * SYNOPSIS
 376  *   detarrayRead (proj, darray, view_num)
 377  *   DETARRAY *darray           Detector array storage location to be filled
 378  *   int      view_num          View number to read
 379  */
 380
 381 bool
 382 Projections::detarrayRead (fnetorderstream& fs, DetectorArray& darray, const int iview)
 383 {
 384   const int detval_bytes = darray.nDet() * sizeof(kfloat32);
 385   const int detheader_bytes = sizeof(kfloat64) /* view_angle */ + sizeof(kint32) /* nDet */;
 386   const int view_bytes = detheader_bytes + detval_bytes;
 387   const off_t start_data = m_headerSize + (iview * view_bytes);
 388   DetectorValue* detval_ptr = darray.detValues();
 389   kfloat64 view_angle;
 390   kuint32 nDet;
 391
 392   fs.seekg (start_data);
 393
 394   fs.readFloat64 (view_angle);
 395   fs.readInt32 (nDet);
 396   darray.setViewAngle (view_angle);
 397   //  darray.setNDet ( nDet);
 398
 399   for (unsigned int i = 0; i < nDet; i++) {
 400       kfloat32 detval;
 401       fs.readFloat32 (detval);
 402       detval_ptr[i] = detval;
 403   }
 404   if (! fs)
 405     return false;
 406
 407   return true;
 408 }
 409
 410
 411 /* NAME
 412  *   detarrayWrite                      Write detector array data to the disk
 413  *
 414  * SYNOPSIS
 415  *   detarrayWrite (darray, view_num)
 416  *   DETARRAY *darray                   Detector array data to be written
 417  *   int      view_num                  View number to write
 418  *
 419  * DESCRIPTION
 420  *       This routine writes the detarray data from the disk sequentially to
 421  *    the file that was opened with open_projections().  Data is written in
 422  *    binary format.
 423  */
 424
 425 bool
 426 Projections::detarrayWrite (fnetorderstream& fs, const DetectorArray& darray, const int iview)
 427 {
 428   const int detval_bytes = darray.nDet() * sizeof(float);
 429   const int detheader_bytes = sizeof(kfloat64) /* view_angle */ + sizeof(kint32) /* nDet */;
 430   const int view_bytes = detheader_bytes + detval_bytes;
 431   const off_t start_data = m_headerSize + (iview * view_bytes);
 432   const DetectorValue* const detval_ptr = darray.detValues();
 433   kfloat64 view_angle = darray.viewAngle();
 434   kuint32 nDet = darray.nDet();
 435
 436   fs.seekp (start_data);
 437   if (! fs) {
 438     sys_error (ERR_SEVERE, "Error seeking detectory array [detarrayWrite]");
 439     return false;
 440   }
 441
 442   fs.writeFloat64 (view_angle);
 443   fs.writeInt32 (nDet);
 444
 445   for (unsigned int i = 0; i < nDet; i++) {
 446     kfloat32 detval = detval_ptr[i];
 447     fs.writeFloat32 (detval);
 448   }
 449
 450   if (! fs)
 451     return (false);
 452
 453   return true;
 454 }
 455
 456 /* NAME
 457  *   printProjectionData                        Print projections data
 458  *
 459  * SYNOPSIS
 460  *   printProjectionData ()
 461  */
 462
 463 void
 464 Projections::printProjectionData (void)
 465 {
 466   printf("Projections Data\n\n");
 467   printf("Description: %s\n", m_remark.c_str());
 468   printf("Geometry: %s\n", Scanner::convertGeometryIDToName (m_geometry));
 469   printf("nView       = %8d           nDet = %8d\n", m_nView, m_nDet);
 470   printf("focalLength = %8.4f  fieldOfView = %8.4f\n", m_focalLength, m_fieldOfView);
 471   printf("rotStart    = %8.4f       rotInc = %8.4f\n", m_rotStart, m_rotInc);
 472   printf("detStart    = %8.4f       detInc = %8.4f\n", m_detStart, m_detInc);
 473   if (m_projData != NULL) {
 474       for (int ir = 0; ir < m_nView; ir++) {
 475           printf("View %d: angle %f\n", ir, m_projData[ir]->viewAngle());
 476           DetectorValue* detval = m_projData[ir]->detValues();
 477           for (int id = 0; id < m_projData[ir]->nDet(); id++)
 478               printf("%8.4f  ", detval[id]);
 479           printf("\n");
 480       }
 481   }
 482 }
 483
 484 void
 485 Projections::printScanInfo (ostringstream& os) const
 486 {
 487   os << "Number of detectors: " << m_nDet << "\n";
 488   os << "    Number of views: " << m_nView<< "\n";
 489   os << "             Remark: " << m_remark.c_str()<< "\n";
 490   os << "           Geometry: " << Scanner::convertGeometryIDToName (m_geometry)<< "\n";
 491   os << "       Focal Length: " << m_focalLength<< "\n";
 492   os << "      Field Of View: " << m_fieldOfView<< "\n";
 493   os << "             phmLen: " << m_phmLen<< "\n";
 494   os << "           detStart: " << m_detStart<< "\n";
 495   os << "             detInc: " << m_detInc<< "\n";
 496   os << "           rotStart: " << m_rotStart<< "\n";
 497   os << "             rotInc: " << m_rotInc<< "\n";
 498 }
 499
 500
 501