** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: projections.cpp,v 1.24 2000/08/31 08:38:58 kevin Exp $
+** $Id: projections.cpp,v 1.29 2000/12/16 02:31:00 kevin Exp $
**
** This program is free software; you can redistribute it and/or modify
** it under the terms of the GNU General Public License (version 2) as
******************************************************************************/
#include "ct.h"
-
+\r
+const kuint16 Projections::m_signature = ('P'*256 + 'J');\r
/* NAME
* Projections Constructor for projections matrix storage
void
Projections::init (const int nView, const int nDet)
{
+ m_label.setLabelType (Array2dFileLabel::L_HISTORY);
m_nView = nView;
m_nDet = nDet;
newProjData ();
void
Projections::initFromScanner (const Scanner& scanner)
{
+ m_label.setLabelType (Array2dFileLabel::L_HISTORY);
deleteProjData();
init (scanner.nView(), scanner.nDet());
return false;
}
- char remarkStorage[_remarksize+1];
- fs.read (remarkStorage, _remarksize);
+ char* pszRemarkStorage = new char [_remarksize+1];
+ fs.read (pszRemarkStorage, _remarksize);
if (! fs) {
sys_error (ERR_SEVERE, "Error reading remark, _remarksize = %d", _remarksize);
return false;
}
- remarkStorage[_remarksize] = 0;
- m_remark = remarkStorage;
+ pszRemarkStorage[_remarksize] = 0;
+ m_remark = pszRemarkStorage;
+ delete pszRemarkStorage;\r
off_t _hsizeread = fs.tellg();
if (!fs || _hsizeread != _hsize) {
}
+bool
+Projections::copyViewData (const string& filename, ostream& os, int startView, int endView)
+{
+ return copyViewData (filename.c_str(), os, startView, endView);
+}
+
+bool
+Projections::copyViewData (const char* const filename, ostream& os, int startView, int endView)
+{
+ frnetorderstream is (filename, ios::in | ios::binary);
+ kuint16 sizeHeader, signature;
+ kuint32 _nView, _nDet;
+
+ is.readInt16 (sizeHeader);
+ is.readInt16 (signature);
+ is.readInt32 (_nView);
+ is.readInt32 (_nDet);
+ int nView = _nView;
+ int nDet = _nDet;
+
+ if (signature != m_signature) {
+ sys_error (ERR_FATAL, "Illegal signature in projection file %s", filename);
+ return false;
+ }
+
+ if (startView < 0)
+ startView = 0;
+ if (startView > nView - 1)
+ startView = nView;
+ if (endView < 0 || endView > nView - 1)
+ endView = nView - 1;
+
+ if (startView > endView) { // swap if start > end
+ int tempView = endView;
+ endView = startView;
+ startView = tempView;
+ }
+
+ int sizeView = 8 /* view_angle */ + 4 /* nDet */ + (4 * nDet);
+ unsigned char* pViewData = new unsigned char [sizeView];
+
+ for (int i = startView; i <= endView; i++) {
+ is.seekg (sizeHeader + i * sizeView);
+ is.read (pViewData, sizeView);
+ os.write (pViewData, sizeView);
+ if (is.fail() || os.fail())
+ break;
+ }
+
+ delete pViewData;
+ if (is.fail())
+ sys_error (ERR_FATAL, "Error reading projection file");
+ if (os.fail())
+ sys_error (ERR_FATAL, "Error writing projection file");
+
+ return (! (is.fail() | os.fail()));
+}
+
+bool
+Projections::copyHeader (const string& filename, ostream& os)
+{
+ return copyHeader (filename.c_str(), os);
+}
+
+bool
+Projections::copyHeader (const char* const filename, ostream& os)
+{
+ frnetorderstream is (filename, ios::in | ios::binary);
+ kuint16 sizeHeader, signature;
+ is.readInt16 (sizeHeader);
+ is.readInt16 (signature);
+ is.seekg (0);
+ if (signature != m_signature) {
+ sys_error (ERR_FATAL, "Illegal signature in projection file %s", filename);
+ return false;
+ }
+
+ unsigned char* pHdrData = new unsigned char [sizeHeader];
+ is.read (pHdrData, sizeHeader);
+ if (is.fail()) {
+ sys_error (ERR_FATAL, "Error reading header");
+ return false;
+ }
+
+ os.write (pHdrData, sizeHeader);
+ if (os.fail()) {
+ sys_error (ERR_FATAL, "Error writing header");
+ return false;
+ }
+
+ return true;
+}
+
bool
Projections::write (const string& filename)
{
return false;
}
+#ifdef HAVE_TIME
time_t t = time(NULL);
tm* lt = localtime(&t);
m_year = lt->tm_year;
m_hour = lt->tm_hour;
m_minute = lt->tm_min;
m_second = lt->tm_sec;
+#endif
if (! headerWrite (fs))
return false;
}
/* NAME
- * prt_projections Print projections data
+ * printProjectionData Print projections data
*
* SYNOPSIS
- * prt_projections (proj)
- * Projections& proj Projection data to be printed
+ * printProjectionData ()
*/
void
-Projections::printProjectionData (void)
+Projections::printProjectionData ()
+{
+ printProjectionData (0, nView() - 1);
+}
+
+void
+Projections::printProjectionData (int startView, int endView)
{
printf("Projections Data\n\n");
printf("Description: %s\n", m_remark.c_str());
printf("rotStart = %8.4f rotInc = %8.4f\n", m_rotStart, m_rotInc);
printf("detStart = %8.4f detInc = %8.4f\n", m_detStart, m_detInc);
if (m_projData != NULL) {
- for (int ir = 0; ir < m_nView; ir++) {
- printf("View %d: angle %f\n", ir, m_projData[ir]->viewAngle());
- DetectorValue* detval = m_projData[ir]->detValues();
- for (int id = 0; id < m_projData[ir]->nDet(); id++)
- printf("%8.4f ", detval[id]);
- printf("\n");
- }
+ if (startView < 0)
+ startView = 0;
+ if (startView > m_nView - 1)
+ startView = m_nView - 1;
+ if (endView > m_nView - 1)
+ endView = m_nView - 1;
+ for (int ir = startView; ir <= endView - 1; ir++) {
+ printf("View %d: angle %f\n", ir, m_projData[ir]->viewAngle());
+ DetectorValue* detval = m_projData[ir]->detValues();
+ for (int id = 0; id < m_projData[ir]->nDet(); id++)
+ printf("%8.4f ", detval[id]);
+ printf("\n");
+ }
}
}
-/* NAME
- * Projections::reconstruct Reconstruct Image from Projections
- *
- * SYNOPSIS
- * im = proj.reconstruct (im, filt_type, filt_param, interp_type)
- * IMAGE *im Output image
- * int filt_type Type of convolution filter to use
- * double filt_param Filter specific parameter
- * Currently, used only with Hamming filters
- * int interp_type Type of interpolation method to use
- *
- * ALGORITHM
- *
- * Calculate one-dimensional filter in spatial domain
- * Allocate & clear (zero) the 2d output image array
- * For each projection view
- * Convolve raysum array with filter
- * Backproject raysums and add (summate) to image array
- * end
- */
-
-bool
-Projections::reconstruct (ImageFile& im, const char* const filterName, double filt_param, const char* const filterMethodName, const int zeropad, const char* filterGenerationName, const char* const interpName, int interpFactor, const char* const backprojectName, const int trace) const
-{
- double detInc = m_detInc;
- int n_filteredProj = m_nDet * interpFactor;
- double filteredProj [n_filteredProj]; // filtered projections
-
-#ifdef HAVE_BSPLINE_INTERP
- int spline_order = 0, zoom_factor = 0;
- if (interp_type == I_BSPLINE) {
- zoom_factor = interpFactor;
- spline_order = 3;
- zoom_factor = 3;
- n_filteredProj = (m_nDet - 1) * (zoom_factor + 1) + 1;
- }
-#endif
-
- double filterBW = 1. / detInc;
- ProcessSignal processSignal (filterName, filterMethodName, filterBW, m_detInc, m_nDet, filt_param, "spatial", filterGenerationName, zeropad, interpFactor, trace, m_geometry, m_focalLength);
-
- if (processSignal.fail()) {
- sys_error (ERR_SEVERE, "%s [Projections::reconstruct]", processSignal.failMessage().c_str());
- return false;
- }
-
- if (trace)
- cout << "Reconstruct: filter="<<filterName<< ", interp="<<interpName<<", backproject="<<backprojectName<<endl;
-
-#if HAVE_SGP
- int nVecFilter = processSignal.getNFilterPoints();
- double plot_xaxis [nVecFilter]; // array for plotting
-
- if (trace > Trace::TRACE_CONSOLE && nVecFilter > 0) {
- int i;
- double f;
- double filterInc = processSignal.getFilterIncrement();
- for (i = 0, f = processSignal.getFilterMin(); i < nVecFilter; i++, f += filterInc)
- plot_xaxis[i] = f;
-
- if (processSignal.getFilter()) {
- SGPDriver sgpDriver ("Filter Function");
- SGP sgp (sgpDriver);
- EZPlot ezplot (sgp);
-
- ezplot.ezset ("title Filter Response");
- ezplot.addCurve (plot_xaxis, processSignal.getFilter(), nVecFilter);
- ezplot.plot();
- cio_put_str ("Press any key to continue");
- cio_kb_getc ();
- }
- }
- if (trace >= Trace::TRACE_CONSOLE) {
- printf ("nview=%d, ndet=%d, det_start=%.4f, detInc=%.4f\n", m_nView, m_nDet, m_detStart, m_detInc);
- }
-#endif //HAVE_SGP
-
- Backprojector bj (*this, im, backprojectName, interpName, interpFactor);
- if (bj.fail()) {
- sys_error (ERR_SEVERE, "%s [Projections::reconstruct]", bj.failMessage().c_str());
- return false;
- }
-
- for (int iview = 0; iview < m_nView; iview++) {
- if (trace >= Trace::TRACE_CONSOLE)
- cout <<"Reconstructing view " << iview << "(last = " << m_nView - 1 << ")\n";
-
- const DetectorArray& darray = getDetectorArray (iview);
- const DetectorValue* detval = darray.detValues();
-
- processSignal.filterSignal (detval, filteredProj);
-
-#ifdef HAVE_BSPLINE_INTERP
- if (interp_type == I_BSPLINE)
- bspline (m_nDet, zoom_factor, spline_order, filteredProj, filteredProj);
-
-#ifdef HAVE_SGP
- if (trace >= Trace::TRACE_PLOT && interp_type == I_BSPLINE) {
- bspline (m_nDet, zoom_factor, spline_order, filteredProj, filteredProj);
- ezplot_1d (filteredProj, n_filteredProj);
- }
-#endif
-#endif
-
- bj.BackprojectView (filteredProj, darray.viewAngle());
-
-#ifdef HAVE_SGP
- if (trace >= Trace::TRACE_PLOT) {
- SGPDriver sgpDriverProj ("Projection");
- SGP sgpProj (sgpDriverProj);
- EZPlot ezplotProj (sgpProj);
-
- ezplotProj.ezset ("clear");
- ezplotProj.ezset ("title Filtered Projection");
- ezplotProj.ezset ("xticks major 5.");
- ezplotProj.ezset ("xlabel ");
- ezplotProj.ezset ("ylabel ");
- ezplotProj.ezset ("yporigin .5.");
- ezplotProj.ezset ("ylength .5.");
- ezplotProj.ezset ("box.");
- ezplotProj.ezset ("grid.");
- ezplotProj.addCurve (plot_xaxis, detval, m_nDet);
- ezplotProj.plot();
- ezplotProj.ezset ("clear");
- ezplotProj.ezset ("xticks major 5.");
- ezplotProj.ezset ("xlabel ");
- ezplotProj.ezset ("ylabel ");
- ezplotProj.ezset ("ylength .5.");
- ezplotProj.ezset ("box");
- ezplotProj.ezset ("grid");
- ezplotProj.addCurve (plot_xaxis, filteredProj, n_filteredProj);
- ezplotProj.plot();
-
- cout << "Do you want to exit with current pic (y/n)? " << flush;
- char str[256];
- fgets(str, sizeof(str), stdin);
- if (tolower(str[0]) == 'y') {
- break;
- }
- }
-#endif //HAVE_SGP
- }
-
- return true;
-}
-
projects / ctsim.git / blobdiff