X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fprojections.cpp;h=f7ebe704934fb736297891c53d721628f03d66b8;hp=8674647146a8c4f64f160c4323ef7be974b304d9;hb=1a050c98763fbbc0662731b0b76953acede6f5d7;hpb=c8b19dfaffba9f06d8b6c40cb1bb83a8964867f7 diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp index 8674647..f7ebe70 100644 --- a/libctsim/projections.cpp +++ b/libctsim/projections.cpp @@ -1,7 +1,7 @@ /***************************************************************************** ** FILE IDENTIFICATION ** -** Name: projections.cpp Projection data classes +** Name: projections.cpp Projection data classes ** Programmer: Kevin Rosenberg ** Date Started: Aug 84 ** @@ -35,14 +35,14 @@ const int Projections::POLAR_INTERP_NEAREST = 0; const int Projections::POLAR_INTERP_BILINEAR = 1; const int Projections::POLAR_INTERP_BICUBIC = 2; -const char* const Projections::s_aszInterpName[] = +const char* const Projections::s_aszInterpName[] = { "nearest", "bilinear", // {"bicubic"}, }; -const char* const Projections::s_aszInterpTitle[] = +const char* const Projections::s_aszInterpTitle[] = { "Nearest", "Bilinear", @@ -54,13 +54,13 @@ const int Projections::s_iInterpCount = sizeof(s_aszInterpName) / sizeof(char*); /* NAME -* Projections Constructor for projections matrix storage +* Projections Constructor for projections matrix storage * * SYNOPSIS * proj = projections_create (filename, nView, nDet) -* Projections& proj Allocated projections structure & matrix -* int nView Number of rotated view -* int nDet Number of detectors +* Projections& proj Allocated projections structure & matrix +* int nView Number of rotated view +* int nDet Number of detectors * */ @@ -92,13 +92,13 @@ int Projections::convertInterpNameToID (const char* const interpName) { int interpID = POLAR_INTERP_INVALID; - + for (int i = 0; i < s_iInterpCount; i++) if (strcasecmp (interpName, s_aszInterpName[i]) == 0) { interpID = i; break; } - + return (interpID); } @@ -106,10 +106,10 @@ const char* Projections::convertInterpIDToName (const int interpID) { static const char *interpName = ""; - + if (interpID >= 0 && interpID < s_iInterpCount) return (s_aszInterpName[interpID]); - + return (interpName); } @@ -117,10 +117,10 @@ const char* Projections::convertInterpIDToTitle (const int interpID) { static const char *interpTitle = ""; - + if (interpID >= 0 && interpID < s_iInterpCount) return (s_aszInterpTitle[interpID]); - + return (interpTitle); } @@ -133,7 +133,7 @@ Projections::init (const int nView, const int nDet) m_nView = nView; m_nDet = nDet; newProjData (); - + time_t t = time (NULL); tm* lt = localtime (&t); m_year = lt->tm_year; @@ -150,7 +150,7 @@ Projections::initFromScanner (const Scanner& scanner) m_label.setLabelType (Array2dFileLabel::L_HISTORY); deleteProjData(); init (scanner.nView(), scanner.nDet()); - + m_rotInc = scanner.rotInc(); m_detInc = scanner.detInc(); m_detStart = scanner.detStart(); @@ -170,56 +170,56 @@ Projections::setNView (int nView) // used by MPI to reduce # of views } // Helical 180 Linear Interpolation. -// This member function takes a set of helical scan projections and -// performs a linear interpolation between pairs of complementary rays +// This member function takes a set of helical scan projections and +// performs a linear interpolation between pairs of complementary rays // to produce a single projection data set approximating what would be // measured at a single axial plane. -// Complementary rays are rays which traverse the same path through the +// Complementary rays are rays which traverse the same path through the // phantom in opposite directions. // // For parallel beam geometry, a ray with a given gantry angle beta and a // detector iDet will have a complementary ray at beta + pi and nDet-iDet // // For equiangular or equilinear beam geometry the complementary ray to -// gantry angle beta and fan-beam angle gamma is at +// gantry angle beta and fan-beam angle gamma is at // beta-hat = beta +2*gamma + pi, and gamma-hat = -gamma. // Note that beta-hat - beta depends on gamma and is not constant. // // The algorithm used here is from Crawford and King, Med. Phys. 17(6) // 1990 p967; what they called method "C", CSH-HH. It uses interpolation only // between pairs of complementary rays on either side of an image plane. -// Input data must sample gantry angles from zero to +// Input data must sample gantry angles from zero to // (2*pi + 2* fan-beam-angle). The data set produced contains gantry // angles from 0 to Pi+fan-beam-angle. This is a "halfscan" data set, -// which still contains redundant data, and can be used with a half scan +// which still contains redundant data, and can be used with a half scan // reconstruction to produce an image. // In this particular implementation a lower triangle from (beta,gamma) = // (0,-fanAngle/2)->(2*fanAngle,-fanAngle/2)->(0,fanAngle/2) contains // zeros, but is actually redundant with data contained in the region // (pi+fanAngle,-fanAngle/2)->(pi+fanAngle, fanAngle/2) ->(pi-fanAngle, -// fanAngle/2). +// fanAngle/2). // -int +int Projections::Helical180LI(int interpolation_view) { - if (m_geometry == Scanner::GEOMETRY_INVALID) + if (m_geometry == Scanner::GEOMETRY_INVALID) { std::cerr << "Invalid geometry " << m_geometry << std::endl; return (2); - } - else if (m_geometry == Scanner::GEOMETRY_PARALLEL) + } + else if (m_geometry == Scanner::GEOMETRY_PARALLEL) { std::cerr << "Helical 180LI not yet implemented for PARALLEL geometry" << std::endl; return (2); } - else if (m_geometry == Scanner::GEOMETRY_EQUILINEAR) + else if (m_geometry == Scanner::GEOMETRY_EQUILINEAR) { std::cerr << "Helical 180LI not yet implemented for EQUILINEAR geometry" << std::endl; return (2); } - else if (m_geometry == Scanner::GEOMETRY_EQUIANGULAR) + else if (m_geometry == Scanner::GEOMETRY_EQUIANGULAR) { return Helical180LI_Equiangular(interpolation_view); } @@ -233,11 +233,11 @@ Projections::Helical180LI(int interpolation_view) int Projections::Helical180LI_Equiangular(int interpView) { - double dbeta = m_rotInc; - double dgamma = m_detInc; + double dbeta = m_rotInc; + double dgamma = m_detInc; double fanAngle = m_dFanBeamAngle; int offsetView=0; - + // is there enough data in the data set? Should have 2(Pi+fanAngle) // coverage minimum if ( m_nView < static_cast((2*( PI + fanAngle ) ) / dbeta) -1 ){ @@ -252,10 +252,10 @@ Projections::Helical180LI_Equiangular(int interpView) } else { - // check if there is PI+fanAngle data on either side of the + // check if there is PI+fanAngle data on either side of the // of the specified image plane - if ( interpView*dbeta < PI+fanAngle || - interpView*dbeta + PI + fanAngle > m_nView*dbeta) + if ( interpView*dbeta < PI+fanAngle || + interpView*dbeta + PI + fanAngle > m_nView*dbeta) { std::cerr << "There isn't PI+fanAngle of data on either side of the requested interpolation view" << std::endl; return(1); @@ -265,7 +265,7 @@ Projections::Helical180LI_Equiangular(int interpView) } int last_interp_view = static_cast ((PI+fanAngle)/dbeta); - + // make a new array for data... class DetectorArray ** newdetarray = new DetectorArray * [last_interp_view+1]; for ( int i=0 ; i <= last_interp_view ; i++ ){ @@ -273,22 +273,22 @@ Projections::Helical180LI_Equiangular(int interpView) newdetarray[i]->setViewAngle((i+offsetView)*dbeta); DetectorValue* newdetval = (newdetarray[i])->detValues(); // and initialize the data to zero - for (int j=0; j < m_nDet; j++) + for (int j=0; j < m_nDet; j++) newdetval[j] = 0.; } int last_acq_view = 2*last_interp_view; for ( int iView = 0 ; iView <= last_acq_view; iView++) { - double beta = iView * dbeta; - + double beta = iView * dbeta; + for ( int iDet = 0; iDet < m_nDet; iDet++) { double gamma = (iDet -(m_nDet-1)/2)* dgamma ; int newiView, newiDet; - if (beta < PI+fanAngle) { //if (PI +fanAngle - beta > dbeta ) - //newbeta = beta; - //newgamma = gamma; - newiDet = iDet; - newiView = iView; + if (beta < PI+fanAngle) { //if (PI +fanAngle - beta > dbeta ) + //newbeta = beta; + //newgamma = gamma; + newiDet = iDet; + newiView = iView; } else // (beta > PI+fanAngle) { @@ -298,55 +298,55 @@ Projections::Helical180LI_Equiangular(int interpView) // newiView = nearest((beta + 2*gamma - PI)/dbeta); //newiView = static_cast(( (iView*dbeta) + 2*(iDet-(m_nDet-1)/2)*dgamma - PI)/dbeta); newiView = nearest(( (iView*dbeta) + 2*(iDet-(m_nDet-1)/2)*dgamma - PI)/dbeta); - } + } #ifdef DEBUG //std::cout << beta << " "<< gamma << " " << newbeta << " " << newgamma <<" " << iView-offsetView << " " << iDet << " " << newiView << " " << newiDet << std::endl; //std::cout << iView-offsetView << " " << iDet << " " << newiView << " " << newiDet << std::endl; #endif - if ( ( beta > fanAngle - 2*gamma) + if ( ( beta > fanAngle - 2*gamma) && ( beta < 2*PI + fanAngle -2*gamma) ) { // not in region 1 or 8 DetectorValue* detval = (m_projData[iView+offsetView])->detValues(); DetectorValue* newdetval = (newdetarray[newiView])->detValues(); - if ( beta > fanAngle - 2*gamma + if ( beta > fanAngle - 2*gamma && beta <= 2*fanAngle ) { // in region 2 - newdetval[newiDet] += + newdetval[newiDet] += (beta +2*gamma - fanAngle)/(PI+2*gamma) * detval[iDet]; - } else if ( beta > 2*fanAngle + } else if ( beta > 2*fanAngle && beta <= PI - 2*gamma) { // in region 3 - newdetval[newiDet] += + newdetval[newiDet] += (beta +2*gamma - fanAngle)/(PI+2*gamma) * detval[iDet]; - } - else if ( beta > PI -2*gamma + } + else if ( beta > PI -2*gamma && beta <= PI + fanAngle ) { // in region 4 - newdetval[newiDet] += + newdetval[newiDet] += (beta +2*gamma - fanAngle)/(PI+2*gamma) * detval[iDet]; - } - else if ( beta > PI + fanAngle + } + else if ( beta > PI + fanAngle && beta <= PI +2*fanAngle -2*gamma) { // in region 5 - newdetval[newiDet] += + newdetval[newiDet] += (2*PI - beta - 2*gamma + fanAngle)/(PI-2*gamma) *detval[iDet]; - } - else if ( beta > PI +2*fanAngle -2*gamma + } + else if ( beta > PI +2*fanAngle -2*gamma && beta <= 2*PI) { // in region 6 - newdetval[newiDet] += + newdetval[newiDet] += (2*PI - beta - 2*gamma + fanAngle)/(PI-2*gamma) *detval[iDet]; - } - else if ( beta > 2*PI + } + else if ( beta > 2*PI && beta <= 2*PI + fanAngle -2*gamma){ // in region 7 - newdetval[newiDet] += + newdetval[newiDet] += (2*PI - beta -2*gamma + fanAngle)/(PI-2*gamma) *detval[iDet]; - } - else - { + } + else + { ; // outside region of interest } } @@ -356,26 +356,26 @@ Projections::Helical180LI_Equiangular(int interpView) m_projData = newdetarray; m_nView = last_interp_view+1; - return (0); + return (0); } // HalfScanFeather: -// A HalfScan Projection Data Set for equiangular geometry, -// covering gantry angles from 0 to pi+fanBeamAngle +// A HalfScan Projection Data Set for equiangular geometry, +// covering gantry angles from 0 to pi+fanBeamAngle // and fan angle gamma from -fanBeamAngle/2 to fanBeamAngle/2 // contains redundant information. If one copy of this data is left as -// zero, (as in the Helical180LI routine above) overweighting is avoided, -// but the discontinuity in the data introduces ringing in the image. +// zero, (as in the Helical180LI routine above) overweighting is avoided, +// but the discontinuity in the data introduces ringing in the image. // This routine makes a copy of the data and applies a weighting to avoid // over-representation, as given in Appendix C of Crawford and King, Med // Phys 17 1990, p967. int Projections::HalfScanFeather(void) { - double dbeta = m_rotInc; - double dgamma = m_detInc; + double dbeta = m_rotInc; + double dgamma = m_detInc; double fanAngle = m_dFanBeamAngle; -// is there enough data? +// is there enough data? if ( m_nView != static_cast(( PI+fanAngle ) / dbeta) +1 ){ std::cerr << "Data set does seem have enough data to be a halfscan data set" << std::endl; return (1); @@ -391,14 +391,14 @@ Projections::HalfScanFeather(void) } for ( int iView2 = 0 ; iView2 < m_nView; iView2++) { - double beta2 = iView2 * dbeta; + double beta2 = iView2 * dbeta; for ( int iDet2 = 0; iDet2 < m_nDet; iDet2++) { double gamma2 = (iDet2 -(m_nDet-1)/2)* dgamma ; - if ( ( beta2 >= PI - 2*gamma2) ) { // in redundant data region + if ( ( beta2 >= PI - 2*gamma2) ) { // in redundant data region int iView1, iDet1; iDet1 = (m_nDet -1) - iDet2; //iView1 = nearest((beta2 + 2*gamma2 - PI)/dbeta); - iView1 = nearest(( (iView2*dbeta) + iView1 = nearest(( (iView2*dbeta) + 2*(iDet2-(m_nDet-1)/2)*dgamma - PI)/dbeta); @@ -408,13 +408,13 @@ Projections::HalfScanFeather(void) detval1[iDet1] = detval2[iDet2] ; double x, w1,w2,beta1, gamma1; - beta1= iView1*dbeta; + beta1= iView1*dbeta; gamma1 = -gamma2; if ( beta1 <= (fanAngle - 2*gamma1) ) x = beta1 / ( fanAngle - 2*gamma1); - else if ( (fanAngle - 2*gamma1 <= beta1 ) && beta1 <= PI - 2*gamma1) - x = 1; - else if ( (PI - 2*gamma1 <= beta1 ) && ( beta1 <=PI + fanAngle) ) + else if ( (fanAngle - 2*gamma1 <= beta1 ) && beta1 <= PI - 2*gamma1) + x = 1; + else if ( (PI - 2*gamma1 <= beta1 ) && ( beta1 <=PI + fanAngle) ) x = (PI +fanAngle - beta1)/(fanAngle + 2*gamma1); else { std::cerr << "Shouldn't be here!"<< std::endl; @@ -422,13 +422,13 @@ Projections::HalfScanFeather(void) } w1 = (3*x - 2*x*x)*x; w2 = 1-w1; - detval1[iDet1] *= w1; + detval1[iDet1] *= w1; detval2[iDet2] *= w2; - } + } } } - // heuristic scaling, why this factor? + // heuristic scaling, why this factor? double scalefactor = m_nView * m_rotInc / PI; for ( int iView = 0 ; iView < m_nView; iView++) { DetectorValue* detval = (m_projData[iView])->detValues(); @@ -437,21 +437,21 @@ Projections::HalfScanFeather(void) } } - return (0); + return (0); } // NAME // newProjData -void +void Projections::newProjData (void) { if (m_projData) sys_error(ERR_WARNING, "m_projData != NULL [newProjData]"); - + if (m_nView > 0 && m_nDet) { m_projData = new DetectorArray* [m_nView]; - + for (int i = 0; i < m_nView; i++) m_projData[i] = new DetectorArray (m_nDet); } @@ -459,20 +459,20 @@ Projections::newProjData (void) /* NAME -* projections_free Free memory allocated to projections +* projections_free Free memory allocated to projections * * SYNOPSIS * projections_free(proj) -* Projections& proj Projectionss to be deallocated +* Projections& proj Projectionss to be deallocated */ -void +void Projections::deleteProjData (void) { if (m_projData != NULL) { for (int i = 0; i < m_nView; i++) delete m_projData[i]; - + delete m_projData; m_projData = NULL; } @@ -484,7 +484,7 @@ Projections::deleteProjData (void) * */ -bool +bool Projections::headerWrite (fnetorderstream& fs) { kuint16 _hsize = m_headerSize; @@ -499,7 +499,7 @@ Projections::headerWrite (fnetorderstream& fs) kuint16 _hour = m_hour; kuint16 _minute = m_minute; kuint16 _second = m_second; - + kfloat64 _calcTime = m_calcTime; kfloat64 _rotStart = m_rotStart; kfloat64 _rotInc = m_rotInc; @@ -513,7 +513,7 @@ Projections::headerWrite (fnetorderstream& fs) fs.seekp(0); if (! fs) return false; - + fs.writeInt16 (_hsize); fs.writeInt16 (_signature); fs.writeInt32 (_nView); @@ -536,14 +536,14 @@ Projections::headerWrite (fnetorderstream& fs) fs.writeInt16 (_second); fs.writeInt16 (_remarksize); fs.write (m_remark.c_str(), _remarksize); - + m_headerSize = fs.tellp(); _hsize = m_headerSize; fs.seekp(0); fs.writeInt16 (_hsize); if (! fs) return false; - + return true; } @@ -557,11 +557,11 @@ Projections::headerRead (fnetorderstream& fs) kuint16 _hsize, _signature, _year, _month, _day, _hour, _minute, _second, _remarksize = 0; kuint32 _nView, _nDet, _geom; kfloat64 _calcTime, _rotStart, _rotInc, _detStart, _detInc, _focalLength, _sourceDetectorLength, _viewDiameter, _fanBeamAngle; - + fs.seekg(0); if (! fs) return false; - + fs.readInt16 (_hsize); fs.readInt16 (_signature); fs.readInt32 (_nView); @@ -583,17 +583,17 @@ Projections::headerRead (fnetorderstream& fs) fs.readInt16 (_minute); fs.readInt16 (_second); fs.readInt16 (_remarksize); - + if (! fs) { sys_error (ERR_SEVERE, "Error reading header information , _remarksize=%d [projections_read_header]", _remarksize); return false; } - + if (_signature != m_signature) { sys_error (ERR_SEVERE, "File %s does not have a valid projection file signature", m_filename.c_str()); return false; } - + char* pszRemarkStorage = new char [_remarksize+1]; fs.read (pszRemarkStorage, _remarksize); if (! fs) { @@ -603,13 +603,13 @@ Projections::headerRead (fnetorderstream& fs) pszRemarkStorage[_remarksize] = 0; m_remark = pszRemarkStorage; delete pszRemarkStorage; - + off_t _hsizeread = fs.tellg(); if (!fs || _hsizeread != _hsize) { 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); return false; } - + m_headerSize = _hsize; m_nView = _nView; m_nDet = _nDet; @@ -629,12 +629,12 @@ Projections::headerRead (fnetorderstream& fs) m_hour = _hour; m_minute = _minute; m_second = _second; - + m_label.setLabelType (Array2dFileLabel::L_HISTORY); m_label.setLabelString (m_remark); m_label.setCalcTime (m_calcTime); m_label.setDateTime (m_year, m_month, m_day, m_hour, m_minute, m_second); - + return true; } @@ -654,39 +654,39 @@ Projections::read (const char* filename) #else frnetorderstream fileRead (m_filename.c_str(), std::ios::in | std::ios::binary); // | std::ios::nocreate); #endif - + if (fileRead.fail()) return false; - + if (! headerRead (fileRead)) return false; - + deleteProjData (); newProjData(); - + for (int i = 0; i < m_nView; i++) { if (! detarrayRead (fileRead, *m_projData[i], i)) break; } - + fileRead.close(); return true; } -bool +bool Projections::copyViewData (const std::string& filename, std::ostream& os, int startView, int endView) { return copyViewData (filename.c_str(), os, startView, endView); } -bool +bool Projections::copyViewData (const char* const filename, std::ostream& os, int startView, int endView) { frnetorderstream is (filename, std::ios::in | std::ios::binary); kuint16 sizeHeader, signature; kuint32 _nView, _nDet; - + is.seekg (0); if (is.fail()) { sys_error (ERR_SEVERE, "Unable to read projection file %s", filename); @@ -699,28 +699,28 @@ Projections::copyViewData (const char* const filename, std::ostream& os, int sta is.readInt32 (_nDet); int nView = _nView; int nDet = _nDet; - + if (signature != m_signature) { sys_error (ERR_SEVERE, "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 (reinterpret_cast(pViewData), sizeView); @@ -728,17 +728,17 @@ Projections::copyViewData (const char* const filename, std::ostream& os, int sta if (is.fail() || os.fail()) break; } - + delete pViewData; - if (is.fail()) + if (is.fail()) sys_error (ERR_SEVERE, "Error reading projection file"); - if (os.fail()) + if (os.fail()) sys_error (ERR_SEVERE, "Error writing projection file"); - + return (! (is.fail() | os.fail())); } -bool +bool Projections::copyHeader (const std::string& filename, std::ostream& os) { return copyHeader (filename.c_str(), os); @@ -756,20 +756,20 @@ Projections::copyHeader (const char* const filename, std::ostream& os) sys_error (ERR_SEVERE, "Illegal signature in projection file %s", filename); return false; } - + unsigned char* pHdrData = new unsigned char [sizeHeader]; is.read (reinterpret_cast(pHdrData), sizeHeader); if (is.fail()) { sys_error (ERR_SEVERE, "Error reading header"); return false; } - + os.write (reinterpret_cast(pHdrData), sizeHeader); if (os.fail()) { sys_error (ERR_SEVERE, "Error writing header"); return false; } - + return true; } @@ -788,10 +788,10 @@ Projections::write (const char* filename) sys_error (ERR_SEVERE, "Error opening file %s for output [projections_create]", filename); return false; } - + if (! headerWrite (fs)) return false; - + if (m_projData != NULL) { for (int i = 0; i < m_nView; i++) { if (! detarrayWrite (fs, *m_projData[i], i)) @@ -800,19 +800,19 @@ Projections::write (const char* filename) } if (! fs) return false; - + fs.close(); - + return true; } /* NAME -* detarrayRead Read a Detector Array structure from the disk +* detarrayRead Read a Detector Array structure from the disk * * SYNOPSIS * detarrayRead (proj, darray, view_num) -* DETARRAY *darray Detector array storage location to be filled -* int view_num View number to read +* DETARRAY *darray Detector array storage location to be filled +* int view_num View number to read */ bool @@ -822,17 +822,17 @@ Projections::detarrayRead (fnetorderstream& fs, DetectorArray& darray, const int const int detheader_bytes = sizeof(kfloat64) /* view_angle */ + sizeof(kint32) /* nDet */; const int view_bytes = detheader_bytes + detval_bytes; const off_t start_data = m_headerSize + (iview * view_bytes); - DetectorValue* detval_ptr = darray.detValues(); + DetectorValue* detval_ptr = darray.detValues(); kfloat64 view_angle; kuint32 nDet; - + fs.seekg (start_data); - + fs.readFloat64 (view_angle); fs.readInt32 (nDet); darray.setViewAngle (view_angle); // darray.setNDet ( nDet); - + for (unsigned int i = 0; i < nDet; i++) { kfloat32 detval; fs.readFloat32 (detval); @@ -840,18 +840,18 @@ Projections::detarrayRead (fnetorderstream& fs, DetectorArray& darray, const int } if (! fs) return false; - + return true; } /* NAME -* detarrayWrite Write detector array data to the disk +* detarrayWrite Write detector array data to the disk * * SYNOPSIS * detarrayWrite (darray, view_num) -* DETARRAY *darray Detector array data to be written -* int view_num View number to write +* DETARRAY *darray Detector array data to be written +* int view_num View number to write * * DESCRIPTION * This routine writes the detarray data from the disk sequentially to @@ -866,32 +866,32 @@ Projections::detarrayWrite (fnetorderstream& fs, const DetectorArray& darray, co const int detheader_bytes = sizeof(kfloat64) /* view_angle */ + sizeof(kint32) /* nDet */; const int view_bytes = detheader_bytes + detval_bytes; const off_t start_data = m_headerSize + (iview * view_bytes); - const DetectorValue* const detval_ptr = darray.detValues(); + const DetectorValue* const detval_ptr = darray.detValues(); kfloat64 view_angle = darray.viewAngle(); kuint32 nDet = darray.nDet(); - + fs.seekp (start_data); if (! fs) { sys_error (ERR_SEVERE, "Error seeking detectory array [detarrayWrite]"); return false; } - + fs.writeFloat64 (view_angle); fs.writeInt32 (nDet); - + for (unsigned int i = 0; i < nDet; i++) { kfloat32 detval = detval_ptr[i]; fs.writeFloat32 (detval); } - + if (! fs) return (false); - + return true; } /* NAME -* printProjectionData Print projections data +* printProjectionData Print projections data * * SYNOPSIS * printProjectionData () @@ -933,7 +933,7 @@ Projections::printProjectionData (int startView, int endView) } } -void +void Projections::printScanInfo (std::ostringstream& os) const { os << "Number of detectors: " << m_nDet << "\n"; @@ -951,7 +951,7 @@ Projections::printScanInfo (std::ostringstream& os) const } -bool +bool Projections::convertPolar (ImageFile& rIF, int iInterpolationID) { unsigned int nx = rIF.nx(); @@ -965,7 +965,7 @@ Projections::convertPolar (ImageFile& rIF, int iInterpolationID) Projections* pProj = this; if (m_geometry == Scanner::GEOMETRY_EQUIANGULAR || m_geometry == Scanner::GEOMETRY_EQUILINEAR) pProj = interpolateToParallel(); - + Array2d adView (nx, ny); Array2d adDet (nx, ny); double** ppdView = adView.getArray(); @@ -982,7 +982,7 @@ Projections::convertPolar (ImageFile& rIF, int iInterpolationID) pProj->calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet, pProj->m_nDet, 1., pProj->m_detInc); - pProj->interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, pProj->m_nView, pProj->m_nDet, + pProj->interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, pProj->m_nView, pProj->m_nDet, pProj->m_nDet, iInterpolationID); for (iView = 0; iView < pProj->m_nView; iView++) @@ -996,7 +996,7 @@ Projections::convertPolar (ImageFile& rIF, int iInterpolationID) } -bool +bool Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad) { #ifndef HAVE_FFTW @@ -1012,7 +1012,7 @@ Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad if (! v || nx == 0 || ny == 0) return false; - + Projections* pProj = this; if (m_geometry == Scanner::GEOMETRY_EQUIANGULAR || m_geometry == Scanner::GEOMETRY_EQUILINEAR) pProj = interpolateToParallel(); @@ -1028,7 +1028,7 @@ Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad std::complex** ppcDetValue = new std::complex* [pProj->m_nView]; //double dInterpScale = (pProj->m_nDet-1) / static_cast(iInterpDet-1); double dInterpScale = pProj->m_nDet / static_cast(iInterpDet); - + double dFFTScale = 1. / static_cast(iInterpDet * iInterpDet); int iMidPoint = iInterpDet / 2; double dMidPoint = static_cast(iInterpDet) / 2.; @@ -1050,7 +1050,7 @@ Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad pcIn[iDet1+iZerosAdded][0] = pcIn[iDet1][0]; pcIn[iDet1+iZerosAdded][1] = pcIn[iDet1][1]; } - for (int iDet2 = iMidPoint; iDet2 < iMidPoint + iZerosAdded; iDet2++) + for (int iDet2 = iMidPoint; iDet2 < iMidPoint + iZerosAdded; iDet2++) pcIn[iDet2][0] = pcIn[iDet2][1] = 0; } @@ -1058,23 +1058,23 @@ Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad ppcDetValue[iView] = new std::complex [iNumInterpDetWithZeros]; for (int iD = 0; iD < iNumInterpDetWithZeros; iD++) { - ppcDetValue[iView][iD] = std::complex (pcIn[iD][0] * dFFTScale, pcIn[iD][1] * dFFTScale); + ppcDetValue[iView][iD] = std::complex (pcIn[iD][0] * dFFTScale, pcIn[iD][1] * dFFTScale); } Fourier::shuffleFourierToNaturalOrder (ppcDetValue[iView], iNumInterpDetWithZeros); } fftw_free(pcIn) ; - fftw_destroy_plan (plan); - + fftw_destroy_plan (plan); + Array2d adView (nx, ny); Array2d adDet (nx, ny); double** ppdView = adView.getArray(); double** ppdDet = adDet.getArray(); - pProj->calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet, iNumInterpDetWithZeros, dZeropadRatio, + pProj->calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet, iNumInterpDetWithZeros, dZeropadRatio, pProj->m_detInc * dInterpScale); - pProj->interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, pProj->m_nView, pProj->m_nDet, + pProj->interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, pProj->m_nView, pProj->m_nDet, iNumInterpDetWithZeros, iInterpolationID); if (m_geometry == Scanner::GEOMETRY_EQUIANGULAR || m_geometry == Scanner::GEOMETRY_EQUILINEAR) @@ -1106,16 +1106,16 @@ Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double yMin = (yMin - yCent) * dZeropadRatio + yCent; yMax = (yMax - yCent) * dZeropadRatio + yCent; - double xInc = (xMax - xMin) / nx; // size of cells + double xInc = (xMax - xMin) / nx; // size of cells double yInc = (yMax - yMin) / ny; - double dDetCenter = (iNumDetWithZeros - 1) / 2.; // index refering to L=0 projection + double dDetCenter = (iNumDetWithZeros - 1) / 2.; // index refering to L=0 projection // +1 is correct for frequency data, ndet-1 is correct for projections // if (isEven (iNumDetWithZeros)) - // dDetCenter = (iNumDetWithZeros + 0) / 2; + // dDetCenter = (iNumDetWithZeros + 0) / 2; // Calculates polar coordinates (view#, det#) for each point on phantom grid - double x = xMin + xInc / 2; // Rectang coords of center of pixel + double x = xMin + xInc / 2; // Rectang coords of center of pixel for (unsigned int ix = 0; ix < nx; x += xInc, ix++) { double y = yMin + yInc / 2; for (unsigned int iy = 0; iy < ny; y += yInc, iy++) { @@ -1128,7 +1128,7 @@ Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double phi -= PI; r = -r; } - + ppdView[ix][iy] = (phi - m_rotStart) / m_rotInc; ppdDet[ix][iy] = (r / dDetInc) + dDetCenter; } @@ -1137,25 +1137,25 @@ Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double void Projections::interpolatePolar (ImageFileArray& v, ImageFileArray& vImag, - unsigned int nx, unsigned int ny, std::complex** ppcDetValue, double** ppdView, + unsigned int nx, unsigned int ny, std::complex** ppcDetValue, double** ppdView, double** ppdDet, unsigned int nView, unsigned int nDet, unsigned int nDetWithZeros, int iInterpolationID) { typedef std::complex complexValue; - BilinearPolarInterpolator* pBilinear = NULL; - BicubicPolyInterpolator* pBicubic = NULL; + BilinearPolarInterpolator* pBilinear = NULL; + BicubicPolyInterpolator* pBicubic = NULL; if (iInterpolationID == POLAR_INTERP_BILINEAR) pBilinear = new BilinearPolarInterpolator (ppcDetValue, nView, nDetWithZeros); else if (iInterpolationID == POLAR_INTERP_BICUBIC) pBicubic = new BicubicPolyInterpolator (ppcDetValue, nView, nDetWithZeros); - + for (unsigned int ix = 0; ix < ny; ix++) { for (unsigned int iy = 0; iy < ny; iy++) { if (iInterpolationID == POLAR_INTERP_NEAREST) { unsigned int iView = nearest (ppdView[ix][iy]); unsigned int iDet = nearest (ppdDet[ix][iy]); if (iView == nView) - iView = 0; + iView = 0; if (iDet >= 0 && iDet < nDetWithZeros && iView >= 0 && iView < nView) { v[ix][iy] = ppcDetValue[iView][iDet].real(); if (vImag) @@ -1164,7 +1164,7 @@ Projections::interpolatePolar (ImageFileArray& v, ImageFileArray& vImag, v[ix][iy] = 0; if (vImag) vImag[ix][iy] = 0; - } + } } else if (iInterpolationID == POLAR_INTERP_BILINEAR) { std::complex vInterp = pBilinear->interpolate (ppdView[ix][iy], ppdDet[ix][iy]); @@ -1195,7 +1195,7 @@ Projections::initFromSomatomAR_STAR (int iNViews, int iNDets, unsigned char* pDa m_rotStart = 0; m_dViewDiameter = sin (m_dFanBeamAngle / 2) * m_dFocalLength * 2; - if (! ((iNViews == 750 && lDataLength == 1560000L) || (iNViews == 950 && lDataLength == 1976000L) + if (! ((iNViews == 750 && lDataLength == 1560000L) || (iNViews == 950 && lDataLength == 1976000L) || (iNViews == 1500 && lDataLength == 3120000))) return false; @@ -1368,12 +1368,12 @@ ParallelRaysums::ParallelRaysums (const Projections* pProjections, int iThetaRan } else if (iGeometry == Scanner::GEOMETRY_EQUILINEAR) { double dFanAngle = atan (dDetPos / pProjections->sourceDetectorLength()); pC->m_dTheta = dViewAngle + dFanAngle; - pC->m_dT = dFocalLength * sin(dFanAngle); + pC->m_dT = dFocalLength * sin(dFanAngle); } else if (iGeometry == Scanner::GEOMETRY_EQUIANGULAR) { // fan angle is same as dDetPos pC->m_dTheta = dViewAngle + dDetPos; - pC->m_dT = dFocalLength * sin (dDetPos); + pC->m_dT = dFocalLength * sin (dDetPos); } if (m_iThetaRange != THETA_RANGE_UNCONSTRAINED) { pC->m_dTheta = normalizeAngle (pC->m_dTheta);