X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fscanner.cpp;h=9f4837da9b6fec2e1beee96f795a03a4af3c4f62;hp=453e913f99feff3a30c17a635880512e89b1bdbc;hb=8a7697ce57b56cdc43698cd1241ad98d49f9b5ac;hpb=ec73e3f9aa11aea5ba6a4255489fdca63d2b4990 diff --git a/libctsim/scanner.cpp b/libctsim/scanner.cpp index 453e913..9f4837d 100644 --- a/libctsim/scanner.cpp +++ b/libctsim/scanner.cpp @@ -7,9 +7,9 @@ ** Date Started: 1984 ** ** This is part of the CTSim program -** Copyright (C) 1983-2000 Kevin Rosenberg +** Copyright (c) 1983-2001 Kevin Rosenberg ** -** $Id: scanner.cpp,v 1.25 2001/01/17 11:00:18 kevin Exp $ +** $Id$ ** ** 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 @@ -32,19 +32,22 @@ const int Scanner::GEOMETRY_INVALID = -1; const int Scanner::GEOMETRY_PARALLEL = 0; const int Scanner::GEOMETRY_EQUIANGULAR = 1; const int Scanner::GEOMETRY_EQUILINEAR = 2; +const int Scanner::GEOMETRY_LINOGRAM = 3; -const char* Scanner::s_aszGeometryName[] = +const char* const Scanner::s_aszGeometryName[] = { - {"parallel"}, - {"equiangular"}, - {"equilinear"}, + "parallel", + "equiangular", + "equilinear", + "linogram", }; -const char* Scanner::s_aszGeometryTitle[] = +const char* const Scanner::s_aszGeometryTitle[] = { - {"Parallel"}, - {"Equiangular"}, - {"Equilinear"}, + "Parallel", + "Equiangular", + "Equilinear", + "Linogram", }; const int Scanner::s_iGeometryCount = sizeof(s_aszGeometryName) / sizeof(const char*); @@ -82,10 +85,13 @@ DetectorArray::~DetectorArray (void) * int nSample Number of rays per detector */ -Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet, int nView, int nSample, const double rot_anglen, const double dFocalLengthRatio, const double dFieldOfViewRatio) +Scanner::Scanner (const Phantom& phm, const char* const geometryName, + int nDet, int nView, int offsetView, + int nSample, const double rot_anglen, + const double dFocalLengthRatio, + const double dCenterDetectorRatio, + const double dViewRatio, const double dScanRatio) { - m_phmLen = phm.maxAxisLength(); // maximal length along an axis - m_fail = false; m_idGeometry = convertGeometryNameToID (geometryName); if (m_idGeometry == GEOMETRY_INVALID) { @@ -105,88 +111,108 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet, m_nDet = nDet; m_nView = nView; + m_iOffsetView = offsetView; m_nSample = nSample; m_dFocalLengthRatio = dFocalLengthRatio; - m_dFieldOfViewRatio = dFieldOfViewRatio; - m_dFocalLength = (m_phmLen * SQRT2 / 2) * dFocalLengthRatio; - m_dFieldOfView = m_phmLen * SQRT2 * dFieldOfViewRatio; + m_dCenterDetectorRatio = dCenterDetectorRatio; + m_dViewRatio = dViewRatio; + m_dScanRatio = dScanRatio; + + m_dViewDiameter = phm.getDiameterBoundaryCircle() * m_dViewRatio; + m_dFocalLength = (m_dViewDiameter / 2) * m_dFocalLengthRatio; + m_dCenterDetectorLength = (m_dViewDiameter / 2) * m_dCenterDetectorRatio; + m_dSourceDetectorLength = m_dFocalLength + m_dCenterDetectorLength; + m_dScanDiameter = m_dViewDiameter * m_dScanRatio; m_dXCenter = phm.xmin() + (phm.xmax() - phm.xmin()) / 2; m_dYCenter = phm.ymin() + (phm.ymax() - phm.ymin()) / 2; m_rotLen = rot_anglen; m_rotInc = m_rotLen / m_nView; if (m_idGeometry == GEOMETRY_PARALLEL) { - m_detLen = m_dFieldOfView; + m_dFanBeamAngle = 0; + m_detLen = m_dScanDiameter; + m_detStart = -m_detLen / 2; m_detInc = m_detLen / m_nDet; - if (m_nDet % 2 == 0) // Adjust for Even number of detectors - m_detInc = m_detLen / (m_nDet - 1); // center detector = (nDet/2)-1 - + double dDetectorArrayEndOffset = 0; + // For even number of detectors, make detInc slightly larger so that center lies + // at nDet/2. Also, extend detector array by one detInc so that all of the phantom is scanned + if (isEven (m_nDet)) { // Adjust for Even number of detectors + m_detInc = m_detLen / (m_nDet - 1); // center detector = (nDet/2) + dDetectorArrayEndOffset = m_detInc; + } + double dHalfDetLen = m_detLen / 2; m_initPos.xs1 = m_dXCenter - dHalfDetLen; m_initPos.ys1 = m_dYCenter + m_dFocalLength; - m_initPos.xs2 = m_dXCenter + dHalfDetLen; + m_initPos.xs2 = m_dXCenter + dHalfDetLen + dDetectorArrayEndOffset; m_initPos.ys2 = m_dYCenter + m_dFocalLength; m_initPos.xd1 = m_dXCenter - dHalfDetLen; - m_initPos.yd1 = m_dYCenter - m_dFocalLength; - m_initPos.xd2 = m_dXCenter + dHalfDetLen; - m_initPos.yd2 = m_dYCenter - m_dFocalLength; - m_initPos.angle = 0.0; + m_initPos.yd1 = m_dYCenter - m_dCenterDetectorLength; + m_initPos.xd2 = m_dXCenter + dHalfDetLen + dDetectorArrayEndOffset; + m_initPos.yd2 = m_dYCenter - m_dCenterDetectorLength; + m_initPos.angle = m_iOffsetView * m_rotInc; + m_detLen += dDetectorArrayEndOffset; } else if (m_idGeometry == GEOMETRY_EQUILINEAR) { -#if 0 - double dAngle = (m_dFieldOfView / 2) / cos (asin (m_dFieldOfView / 2 / m_dFocalLength)); -#else - double dHalfSquare = m_dFieldOfView / SQRT2 / 2; - double dFocalPastPhm = m_dFocalLength - dHalfSquare; - if (dFocalPastPhm <= 0.) { + if (m_dScanDiameter / 2 >= m_dFocalLength) { m_fail = true; - m_failMessage = "Focal Point inside of phantom"; + m_failMessage = "Invalid geometry: Focal length must be larger than scan length"; return; } - double dAngle = atan( dHalfSquare / dFocalPastPhm ); -#endif - double dHalfDetLen = 2 * m_dFocalLength * tan (dAngle); + + const double dAngle = asin ((m_dScanDiameter / 2) / m_dFocalLength); + const double dHalfDetLen = m_dSourceDetectorLength * tan (dAngle); m_detLen = dHalfDetLen * 2; + m_detStart = -dHalfDetLen; m_detInc = m_detLen / m_nDet; - if (m_nDet % 2 == 0) // Adjust for Even number of detectors - m_detInc = m_detLen / (m_nDet - 1); // center detector = (nDet/2)-1 - - m_initPos.angle = 0.0; + double dDetectorArrayEndOffset = 0; + if (isEven (m_nDet)) { // Adjust for Even number of detectors + m_detInc = m_detLen / (m_nDet - 1); // center detector = (nDet/2) + dDetectorArrayEndOffset = m_detInc; + m_detLen += dDetectorArrayEndOffset; + } + + m_dFanBeamAngle = dAngle * 2; m_initPos.xs1 = m_dXCenter; m_initPos.ys1 = m_dYCenter + m_dFocalLength; m_initPos.xs2 = m_dXCenter; m_initPos.ys2 = m_dYCenter + m_dFocalLength; m_initPos.xd1 = m_dXCenter - dHalfDetLen; - m_initPos.yd1 = m_dYCenter - m_dFocalLength; - m_initPos.xd2 = m_dXCenter + dHalfDetLen; - m_initPos.yd2 = m_dYCenter - m_dFocalLength; - m_initPos.angle = 0.0; + m_initPos.yd1 = m_dYCenter - m_dCenterDetectorLength; + m_initPos.xd2 = m_dXCenter + dHalfDetLen + dDetectorArrayEndOffset; + m_initPos.yd2 = m_dYCenter - m_dCenterDetectorLength; + m_initPos.angle = m_iOffsetView * m_rotInc; } else if (m_idGeometry == GEOMETRY_EQUIANGULAR) { -#if 0 - double dAngle = atan ((m_dFieldOfView / 2) / m_dFocalLength); -#else - double dHalfSquare = m_dFieldOfView / SQRT2 / 2; - double dFocalPastPhm = m_dFocalLength - dHalfSquare; - if (dFocalPastPhm <= 0.) { + if (m_dScanDiameter / 2 > m_dFocalLength) { m_fail = true; - m_failMessage = "Focal Point inside of phantom"; + m_failMessage = "Invalid geometry: Focal length must be larger than scan length"; return; } - double dAngle = atan ( dHalfSquare / dFocalPastPhm ); -#endif + const double dAngle = asin ((m_dScanDiameter / 2) / m_dFocalLength); + m_detLen = 2 * dAngle; + m_detStart = -dAngle; m_detInc = m_detLen / m_nDet; - if (m_nDet % 2 == 0) // Adjust for Even number of detectors - m_detInc = m_detLen / (m_nDet - 1); // center detector = (nDet/2)-1 - m_dAngularDetIncrement = m_detInc * 2; // Angular Position 2x gamma angle - m_dAngularDetLen = m_detLen * 2; + double dDetectorArrayEndOffset = 0; + if (isEven (m_nDet)) { // Adjust for Even number of detectors + m_detInc = m_detLen / (m_nDet - 1); // center detector = (nDet/2) + dDetectorArrayEndOffset = m_detInc; + } + // adjust for center-detector length + double dA1 = acos ((m_dScanDiameter / 2) / m_dCenterDetectorLength); + double dAngularScale = 2 * (HALFPI + dAngle - dA1) / m_detLen; + + m_dAngularDetLen = dAngularScale * (m_detLen + dDetectorArrayEndOffset); + m_dAngularDetIncrement = dAngularScale * m_detInc; m_initPos.dAngularDet = -m_dAngularDetLen / 2; - m_initPos.angle = 0; + m_dFanBeamAngle = dAngle * 2; + m_initPos.angle = m_iOffsetView * m_rotInc; m_initPos.xs1 = m_dXCenter; m_initPos.ys1 = m_dYCenter + m_dFocalLength;; m_initPos.xs2 = m_dXCenter; m_initPos.ys2 = m_dYCenter + m_dFocalLength; + m_detLen += dDetectorArrayEndOffset; } // Calculate incrementatal rotation matrix @@ -256,14 +282,25 @@ Scanner::convertGeometryNameToID (const char* const geomName) void Scanner::collectProjections (Projections& proj, const Phantom& phm, const int trace, SGP* pSGP) { - collectProjections (proj, phm, 0, proj.nView(), true, trace, pSGP); + collectProjections (proj, phm, m_startView, proj.nView(), m_iOffsetView, true, trace, pSGP); } void -Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iStartView, const int iNumViews, bool bStoreAtViewPosition, const int trace, SGP* pSGP) +Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iStartView, + const int iNumViews, const int iOffsetView, bool bStoreAtViewPosition, + const int trace, SGP* pSGP) +{ + int iStorageOffset = (bStoreAtViewPosition ? iStartView : 0); + collectProjections (proj, phm, iStartView, iNumViews, iOffsetView, iStorageOffset, trace, pSGP); +} + +void +Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iStartView, + const int iNumViews, const int iOffsetView, int iStorageOffset, + const int trace, SGP* pSGP) { m_trace = trace; - double start_angle = iStartView * proj.rotInc(); + double start_angle = (iStartView + iOffsetView) * proj.rotInc(); // Calculate initial rotation matrix GRFMTX_2D rotmtx_initial, temp; @@ -293,43 +330,48 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS int iView; double viewAngle; for (iView = 0, viewAngle = start_angle; iView < iNumViews; iView++, viewAngle += proj.rotInc()) { - int iStoragePosition = iView; - if (bStoreAtViewPosition) - iStoragePosition += iStartView; - + int iStoragePosition = iView + iStorageOffset; + DetectorArray& detArray = proj.getDetectorArray( iStoragePosition ); #ifdef HAVE_SGP if (pSGP && m_trace >= Trace::TRACE_PHANTOM) { m_pSGP = pSGP; - double dWindowSize = dmax (m_detLen, m_dFocalLength * 2) * SQRT2; + double dWindowSize = dmax (m_detLen, m_dSourceDetectorLength) * 2; double dHalfWindowSize = dWindowSize / 2; m_dXMinWin = m_dXCenter - dHalfWindowSize; m_dXMaxWin = m_dXCenter + dHalfWindowSize; m_dYMinWin = m_dYCenter - dHalfWindowSize; m_dYMaxWin = m_dYCenter + dHalfWindowSize; - double dHalfPhmLen = m_phmLen / 2; m_pSGP->setWindow (m_dXMinWin, m_dYMinWin, m_dXMaxWin, m_dYMaxWin); m_pSGP->setRasterOp (RO_COPY); + m_pSGP->setColor (C_RED); m_pSGP->moveAbs (0., 0.); - m_pSGP->drawRect (m_dXCenter - dHalfPhmLen, m_dYCenter - dHalfPhmLen, m_dXCenter + dHalfPhmLen, m_dYCenter + dHalfPhmLen); + m_pSGP->drawCircle (m_dViewDiameter / 2); + m_pSGP->moveAbs (0., 0.); + m_pSGP->setColor (C_GREEN); m_pSGP->drawCircle (m_dFocalLength); m_pSGP->setColor (C_BLUE); +#if MSVC + m_pSGP->setTextPointSize (9); +#else + m_pSGP->setTextPointSize (14); +#endif phm.draw (*m_pSGP); m_dTextHeight = m_pSGP->getCharHeight (); traceShowParam ("Phantom:", "%s", PROJECTION_TRACE_ROW_PHANT_ID, C_BLACK, phm.name().c_str()); traceShowParam ("Geometry:", "%s", PROJECTION_TRACE_ROW_GEOMETRY, C_BLUE, convertGeometryIDToName(m_idGeometry)); traceShowParam ("Focal Length Ratio:", "%.2f", PROJECTION_TRACE_ROW_FOCAL_LENGTH, C_BLUE, m_dFocalLengthRatio); - traceShowParam ("Field Of View Ratio:", "%.2f", PROJECTION_TRACE_ROW_FIELD_OF_VIEW, C_BLUE, m_dFieldOfViewRatio); + // traceShowParam ("Field Of View Ratio:", "%.2f", PROJECTION_TRACE_ROW_FIELD_OF_VIEW, C_BLUE, m_dFieldOfViewRatio); traceShowParam ("Num Detectors:", "%d", PROJECTION_TRACE_ROW_NDET, C_BLUE, proj.nDet()); traceShowParam ("Num Views:", "%d", PROJECTION_TRACE_ROW_NVIEW, C_BLUE, proj.nView()); traceShowParam ("Samples / Ray:", "%d", PROJECTION_TRACE_ROW_SAMPLES, C_BLUE, m_nSample); - m_pSGP->setMarker (SGP::MARK_BDIAMOND, C_LTGREEN); + m_pSGP->setMarker (SGP::MARKER_BDIAMOND); } #endif @@ -355,7 +397,7 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS m_pSGP->lineAbs (xs2, ys2); m_pSGP->setPenWidth (2); m_pSGP->moveAbs (0., 0.); - m_pSGP->drawArc (m_dFocalLength, viewAngle + 3 * HALFPI - (m_dAngularDetLen/2), viewAngle + 3 * HALFPI + (m_dAngularDetLen/2)); + m_pSGP->drawArc (m_dCenterDetectorLength, viewAngle + 3 * HALFPI - (m_dAngularDetLen/2), viewAngle + 3 * HALFPI + (m_dAngularDetLen/2)); } m_pSGP->setPenWidth (1); } @@ -438,10 +480,8 @@ Scanner::projectSingleView (const Phantom& phm, DetectorArray& detArray, const d if (phm.getComposition() == P_UNIT_PULSE) { // put unit pulse in center of view for (int d = 0; d < detArray.nDet(); d++) - if (detArray.nDet() / 2 == d && (d % 2) == 1) - detval[d] = 1; - else detval[d] = 0; + detval[ detArray.nDet() / 2 ] = 1; } else { for (int d = 0; d < detArray.nDet(); d++) { double xs = xs_maj; @@ -456,8 +496,8 @@ Scanner::projectSingleView (const Phantom& phm, DetectorArray& detArray, const d double sum = 0.0; for (unsigned int i = 0; i < m_nSample; i++) { if (m_idGeometry == GEOMETRY_EQUIANGULAR) { - xd = m_dFocalLength * cos (dAngle); - yd = m_dFocalLength * sin (dAngle); + xd = m_dCenterDetectorLength * cos (dAngle); + yd = m_dCenterDetectorLength * sin (dAngle); } #ifdef HAVE_SGP @@ -535,13 +575,24 @@ Scanner::traceShowParamRasterOp (int iRasterOp, const char *szLabel, const char #ifdef HAVE_SGP if (m_pSGP) { m_pSGP->setRasterOp (iRasterOp); - double dYPos = m_dYMaxWin - (row * m_dTextHeight); - m_pSGP->moveAbs (m_dXMinWin, dYPos); m_pSGP->setTextColor (color, -1); - m_pSGP->drawText (szLabel); double dValueOffset = (m_dXMaxWin - m_dXMinWin) / 4; - m_pSGP->moveAbs (m_dXMinWin + dValueOffset, dYPos); - m_pSGP->drawText (szValue); + if (row < 4) { + double dYPos = m_dYMaxWin - (row * m_dTextHeight); + double dXPos = m_dXMinWin; + m_pSGP->moveAbs (dXPos, dYPos); + m_pSGP->drawText (szLabel); + m_pSGP->moveAbs (dXPos + dValueOffset, dYPos); + m_pSGP->drawText (szValue); + } else { + row -= 4; + double dYPos = m_dYMaxWin - (row * m_dTextHeight); + double dXPos = m_dXMinWin + (m_dXMaxWin - m_dXMinWin) * 0.5; + m_pSGP->moveAbs (dXPos, dYPos); + m_pSGP->drawText (szLabel); + m_pSGP->moveAbs (dXPos + dValueOffset, dYPos); + m_pSGP->drawText (szValue); + } } else #endif {