X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fscanner.cpp;h=c4fc4d06a2b8041c246ba1a47aff538412d54c12;hp=709dc0f5698016f515758a031a83c75f21e2c1c5;hb=1e88cf0f7fa4f690ea9f110e8ed3f2b5338d0a10;hpb=74a34e63a9a356e1467acdba65497ab15190dde0 diff --git a/libctsim/scanner.cpp b/libctsim/scanner.cpp index 709dc0f..c4fc4d0 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.10 2000/08/03 09:21:12 kevin Exp $ +** $Id: scanner.cpp,v 1.11 2000/08/25 15:59:13 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 @@ -82,7 +82,7 @@ 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) +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) { m_phmLen = phm.maxAxisLength(); // maximal length along an axis @@ -101,29 +101,67 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet, m_nSample = 1; if (nDet < 1) nDet = 1; - if ((nDet % 2) == 0) - ++nDet; // ensure odd number of detectors + // if ((nDet % 2) == 0) + // ++nDet; // ensure odd number of detectors m_nDet = nDet; m_nView = nView; m_nSample = nSample; - m_detLen = SQRT2 * m_phmLen * ((m_nDet + N_EXTRA_DETECTORS) / static_cast(m_nDet)); - - m_rotLen = rot_anglen; - - m_radius = m_detLen / 2; - m_detInc = m_detLen / m_nDet; - m_rotInc = m_rotLen / m_nView; - - m_initPos.xd1 = m_detLen / 2; - m_initPos.yd1 = -m_detLen / 2; - m_initPos.xd2 = m_detLen / 2; - m_initPos.yd2 = m_detLen / 2; - m_initPos.xs1 = -m_detLen / 2; - m_initPos.ys1 = -m_detLen / 2; - m_initPos.xs2 = -m_detLen / 2; - m_initPos.ys2 = m_detLen / 2; - m_initPos.angle = 0.0; + m_dFocalLength = (m_phmLen * SQRT2 / 2) * dFocalLengthRatio; + m_dFieldOfView = m_phmLen * SQRT2 * dFieldOfViewRatio; + + 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; + + double dHalfDetLen = m_detLen / 2; + m_initPos.xs1 = -m_dFocalLength; + m_initPos.ys1 = -dHalfDetLen; + m_initPos.xs2 = -m_dFocalLength; + m_initPos.ys2 = dHalfDetLen; + m_initPos.xd1 = m_dFocalLength; + m_initPos.yd1 = -dHalfDetLen; + m_initPos.xd2 = m_dFocalLength; + m_initPos.yd2 = dHalfDetLen; + m_initPos.angle = 0.0; + } else if (m_idGeometry == GEOMETRY_EQUILINEAR) { + double dHalfSquare = m_phmLen / 2; + double dFocalPastPhm = m_dFocalLength - dHalfSquare; + if (dFocalPastPhm <= 0.) { + m_fail = true; + m_failMessage = "Focal Point inside of phantom"; + return; + } + double dAngle = atan( dHalfSquare / dFocalPastPhm ); + 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; + + m_initPos.xs1 = -m_dFocalLength; + m_initPos.ys1 = 0; + m_initPos.xs2 = -m_dFocalLength; + m_initPos.ys2 = 0; + m_initPos.xd1 = m_dFocalLength; + m_initPos.yd1 = -dHalfDetLen; + m_initPos.xd2 = m_dFocalLength; + m_initPos.yd2 = dHalfDetLen; + m_initPos.angle = 0.0; + } else if (m_idGeometry == GEOMETRY_EQUIANGULAR) { + double dHalfSquare = m_phmLen / 2; + double dFocalPastPhm = m_dFocalLength - dHalfSquare; + if (dFocalPastPhm <= 0.) { + m_fail = true; + m_failMessage = "Focal Point inside of phantom"; + return; + } + double dAngle = atan( dHalfSquare / dFocalPastPhm ); + m_detLen = dAngle * 2; + } } Scanner::~Scanner (void) @@ -169,23 +207,30 @@ Scanner::convertGeometryNameToID (const char* const geomName) /* NAME - * raysum_collect Calculate ray sums for a Phantom + * collectProjections Calculate projections for a Phantom * * SYNOPSIS - * rs = raysum_collect (det, phm, start_view, trace, unit_pulse) - * Scanner& det Scanner specifications** - * RAYSUM *rs Calculated ray sum matrix - * Phantom& phm Phantom we are taking ray sums of - * int trace Boolean flag to signal ray sum tracing + * collectProjections (proj, phm, start_view, nView, bStoreViewPos, trace) + * Projectrions& proj Projection storage + * Phantom& phm Phantom for which we collect projections + * bool bStoreViewPos TRUE then storage proj at normal view position + * int trace Trace level */ + void -Scanner::collectProjections (Projections& proj, const Phantom& phm, const int start_view = 0, const int trace = TRACE_NONE, SGP* pSGP = NULL) +Scanner::collectProjections (Projections& proj, const Phantom& phm, const int trace, SGP* pSGP) +{ + collectProjections (proj, phm, 0, proj.nView(), true, trace, pSGP); +} + +void +Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iStartView, const int iNumViews, bool bStoreAtViewPosition, const int trace, SGP* pSGP) { GRFMTX_2D rotmtx_initial, temp; GRFMTX_2D rotmtx_incr; - double start_angle = start_view * proj.rotInc(); + double start_angle = iStartView * proj.rotInc(); double xcent = phm.xmin() + (phm.xmax() - phm.xmin()) / 2; double ycent = phm.ymin() + (phm.ymax() - phm.ymin()) / 2; @@ -200,38 +245,6 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int st m_trace = trace; -#ifdef HAVE_SGP - if (pSGP && m_trace >= TRACE_PHM) { - double halfPhmLen = m_phmLen / 2; - double wsize = SQRT2 * halfPhmLen; - - pSGP->setRasterOp (RO_SET); - pSGP->eraseWindow (); - pSGP->setColor (C_LTBLUE); - pSGP->setWindow (xcent - wsize, ycent - wsize, xcent + wsize, ycent + wsize); - pSGP->setColor (C_BROWN); - pSGP->drawRect (xcent - halfPhmLen, ycent - halfPhmLen, xcent + halfPhmLen, ycent + halfPhmLen); - pSGP->setColor (C_BROWN); - pSGP->moveAbs (0., 0.); - pSGP->drawCircle (wsize); - - traceShowParam ("X-Ray Simulator", "%s", RAYSUM_TRACE_ROW_TITLE, C_BLACK, " "); - traceShowParam ("---------------", "%s", RAYSUM_TRACE_ROW_TITLE2, C_BLACK, " "); - traceShowParam ("Phantom:", "%s", RAYSUM_TRACE_ROW_PHANT_ID, C_YELLOW, " Herman"); - traceShowParam ("Chomaticity :", "%s", RAYSUM_TRACE_ROW_CHROMATIC, C_LTGREEN, "Mono"); - traceShowParam ("Scatter :", "%5.1f", RAYSUM_TRACE_ROW_SCATTER, C_LTGREEN, 0.); - traceShowParam ("Photon Uncert:", "%5.1f", RAYSUM_TRACE_ROW_PHOT_STAT, C_LTGREEN, 0.); - traceShowParam ("Num Detectors:", "%5d", RAYSUM_TRACE_ROW_NDET, C_LTRED, proj.nDet()); - traceShowParam ("Num Views :", "%5d", RAYSUM_TRACE_ROW_NVIEW, C_LTRED, proj.nView()); - traceShowParam ("Samples / Ray:", "%5d", RAYSUM_TRACE_ROW_SAMPLES, C_LTRED, m_nSample); - - pSGP->setColor (C_LTGREEN); - phm.draw (*pSGP); - - pSGP->setMarker (SGP::MARK_BDIAMOND, C_LTGREEN); - } -#endif - /* Calculate initial rotation matrix */ xlat_mtx2 (rotmtx_initial, -xcent, -ycent); rot_mtx2 (temp, start_angle); @@ -253,8 +266,45 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int st int iview; double viewAngle; - for (iview = 0, viewAngle = start_angle; iview < proj.nView(); iview++, viewAngle += proj.rotInc()) { - DetectorArray& detArray = proj.getDetectorArray( iview ); + for (iview = 0, viewAngle = start_angle; iview < iNumViews; iview++, viewAngle += proj.rotInc()) { + int iStoragePosition = iview; + if (bStoreAtViewPosition) + iStoragePosition += iStartView; + + DetectorArray& detArray = proj.getDetectorArray( iStoragePosition ); + +#ifdef HAVE_SGP + if (pSGP && m_trace >= TRACE_PHM) { + pSGP->eraseWindow(); + double dWindowSize = max(m_detLen, m_dFocalLength * 2) * SQRT2; + double dHalfWindowSize = dWindowSize / 2; + double dHalfPhmLen = m_phmLen / 2; + + pSGP->setRasterOp (RO_SET); + pSGP->eraseWindow (); + pSGP->setWindow (xcent - dHalfWindowSize, ycent - dHalfWindowSize, xcent + dHalfWindowSize, ycent + dHalfWindowSize); + pSGP->setColor (C_BROWN); + pSGP->moveAbs (0., 0.); + pSGP->drawRect (xcent - dHalfPhmLen, ycent - dHalfPhmLen, xcent + dHalfPhmLen, ycent + dHalfPhmLen); + +#if 0 + traceShowParam (pSGP, "X-Ray Simulator", "%s", RAYSUM_TRACE_ROW_TITLE, C_BLACK, " "); + traceShowParam (pSGP, "---------------", "%s", RAYSUM_TRACE_ROW_TITLE2, C_BLACK, " "); + traceShowParam (pSGP, "Phantom:", "%s", RAYSUM_TRACE_ROW_PHANT_ID, C_YELLOW, " Herman"); + traceShowParam (pSGP, "Chomaticity :", "%s", RAYSUM_TRACE_ROW_CHROMATIC, C_LTGREEN, "Mono"); + traceShowParam (pSGP, "Scatter :", "%5.1f", RAYSUM_TRACE_ROW_SCATTER, C_LTGREEN, 0.); + traceShowParam (pSGP, "Photon Uncert:", "%5.1f", RAYSUM_TRACE_ROW_PHOT_STAT, C_LTGREEN, 0.); + traceShowParam (pSGP, "Num Detectors:", "%5d", RAYSUM_TRACE_ROW_NDET, C_LTRED, proj.nDet()); + traceShowParam (pSGP, "Num Views :", "%5d", RAYSUM_TRACE_ROW_NVIEW, C_LTRED, proj.nViews()); + traceShowParam (pSGP, "Samples / Ray:", "%5d", RAYSUM_TRACE_ROW_SAMPLES, C_LTRED, m_nSample); +#endif + + pSGP->setColor (C_RED); + phm.draw (*pSGP); + + pSGP->setMarker (SGP::MARK_BDIAMOND, C_LTGREEN); + } +#endif #ifdef HAVE_SGP if (pSGP && m_trace >= TRACE_PHM) { @@ -264,9 +314,10 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int st pSGP->lineAbs (xd2, yd2); pSGP->moveAbs (xs1, ys1); pSGP->lineAbs (xs2, ys2); + pSGP->setRasterOp (RO_SET); } - if (m_trace) - traceShowParam ("Current View :", "%5d", RAYSUM_TRACE_ROW_CURR_VIEW, C_LTMAGENTA, iview); + if (m_trace >= TRACE_TEXT) + traceShowParam (pSGP, "Current View :", "%5d", RAYSUM_TRACE_ROW_CURR_VIEW, C_LTMAGENTA, iview); #endif projectSingleView (phm, detArray, xd1, yd1, xd2, yd2, xs1, ys1, xs2, ys2, pSGP); @@ -276,10 +327,12 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int st if (pSGP && m_trace >= TRACE_PHM) { // rs_plot (detArray, xd1, yd1, xcent, ycent, theta); pSGP->setColor (C_RED); + pSGP->setRasterOp (RO_XOR); pSGP->moveAbs (xd1, yd1); pSGP->lineAbs (xd2, yd2); pSGP->moveAbs (xs1, ys1); pSGP->lineAbs (xs2, ys2); + pSGP->setRasterOp (RO_SET); } #endif xform_mtx2 (rotmtx_incr, xd1, yd1); // rotate detector endpoints @@ -350,20 +403,25 @@ Scanner::projectSingleView (const Phantom& phm, DetectorArray& detArray, const d #ifdef HAVE_SGP if (pSGP && m_trace >= TRACE_RAYS) { pSGP->setColor (C_LTBLUE); + pSGP->setRasterOp (RO_XOR); pSGP->moveAbs (xs, ys); pSGP->lineAbs (xd, yd); + pSGP->setRasterOp (RO_SET); } #endif sum += projectSingleLine (phm, xd, yd, xs, ys, pSGP); #ifdef HAVE_SGP - if (m_trace >= TRACE_RAYS) - traceShowParam ("Attenuation :", "%5.2f", RAYSUM_TRACE_ROW_ATTEN, C_LTMAGENTA, "sum"); + if (m_trace >= TRACE_RAYS) { + traceShowParam (pSGP, "Attenuation :", "%5.2f", RAYSUM_TRACE_ROW_ATTEN, C_LTMAGENTA, sum); + } if (pSGP && m_trace >= TRACE_RAYS) { pSGP->setColor (C_LTBLUE); + pSGP->setRasterOp (RO_XOR); pSGP->moveAbs (xs, ys); pSGP->lineAbs (xd, yd); + pSGP->setRasterOp (RO_SET); } #endif xd += ddx2; @@ -381,20 +439,30 @@ Scanner::projectSingleView (const Phantom& phm, DetectorArray& detArray, const d void -Scanner::traceShowParam (const char *label, const char *fmt, int row, int color, ...) +Scanner::traceShowParam (SGP* pSGP, const char *label, const char *fmt, int row, int color, ...) { - char s[80]; + char s[256]; va_list arg; va_start(arg, color); - // cio_set_cpos (raysum_trace_menu_column, row); snprintf (s, sizeof(s), label, "%s"); - // cio_set_text_clr (color - 8, 0); - cio_put_str (s); + string strOut(s); vsnprintf (s, sizeof(s), fmt, arg); + strOut += s; + + // cio_set_cpos (raysum_trace_menu_column, row); + // cio_set_text_clr (color - 8, 0); // cio_set_text_clr (color, 0); - cio_put_str (s); - cout << "\n"; + + if (pSGP) { + pSGP->moveAbs (0., row * 0.04); + pSGP->setTextColor (color, -1); + pSGP->drawText (strOut); + } else { + cio_put_str (strOut.c_str()); + cio_put_str ("\n"); + } + va_end(arg); } @@ -442,11 +510,13 @@ Scanner::projectLineAgainstPElem (const PhantomElement& pelem, double x1, double #ifdef HAVE_SGP if (pSGP && m_trace == TRACE_CLIPPING) { + pSGP->setRasterOp (RO_XOR); pSGP->moveAbs (x1, y1); pSGP->lineAbs (x2, y2); cio_tone (8000., 0.05); pSGP->moveAbs (x1, y1); pSGP->lineAbs (x2, y2); + pSGP->setRasterOp (RO_SET); } #endif