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.26 2000/09/07 05:51:52 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 ******************************************************************************/
31 * Projections Constructor for projections matrix storage
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
41 Projections::Projections (const Scanner& scanner)
44 initFromScanner (scanner);
48 Projections::Projections (const int nView, const int nDet)
54 Projections::Projections (void)
60 Projections::~Projections (void)
67 Projections::init (const int nView, const int nDet)
69 m_label.setLabelType (Array2dFileLabel::L_HISTORY);
76 Projections::initFromScanner (const Scanner& scanner)
78 m_label.setLabelType (Array2dFileLabel::L_HISTORY);
80 init (scanner.nView(), scanner.nDet());
82 m_phmLen = scanner.phmLen();
83 m_rotInc = scanner.rotInc();
84 m_detInc = scanner.detInc();
85 m_geometry = scanner.geometry();
86 m_focalLength = scanner.focalLength();
87 m_fieldOfView = scanner.fieldOfView();
89 m_detStart = -(scanner.detLen() / 2);
91 if (m_geometry == Scanner::GEOMETRY_EQUILINEAR) {
93 cout << "Kludge: detInc /= 2 in Projections::initFromScanner" << endl;
99 Projections::setNView (int nView) // used by MPI to reduce # of views
102 init (nView, m_nDet);
109 Projections::newProjData (void)
112 sys_error(ERR_WARNING, "m_projData != NULL [newProjData]");
114 if (m_nView > 0 && m_nDet) {
115 m_projData = new DetectorArray* [m_nView];
117 for (int i = 0; i < m_nView; i++)
118 m_projData[i] = new DetectorArray (m_nDet);
124 * projections_free Free memory allocated to projections
127 * projections_free(proj)
128 * Projections& proj Projectionss to be deallocated
132 Projections::deleteProjData (void)
134 if (m_projData != NULL) {
135 for (int i = 0; i < m_nView; i++)
136 delete m_projData[i];
145 * Projections::headerWwrite Write data header for projections file
150 Projections::headerWrite (fnetorderstream& fs)
152 kuint16 _hsize = m_headerSize;
153 kuint16 _signature = m_signature;
154 kuint32 _nView = m_nView;
155 kuint32 _nDet = m_nDet;
156 kuint32 _geom = m_geometry;
157 kuint16 _remarksize = m_remark.length();
158 kuint16 _year = m_year;
159 kuint16 _month = m_month;
160 kuint16 _day = m_day;
161 kuint16 _hour = m_hour;
162 kuint16 _minute = m_minute;
163 kuint16 _second = m_second;
165 kfloat64 _calcTime = m_calcTime;
166 kfloat64 _rotStart = m_rotStart;
167 kfloat64 _rotInc = m_rotInc;
168 kfloat64 _detStart = m_detStart;
169 kfloat64 _detInc = m_detInc;
170 kfloat64 _phmLen = m_phmLen;
171 kfloat64 _fieldOfView = m_fieldOfView;
172 kfloat64 _focalLength = m_focalLength;
178 fs.writeInt16 (_hsize);
179 fs.writeInt16 (_signature);
180 fs.writeInt32 (_nView);
181 fs.writeInt32 (_nDet);
182 fs.writeInt32 (_geom);
183 fs.writeFloat64 (_calcTime);
184 fs.writeFloat64 (_rotStart);
185 fs.writeFloat64 (_rotInc);
186 fs.writeFloat64 (_detStart);
187 fs.writeFloat64 (_detInc);
188 fs.writeFloat64 (_phmLen);
189 fs.writeFloat64 (_focalLength);
190 fs.writeFloat64 (_fieldOfView);
191 fs.writeInt16 (_year);
192 fs.writeInt16 (_month);
193 fs.writeInt16 (_day);
194 fs.writeInt16 (_hour);
195 fs.writeInt16 (_minute);
196 fs.writeInt16 (_second);
197 fs.writeInt16 (_remarksize);
198 fs.write (m_remark.c_str(), _remarksize);
200 m_headerSize = fs.tellp();
201 _hsize = m_headerSize;
203 fs.writeInt16 (_hsize);
211 * projections_read_header Read data header for projections file
215 Projections::headerRead (fnetorderstream& fs)
217 kuint16 _hsize, _signature, _year, _month, _day, _hour, _minute, _second, _remarksize = 0;
218 kuint32 _nView, _nDet, _geom;
219 kfloat64 _calcTime, _rotStart, _rotInc, _detStart, _detInc, _phmLen, _focalLength, _fieldOfView;
225 fs.readInt16 (_hsize);
226 fs.readInt16 (_signature);
227 fs.readInt32 (_nView);
228 fs.readInt32 (_nDet);
229 fs.readInt32 (_geom);
230 fs.readFloat64 (_calcTime);
231 fs.readFloat64 (_rotStart);
232 fs.readFloat64 (_rotInc);
233 fs.readFloat64 (_detStart);
234 fs.readFloat64 (_detInc);
235 fs.readFloat64 (_phmLen);
236 fs.readFloat64 (_focalLength);
237 fs.readFloat64 (_fieldOfView);
238 fs.readInt16 (_year);
239 fs.readInt16 (_month);
241 fs.readInt16 (_hour);
242 fs.readInt16 (_minute);
243 fs.readInt16 (_second);
244 fs.readInt16 (_remarksize);
247 sys_error (ERR_SEVERE, "Error reading header information , _remarksize=%d [projections_read_header]", _remarksize);
251 if (_signature != m_signature) {
252 sys_error (ERR_SEVERE, "File %s does not have a valid projection file signature", m_filename.c_str());
256 char remarkStorage[_remarksize+1];
257 fs.read (remarkStorage, _remarksize);
259 sys_error (ERR_SEVERE, "Error reading remark, _remarksize = %d", _remarksize);
262 remarkStorage[_remarksize] = 0;
263 m_remark = remarkStorage;
265 off_t _hsizeread = fs.tellg();
266 if (!fs || _hsizeread != _hsize) {
267 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);
271 m_headerSize = _hsize;
275 m_calcTime = _calcTime;
276 m_rotStart = _rotStart;
278 m_detStart = _detStart;
281 m_focalLength = _focalLength;
282 m_fieldOfView = _fieldOfView;
290 m_label.setLabelType (Array2dFileLabel::L_HISTORY);
291 m_label.setLabelString (m_remark);
292 m_label.setCalcTime (m_calcTime);
293 m_label.setDateTime (m_year, m_month, m_day, m_hour, m_minute, m_second);
299 Projections::read (const string& filename)
301 return read (filename.c_str());
306 Projections::read (const char* filename)
308 frnetorderstream fileRead (filename, ios::in | ios::binary);
309 m_filename = filename;
314 if (! headerRead (fileRead))
320 for (int i = 0; i < m_nView; i++) {
321 if (! detarrayRead (fileRead, *m_projData[i], i))
331 Projections::write (const string& filename)
333 return write (filename.c_str());
337 Projections::write (const char* filename)
339 frnetorderstream fs (filename, ios::out | ios::binary | ios::trunc | ios::ate);
340 m_filename = filename;
342 sys_error (ERR_SEVERE, "Error opening file %s for output [projections_create]", filename);
347 time_t t = time(NULL);
348 tm* lt = localtime(&t);
349 m_year = lt->tm_year;
350 m_month = lt->tm_mon;
352 m_hour = lt->tm_hour;
353 m_minute = lt->tm_min;
354 m_second = lt->tm_sec;
357 if (! headerWrite (fs))
360 if (m_projData != NULL) {
361 for (int i = 0; i < m_nView; i++) {
362 if (! detarrayWrite (fs, *m_projData[i], i))
375 * detarrayRead Read a Detector Array structure from the disk
378 * detarrayRead (proj, darray, view_num)
379 * DETARRAY *darray Detector array storage location to be filled
380 * int view_num View number to read
384 Projections::detarrayRead (fnetorderstream& fs, DetectorArray& darray, const int iview)
386 const int detval_bytes = darray.nDet() * sizeof(kfloat32);
387 const int detheader_bytes = sizeof(kfloat64) /* view_angle */ + sizeof(kint32) /* nDet */;
388 const int view_bytes = detheader_bytes + detval_bytes;
389 const off_t start_data = m_headerSize + (iview * view_bytes);
390 DetectorValue* detval_ptr = darray.detValues();
394 fs.seekg (start_data);
396 fs.readFloat64 (view_angle);
398 darray.setViewAngle (view_angle);
399 // darray.setNDet ( nDet);
401 for (unsigned int i = 0; i < nDet; i++) {
403 fs.readFloat32 (detval);
404 detval_ptr[i] = detval;
414 * detarrayWrite Write detector array data to the disk
417 * detarrayWrite (darray, view_num)
418 * DETARRAY *darray Detector array data to be written
419 * int view_num View number to write
422 * This routine writes the detarray data from the disk sequentially to
423 * the file that was opened with open_projections(). Data is written in
428 Projections::detarrayWrite (fnetorderstream& fs, const DetectorArray& darray, const int iview)
430 const int detval_bytes = darray.nDet() * sizeof(float);
431 const int detheader_bytes = sizeof(kfloat64) /* view_angle */ + sizeof(kint32) /* nDet */;
432 const int view_bytes = detheader_bytes + detval_bytes;
433 const off_t start_data = m_headerSize + (iview * view_bytes);
434 const DetectorValue* const detval_ptr = darray.detValues();
435 kfloat64 view_angle = darray.viewAngle();
436 kuint32 nDet = darray.nDet();
438 fs.seekp (start_data);
440 sys_error (ERR_SEVERE, "Error seeking detectory array [detarrayWrite]");
444 fs.writeFloat64 (view_angle);
445 fs.writeInt32 (nDet);
447 for (unsigned int i = 0; i < nDet; i++) {
448 kfloat32 detval = detval_ptr[i];
449 fs.writeFloat32 (detval);
459 * printProjectionData Print projections data
462 * printProjectionData ()
466 Projections::printProjectionData (void)
468 printf("Projections Data\n\n");
469 printf("Description: %s\n", m_remark.c_str());
470 printf("Geometry: %s\n", Scanner::convertGeometryIDToName (m_geometry));
471 printf("nView = %8d nDet = %8d\n", m_nView, m_nDet);
472 printf("focalLength = %8.4f fieldOfView = %8.4f\n", m_focalLength, m_fieldOfView);
473 printf("rotStart = %8.4f rotInc = %8.4f\n", m_rotStart, m_rotInc);
474 printf("detStart = %8.4f detInc = %8.4f\n", m_detStart, m_detInc);
475 if (m_projData != NULL) {
476 for (int ir = 0; ir < m_nView; ir++) {
477 printf("View %d: angle %f\n", ir, m_projData[ir]->viewAngle());
478 DetectorValue* detval = m_projData[ir]->detValues();
479 for (int id = 0; id < m_projData[ir]->nDet(); id++)
480 printf("%8.4f ", detval[id]);
487 Projections::printScanInfo (ostringstream& os) const
489 os << "Number of detectors: " << m_nDet << "\n";
490 os << " Number of views: " << m_nView<< "\n";
491 os << " Remark: " << m_remark.c_str()<< "\n";
492 os << " Geometry: " << Scanner::convertGeometryIDToName (m_geometry)<< "\n";
493 os << " Focal Length: " << m_focalLength<< "\n";
494 os << " Field Of View: " << m_fieldOfView<< "\n";
495 os << " phmLen: " << m_phmLen<< "\n";
496 os << " detStart: " << m_detStart<< "\n";
497 os << " detInc: " << m_detInc<< "\n";
498 os << " rotStart: " << m_rotStart<< "\n";
499 os << " rotInc: " << m_rotInc<< "\n";