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.31 2000/12/16 02:54:47 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 \r
  29 const kuint16 Projections::m_signature = ('P'*256 + 'J');\r
  30
  31 /* NAME
  32  *    Projections               Constructor for projections matrix storage
  33  *
  34  * SYNOPSIS
  35  *    proj = projections_create (filename, nView, nDet)
  36  *    Projections& proj         Allocated projections structure & matrix
  37  *    int nView                 Number of rotated view
  38  *    int nDet                  Number of detectors
  39  *
  40  */
  41
  42 Projections::Projections (const Scanner& scanner)
  43   : m_projData(0)
  44 {
  45   initFromScanner (scanner);
  46 }
  47
  48
  49 Projections::Projections (const int nView, const int nDet)
  50   : m_projData(0)
  51 {
  52   init (nView, nDet);
  53 }
  54
  55 Projections::Projections (void)
  56   : m_projData(0)
  57 {
  58   init (0, 0);
  59 }
  60
  61 Projections::~Projections (void)
  62 {
  63   deleteProjData();
  64 }
  65
  66
  67 void
  68 Projections::init (const int nView, const int nDet)
  69 {
  70   m_label.setLabelType (Array2dFileLabel::L_HISTORY);
  71   m_nView = nView;
  72   m_nDet = nDet;
  73   newProjData ();
  74 }
  75
  76 void
  77 Projections::initFromScanner (const Scanner& scanner)
  78 {
  79   m_label.setLabelType (Array2dFileLabel::L_HISTORY);
  80   deleteProjData();
  81   init (scanner.nView(), scanner.nDet());
  82
  83   m_phmLen = scanner.phmLen();
  84   m_rotInc = scanner.rotInc();
  85   m_detInc = scanner.detInc();
  86   m_geometry = scanner.geometry();
  87   m_focalLength = scanner.focalLength();
  88   m_fieldOfView = scanner.fieldOfView();
  89   m_rotStart = 0;
  90   m_detStart =  -(scanner.detLen() / 2);
  91 #if 0
  92   if (m_geometry == Scanner::GEOMETRY_EQUILINEAR) {
  93     m_detInc /= 2;
  94     cout << "Kludge: detInc /= 2 in Projections::initFromScanner" << endl;
  95   }
  96 #endif
  97 }
  98
  99 void
 100 Projections::setNView (int nView)  // used by MPI to reduce # of views
 101 {
 102   deleteProjData();
 103   init (nView, m_nDet);
 104 }
 105
 106 // NAME
 107 // newProjData
 108
 109 void
 110 Projections::newProjData (void)
 111 {
 112   if (m_projData)
 113     sys_error(ERR_WARNING, "m_projData != NULL [newProjData]");
 114
 115   if (m_nView > 0 && m_nDet) {
 116     m_projData = new DetectorArray* [m_nView];
 117
 118     for (int i = 0; i < m_nView; i++)
 119       m_projData[i] = new DetectorArray (m_nDet);
 120   }
 121 }
 122
 123
 124 /* NAME
 125  *   projections_free                   Free memory allocated to projections
 126  *
 127  * SYNOPSIS
 128  *   projections_free(proj)
 129  *   Projections& proj                          Projectionss to be deallocated
 130  */
 131
 132 void
 133 Projections::deleteProjData (void)
 134 {
 135   if (m_projData != NULL) {
 136     for (int i = 0; i < m_nView; i++)
 137       delete m_projData[i];
 138
 139     delete m_projData;
 140     m_projData = NULL;
 141   }
 142 }
 143
 144
 145 /* NAME
 146  *    Projections::headerWwrite         Write data header for projections file
 147  *
 148  */
 149
 150 bool
 151 Projections::headerWrite (fnetorderstream& fs)
 152 {
 153   kuint16 _hsize = m_headerSize;
 154   kuint16 _signature = m_signature;
 155   kuint32 _nView = m_nView;
 156   kuint32 _nDet = m_nDet;
 157   kuint32 _geom = m_geometry;
 158   kuint16 _remarksize = m_remark.length();
 159   kuint16 _year = m_year;
 160   kuint16 _month = m_month;
 161   kuint16 _day = m_day;
 162   kuint16 _hour = m_hour;
 163   kuint16 _minute = m_minute;
 164   kuint16 _second = m_second;
 165
 166   kfloat64 _calcTime = m_calcTime;
 167   kfloat64 _rotStart = m_rotStart;
 168   kfloat64 _rotInc = m_rotInc;
 169   kfloat64 _detStart = m_detStart;
 170   kfloat64 _detInc = m_detInc;
 171   kfloat64 _phmLen = m_phmLen;
 172   kfloat64 _fieldOfView = m_fieldOfView;
 173   kfloat64 _focalLength = m_focalLength;
 174
 175   fs.seekp(0);
 176   if (! fs)
 177     return false;
 178
 179   fs.writeInt16 (_hsize);
 180   fs.writeInt16 (_signature);
 181   fs.writeInt32 (_nView);
 182   fs.writeInt32 (_nDet);
 183   fs.writeInt32 (_geom);
 184   fs.writeFloat64 (_calcTime);
 185   fs.writeFloat64 (_rotStart);
 186   fs.writeFloat64 (_rotInc);
 187   fs.writeFloat64 (_detStart);
 188   fs.writeFloat64 (_detInc);
 189   fs.writeFloat64 (_phmLen);
 190   fs.writeFloat64 (_focalLength);
 191   fs.writeFloat64 (_fieldOfView);
 192   fs.writeInt16 (_year);
 193   fs.writeInt16 (_month);
 194   fs.writeInt16 (_day);
 195   fs.writeInt16 (_hour);
 196   fs.writeInt16 (_minute);
 197   fs.writeInt16 (_second);
 198   fs.writeInt16 (_remarksize);
 199   fs.write (m_remark.c_str(), _remarksize);
 200
 201   m_headerSize = fs.tellp();
 202   _hsize = m_headerSize;
 203   fs.seekp(0);
 204   fs.writeInt16 (_hsize);
 205   if (! fs)
 206       return false;
 207
 208   return true;
 209 }
 210
 211 /* NAME
 212  *    projections_read_header         Read data header for projections file
 213  *
 214  */
 215 bool
 216 Projections::headerRead (fnetorderstream& fs)
 217 {
 218   kuint16 _hsize, _signature, _year, _month, _day, _hour, _minute, _second, _remarksize = 0;
 219   kuint32 _nView, _nDet, _geom;
 220   kfloat64 _calcTime, _rotStart, _rotInc, _detStart, _detInc, _phmLen, _focalLength, _fieldOfView;
 221
 222   fs.seekg(0);
 223   if (! fs)
 224       return false;
 225
 226   fs.readInt16 (_hsize);
 227   fs.readInt16 (_signature);
 228   fs.readInt32 (_nView);
 229   fs.readInt32 (_nDet);
 230   fs.readInt32 (_geom);
 231   fs.readFloat64 (_calcTime);
 232   fs.readFloat64 (_rotStart);
 233   fs.readFloat64 (_rotInc);
 234   fs.readFloat64 (_detStart);
 235   fs.readFloat64 (_detInc);
 236   fs.readFloat64 (_phmLen);
 237   fs.readFloat64 (_focalLength);
 238   fs.readFloat64 (_fieldOfView);
 239   fs.readInt16 (_year);
 240   fs.readInt16 (_month);
 241   fs.readInt16 (_day);
 242   fs.readInt16 (_hour);
 243   fs.readInt16 (_minute);
 244   fs.readInt16 (_second);
 245   fs.readInt16 (_remarksize);
 246
 247   if (! fs) {
 248       sys_error (ERR_SEVERE, "Error reading header information , _remarksize=%d [projections_read_header]", _remarksize);
 249       return false;
 250   }
 251
 252   if (_signature != m_signature) {
 253     sys_error (ERR_SEVERE, "File %s does not have a valid projection file signature", m_filename.c_str());
 254     return false;
 255   }
 256
 257   char* pszRemarkStorage = new char [_remarksize+1];
 258   fs.read (pszRemarkStorage, _remarksize);
 259   if (! fs) {
 260     sys_error (ERR_SEVERE, "Error reading remark, _remarksize = %d", _remarksize);
 261     return false;
 262   }
 263   pszRemarkStorage[_remarksize] = 0;
 264   m_remark = pszRemarkStorage;
 265   delete pszRemarkStorage;\r
 266
 267   off_t _hsizeread = fs.tellg();
 268   if (!fs || _hsizeread != _hsize) {
 269       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);
 270       return false;
 271   }
 272
 273   m_headerSize = _hsize;
 274   m_nView = _nView;
 275   m_nDet = _nDet;
 276   m_geometry = _geom;
 277   m_calcTime = _calcTime;
 278   m_rotStart = _rotStart;
 279   m_rotInc = _rotInc;
 280   m_detStart = _detStart;
 281   m_detInc = _detInc;
 282   m_phmLen = _phmLen;
 283   m_focalLength = _focalLength;
 284   m_fieldOfView = _fieldOfView;
 285   m_year = _year;
 286   m_month = _month;
 287   m_day = _day;
 288   m_hour = _hour;
 289   m_minute = _minute;
 290   m_second = _second;
 291
 292   m_label.setLabelType (Array2dFileLabel::L_HISTORY);
 293   m_label.setLabelString (m_remark);
 294   m_label.setCalcTime (m_calcTime);
 295   m_label.setDateTime (m_year, m_month, m_day, m_hour, m_minute, m_second);
 296
 297   return true;
 298 }
 299
 300 bool
 301 Projections::read (const string& filename)
 302 {
 303   return read (filename.c_str());
 304 }
 305
 306
 307 bool
 308 Projections::read (const char* filename)
 309 {
 310   frnetorderstream fileRead (filename, ios::in | ios::binary);
 311   m_filename = filename;
 312
 313   if (! fileRead)
 314     return false;
 315
 316   if (! headerRead (fileRead))
 317     return false;
 318
 319   deleteProjData ();
 320   newProjData();
 321
 322   for (int i = 0; i < m_nView; i++) {
 323     if (! detarrayRead (fileRead, *m_projData[i], i))
 324       break;
 325   }
 326
 327   fileRead.close();
 328   return true;
 329 }
 330
 331
 332 bool
 333 Projections::copyViewData (const string& filename, ostream& os, int startView, int endView)
 334 {
 335     return copyViewData (filename.c_str(), os, startView, endView);
 336 }
 337
 338 bool
 339 Projections::copyViewData (const char* const filename, ostream& os, int startView, int endView)
 340 {
 341   frnetorderstream is (filename, ios::in | ios::binary);
 342   kuint16 sizeHeader, signature;
 343   kuint32 _nView, _nDet;
 344
 345   is.readInt16 (sizeHeader);
 346   is.readInt16 (signature);
 347   is.readInt32 (_nView);
 348   is.readInt32 (_nDet);
 349   int nView = _nView;
 350   int nDet = _nDet;
 351
 352   if (signature != m_signature) {
 353     sys_error (ERR_FATAL, "Illegal signature in projection file %s", filename);
 354     return false;
 355   }
 356
 357   if (startView < 0)
 358       startView = 0;
 359   if (startView > nView - 1)
 360       startView = nView;
 361   if (endView < 0 || endView > nView - 1)
 362       endView = nView - 1;
 363
 364   if (startView > endView) { // swap if start > end
 365       int tempView = endView;
 366       endView = startView;
 367       startView = tempView;
 368   }
 369
 370   int sizeView = 8 /* view_angle */ + 4 /* nDet */ + (4 * nDet);
 371   unsigned char* pViewData = new unsigned char [sizeView];
 372
 373   for (int i = startView; i <= endView; i++) {
 374       is.seekg (sizeHeader + i * sizeView);
 375       is.read (reinterpret_cast<char*>(pViewData), sizeView);
 376       os.write (reinterpret_cast<char*>(pViewData), sizeView);
 377       if (is.fail() || os.fail())
 378           break;
 379   }
 380
 381   delete pViewData;
 382   if (is.fail())
 383     sys_error (ERR_FATAL, "Error reading projection file");
 384   if (os.fail())
 385     sys_error (ERR_FATAL, "Error writing projection file");
 386
 387   return (! (is.fail() | os.fail()));
 388 }
 389
 390 bool
 391 Projections::copyHeader (const string& filename, ostream& os)
 392 {
 393     return copyHeader (filename.c_str(), os);
 394 }
 395
 396 bool
 397 Projections::copyHeader (const char* const filename, ostream& os)
 398 {
 399   frnetorderstream is (filename, ios::in | ios::binary);
 400   kuint16 sizeHeader, signature;
 401   is.readInt16 (sizeHeader);
 402   is.readInt16 (signature);
 403   is.seekg (0);
 404   if (signature != m_signature) {
 405     sys_error (ERR_FATAL, "Illegal signature in projection file %s", filename);
 406     return false;
 407   }
 408
 409   unsigned char* pHdrData = new unsigned char [sizeHeader];
 410   is.read (reinterpret_cast<char*>(pHdrData), sizeHeader);
 411   if (is.fail()) {
 412       sys_error (ERR_FATAL, "Error reading header");
 413       return false;
 414   }
 415
 416   os.write (reinterpret_cast<char*>(pHdrData), sizeHeader);
 417   if (os.fail()) {
 418       sys_error (ERR_FATAL, "Error writing header");
 419       return false;
 420   }
 421
 422   return true;
 423 }
 424
 425 bool
 426 Projections::write (const string& filename)
 427 {
 428   return write (filename.c_str());
 429 }
 430
 431 bool
 432 Projections::write (const char* filename)
 433 {
 434   frnetorderstream fs (filename, ios::out | ios::binary | ios::trunc | ios::ate);
 435   m_filename = filename;
 436   if (! fs) {
 437     sys_error (ERR_SEVERE, "Error opening file %s for output [projections_create]", filename);
 438     return false;
 439   }
 440
 441 #ifdef HAVE_TIME
 442   time_t t = time(NULL);
 443   tm* lt = localtime(&t);
 444   m_year = lt->tm_year;
 445   m_month = lt->tm_mon;
 446   m_day = lt->tm_mday;
 447   m_hour = lt->tm_hour;
 448   m_minute = lt->tm_min;
 449   m_second = lt->tm_sec;
 450 #endif
 451
 452   if (! headerWrite (fs))
 453       return false;
 454
 455   if (m_projData != NULL) {
 456     for (int i = 0; i < m_nView; i++) {
 457       if (! detarrayWrite (fs, *m_projData[i], i))
 458         break;
 459     }
 460   }
 461   if (! fs)
 462     return false;
 463
 464  fs.close();
 465
 466   return true;
 467 }
 468
 469 /* NAME
 470  *   detarrayRead               Read a Detector Array structure from the disk
 471  *
 472  * SYNOPSIS
 473  *   detarrayRead (proj, darray, view_num)
 474  *   DETARRAY *darray           Detector array storage location to be filled
 475  *   int      view_num          View number to read
 476  */
 477
 478 bool
 479 Projections::detarrayRead (fnetorderstream& fs, DetectorArray& darray, const int iview)
 480 {
 481   const int detval_bytes = darray.nDet() * sizeof(kfloat32);
 482   const int detheader_bytes = sizeof(kfloat64) /* view_angle */ + sizeof(kint32) /* nDet */;
 483   const int view_bytes = detheader_bytes + detval_bytes;
 484   const off_t start_data = m_headerSize + (iview * view_bytes);
 485   DetectorValue* detval_ptr = darray.detValues();
 486   kfloat64 view_angle;
 487   kuint32 nDet;
 488
 489   fs.seekg (start_data);
 490
 491   fs.readFloat64 (view_angle);
 492   fs.readInt32 (nDet);
 493   darray.setViewAngle (view_angle);
 494   //  darray.setNDet ( nDet);
 495
 496   for (unsigned int i = 0; i < nDet; i++) {
 497       kfloat32 detval;
 498       fs.readFloat32 (detval);
 499       detval_ptr[i] = detval;
 500   }
 501   if (! fs)
 502     return false;
 503
 504   return true;
 505 }
 506
 507
 508 /* NAME
 509  *   detarrayWrite                      Write detector array data to the disk
 510  *
 511  * SYNOPSIS
 512  *   detarrayWrite (darray, view_num)
 513  *   DETARRAY *darray                   Detector array data to be written
 514  *   int      view_num                  View number to write
 515  *
 516  * DESCRIPTION
 517  *       This routine writes the detarray data from the disk sequentially to
 518  *    the file that was opened with open_projections().  Data is written in
 519  *    binary format.
 520  */
 521
 522 bool
 523 Projections::detarrayWrite (fnetorderstream& fs, const DetectorArray& darray, const int iview)
 524 {
 525   const int detval_bytes = darray.nDet() * sizeof(float);
 526   const int detheader_bytes = sizeof(kfloat64) /* view_angle */ + sizeof(kint32) /* nDet */;
 527   const int view_bytes = detheader_bytes + detval_bytes;
 528   const off_t start_data = m_headerSize + (iview * view_bytes);
 529   const DetectorValue* const detval_ptr = darray.detValues();
 530   kfloat64 view_angle = darray.viewAngle();
 531   kuint32 nDet = darray.nDet();
 532
 533   fs.seekp (start_data);
 534   if (! fs) {
 535     sys_error (ERR_SEVERE, "Error seeking detectory array [detarrayWrite]");
 536     return false;
 537   }
 538
 539   fs.writeFloat64 (view_angle);
 540   fs.writeInt32 (nDet);
 541
 542   for (unsigned int i = 0; i < nDet; i++) {
 543     kfloat32 detval = detval_ptr[i];
 544     fs.writeFloat32 (detval);
 545   }
 546
 547   if (! fs)
 548     return (false);
 549
 550   return true;
 551 }
 552
 553 /* NAME
 554  *   printProjectionData                        Print projections data
 555  *
 556  * SYNOPSIS
 557  *   printProjectionData ()
 558  */
 559
 560 void
 561 Projections::printProjectionData ()
 562 {
 563   printProjectionData (0, nView() - 1);
 564 }
 565
 566 void
 567 Projections::printProjectionData (int startView, int endView)
 568 {
 569   printf("Projections Data\n\n");
 570   printf("Description: %s\n", m_remark.c_str());
 571   printf("Geometry: %s\n", Scanner::convertGeometryIDToName (m_geometry));
 572   printf("nView       = %8d           nDet = %8d\n", m_nView, m_nDet);
 573   printf("focalLength = %8.4f  fieldOfView = %8.4f\n", m_focalLength, m_fieldOfView);
 574   printf("rotStart    = %8.4f       rotInc = %8.4f\n", m_rotStart, m_rotInc);
 575   printf("detStart    = %8.4f       detInc = %8.4f\n", m_detStart, m_detInc);
 576   if (m_projData != NULL) {
 577     if (startView < 0)
 578       startView = 0;
 579     if (startView > m_nView - 1)
 580       startView = m_nView - 1;
 581     if (endView > m_nView - 1)
 582       endView = m_nView - 1;
 583     for (int ir = startView; ir <= endView - 1; ir++) {
 584       printf("View %d: angle %f\n", ir, m_projData[ir]->viewAngle());
 585       DetectorValue* detval = m_projData[ir]->detValues();
 586       for (int id = 0; id < m_projData[ir]->nDet(); id++)
 587         printf("%8.4f  ", detval[id]);
 588       printf("\n");
 589     }
 590   }
 591 }
 592
 593 void
 594 Projections::printScanInfo (ostringstream& os) const
 595 {
 596   os << "Number of detectors: " << m_nDet << "\n";
 597   os << "    Number of views: " << m_nView<< "\n";
 598   os << "             Remark: " << m_remark.c_str()<< "\n";
 599   os << "           Geometry: " << Scanner::convertGeometryIDToName (m_geometry)<< "\n";
 600   os << "       Focal Length: " << m_focalLength<< "\n";
 601   os << "      Field Of View: " << m_fieldOfView<< "\n";
 602   os << "             phmLen: " << m_phmLen<< "\n";
 603   os << "           detStart: " << m_detStart<< "\n";
 604   os << "             detInc: " << m_detInc<< "\n";
 605   os << "           rotStart: " << m_rotStart<< "\n";
 606   os << "             rotInc: " << m_rotInc<< "\n";
 607 }
 608
 609
 610