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.33 2000/12/17 22:30:09 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 m_filename = filename;
\r
312 frnetorderstream fileRead (m_filename.c_str(), std::ios::in | std::ios::binary);
\r
314 frnetorderstream fileRead (m_filename.c_str(), std::ios::in | std::ios::binary | std::ios::nocreate);
\r
320 if (! headerRead (fileRead))
326 for (int i = 0; i < m_nView; i++) {
327 if (! detarrayRead (fileRead, *m_projData[i], i))
337 Projections::copyViewData (const std::string& filename, std::ostream& os, int startView, int endView)
339 return copyViewData (filename.c_str(), os, startView, endView);
343 Projections::copyViewData (const char* const filename, std::ostream& os, int startView, int endView)
345 frnetorderstream is (filename, ios::in | ios::binary);
346 kuint16 sizeHeader, signature;
347 kuint32 _nView, _nDet;
349 is.readInt16 (sizeHeader);
350 is.readInt16 (signature);
351 is.readInt32 (_nView);
352 is.readInt32 (_nDet);
356 if (signature != m_signature) {
357 sys_error (ERR_FATAL, "Illegal signature in projection file %s", filename);
363 if (startView > nView - 1)
365 if (endView < 0 || endView > nView - 1)
368 if (startView > endView) { // swap if start > end
369 int tempView = endView;
371 startView = tempView;
374 int sizeView = 8 /* view_angle */ + 4 /* nDet */ + (4 * nDet);
375 unsigned char* pViewData = new unsigned char [sizeView];
377 for (int i = startView; i <= endView; i++) {
378 is.seekg (sizeHeader + i * sizeView);
379 is.read (reinterpret_cast<char*>(pViewData), sizeView);
380 os.write (reinterpret_cast<char*>(pViewData), sizeView);
381 if (is.fail() || os.fail())
387 sys_error (ERR_FATAL, "Error reading projection file");
389 sys_error (ERR_FATAL, "Error writing projection file");
391 return (! (is.fail() | os.fail()));
395 Projections::copyHeader (const std::string& filename, std::ostream& os)
397 return copyHeader (filename.c_str(), os);
401 Projections::copyHeader (const char* const filename, std::ostream& os)
403 frnetorderstream is (filename, std::ios::in | std::ios::binary);
404 kuint16 sizeHeader, signature;
405 is.readInt16 (sizeHeader);
406 is.readInt16 (signature);
408 if (signature != m_signature) {
409 sys_error (ERR_FATAL, "Illegal signature in projection file %s", filename);
413 unsigned char* pHdrData = new unsigned char [sizeHeader];
414 is.read (reinterpret_cast<char*>(pHdrData), sizeHeader);
416 sys_error (ERR_FATAL, "Error reading header");
420 os.write (reinterpret_cast<char*>(pHdrData), sizeHeader);
422 sys_error (ERR_FATAL, "Error writing header");
430 Projections::write (const std::string& filename)
432 return write (filename.c_str());
436 Projections::write (const char* filename)
438 frnetorderstream fs (filename, std::ios::out | std::ios::binary | std::ios::trunc | std::ios::ate);
439 m_filename = filename;
441 sys_error (ERR_SEVERE, "Error opening file %s for output [projections_create]", filename);
446 time_t t = time(NULL);
447 tm* lt = localtime(&t);
448 m_year = lt->tm_year;
449 m_month = lt->tm_mon;
451 m_hour = lt->tm_hour;
452 m_minute = lt->tm_min;
453 m_second = lt->tm_sec;
456 if (! headerWrite (fs))
459 if (m_projData != NULL) {
460 for (int i = 0; i < m_nView; i++) {
461 if (! detarrayWrite (fs, *m_projData[i], i))
474 * detarrayRead Read a Detector Array structure from the disk
477 * detarrayRead (proj, darray, view_num)
478 * DETARRAY *darray Detector array storage location to be filled
479 * int view_num View number to read
483 Projections::detarrayRead (fnetorderstream& fs, DetectorArray& darray, const int iview)
485 const int detval_bytes = darray.nDet() * sizeof(kfloat32);
486 const int detheader_bytes = sizeof(kfloat64) /* view_angle */ + sizeof(kint32) /* nDet */;
487 const int view_bytes = detheader_bytes + detval_bytes;
488 const off_t start_data = m_headerSize + (iview * view_bytes);
489 DetectorValue* detval_ptr = darray.detValues();
493 fs.seekg (start_data);
495 fs.readFloat64 (view_angle);
497 darray.setViewAngle (view_angle);
498 // darray.setNDet ( nDet);
500 for (unsigned int i = 0; i < nDet; i++) {
502 fs.readFloat32 (detval);
503 detval_ptr[i] = detval;
513 * detarrayWrite Write detector array data to the disk
516 * detarrayWrite (darray, view_num)
517 * DETARRAY *darray Detector array data to be written
518 * int view_num View number to write
521 * This routine writes the detarray data from the disk sequentially to
522 * the file that was opened with open_projections(). Data is written in
527 Projections::detarrayWrite (fnetorderstream& fs, const DetectorArray& darray, const int iview)
529 const int detval_bytes = darray.nDet() * sizeof(float);
530 const int detheader_bytes = sizeof(kfloat64) /* view_angle */ + sizeof(kint32) /* nDet */;
531 const int view_bytes = detheader_bytes + detval_bytes;
532 const off_t start_data = m_headerSize + (iview * view_bytes);
533 const DetectorValue* const detval_ptr = darray.detValues();
534 kfloat64 view_angle = darray.viewAngle();
535 kuint32 nDet = darray.nDet();
537 fs.seekp (start_data);
539 sys_error (ERR_SEVERE, "Error seeking detectory array [detarrayWrite]");
543 fs.writeFloat64 (view_angle);
544 fs.writeInt32 (nDet);
546 for (unsigned int i = 0; i < nDet; i++) {
547 kfloat32 detval = detval_ptr[i];
548 fs.writeFloat32 (detval);
558 * printProjectionData Print projections data
561 * printProjectionData ()
565 Projections::printProjectionData ()
567 printProjectionData (0, nView() - 1);
571 Projections::printProjectionData (int startView, int endView)
573 printf("Projections Data\n\n");
574 printf("Description: %s\n", m_remark.c_str());
575 printf("Geometry: %s\n", Scanner::convertGeometryIDToName (m_geometry));
576 printf("nView = %8d nDet = %8d\n", m_nView, m_nDet);
577 printf("focalLength = %8.4f fieldOfView = %8.4f\n", m_focalLength, m_fieldOfView);
578 printf("rotStart = %8.4f rotInc = %8.4f\n", m_rotStart, m_rotInc);
579 printf("detStart = %8.4f detInc = %8.4f\n", m_detStart, m_detInc);
580 if (m_projData != NULL) {
583 if (startView > m_nView - 1)
584 startView = m_nView - 1;
585 if (endView > m_nView - 1)
586 endView = m_nView - 1;
587 for (int ir = startView; ir <= endView - 1; ir++) {
588 printf("View %d: angle %f\n", ir, m_projData[ir]->viewAngle());
589 DetectorValue* detval = m_projData[ir]->detValues();
590 for (int id = 0; id < m_projData[ir]->nDet(); id++)
591 printf("%8.4f ", detval[id]);
598 Projections::printScanInfo (std::ostringstream& os) const
600 os << "Number of detectors: " << m_nDet << "\n";
601 os << " Number of views: " << m_nView<< "\n";
602 os << " Remark: " << m_remark.c_str()<< "\n";
603 os << " Geometry: " << Scanner::convertGeometryIDToName (m_geometry)<< "\n";
604 os << " Focal Length: " << m_focalLength<< "\n";
605 os << " Field Of View: " << m_fieldOfView<< "\n";
606 os << " phmLen: " << m_phmLen<< "\n";
607 os << " detStart: " << m_detStart<< "\n";
608 os << " detInc: " << m_detInc<< "\n";
609 os << " rotStart: " << m_rotStart<< "\n";
610 os << " rotInc: " << m_rotInc<< "\n";