X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fscanner.cpp;h=204c9cb13b86b30a3c8130ad0eea0d1a05af3c61;hp=0319912339cbee2c8755852cc0bcdbbca561aee5;hb=dca8d386ab6a4e33050b456cff4547d6adfdf249;hpb=ac39e654e1385006647c73b253d4d79ccf09a919 diff --git a/libctsim/scanner.cpp b/libctsim/scanner.cpp index 0319912..204c9cb 100644 --- a/libctsim/scanner.cpp +++ b/libctsim/scanner.cpp @@ -1,15 +1,13 @@ /***************************************************************************** ** FILE IDENTIFICATION ** -** Name: scanner.cpp +** Name: scanner.cpp ** Purpose: Classes for CT scanner ** Programmer: Kevin Rosenberg ** Date Started: 1984 ** ** This is part of the CTSim program -** Copyright (c) 1983-2001 Kevin Rosenberg -** -** $Id: scanner.cpp,v 1.38 2001/04/02 03:49:52 kevin Exp $ +** Copyright (c) 1983-2009 Kevin Rosenberg ** ** 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,26 +30,29 @@ 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*); // NAME -// DetectorArray Construct a DetectorArray +// DetectorArray Construct a DetectorArray DetectorArray::DetectorArray (const int nDet) { @@ -61,7 +62,7 @@ DetectorArray::DetectorArray (const int nDet) // NAME -// ~DetectorArray Free memory allocated to a detector array +// ~DetectorArray Free memory allocated to a detector array DetectorArray::~DetectorArray (void) { @@ -71,20 +72,22 @@ DetectorArray::~DetectorArray (void) /* NAME -* Scanner::Scanner Construct a user specified detector structure +* Scanner::Scanner Construct a user specified detector structure * * SYNOPSIS * Scanner (phm, nDet, nView, nSample) -* Phantom& phm PHANTOM that we are making detector for +* Phantom& phm PHANTOM that we are making detector for * int geomety Geometry of detector -* int nDet Number of detector along detector array -* int nView Number of rotated views -* int nSample Number of rays per detector +* int nDet Number of detector along detector array +* int nView Number of rotated views +* 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 dCenterDetectorRatio, +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_fail = false; @@ -95,7 +98,7 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, m_failMessage += geometryName; return; } - + if (nView < 1 || nDet < 1) { m_fail = true; m_failMessage = "nView & nDet must be greater than 0"; @@ -103,9 +106,10 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, } if (nSample < 1) m_nSample = 1; - + m_nDet = nDet; m_nView = nView; + m_iOffsetView = offsetView; m_nSample = nSample; m_dFocalLengthRatio = dFocalLengthRatio; m_dCenterDetectorRatio = dCenterDetectorRatio; @@ -117,7 +121,7 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, 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; @@ -144,7 +148,7 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, m_initPos.yd1 = m_dYCenter - m_dCenterDetectorLength; m_initPos.xd2 = m_dXCenter + dHalfDetLen + dDetectorArrayEndOffset; m_initPos.yd2 = m_dYCenter - m_dCenterDetectorLength; - m_initPos.angle = 0.0; + m_initPos.angle = m_iOffsetView * m_rotInc; m_detLen += dDetectorArrayEndOffset; } else if (m_idGeometry == GEOMETRY_EQUILINEAR) { if (m_dScanDiameter / 2 >= m_dFocalLength) { @@ -155,7 +159,7 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, 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; @@ -165,9 +169,8 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, dDetectorArrayEndOffset = m_detInc; m_detLen += dDetectorArrayEndOffset; } - + m_dFanBeamAngle = dAngle * 2; - m_initPos.angle = 0.0; m_initPos.xs1 = m_dXCenter; m_initPos.ys1 = m_dYCenter + m_dFocalLength; m_initPos.xs2 = m_dXCenter; @@ -176,7 +179,7 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, m_initPos.yd1 = m_dYCenter - m_dCenterDetectorLength; m_initPos.xd2 = m_dXCenter + dHalfDetLen + dDetectorArrayEndOffset; m_initPos.yd2 = m_dYCenter - m_dCenterDetectorLength; - m_initPos.angle = 0.0; + m_initPos.angle = m_iOffsetView * m_rotInc; } else if (m_idGeometry == GEOMETRY_EQUIANGULAR) { if (m_dScanDiameter / 2 > m_dFocalLength) { m_fail = true; @@ -200,24 +203,24 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, m_dAngularDetLen = dAngularScale * (m_detLen + dDetectorArrayEndOffset); m_dAngularDetIncrement = dAngularScale * m_detInc; m_initPos.dAngularDet = -m_dAngularDetLen / 2; - + m_dFanBeamAngle = dAngle * 2; - m_initPos.angle = 0; + 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 + + // Calculate incrementatal rotation matrix GRFMTX_2D temp; xlat_mtx2 (m_rotmtxIncrement, -m_dXCenter, -m_dYCenter); rot_mtx2 (temp, m_rotInc); mult_mtx2 (m_rotmtxIncrement, temp, m_rotmtxIncrement); xlat_mtx2 (temp, m_dXCenter, m_dYCenter); mult_mtx2 (m_rotmtxIncrement, temp, m_rotmtxIncrement); - + } Scanner::~Scanner (void) @@ -229,10 +232,10 @@ const char* Scanner::convertGeometryIDToName (const int geomID) { const char *name = ""; - + if (geomID >= 0 && geomID < s_iGeometryCount) return (s_aszGeometryName[geomID]); - + return (name); } @@ -240,94 +243,105 @@ const char* Scanner::convertGeometryIDToTitle (const int geomID) { const char *title = ""; - + if (geomID >= 0 && geomID < s_iGeometryCount) return (s_aszGeometryName[geomID]); - + return (title); } int -Scanner::convertGeometryNameToID (const char* const geomName) +Scanner::convertGeometryNameToID (const char* const geomName) { int id = GEOMETRY_INVALID; - - for (int i = 0; i < s_iGeometryCount; i++) + + for (int i = 0; i < s_iGeometryCount; i++) { if (strcasecmp (geomName, s_aszGeometryName[i]) == 0) { id = i; break; } - - return (id); + } + return (id); } /* NAME -* collectProjections Calculate projections for a Phantom +* collectProjections Calculate projections for a Phantom * * SYNOPSIS * collectProjections (proj, phm, start_view, nView, bStoreViewPos, trace) * Projectrions& proj Projection storage -* Phantom& phm Phantom for which we collect projections +* Phantom& phm Phantom for which we collect projections * bool bStoreViewPos TRUE then storage proj at normal view position -* int trace Trace level +* int trace Trace level */ 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, iStorageOffset, trace, pSGP); + collectProjections (proj, phm, iStartView, iNumViews, iOffsetView, iStorageOffset, trace, pSGP); +} + +static void mtx2_offset_rot (GRFMTX_2D m, double angle, double x, double y) { + GRFMTX_2D temp; + xlat_mtx2 (m, -x, -y); + rot_mtx2 (temp, angle); + mult_mtx2 (m, temp, m); + xlat_mtx2 (temp, x, y); + mult_mtx2 (m, temp, m); } void -Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iStartView, const int iNumViews, - int iStorageOffset, const int trace, SGP* pSGP) +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(); - - // Calculate initial rotation matrix - GRFMTX_2D rotmtx_initial, temp; - xlat_mtx2 (rotmtx_initial, -m_dXCenter, -m_dYCenter); - rot_mtx2 (temp, start_angle); - mult_mtx2 (rotmtx_initial, temp, rotmtx_initial); - xlat_mtx2 (temp, m_dXCenter, m_dYCenter); - mult_mtx2 (rotmtx_initial, temp, rotmtx_initial); - - 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; - yd2 = m_initPos.yd2; - xform_mtx2 (rotmtx_initial, xd1, yd1); // rotate detector endpoints - xform_mtx2 (rotmtx_initial, xd2, yd2); // to initial view_angle - } - - double xs1 = m_initPos.xs1; - double ys1 = m_initPos.ys1; - double xs2 = m_initPos.xs2; - double ys2 = m_initPos.ys2; - xform_mtx2 (rotmtx_initial, xs1, ys1); // rotate source endpoints to - xform_mtx2 (rotmtx_initial, xs2, ys2); // initial view angle - - int iView; - double viewAngle; - for (iView = 0, viewAngle = start_angle; iView < iNumViews; iView++, viewAngle += proj.rotInc()) { - int iStoragePosition = iView + iStorageOffset; + double start_angle = (iStartView + iOffsetView) * proj.rotInc(); + int parallel_enabled = 1; + UNUSED(parallel_enabled); + +#if HAVE_SGP + if (pSGP && (m_trace >= Trace::TRACE_PHANTOM)) + parallel_enabled = 0; +#endif + +#if HAVE_OPENMP + #pragma omp parallel for if (parallel_enabled) +#endif + for (int iView = 0; iView < iNumViews; iView++) { + double viewAngle = start_angle + (iView * proj.rotInc()); + + // With OpenMP, need to calculate source and detector positions at each view + GRFMTX_2D rotmtx; + mtx2_offset_rot (rotmtx, viewAngle, m_dXCenter, m_dYCenter); + 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; yd2 = m_initPos.yd2; + xform_mtx2 (rotmtx, xd1, yd1); // rotate detector endpoints + xform_mtx2 (rotmtx, xd2, yd2); // to initial view_angle + } + + double xs1 = m_initPos.xs1, ys1 = m_initPos.ys1; + double xs2 = m_initPos.xs2, ys2 = m_initPos.ys2; + xform_mtx2 (rotmtx, xs1, ys1); // rotate source endpoints to + xform_mtx2 (rotmtx, xs2, ys2); // initial view angle + int iStoragePosition = iView + iStorageOffset; DetectorArray& detArray = proj.getDetectorArray( iStoragePosition ); - -#ifdef HAVE_SGP + +#ifdef HAVE_SGP if (pSGP && m_trace >= Trace::TRACE_PHANTOM) { m_pSGP = pSGP; double dWindowSize = dmax (m_detLen, m_dSourceDetectorLength) * 2; @@ -336,7 +350,7 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS m_dXMaxWin = m_dXCenter + dHalfWindowSize; m_dYMinWin = m_dYCenter - dHalfWindowSize; m_dYMaxWin = m_dYCenter + dHalfWindowSize; - + m_pSGP->setWindow (m_dXMinWin, m_dYMinWin, m_dXMaxWin, m_dYMaxWin); m_pSGP->setRasterOp (RO_COPY); @@ -351,21 +365,17 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS m_pSGP->setTextPointSize (9); 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::MARKER_BDIAMOND); - } -#endif - -#ifdef HAVE_SGP - if (m_pSGP && m_trace >= Trace::TRACE_PHANTOM) { + m_pSGP->setColor (C_BLACK); m_pSGP->setPenWidth (2); if (m_idGeometry == GEOMETRY_PARALLEL) { @@ -373,7 +383,7 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS m_pSGP->lineAbs (xs2, ys2); m_pSGP->moveAbs (xd1, yd1); m_pSGP->lineAbs (xd2, yd2); - } else if (m_idGeometry == GEOMETRY_EQUILINEAR) { + } else if (m_idGeometry == GEOMETRY_EQUILINEAR) { m_pSGP->setPenWidth (4); m_pSGP->moveAbs (xs1, ys1); m_pSGP->lineAbs (xs2, ys2); @@ -393,37 +403,33 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS 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 + if (m_trace == Trace::TRACE_CONSOLE) std::cout << "Current View: " << iView+iStartView << std::endl; - + projectSingleView (phm, detArray, xd1, yd1, xd2, yd2, xs1, ys1, xs2, ys2, viewAngle + 3 * HALFPI); detArray.setViewAngle (viewAngle); - + #ifdef HAVE_SGP if (m_pSGP && m_trace >= Trace::TRACE_PHANTOM) { - // rs_plot (detArray, xd1, yd1, dXCenter, dYCenter, theta); + // rs_plot (detArray, xd1, yd1, dXCenter, dYCenter, theta); } #endif - xform_mtx2 (m_rotmtxIncrement, xs1, ys1); - xform_mtx2 (m_rotmtxIncrement, xs2, ys2); - if (m_idGeometry != GEOMETRY_EQUIANGULAR) { - xform_mtx2 (m_rotmtxIncrement, xd1, yd1); // rotate detector endpoints - xform_mtx2 (m_rotmtxIncrement, xd2, yd2); - } + } /* for each iView */ } /* NAME -* rayview Calculate raysums for a view at any angle +* rayview Calculate raysums for a view at any angle * * SYNOPSIS * rayview (phm, detArray, xd1, nSample, yd1, xd2, yd2, xs1, ys1, xs2, ys2) -* Phantom& phm Phantom to scan -* DETARRAY *detArray Storage of values for detector array -* Scanner& det Scanner parameters -* double xd1, yd1, xd2, yd2 Beginning & ending detector positions -* double xs1, ys1, xs2, ys2 Beginning & ending source positions +* Phantom& phm Phantom to scan +* DETARRAY *detArray Storage of values for detector array +* Scanner& det Scanner parameters +* double xd1, yd1, xd2, yd2 Beginning & ending detector positions +* double xs1, ys1, xs2, ys2 Beginning & ending source positions * * RAY POSITIONING * For each detector, have there are a variable number of rays traced. @@ -437,15 +443,15 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS * In that case, the detector position is the center of the detector cell. */ -void +void Scanner::projectSingleView (const Phantom& phm, DetectorArray& detArray, const double xd1, const double yd1, const double xd2, const double yd2, const double xs1, const double ys1, const double xs2, const double ys2, const double dDetAngle) { - - double sdx = (xs2 - xs1) / detArray.nDet(); // change in coords + + double sdx = (xs2 - xs1) / detArray.nDet(); // change in coords double sdy = (ys2 - ys1) / detArray.nDet(); // between source - double xs_maj = xs1 + (sdx / 2); // put ray source in center of cell + double xs_maj = xs1 + (sdx / 2); // put ray source in center of cell double ys_maj = ys1 + (sdy / 2); - + 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) { @@ -454,19 +460,19 @@ Scanner::projectSingleView (const Phantom& phm, DetectorArray& detArray, const d dAngleSampleOffset = dAngleSampleInc / 2; dAngleMajor = dDetAngle - (m_dAngularDetLen/2) + dAngleSampleOffset; } else { - ddx = (xd2 - xd1) / detArray.nDet(); // change in coords + ddx = (xd2 - xd1) / detArray.nDet(); // change in coords ddy = (yd2 - yd1) / detArray.nDet(); // between detectors - ddx2 = ddx / m_nSample; // Incr. between rays with detector cell - ddy2 = ddy / m_nSample; // Doesn't include detector endpoints + ddx2 = ddx / m_nSample; // Incr. between rays with detector cell + ddy2 = ddy / m_nSample; // Doesn't include detector endpoints ddx2_ofs = ddx2 / 2; // offset of 1st ray from start of detector cell ddy2_ofs = ddy2 / 2; - + xd_maj = xd1 + ddx2_ofs; // Incr. between detector cells yd_maj = yd1 + ddy2_ofs; } - + DetectorValue* detval = detArray.detValues(); - + if (phm.getComposition() == P_UNIT_PULSE) { // put unit pulse in center of view for (int d = 0; d < detArray.nDet(); d++) detval[d] = 0; @@ -488,7 +494,7 @@ Scanner::projectSingleView (const Phantom& phm, DetectorArray& detArray, const d xd = m_dCenterDetectorLength * cos (dAngle); yd = m_dCenterDetectorLength * sin (dAngle); } - + #ifdef HAVE_SGP if (m_pSGP && m_trace >= Trace::TRACE_PROJECTIONS) { m_pSGP->setColor (C_YELLOW); @@ -497,14 +503,14 @@ Scanner::projectSingleView (const Phantom& phm, DetectorArray& detArray, const d m_pSGP->lineAbs (xd, yd); } #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_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) dAngle += dAngleSampleInc; @@ -513,7 +519,7 @@ Scanner::projectSingleView (const Phantom& phm, DetectorArray& detArray, const d yd += ddy2; } } // for each sample in detector - + detval[d] = sum / m_nSample; xs_maj += sdx; ys_maj += sdy; @@ -528,22 +534,22 @@ Scanner::projectSingleView (const Phantom& phm, DetectorArray& detArray, const d } -void +void Scanner::traceShowParam (const char *szLabel, const char *fmt, int row, int color, ...) -{ +{ 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 +#endif va_end(arg); } -void +void Scanner::traceShowParamXOR (const char *szLabel, const char *fmt, int row, int color, ...) -{ +{ va_list arg; va_start(arg, color); #ifdef HAVE_SGP @@ -554,13 +560,13 @@ Scanner::traceShowParamXOR (const char *szLabel, const char *fmt, int row, int c va_end(arg); } -void +void Scanner::traceShowParamRasterOp (int iRasterOp, const char *szLabel, const char *fmt, int row, int color, va_list args) -{ +{ char szValue[256]; - + vsnprintf (szValue, sizeof(szValue), fmt, args); - + #ifdef HAVE_SGP if (m_pSGP) { m_pSGP->setRasterOp (iRasterOp); @@ -582,7 +588,7 @@ Scanner::traceShowParamRasterOp (int iRasterOp, const char *szLabel, const char m_pSGP->moveAbs (dXPos + dValueOffset, dYPos); m_pSGP->drawText (szValue); } - } else + } else #endif { cio_put_str (szLabel); @@ -591,41 +597,206 @@ Scanner::traceShowParamRasterOp (int iRasterOp, const char *szLabel, const char } } +void swap_xy_points (double& x1, double& y1, double& x2, double& y2) +{ + double temp = x1; x1 = x2; x2 = temp; + temp = y1; y1 = y2; y2 = temp; +} + +class WeightedPoint { +public: + int x, y; + double weight; + WeightedPoint (int _x, int _y, double _weight) + : x(_x), y(_y), weight(_weight) + {} +}; + + +/* FUNCTION + * Name: projection_pixel_weights + * Purpose: Returns a vector of WeightedPoint with the length of + * line that intersects with each pixel + */ + +void +projection_pixel_weights (std::vector& wp, const int nx, const int ny, + double x1, double y1, double x2, double y2) +{ + double ylen = fabs(y2-y1); + double xlen = fabs(x2-x1); + bool swap_xy = false, invert_slope = false; + double slope; + + if (ylen > xlen) { + swap_xy = true; + slope = xlen / ylen; + if (y2 < y1) // swap start/end so always moving from bottom to top + swap_xy_points (x1, y1, x2, y2); + if (x2 < x1) { + invert_slope = true; + } + } else { +#if DEBUG + if (ylen == xlen) + sys_error(ERR_WARNING, "Slope == 1"); +#endif + slope = ylen / xlen; + if (x2 < x1) // swap start/end so always moving from left to right in image + swap_xy_points (x1, y1, x2, y2); + if (y2 < y1) { + invert_slope = true; + } + } + double angle = atan(fabs(slope)); + double minor_dist = sin(angle); // distance along minor axis + double pixel_len = 1 / cos(angle); + + int minor_dir = 1; + if (invert_slope) { + minor_dir = -1; + slope = -slope; + } + + double x = x1, y = y1; + int ix = floor(x); + int iy = floor(y); + double ydelta = y - iy; + double xdelta = x - ix; + + double min_delta; + int *imaj, *imin; + int max_maj, max_min; + if (swap_xy) { + min_delta = xdelta; + imaj = &iy; + imin = &ix; + max_maj = ny; + max_min = nx; + } else { + min_delta = ydelta; + imaj = &ix; + imin = &iy; + max_maj = nx; + max_min = ny; + } + +#if DEBUG + sys_error(ERR_TRACE, "m=%6.3f swap_xy=%d invert=%d len=%8.6f min_delta=%.4g minor_dist=%6.3f (%.3f,%.3f)-(%.3f,%.3f)", + slope, swap_xy, invert_slope, pixel_len, min_delta, minor_dist, x1, y1, x2, y2); +#endif + + // if position of minor axis is at edge of image, but will be moving into pixel within image + if (*imin == max_min && invert_slope) { + (*imin)--; // select the pixel within image +#if DEBUG + sys_error(ERR_TRACE, "Moving pixel inside image, adding %f to min_delta", (1+slope)); +#endif + min_delta += (1+slope); + } + + while (*imaj < max_maj && *imin < max_min && *imin >= 0) { + double next_min_delta = min_delta + slope; + + if (((!invert_slope) && (next_min_delta < 1)) || + (invert_slope && (next_min_delta > 0))) { + // stay within same pixel + double w = pixel_len; + WeightedPoint p (ix, iy, w); + wp.push_back(p); +#if DEBUG + sys_error(ERR_TRACE, " Full pixel: (%3d,%3d)=%.4g, min_delta=%.4g", ix, iy, w, min_delta); +#endif + min_delta = next_min_delta; + } else { + // Scale partial pixel_len into pixel + double norm_delta = invert_slope ? min_delta : (1 - min_delta); + double p1_line = norm_delta * pixel_len; + WeightedPoint p1 (ix, iy, p1_line); + wp.push_back (p1); +#if DEBUG + sys_error(ERR_TRACE, " Part pixel: (%3d,%3d)=%.4g, min_delta=%.4g", ix, iy, p1_line, min_delta); +#endif + (*imin) += minor_dir; + min_delta = next_min_delta - minor_dir; + } + (*imaj)++; + } + +} /* NAME -* projectSingleLine INTERNAL: Calculates raysum along a line for a Phantom +* projectSingleLine INTERNAL: Calculates raysum along a line for a Phantom * * SYNOPSIS * rsum = phm_ray_attenuation (phm, x1, y1, x2, y2) -* double rsum Ray sum of Phantom along given line -* Phantom& phm; Phantom from which to calculate raysum -* double *x1, *y1, *x2, y2 Endpoints of ray path (in Phantom coords) +* double rsum Ray sum of Phantom along given line +* Phantom& phm; Phantom from which to calculate raysum +* double *x1, *y1, *x2, y2 Endpoints of ray path (in Phantom coords) */ -double -Scanner::projectSingleLine (const Phantom& phm, const double x1, const double y1, const double x2, const double y2) +double +Scanner::projectSingleLine (const Phantom& phm, double x1, double y1, double x2, double y2) { - // check ray against each pelem in Phantom double rsum = 0.0; - for (PElemConstIterator i = phm.listPElem().begin(); i != phm.listPElem().end(); i++) - rsum += projectLineAgainstPElem (**i, x1, y1, x2, y2); - + + if (phm.isImagefile()) { + // Project through an imagefile + + const ImageFile* im = phm.getImagefile(); + const ImageFileArray v = im->getArray(); + + // Get image axis extents + int nx = im->nx(), ny = im->ny(); + double xmin=0, xmax=nx, ymin=0, ymax=ny; // default coordinate + if (! im->getAxisExtent (xmin, xmax, ymin, ymax)) { + sys_error(ERR_WARNING, "Axis extent not available [Scanner::projectSingleLine]"); + } + + // Clip line in image object coordinates + double rect[4]; + rect[0] = xmin; rect[1] = ymin; + rect[2] = xmax; rect[3] = ymax; + bool accept = clip_rect (x1, y1, x2, y2, rect); + if (! accept) + return (0.0); + + // Convert to pixel coordinates + double xlen = xmax - xmin; + double ylen = ymax - ymin; + double px1 = nx * (x1 - xmin) / xlen; + double px2 = nx * (x2 - xmin) / xlen; + double py1 = ny * (y1 - ymin) / ylen; + double py2 = ny * (y2 - ymin) / ylen; + + std::vector wp; + projection_pixel_weights (wp, nx, ny, px1, py1, px2, py2); + for (unsigned int i = 0; i < wp.size(); i++) { + WeightedPoint& p = wp[i]; + rsum += v[p.x][p.y] * p.weight; + } + } else { + + // Project through each pelem in Phantom + for (PElemConstIterator i = phm.listPElem().begin(); i != phm.listPElem().end(); i++) + rsum += projectLineAgainstPElem (**i, x1, y1, x2, y2); + } return (rsum); } /* NAME -* pelem_ray_attenuation Calculate raysum of an pelem along one line +* pelem_ray_attenuation Calculate raysum of an pelem along one line * * SYNOPSIS * rsum = pelem_ray_attenuation (pelem, x1, y1, x2, y2) -* double rsum Computed raysum -* PhantomElement& pelem Pelem to scan -* double x1, y1, x2, y2 Endpoints of raysum line +* double rsum Computed raysum +* PhantomElement& pelem Pelem to scan +* double x1, y1, x2, y2 Endpoints of raysum line */ -double +double Scanner::projectLineAgainstPElem (const PhantomElement& pelem, double x1, double y1, double x2, double y2) { if (! pelem.clipLineWorldCoords (x1, y1, x2, y2)) { @@ -633,7 +804,7 @@ Scanner::projectLineAgainstPElem (const PhantomElement& pelem, double x1, double cio_tone (1000., 0.05); return (0.0); } - + #ifdef HAVE_SGP if (m_pSGP && m_trace == Trace::TRACE_CLIPPING) { m_pSGP->setRasterOp (RO_XOR); @@ -645,7 +816,7 @@ Scanner::projectLineAgainstPElem (const PhantomElement& pelem, double x1, double m_pSGP->setRasterOp (RO_SET); } #endif - + double len = lineLength (x1, y1, x2, y2); return (len * pelem.atten()); }