X-Git-Url: http://git.kpe.io/?a=blobdiff_plain;f=libctsim%2Fscanner.cpp;h=34864fb018a6ce94a1ba63bd1cb2904777226614;hb=0219ede69e1c3afc6e160b8f276bfd4617acbc08;hp=89745061240059b7aa3aaeece7dad440afcad907;hpb=e8462f7431582627e44906239077f1c696eefba1;p=ctsim.git diff --git a/libctsim/scanner.cpp b/libctsim/scanner.cpp index 8974506..34864fb 100644 --- a/libctsim/scanner.cpp +++ b/libctsim/scanner.cpp @@ -292,6 +292,15 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS 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, const int iOffsetView, int iStorageOffset, @@ -300,36 +309,30 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS m_trace = trace; double start_angle = (iStartView + iOffsetView) * 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 - } +#ifdef HAVE_OPENMP + #pragma omp parallel for +#endif - 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 + for (int iView = 0; iView < iNumViews; iView++) { + double viewAngle = start_angle + (iView * proj.rotInc()); - int iView; - double viewAngle; - for (iView = 0, viewAngle = start_angle; iView < iNumViews; iView++, viewAngle += proj.rotInc()) { - int iStoragePosition = iView + iStorageOffset; + // 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 @@ -366,11 +369,7 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS 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) { @@ -398,6 +397,7 @@ 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; @@ -409,12 +409,7 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS // 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 */ } @@ -609,13 +604,45 @@ Scanner::traceShowParamRasterOp (int iRasterOp, const char *szLabel, const char */ double -Scanner::projectSingleLine (const Phantom& phm, const double x1, const double y1, const double x2, const double y2) +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(); + int nx = im->nx(), ny = im->ny(); + + // Convert endpoints into image pixel coordinates + 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]"); + } + 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); + + double xlen = xmax - xmin, 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; + + // Use Bresenham integer line stepping to step to each pixel, and walked through image + + sys_error(ERR_WARNING, "Not yet able to project through imagefile. Line (%.3f,%.3f) - (%.3f,%.3f)", px1, py1, px2, py2); + + } 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); }