X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fscanner.cpp;h=48f23943a739ec654d361fba8fdcc35f1fad92d5;hp=c45d4fd3234c1aa71561728d3bfb952b0840eeba;hb=9703ace3c451ce079967284bf191783736dbc77f;hpb=6bfb747f8163381d94a02c03a0351e9ca6815d22 diff --git a/libctsim/scanner.cpp b/libctsim/scanner.cpp index c45d4fd..48f2394 100644 --- a/libctsim/scanner.cpp +++ b/libctsim/scanner.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: scanner.cpp,v 1.12 2000/08/27 20:32:55 kevin Exp $ +** $Id: scanner.cpp,v 1.18 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 @@ -30,21 +30,21 @@ const int Scanner::GEOMETRY_INVALID = -1; const int Scanner::GEOMETRY_PARALLEL = 0; -const int Scanner::GEOMETRY_EQUILINEAR = 1; -const int Scanner::GEOMETRY_EQUIANGULAR = 2; +const int Scanner::GEOMETRY_EQUIANGULAR = 1; +const int Scanner::GEOMETRY_EQUILINEAR = 2; const char* Scanner::s_aszGeometryName[] = { {"parallel"}, - {"equilinear"}, {"equiangular"}, + {"equilinear"}, }; const char* Scanner::s_aszGeometryTitle[] = { {"Parallel"}, - {"Equilinear"}, {"Equiangular"}, + {"Equilinear"}, }; const int Scanner::s_iGeometryCount = sizeof(s_aszGeometryName) / sizeof(const char*); @@ -95,14 +95,13 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet, return; } - if (nView < 1) - nView = 1; + if (nView < 1 || nDet < 1) { + m_fail = true; + m_failMessage = "nView & nDet must be greater than 0"; + return; + } if (nSample < 1) m_nSample = 1; - if (nDet < 1) - nDet = 1; - // if ((nDet % 2) == 0) - // ++nDet; // ensure odd number of detectors m_nDet = nDet; m_nView = nView; @@ -114,23 +113,28 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet, 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_detInc = m_detLen / m_nDet; - m_rotLen = rot_anglen; - m_rotInc = m_rotLen / m_nView; + if (m_nDet % 2 == 0) // Adjust for Even number of detectors + m_detInc = m_detLen / (m_nDet - 1); // center detector = (nDet/2)-1 double dHalfDetLen = m_detLen / 2; - m_initPos.xs1 = m_dXCenter - m_dFocalLength; - m_initPos.ys1 = m_dYCenter - dHalfDetLen; - m_initPos.xs2 = m_dXCenter - m_dFocalLength; - m_initPos.ys2 = m_dYCenter + dHalfDetLen; - m_initPos.xd1 = m_dXCenter + m_dFocalLength; - m_initPos.yd1 = m_dYCenter - dHalfDetLen; - m_initPos.xd2 = m_dXCenter + m_dFocalLength; - m_initPos.yd2 = m_dYCenter + dHalfDetLen; + m_initPos.xs1 = m_dXCenter - dHalfDetLen; + m_initPos.ys1 = m_dYCenter + m_dFocalLength; + m_initPos.xs2 = m_dXCenter + dHalfDetLen; + 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; } 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.) { @@ -139,23 +143,28 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet, return; } double dAngle = atan( dHalfSquare / dFocalPastPhm ); +#endif double dHalfDetLen = 2 * m_dFocalLength * tan (dAngle); m_detLen = dHalfDetLen * 2; m_detInc = m_detLen / m_nDet; - m_rotLen = rot_anglen; - m_rotInc = m_rotLen / m_nView; + 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.xs1 = m_dXCenter - m_dFocalLength; - m_initPos.ys1 = m_dYCenter; - m_initPos.xs2 = m_dXCenter - m_dFocalLength; - m_initPos.ys2 = m_dYCenter; - m_initPos.xd1 = m_dXCenter + m_dFocalLength; - m_initPos.yd1 = m_dYCenter - dHalfDetLen; - m_initPos.xd2 = m_dXCenter + m_dFocalLength; - m_initPos.yd2 = m_dYCenter + dHalfDetLen; + m_initPos.angle = 0.0; + 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; } 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.) { @@ -163,17 +172,21 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet, m_failMessage = "Focal Point inside of phantom"; return; } - double dAngle = 2 * atan( dHalfSquare / dFocalPastPhm ); + double dAngle = atan ( dHalfSquare / dFocalPastPhm ); +#endif m_detLen = 2 * dAngle; m_detInc = m_detLen / m_nDet; - m_rotLen = rot_anglen; - m_rotInc = m_rotLen / m_nView; - - m_initPos.xs1 = m_dXCenter - m_dFocalLength; - m_initPos.ys1 = m_dYCenter; - m_initPos.xs2 = m_dXCenter - m_dFocalLength; - m_initPos.ys2 = m_dYCenter; - m_initPos.angle = -dAngle; + 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; + m_initPos.dAngularDet = -m_dAngularDetLen / 2; + + m_initPos.angle = 0; + m_initPos.xs1 = m_dXCenter; + m_initPos.ys1 = m_dYCenter + m_dFocalLength;; + m_initPos.xs2 = m_dXCenter; + m_initPos.ys2 = m_dYCenter + m_dFocalLength; } // Calculate incrementatal rotation matrix @@ -260,10 +273,8 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS xlat_mtx2 (temp, m_dXCenter, m_dYCenter); mult_mtx2 (rotmtx_initial, temp, rotmtx_initial); - double xd1=0, yd1=0, xd2=0, yd2=0, dDetAngle=0; - if (m_idGeometry == GEOMETRY_EQUIANGULAR) - dDetAngle = m_initPos.angle + start_angle; - else { + double xd1=0, yd1=0, xd2=0, yd2=0; + if (m_idGeometry != GEOMETRY_EQUIANGULAR) { xd1 = m_initPos.xd1; yd1 = m_initPos.yd1; xd2 = m_initPos.xd2; @@ -303,17 +314,15 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS m_pSGP->eraseWindow (); m_pSGP->setWindow (m_dXMinWin, m_dYMinWin, m_dXMaxWin, m_dYMaxWin); m_pSGP->setRasterOp (RO_COPY); - m_pSGP->setColor (C_BLUE); + 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->moveAbs (0., 0.); m_pSGP->drawCircle (m_dFocalLength); - m_pSGP->setColor (C_RED); + m_pSGP->setColor (C_BLUE); phm.draw (*m_pSGP); m_dTextHeight = m_pSGP->getCharHeight (); - traceShowParam ("Projection Collector", "%s", PROJECTION_TRACE_ROW_TITLE, C_BLACK, " "); - traceShowParam ("________________", "%s", PROJECTION_TRACE_ROW_TITLE2, C_LTGRAY, " "); 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); @@ -348,15 +357,17 @@ 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, start_angle + m_initPos.angle, start_angle - m_initPos.angle); + m_pSGP->drawArc (m_dFocalLength, viewAngle + 3 * HALFPI - (m_dAngularDetLen/2), viewAngle + 3 * HALFPI + (m_dAngularDetLen/2)); } m_pSGP->setPenWidth (1); } - if (m_trace >= Trace::TRACE_CONSOLE) - traceShowParam ("Current View:", "%d (%.0f%%)", PROJECTION_TRACE_ROW_CURR_VIEW, C_RED, iView + iStartView, (iView + iStartView) / m_nView * 100.); + if (m_trace > Trace::TRACE_CONSOLE) + traceShowParam ("Current View:", "%d (%.0f%%)", PROJECTION_TRACE_ROW_CURR_VIEW, C_RED, iView + iStartView, (iView + iStartView) / static_cast(m_nView) * 100.); #endif - - projectSingleView (phm, detArray, xd1, yd1, xd2, yd2, xs1, ys1, xs2, ys2, dDetAngle); + if (m_trace == Trace::TRACE_CONSOLE) + cout << "Current View: " << iView+iStartView << endl; + + projectSingleView (phm, detArray, xd1, yd1, xd2, yd2, xs1, ys1, xs2, ys2, viewAngle + 3 * HALFPI); detArray.setViewAngle (viewAngle); #ifdef HAVE_SGP @@ -366,9 +377,7 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS #endif xform_mtx2 (m_rotmtxIncrement, xs1, ys1); xform_mtx2 (m_rotmtxIncrement, xs2, ys2); - if (m_idGeometry == GEOMETRY_EQUIANGULAR) - dDetAngle += m_detInc; - else { + if (m_idGeometry != GEOMETRY_EQUIANGULAR) { xform_mtx2 (m_rotmtxIncrement, xd1, yd1); // rotate detector endpoints xform_mtx2 (m_rotmtxIncrement, xd2, yd2); } @@ -411,10 +420,10 @@ Scanner::projectSingleView (const Phantom& phm, DetectorArray& detArray, const d double ddx=0, ddy=0, ddx2=0, ddy2=0, ddx2_ofs=0, ddy2_ofs=0, xd_maj=0, yd_maj=0; double dAngleInc=0, dAngleSampleInc=0, dAngleSampleOffset=0, dAngleMajor=0; if (m_idGeometry == GEOMETRY_EQUIANGULAR) { - dAngleInc = m_detInc; + dAngleInc = m_dAngularDetIncrement; dAngleSampleInc = dAngleInc / m_nSample; dAngleSampleOffset = dAngleSampleInc / 2; - dAngleMajor = dDetAngle + dAngleSampleOffset; + dAngleMajor = dDetAngle - (m_dAngularDetLen/2) + dAngleSampleOffset; } else { ddx = (xd2 - xd1) / detArray.nDet(); // change in coords ddy = (yd2 - yd1) / detArray.nDet(); // between detectors @@ -456,27 +465,18 @@ Scanner::projectSingleView (const Phantom& phm, DetectorArray& detArray, const d #ifdef HAVE_SGP if (m_pSGP && m_trace >= Trace::TRACE_PROJECTIONS) { m_pSGP->setColor (C_YELLOW); - m_pSGP->setRasterOp (RO_OR_REVERSE); + m_pSGP->setRasterOp (RO_AND); m_pSGP->moveAbs (xs, ys); m_pSGP->lineAbs (xd, yd); - m_pSGP->setRasterOp (RO_SET); } #endif sum += projectSingleLine (phm, xd, yd, xs, ys); #ifdef HAVE_SGP - if (m_trace >= Trace::TRACE_CLIPPING) { - traceShowParam ("Attenuation:", "%s", PROJECTION_TRACE_ROW_ATTEN, C_LTMAGENTA, " "); - traceShowParam ("Attenuation:", "%.3f", PROJECTION_TRACE_ROW_ATTEN, C_LTMAGENTA, sum); - } - - // if (m_pSGP && m_trace >= Trace::TRACE_PROJECTIONS) { - // m_pSGP->setColor (C_YELLOW); - // m_pSGP->setRasterOp (RO_XOR); - // m_pSGP->moveAbs (xs, ys); - // m_pSGP->lineAbs (xd, yd); - // m_pSGP->setRasterOp (RO_SET); + // if (m_trace >= Trace::TRACE_CLIPPING) { + // traceShowParam ("Attenuation:", "%s", PROJECTION_TRACE_ROW_ATTEN, C_LTMAGENTA, " "); + // traceShowParam ("Attenuation:", "%.3f", PROJECTION_TRACE_ROW_ATTEN, C_LTMAGENTA, sum); // } #endif if (m_idGeometry == GEOMETRY_EQUIANGULAR) @@ -506,7 +506,11 @@ Scanner::traceShowParam (const char *szLabel, const char *fmt, int row, int colo { va_list arg; va_start(arg, color); +#ifdef HAVE_SGP traceShowParamRasterOp (RO_COPY, szLabel, fmt, row, color, arg); +#else + traceShowParamRasterOp (0, szLabel, fmt, row, color, arg); +#endif va_end(arg); } @@ -515,7 +519,11 @@ Scanner::traceShowParamXOR (const char *szLabel, const char *fmt, int row, int c { va_list arg; va_start(arg, color); +#ifdef HAVE_SGP traceShowParamRasterOp (RO_XOR, szLabel, fmt, row, color, arg); +#else + traceShowParamRasterOp (0, szLabel, fmt, row, color, arg); +#endif va_end(arg); } @@ -524,26 +532,33 @@ Scanner::traceShowParamRasterOp (int iRasterOp, const char *szLabel, const char { char szValue[256]; +#ifdef MSVC + vsprintf (szValue, fmt, args); +#else vsnprintf (szValue, sizeof(szValue), fmt, args); +#endif // cio_set_cpos (raysum_trace_menu_column, row); // cio_set_text_clr (color - 8, 0); // cio_set_text_clr (color, 0); +#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) / 5; + double dValueOffset = (m_dXMaxWin - m_dXMinWin) / 4; m_pSGP->moveAbs (m_dXMinWin + dValueOffset, dYPos); m_pSGP->drawText (szValue); - } else { - cio_put_str (szLabel); - cio_put_str (szValue); - cio_put_str ("\n"); - } + } else +#endif + { + cio_put_str (szLabel); + cio_put_str (szValue); + cio_put_str ("\n"); + } }