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.28 2000/12/06 01:46:43 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 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 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::write (const string& filename)
335 return write (filename.c_str());
339 Projections::write (const char* filename)
341 frnetorderstream fs (filename, ios::out | ios::binary | ios::trunc | ios::ate);
342 m_filename = filename;
344 sys_error (ERR_SEVERE, "Error opening file %s for output [projections_create]", filename);
349 time_t t = time(NULL);
350 tm* lt = localtime(&t);
351 m_year = lt->tm_year;
352 m_month = lt->tm_mon;
354 m_hour = lt->tm_hour;
355 m_minute = lt->tm_min;
356 m_second = lt->tm_sec;
359 if (! headerWrite (fs))
362 if (m_projData != NULL) {
363 for (int i = 0; i < m_nView; i++) {
364 if (! detarrayWrite (fs, *m_projData[i], i))
377 * detarrayRead Read a Detector Array structure from the disk
380 * detarrayRead (proj, darray, view_num)
381 * DETARRAY *darray Detector array storage location to be filled
382 * int view_num View number to read
386 Projections::detarrayRead (fnetorderstream& fs, DetectorArray& darray, const int iview)
388 const int detval_bytes = darray.nDet() * sizeof(kfloat32);
389 const int detheader_bytes = sizeof(kfloat64) /* view_angle */ + sizeof(kint32) /* nDet */;
390 const int view_bytes = detheader_bytes + detval_bytes;
391 const off_t start_data = m_headerSize + (iview * view_bytes);
392 DetectorValue* detval_ptr = darray.detValues();
396 fs.seekg (start_data);
398 fs.readFloat64 (view_angle);
400 darray.setViewAngle (view_angle);
401 // darray.setNDet ( nDet);
403 for (unsigned int i = 0; i < nDet; i++) {
405 fs.readFloat32 (detval);
406 detval_ptr[i] = detval;
416 * detarrayWrite Write detector array data to the disk
419 * detarrayWrite (darray, view_num)
420 * DETARRAY *darray Detector array data to be written
421 * int view_num View number to write
424 * This routine writes the detarray data from the disk sequentially to
425 * the file that was opened with open_projections(). Data is written in
430 Projections::detarrayWrite (fnetorderstream& fs, const DetectorArray& darray, const int iview)
432 const int detval_bytes = darray.nDet() * sizeof(float);
433 const int detheader_bytes = sizeof(kfloat64) /* view_angle */ + sizeof(kint32) /* nDet */;
434 const int view_bytes = detheader_bytes + detval_bytes;
435 const off_t start_data = m_headerSize + (iview * view_bytes);
436 const DetectorValue* const detval_ptr = darray.detValues();
437 kfloat64 view_angle = darray.viewAngle();
438 kuint32 nDet = darray.nDet();
440 fs.seekp (start_data);
442 sys_error (ERR_SEVERE, "Error seeking detectory array [detarrayWrite]");
446 fs.writeFloat64 (view_angle);
447 fs.writeInt32 (nDet);
449 for (unsigned int i = 0; i < nDet; i++) {
450 kfloat32 detval = detval_ptr[i];
451 fs.writeFloat32 (detval);
461 * printProjectionData Print projections data
464 * printProjectionData ()
468 Projections::printProjectionData (void)
470 printf("Projections Data\n\n");
471 printf("Description: %s\n", m_remark.c_str());
472 printf("Geometry: %s\n", Scanner::convertGeometryIDToName (m_geometry));
473 printf("nView = %8d nDet = %8d\n", m_nView, m_nDet);
474 printf("focalLength = %8.4f fieldOfView = %8.4f\n", m_focalLength, m_fieldOfView);
475 printf("rotStart = %8.4f rotInc = %8.4f\n", m_rotStart, m_rotInc);
476 printf("detStart = %8.4f detInc = %8.4f\n", m_detStart, m_detInc);
477 if (m_projData != NULL) {
478 for (int ir = 0; ir < m_nView; ir++) {
479 printf("View %d: angle %f\n", ir, m_projData[ir]->viewAngle());
480 DetectorValue* detval = m_projData[ir]->detValues();
481 for (int id = 0; id < m_projData[ir]->nDet(); id++)
482 printf("%8.4f ", detval[id]);
489 Projections::printScanInfo (ostringstream& os) const
491 os << "Number of detectors: " << m_nDet << "\n";
492 os << " Number of views: " << m_nView<< "\n";
493 os << " Remark: " << m_remark.c_str()<< "\n";
494 os << " Geometry: " << Scanner::convertGeometryIDToName (m_geometry)<< "\n";
495 os << " Focal Length: " << m_focalLength<< "\n";
496 os << " Field Of View: " << m_fieldOfView<< "\n";
497 os << " phmLen: " << m_phmLen<< "\n";
498 os << " detStart: " << m_detStart<< "\n";
499 os << " detInc: " << m_detInc<< "\n";
500 os << " rotStart: " << m_rotStart<< "\n";
501 os << " rotInc: " << m_rotInc<< "\n";