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