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.2 2000/06/19 17:58:20 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 (void)
 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   if (m_fd < 0) {
 155     sys_error (ERR_SEVERE, "m_fd < 0 [projections_write_header]");
 156     return false;
 157   }
 158
 159   if (lseek (m_fd, (off_t) 0, SEEK_SET) != (off_t) 0) {
 160     sys_error (ERR_WARNING, "Error seeking to beginning of fiel [projections_write_header]");
 161     return false;
 162   }
 163
 164   if (! write_nint32 (&_hsize, m_fd) ||
 165       ! write_nint32 (&_nView, m_fd) ||
 166       ! write_nint32 (&_nDet, m_fd) ||
 167       ! write_nint32 (&_geom, m_fd) ||
 168       ! write_nfloat64 (&_calcTime, m_fd) ||
 169       ! write_nfloat64 (&_rotStart, m_fd) ||
 170       ! write_nfloat64 (&_rotInc, m_fd) ||
 171       ! write_nfloat64 (&_detStart, m_fd) ||
 172       ! write_nfloat64 (&_detInc, m_fd) ||
 173       ! write_nfloat64 (&_phmLen, m_fd) ||
 174       ! write_nint32 (&_remarksize, m_fd) ||
 175       (::write (m_fd, m_remark.c_str(), _remarksize) != _remarksize)) {
 176     sys_error (ERR_SEVERE, "Error writing header information [projections_write_header] %ld", _remarksize);
 177     return false;
 178   }
 179   m_headerSize = _hsize = lseek(m_fd, (off_t) 0, SEEK_CUR);
 180   if ((lseek(m_fd, 0, SEEK_SET) != 0) || ! write_nint32 (&_hsize, m_fd)) {
 181     sys_error (ERR_SEVERE, "Error writing header information [projections_write_header]");
 182     return false;
 183   }
 184
 185   return true;
 186 }
 187
 188
 189 /* NAME
 190  *    projections_read_header         Read data header for projections file
 191  *
 192  */
 193
 194 bool
 195 Projections::headerRead (void)
 196 {
 197   kuint32 _hsize;
 198   kuint32 _nView;
 199   kuint32 _nDet;
 200   kuint32 _geom;
 201   kuint32 _remarksize = 0;
 202   kfloat64 _calcTime;
 203   kfloat64 _rotStart;
 204   kfloat64 _rotInc;
 205   kfloat64 _detStart;
 206   kfloat64 _detInc;
 207   kfloat64 _phmLen;
 208
 209   if (m_fd < 0) {
 210     sys_error (ERR_SEVERE, "m_fd < 0 [projections_read_header]");
 211     return false;
 212   }
 213
 214   if (lseek (m_fd, (off_t) 0, SEEK_SET) != (off_t) 0) {
 215     sys_error (ERR_WARNING, "Error seeking to beginning of field [projections_read_header]");
 216     return false;
 217   }
 218
 219   off_t testPos;
 220   if (! read_nint32 (&_hsize, m_fd)) {
 221     sys_error (ERR_SEVERE, "Error reading header size , _hsize=%d [projections_read_header]", _hsize);
 222     return false;
 223   }
 224   if ( ! read_nint32 (&_nView, m_fd) ||
 225       ! read_nint32 (&_nDet, m_fd) ||
 226       ! read_nint32 (&_geom, m_fd) ||
 227       ! read_nfloat64 (&_calcTime, m_fd) ||
 228       ! read_nfloat64 (&_rotStart, m_fd) ||
 229       ! read_nfloat64 (&_rotInc, m_fd) ||
 230       ! read_nfloat64 (&_detStart, m_fd) ||
 231       ! read_nfloat64 (&_detInc, m_fd) ||
 232       ! read_nfloat64 (&_phmLen, m_fd) ||
 233        ! read_nint32 (&_remarksize, m_fd))
 234       {
 235           sys_error (ERR_SEVERE, "Error reading header information , _remarksize=%d [projections_read_header]", _remarksize);
 236           return false;
 237       }
 238
 239   char remarkStorage[_remarksize+1];
 240   if (::read (m_fd, &remarkStorage, _remarksize) != _remarksize) {
 241     sys_error (ERR_SEVERE, "Error reading remark, _remarksize = %d", _remarksize);
 242     return false;
 243   }
 244   remarkStorage[_remarksize] = 0;
 245   m_remark = remarkStorage;
 246
 247   off_t _hsizeread = lseek(m_fd, 0, SEEK_CUR);
 248   if (_hsizeread != _hsize) {
 249       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);
 250       return false;
 251   }
 252
 253   m_headerSize = _hsize;
 254   m_nView = _nView;
 255   m_nDet = _nDet;
 256   m_geometry = _geom;
 257   m_calcTime = _calcTime;
 258   m_rotStart = _rotStart;
 259   m_rotInc = _rotInc;
 260   m_detStart = _detStart;
 261   m_detInc = _detInc;
 262   m_phmLen = _phmLen;
 263
 264   return true;
 265 }
 266
 267
 268 bool
 269 Projections::read (const char* filename)
 270 {
 271   if ((m_fd = open(filename, OPEN_RDONLY | O_BINARY)) < 0) {
 272     sys_error (ERR_SEVERE, "Error opening projections file %s for input [open_projections]", filename);
 273     return false;
 274   }
 275
 276   if (! headerRead()) {
 277     sys_error (ERR_SEVERE, "Error reading projections file %s [open_projections]", filename);
 278     return false;
 279   }
 280
 281   deleteProjData ();
 282   newProjData ();
 283
 284   for (int i = 0; i < m_nView; i++) {
 285     if (! detarrayRead (*m_projData[i], i))
 286       break;
 287   }
 288
 289   close (m_fd);
 290   m_fd = -1;
 291
 292   return true;
 293 }
 294
 295
 296 bool
 297 Projections::write (const char* filename)
 298 {
 299     if ((m_fd = open (filename, O_WRONLY | O_CREAT | O_TRUNC | O_BINARY
 300 #ifndef _WIN32
 301                         , S_IWRITE | S_IREAD | S_IRGRP | S_IROTH
 302 #endif
 303                         )) < 0) {
 304       sys_error (ERR_SEVERE, "Error opening file %s for output [projections_create]", filename);
 305       return false;
 306     }
 307     if (! headerWrite ()) {
 308       sys_error(ERR_SEVERE, "Error writing to file %s [create_projections]", filename);
 309       return false;
 310     }
 311
 312     headerWrite ();
 313
 314   if (m_projData != NULL) {
 315     for (int i = 0; i < m_nView; i++) {
 316       if (! detarrayWrite (*m_projData[i], i))
 317         break;
 318     }
 319   }
 320
 321   close (m_fd);
 322   m_fd = -1;
 323
 324   return true;
 325 }
 326
 327 /* NAME
 328  *   detarrayRead               Read a Detector Array structure from the disk
 329  *
 330  * SYNOPSIS
 331  *   detarrayRead (proj, darray, view_num)
 332  *   DETARRAY *darray           Detector array storage location to be filled
 333  *   int      view_num          View number to read
 334  */
 335
 336 bool
 337 Projections::detarrayRead (DetectorArray& darray, const int iview)
 338 {
 339   const int detval_bytes = darray.nDet() * sizeof(kfloat32);
 340   const int detheader_bytes = sizeof(kfloat64) /* view_angle */ + sizeof(kint32) /* nDet */;
 341   const int view_bytes = detheader_bytes + detval_bytes;
 342   const off_t start_data = m_headerSize + (iview * view_bytes);
 343   DetectorValue* detval_ptr = darray.detValues();
 344   kfloat64 view_angle;
 345   kuint32 nDet;
 346
 347   if (m_fd < 0) {
 348     sys_error (ERR_SEVERE, "Tried to read detector array on closed file [detarrayRead]");
 349     return false;
 350   }
 351
 352   if (lseek(m_fd, start_data, SEEK_SET) != start_data) {
 353     sys_error (ERR_SEVERE, "Error seeking detectory array [detarrayRead]");
 354     return false;
 355   }
 356
 357   if (! read_nfloat64 (&view_angle, m_fd) ||
 358       ! read_nint32 (&nDet, m_fd)) {
 359       sys_error (ERR_SEVERE, "Error reading detector array [detarrayRead]");
 360       return false;
 361   }
 362   darray.setViewAngle (view_angle);
 363   //  darray.setNDet ( nDet);
 364
 365   for (int i = 0; i < nDet; i++) {
 366       kfloat32 detval;
 367       if (! read_nfloat32 (&detval, m_fd)) {
 368           sys_error (ERR_SEVERE, "Error reading detector array [detarrayRead]");
 369           return false;
 370       }
 371       detval_ptr[i] = detval;
 372   }
 373
 374   return true;
 375 }
 376
 377
 378 /* NAME
 379  *   detarrayWrite                      Write detector array data to the disk
 380  *
 381  * SYNOPSIS
 382  *   detarrayWrite (darray, view_num)
 383  *   DETARRAY *darray                   Detector array data to be written
 384  *   int      view_num                  View number to write
 385  *
 386  * DESCRIPTION
 387  *       This routine writes the detarray data from the disk sequentially to
 388  *    the file that was opened with open_projections().  Data is written in
 389  *    binary format.
 390  */
 391
 392 bool
 393 Projections::detarrayWrite (const DetectorArray& darray, const int iview)
 394 {
 395   const int detval_bytes = darray.nDet() * sizeof(float);
 396   const int detheader_bytes = sizeof(kfloat64) /* view_angle */ + sizeof(kint32) /* nDet */;
 397   const int view_bytes = detheader_bytes + detval_bytes;
 398   const off_t start_data = m_headerSize + (iview * view_bytes);
 399   const DetectorValue* const detval_ptr = darray.detValues();
 400   kfloat64 view_angle = darray.viewAngle();
 401   kuint32 nDet = darray.nDet();
 402
 403   if (m_fd < 0) {
 404     sys_error (ERR_SEVERE, "Tried to write detector array with no file open [detarrayWrite]");
 405     return false;
 406   }
 407
 408   if (lseek(m_fd, start_data, SEEK_SET) != start_data) {
 409       sys_error (ERR_SEVERE, "Error seeking detectory array [detarrayWrite]");
 410       return false;
 411   }
 412
 413   if ( ! write_nfloat64 (&view_angle, m_fd)) {
 414       sys_error (ERR_SEVERE, "Error writing detector array [detarrayWrite]");
 415       return false;
 416   }
 417   if ( ! write_nint32 (&nDet, m_fd)) {
 418       sys_error (ERR_SEVERE, "Error writing detector array [detarrayWrite]");
 419       return false;
 420   }
 421   for (int i = 0; i < nDet; i++) {
 422       kfloat32 detval = detval_ptr[i];
 423       if ( ! write_nfloat32 (&detval, m_fd)) {
 424           sys_error (ERR_SEVERE, "Error writing detector array [detarrayWrite]");
 425           return false;
 426       }
 427   }
 428
 429   return true;
 430 }
 431
 432
 433 /* NAME
 434  *   prt_projections                    Print projections data
 435  *
 436  * SYNOPSIS
 437  *   prt_projections (proj)
 438  *   Projections& proj                  Projection data to be printed
 439  */
 440
 441 void
 442 Projections::printProjectionData (void)
 443 {
 444   printf("Projections Print\n\n");
 445   printf("Description: %s\n", m_remark.c_str());
 446   printf("nView = %d  nDet = %d\n", m_nView, m_nDet);
 447   printf("rotStart = %8.4f   rotInc = %8.4f\n", m_rotStart, m_rotInc);
 448   printf("detStart = %8.4f   detInc = %8.4f\n", m_detStart, m_detInc);
 449
 450   if (m_projData != NULL) {
 451       for (int ir = 0; ir < m_nView; ir++) {
 452           DetectorValue* detval = m_projData[ir]->detValues();
 453           for (int id = 0; id < m_projData[ir]->nDet(); id++)
 454               printf("%8.4f  ", detval[id]);
 455           printf("\n");
 456       }
 457   }
 458 }
 459
 460 void
 461 Projections::printScanInfo (void) const
 462 {
 463   printf ("Number of detectors: %d\n", m_nDet);
 464   printf ("    Number of views: %d\n", m_nView);
 465   printf ("             Remark: %s\n", m_remark.c_str());
 466   printf ("             phmLen: %f\n", m_phmLen);
 467   printf ("           detStart: %f\n", m_detStart);
 468   printf ("             detInc: %f\n", m_detInc);
 469   printf ("           rotStart: %f\n", m_rotStart);
 470   printf ("             rotInc: %f\n", m_rotInc);
 471 }
 472
 473