1 /*****************************************************************************
4 ** Name: projections.cpp Projection data classes
5 ** Programmer: Kevin Rosenberg
6 ** Date Started: Aug 84
8 ** This is part of the CTSim program
9 ** Copyright (C) 1983-2000 Kevin Rosenberg
11 ** $Id: projections.cpp,v 1.32 2000/12/16 06:12:47 kevin Exp $
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.
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.
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 ******************************************************************************/
29 const kuint16 Projections::m_signature = ('P'*256 + 'J');
\r
32 * Projections Constructor for projections matrix storage
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
42 Projections::Projections (const Scanner& scanner)
45 initFromScanner (scanner);
49 Projections::Projections (const int nView, const int nDet)
55 Projections::Projections (void)
61 Projections::~Projections (void)
68 Projections::init (const int nView, const int nDet)
70 m_label.setLabelType (Array2dFileLabel::L_HISTORY);
77 Projections::initFromScanner (const Scanner& scanner)
79 m_label.setLabelType (Array2dFileLabel::L_HISTORY);
81 init (scanner.nView(), scanner.nDet());
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();
90 m_detStart = -(scanner.detLen() / 2);
92 if (m_geometry == Scanner::GEOMETRY_EQUILINEAR) {
94 std::cout << "Kludge: detInc /= 2 in Projections::initFromScanner" << endl;
100 Projections::setNView (int nView) // used by MPI to reduce # of views
103 init (nView, m_nDet);
110 Projections::newProjData (void)
113 sys_error(ERR_WARNING, "m_projData != NULL [newProjData]");
115 if (m_nView > 0 && m_nDet) {
116 m_projData = new DetectorArray* [m_nView];
118 for (int i = 0; i < m_nView; i++)
119 m_projData[i] = new DetectorArray (m_nDet);
125 * projections_free Free memory allocated to projections
128 * projections_free(proj)
129 * Projections& proj Projectionss to be deallocated
133 Projections::deleteProjData (void)
135 if (m_projData != NULL) {
136 for (int i = 0; i < m_nView; i++)
137 delete m_projData[i];
146 * Projections::headerWwrite Write data header for projections file
151 Projections::headerWrite (fnetorderstream& fs)
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;
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;
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);
201 m_headerSize = fs.tellp();
202 _hsize = m_headerSize;
204 fs.writeInt16 (_hsize);
212 * projections_read_header Read data header for projections file
216 Projections::headerRead (fnetorderstream& fs)
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;
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);
242 fs.readInt16 (_hour);
243 fs.readInt16 (_minute);
244 fs.readInt16 (_second);
245 fs.readInt16 (_remarksize);
248 sys_error (ERR_SEVERE, "Error reading header information , _remarksize=%d [projections_read_header]", _remarksize);
252 if (_signature != m_signature) {
253 sys_error (ERR_SEVERE, "File %s does not have a valid projection file signature", m_filename.c_str());
257 char* pszRemarkStorage = new char [_remarksize+1];
258 fs.read (pszRemarkStorage, _remarksize);
260 sys_error (ERR_SEVERE, "Error reading remark, _remarksize = %d", _remarksize);
263 pszRemarkStorage[_remarksize] = 0;
264 m_remark = pszRemarkStorage;
265 delete pszRemarkStorage;
\r
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);
273 m_headerSize = _hsize;
277 m_calcTime = _calcTime;
278 m_rotStart = _rotStart;
280 m_detStart = _detStart;
283 m_focalLength = _focalLength;
284 m_fieldOfView = _fieldOfView;
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);
301 Projections::read (const std::string& filename)
303 return read (filename.c_str());
308 Projections::read (const char* filename)
310 frnetorderstream fileRead (filename, ios::in | ios::binary);
311 m_filename = filename;
316 if (! headerRead (fileRead))
322 for (int i = 0; i < m_nView; i++) {
323 if (! detarrayRead (fileRead, *m_projData[i], i))
333 Projections::copyViewData (const std::string& filename, std::ostream& os, int startView, int endView)
335 return copyViewData (filename.c_str(), os, startView, endView);
339 Projections::copyViewData (const char* const filename, std::ostream& os, int startView, int endView)
341 frnetorderstream is (filename, ios::in | ios::binary);
342 kuint16 sizeHeader, signature;
343 kuint32 _nView, _nDet;
345 is.readInt16 (sizeHeader);
346 is.readInt16 (signature);
347 is.readInt32 (_nView);
348 is.readInt32 (_nDet);
352 if (signature != m_signature) {
353 sys_error (ERR_FATAL, "Illegal signature in projection file %s", filename);
359 if (startView > nView - 1)
361 if (endView < 0 || endView > nView - 1)
364 if (startView > endView) { // swap if start > end
365 int tempView = endView;
367 startView = tempView;
370 int sizeView = 8 /* view_angle */ + 4 /* nDet */ + (4 * nDet);
371 unsigned char* pViewData = new unsigned char [sizeView];
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())
383 sys_error (ERR_FATAL, "Error reading projection file");
385 sys_error (ERR_FATAL, "Error writing projection file");
387 return (! (is.fail() | os.fail()));
391 Projections::copyHeader (const std::string& filename, std::ostream& os)
393 return copyHeader (filename.c_str(), os);
397 Projections::copyHeader (const char* const filename, std::ostream& os)
399 frnetorderstream is (filename, std::ios::in | std::ios::binary);
400 kuint16 sizeHeader, signature;
401 is.readInt16 (sizeHeader);
402 is.readInt16 (signature);
404 if (signature != m_signature) {
405 sys_error (ERR_FATAL, "Illegal signature in projection file %s", filename);
409 unsigned char* pHdrData = new unsigned char [sizeHeader];
410 is.read (reinterpret_cast<char*>(pHdrData), sizeHeader);
412 sys_error (ERR_FATAL, "Error reading header");
416 os.write (reinterpret_cast<char*>(pHdrData), sizeHeader);
418 sys_error (ERR_FATAL, "Error writing header");
426 Projections::write (const std::string& filename)
428 return write (filename.c_str());
432 Projections::write (const char* filename)
434 frnetorderstream fs (filename, std::ios::out | std::ios::binary | std::ios::trunc | std::ios::ate);
435 m_filename = filename;
437 sys_error (ERR_SEVERE, "Error opening file %s for output [projections_create]", filename);
442 time_t t = time(NULL);
443 tm* lt = localtime(&t);
444 m_year = lt->tm_year;
445 m_month = lt->tm_mon;
447 m_hour = lt->tm_hour;
448 m_minute = lt->tm_min;
449 m_second = lt->tm_sec;
452 if (! headerWrite (fs))
455 if (m_projData != NULL) {
456 for (int i = 0; i < m_nView; i++) {
457 if (! detarrayWrite (fs, *m_projData[i], i))
470 * detarrayRead Read a Detector Array structure from the disk
473 * detarrayRead (proj, darray, view_num)
474 * DETARRAY *darray Detector array storage location to be filled
475 * int view_num View number to read
479 Projections::detarrayRead (fnetorderstream& fs, DetectorArray& darray, const int iview)
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();
489 fs.seekg (start_data);
491 fs.readFloat64 (view_angle);
493 darray.setViewAngle (view_angle);
494 // darray.setNDet ( nDet);
496 for (unsigned int i = 0; i < nDet; i++) {
498 fs.readFloat32 (detval);
499 detval_ptr[i] = detval;
509 * detarrayWrite Write detector array data to the disk
512 * detarrayWrite (darray, view_num)
513 * DETARRAY *darray Detector array data to be written
514 * int view_num View number to write
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
523 Projections::detarrayWrite (fnetorderstream& fs, const DetectorArray& darray, const int iview)
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();
533 fs.seekp (start_data);
535 sys_error (ERR_SEVERE, "Error seeking detectory array [detarrayWrite]");
539 fs.writeFloat64 (view_angle);
540 fs.writeInt32 (nDet);
542 for (unsigned int i = 0; i < nDet; i++) {
543 kfloat32 detval = detval_ptr[i];
544 fs.writeFloat32 (detval);
554 * printProjectionData Print projections data
557 * printProjectionData ()
561 Projections::printProjectionData ()
563 printProjectionData (0, nView() - 1);
567 Projections::printProjectionData (int startView, int endView)
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) {
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]);
594 Projections::printScanInfo (std::ostringstream& os) const
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";