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.1 2000/06/19 02:59:34 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 {
  43   init (scanner.nView(), scanner.nDet());
  44
  45   m_phmLen = scanner.phmLen();
  46   m_rotInc = scanner.rotInc();
  47   m_detInc = scanner.detInc();
  48   m_rotStart = 0;
  49   m_detStart =  -scanner.radius() + (scanner.detInc() / 2);
  50   m_phmLen = scanner.phmLen();
  51 }
  52
  53
  54 Projections::Projections (const int nView, const int nDet)
  55 {
  56   init (nView, nDet);
  57 }
  58
  59 Projections::Projections (void)
  60 {
  61   init (0, 0);
  62 }
  63
  64 Projections::~Projections (void)
  65 {
  66   deleteProjData();
  67 }
  68
  69
  70 void
  71 Projections::init (const int nView, const int nDet)
  72 {
  73   m_nView = nView;
  74   m_nDet = nDet;
  75   m_projData = NULL;
  76   newProjData ();
  77   m_fd = -1;
  78 }
  79
  80
  81 // NAME
  82 // newProjData
  83
  84 void
  85 Projections::newProjData (void)
  86 {
  87   if (m_projData)
  88     sys_error(ERR_WARNING, "m_projData != NULL [newProjData]");
  89
  90   if (m_nView > 0 && m_nDet) {
  91     m_projData = new DetectorArray* [m_nView];
  92
  93     for (int i = 0; i < m_nView; i++)
  94       m_projData[i] = new DetectorArray (m_nDet);
  95   }
  96 }
  97
  98
  99 /* NAME
 100  *   projections_free                   Free memory allocated to projections
 101  *
 102  * SYNOPSIS
 103  *   projections_free(proj)
 104  *   Projections& proj                          Projectionss to be deallocated
 105  */
 106
 107 void
 108 Projections::deleteProjData (void)
 109 {
 110   if (m_projData != NULL) {
 111     for (int i = 0; i < m_nView; i++)
 112       delete m_projData[i];
 113
 114     delete m_projData;
 115     m_projData = NULL;
 116   }
 117 }
 118
 119
 120 /* NAME
 121  *    projections_write_header         Write data header for projections file
 122  *
 123  */
 124
 125 bool
 126 Projections::headerWrite (void)
 127 {
 128   kuint32 _hsize = m_headerSize;
 129   kuint32 _nView = m_nView;
 130   kuint32 _nDet = m_nDet;
 131   kuint32 _geom = m_geometry;
 132   kuint32 _remarksize = m_remark.length();
 133   kfloat64 _calcTime = m_calcTime;
 134   kfloat64 _rotStart = m_rotStart;
 135   kfloat64 _rotInc = m_rotInc;
 136   kfloat64 _detStart = m_detStart;
 137   kfloat64 _detInc = m_detInc;
 138   kfloat64 _phmLen = m_phmLen;
 139
 140   if (m_fd < 0) {
 141     sys_error (ERR_SEVERE, "m_fd < 0 [projections_write_header]");
 142     return false;
 143   }
 144
 145   if (lseek (m_fd, (off_t) 0, SEEK_SET) != (off_t) 0) {
 146     sys_error (ERR_WARNING, "Error seeking to beginning of fiel [projections_write_header]");
 147     return false;
 148   }
 149
 150   if (! write_nint32 (&_hsize, m_fd) ||
 151       ! write_nint32 (&_nView, m_fd) ||
 152       ! write_nint32 (&_nDet, m_fd) ||
 153       ! write_nint32 (&_geom, m_fd) ||
 154       ! write_nfloat64 (&_calcTime, m_fd) ||
 155       ! write_nfloat64 (&_rotStart, m_fd) ||
 156       ! write_nfloat64 (&_rotInc, m_fd) ||
 157       ! write_nfloat64 (&_detStart, m_fd) ||
 158       ! write_nfloat64 (&_detInc, m_fd) ||
 159       ! write_nfloat64 (&_phmLen, m_fd) ||
 160       ! write_nint32 (&_remarksize, m_fd) ||
 161       (::write (m_fd, m_remark.c_str(), _remarksize) != _remarksize)) {
 162     sys_error (ERR_SEVERE, "Error writing header information [projections_write_header] %ld", _remarksize);
 163     return false;
 164   }
 165   m_headerSize = _hsize = lseek(m_fd, (off_t) 0, SEEK_CUR);
 166   if ((lseek(m_fd, 0, SEEK_SET) != 0) || ! write_nint32 (&_hsize, m_fd)) {
 167     sys_error (ERR_SEVERE, "Error writing header information [projections_write_header]");
 168     return false;
 169   }
 170
 171   return true;
 172 }
 173
 174
 175 /* NAME
 176  *    projections_read_header         Read data header for projections file
 177  *
 178  */
 179
 180 bool
 181 Projections::headerRead (void)
 182 {
 183   kuint32 _hsize;
 184   kuint32 _nView;
 185   kuint32 _nDet;
 186   kuint32 _geom;
 187   kuint32 _remarksize = 0;
 188   kfloat64 _calcTime;
 189   kfloat64 _rotStart;
 190   kfloat64 _rotInc;
 191   kfloat64 _detStart;
 192   kfloat64 _detInc;
 193   kfloat64 _phmLen;
 194
 195   if (m_fd < 0) {
 196     sys_error (ERR_SEVERE, "m_fd < 0 [projections_read_header]");
 197     return false;
 198   }
 199
 200   if (lseek (m_fd, (off_t) 0, SEEK_SET) != (off_t) 0) {
 201     sys_error (ERR_WARNING, "Error seeking to beginning of field [projections_read_header]");
 202     return false;
 203   }
 204
 205   off_t testPos;
 206   if (! read_nint32 (&_hsize, m_fd)) {
 207     sys_error (ERR_SEVERE, "Error reading header size , _hsize=%d [projections_read_header]", _hsize);
 208     return false;
 209   }
 210   if ( ! read_nint32 (&_nView, m_fd) ||
 211       ! read_nint32 (&_nDet, m_fd) ||
 212       ! read_nint32 (&_geom, m_fd) ||
 213       ! read_nfloat64 (&_calcTime, m_fd) ||
 214       ! read_nfloat64 (&_rotStart, m_fd) ||
 215       ! read_nfloat64 (&_rotInc, m_fd) ||
 216       ! read_nfloat64 (&_detStart, m_fd) ||
 217       ! read_nfloat64 (&_detInc, m_fd) ||
 218       ! read_nfloat64 (&_phmLen, m_fd) ||
 219        ! read_nint32 (&_remarksize, m_fd))
 220       {
 221           sys_error (ERR_SEVERE, "Error reading header information , _remarksize=%d [projections_read_header]", _remarksize);
 222           return false;
 223       }
 224
 225   char remarkStorage[_remarksize+1];
 226   if (::read (m_fd, &remarkStorage, _remarksize) != _remarksize) {
 227     sys_error (ERR_SEVERE, "Error reading remark, _remarksize = %d", _remarksize);
 228     return false;
 229   }
 230   remarkStorage[_remarksize] = 0;
 231   m_remark = remarkStorage;
 232
 233   off_t _hsizeread = lseek(m_fd, 0, SEEK_CUR);
 234   if (_hsizeread != _hsize) {
 235       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);
 236       return false;
 237   }
 238
 239   m_headerSize = _hsize;
 240   m_nView = _nView;
 241   m_nDet = _nDet;
 242   m_geometry = _geom;
 243   m_calcTime = _calcTime;
 244   m_rotStart = _rotStart;
 245   m_rotInc = _rotInc;
 246   m_detStart = _detStart;
 247   m_detInc = _detInc;
 248   m_phmLen = _phmLen;
 249
 250   return true;
 251 }
 252
 253
 254 bool
 255 Projections::read (const char* filename)
 256 {
 257   if ((m_fd = open(filename, OPEN_RDONLY | O_BINARY)) < 0) {
 258     sys_error (ERR_SEVERE, "Error opening projections file %s for input [open_projections]", filename);
 259     return false;
 260   }
 261
 262   if (! headerRead()) {
 263     sys_error (ERR_SEVERE, "Error reading projections file %s [open_projections]", filename);
 264     return false;
 265   }
 266
 267   deleteProjData ();
 268   newProjData ();
 269
 270   for (int i = 0; i < m_nView; i++) {
 271     if (! detarrayRead (*m_projData[i], i))
 272       break;
 273   }
 274
 275   close (m_fd);
 276   m_fd = -1;
 277
 278   return true;
 279 }
 280
 281
 282 bool
 283 Projections::write (const char* filename)
 284 {
 285     if ((m_fd = open (filename, O_WRONLY | O_CREAT | O_TRUNC | O_BINARY
 286 #ifndef _WIN32
 287                         , S_IWRITE | S_IREAD | S_IRGRP | S_IROTH
 288 #endif
 289                         )) < 0) {
 290       sys_error (ERR_SEVERE, "Error opening file %s for output [projections_create]", filename);
 291       return false;
 292     }
 293     if (! headerWrite ()) {
 294       sys_error(ERR_SEVERE, "Error writing to file %s [create_projections]", filename);
 295       return false;
 296     }
 297
 298     headerWrite ();
 299
 300   if (m_projData != NULL) {
 301     for (int i = 0; i < m_nView; i++) {
 302       if (! detarrayWrite (*m_projData[i], i))
 303         break;
 304     }
 305   }
 306
 307   close (m_fd);
 308   m_fd = -1;
 309
 310   return true;
 311 }
 312
 313 /* NAME
 314  *   detarrayRead               Read a Detector Array structure from the disk
 315  *
 316  * SYNOPSIS
 317  *   detarrayRead (proj, darray, view_num)
 318  *   DETARRAY *darray           Detector array storage location to be filled
 319  *   int      view_num          View number to read
 320  */
 321
 322 bool
 323 Projections::detarrayRead (DetectorArray& darray, const int iview)
 324 {
 325   const int detval_bytes = darray.nDet() * sizeof(kfloat32);
 326   const int detheader_bytes = sizeof(kfloat64) /* view_angle */ + sizeof(kint32) /* nDet */;
 327   const int view_bytes = detheader_bytes + detval_bytes;
 328   const off_t start_data = m_headerSize + (iview * view_bytes);
 329   DetectorValue* detval_ptr = darray.detValues();
 330   kfloat64 view_angle;
 331   kuint32 nDet;
 332
 333   if (m_fd < 0) {
 334     sys_error (ERR_SEVERE, "Tried to read detector array on closed file [detarrayRead]");
 335     return false;
 336   }
 337
 338   if (lseek(m_fd, start_data, SEEK_SET) != start_data) {
 339     sys_error (ERR_SEVERE, "Error seeking detectory array [detarrayRead]");
 340     return false;
 341   }
 342
 343   if (! read_nfloat64 (&view_angle, m_fd) ||
 344       ! read_nint32 (&nDet, m_fd)) {
 345       sys_error (ERR_SEVERE, "Error reading detector array [detarrayRead]");
 346       return false;
 347   }
 348   darray.setViewAngle (view_angle);
 349   //  darray.setNDet ( nDet);
 350
 351   for (int i = 0; i < nDet; i++) {
 352       kfloat32 detval;
 353       if (! read_nfloat32 (&detval, m_fd)) {
 354           sys_error (ERR_SEVERE, "Error reading detector array [detarrayRead]");
 355           return false;
 356       }
 357       detval_ptr[i] = detval;
 358   }
 359
 360   return true;
 361 }
 362
 363
 364 /* NAME
 365  *   detarrayWrite                      Write detector array data to the disk
 366  *
 367  * SYNOPSIS
 368  *   detarrayWrite (darray, view_num)
 369  *   DETARRAY *darray                   Detector array data to be written
 370  *   int      view_num                  View number to write
 371  *
 372  * DESCRIPTION
 373  *       This routine writes the detarray data from the disk sequentially to
 374  *    the file that was opened with open_projections().  Data is written in
 375  *    binary format.
 376  */
 377
 378 bool
 379 Projections::detarrayWrite (const DetectorArray& darray, const int iview)
 380 {
 381   const int detval_bytes = darray.nDet() * sizeof(float);
 382   const int detheader_bytes = sizeof(kfloat64) /* view_angle */ + sizeof(kint32) /* nDet */;
 383   const int view_bytes = detheader_bytes + detval_bytes;
 384   const off_t start_data = m_headerSize + (iview * view_bytes);
 385   const DetectorValue* const detval_ptr = darray.detValues();
 386   kfloat64 view_angle = darray.viewAngle();
 387   kuint32 nDet = darray.nDet();
 388
 389   if (m_fd < 0) {
 390     sys_error (ERR_SEVERE, "Tried to write detector array with no file open [detarrayWrite]");
 391     return false;
 392   }
 393
 394   if (lseek(m_fd, start_data, SEEK_SET) != start_data) {
 395       sys_error (ERR_SEVERE, "Error seeking detectory array [detarrayWrite]");
 396       return false;
 397   }
 398
 399   if ( ! write_nfloat64 (&view_angle, m_fd)) {
 400       sys_error (ERR_SEVERE, "Error writing detector array [detarrayWrite]");
 401       return false;
 402   }
 403   if ( ! write_nint32 (&nDet, m_fd)) {
 404       sys_error (ERR_SEVERE, "Error writing detector array [detarrayWrite]");
 405       return false;
 406   }
 407   for (int i = 0; i < nDet; i++) {
 408       kfloat32 detval = detval_ptr[i];
 409       if ( ! write_nfloat32 (&detval, m_fd)) {
 410           sys_error (ERR_SEVERE, "Error writing detector array [detarrayWrite]");
 411           return false;
 412       }
 413   }
 414
 415   return true;
 416 }
 417
 418
 419 /* NAME
 420  *   prt_projections                    Print projections data
 421  *
 422  * SYNOPSIS
 423  *   prt_projections (proj)
 424  *   Projections& proj                  Projection data to be printed
 425  */
 426
 427 void
 428 Projections::printProjectionData (void)
 429 {
 430   printf("Projections Print\n\n");
 431   printf("Description: %s\n", m_remark.c_str());
 432   printf("nView = %d  nDet = %d\n", m_nView, m_nDet);
 433   printf("rotStart = %8.4f   rotInc = %8.4f\n", m_rotStart, m_rotInc);
 434   printf("detStart = %8.4f   detInc = %8.4f\n", m_detStart, m_detInc);
 435
 436   if (m_projData != NULL) {
 437       for (int ir = 0; ir < m_nView; ir++) {
 438           DetectorValue* detval = m_projData[ir]->detValues();
 439           for (int id = 0; id < m_projData[ir]->nDet(); id++)
 440               printf("%8.4f  ", detval[id]);
 441           printf("\n");
 442       }
 443   }
 444 }
 445
 446 void
 447 Projections::printScanInfo (void) const
 448 {
 449   printf ("Number of detectors: %d\n", m_nDet);
 450   printf ("    Number of views: %d\n", m_nView);
 451   printf ("             Remark: %s\n", m_remark.c_str());
 452   printf ("             phmLen: %f\n", m_phmLen);
 453   printf ("           detStart: %f\n", m_detStart);
 454   printf ("             detInc: %f\n", m_detInc);
 455   printf ("           rotStart: %f\n", m_rotStart);
 456   printf ("             rotInc: %f\n", m_rotInc);
 457 }
 458
 459