X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fscanner.cpp;h=34864fb018a6ce94a1ba63bd1cb2904777226614;hp=11cd0926ce68874d353bf3454c202b75298c2d1f;hb=747a2ec9e0f3c49723b36da0cc77270fbecc9dfe;hpb=1a050c98763fbbc0662731b0b76953acede6f5d7

diff --git a/libctsim/scanner.cpp b/libctsim/scanner.cpp
index 11cd092..34864fb 100644
--- a/libctsim/scanner.cpp
+++ b/libctsim/scanner.cpp
@@ -7,9 +7,7 @@
 **   Date Started:  1984
 **
 **  This is part of the CTSim program
-**  Copyright (c) 1983-2001 Kevin Rosenberg
-**
-**  $Id$
+**  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
@@ -257,13 +255,13 @@ 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);
 }
 
 
@@ -294,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,
@@ -302,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
@@ -355,11 +356,7 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS
       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 ();
 
@@ -372,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) {
@@ -404,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<double>(m_nView) * 100.);
 #endif
+
     if (m_trace == Trace::TRACE_CONSOLE)
       std::cout << "Current View: " << iView+iStartView << std::endl;
 
@@ -415,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 */
 }
 
@@ -615,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);
 }