X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fphantom.cpp;h=9125e36be718af7ea36959b42bf5461cd673f0ab;hp=2063683d0675d8532a2a31feb635c4ecf8e3a82a;hb=f7d2b7144f32a7bd157b7689022e62944b82fcc1;hpb=286d655a25df2668bd65ad365676c6ecc94415a1

diff --git a/libctsim/phantom.cpp b/libctsim/phantom.cpp
index 2063683..9125e36 100644
--- a/libctsim/phantom.cpp
+++ b/libctsim/phantom.cpp
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: phantom.cpp,v 1.10 2000/07/22 15:45:33 kevin Exp $
+**  $Id: phantom.cpp,v 1.23 2000/12/21 03:40:58 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
@@ -27,29 +27,31 @@
 
 #include "ct.h"
 
+const int PhantomElement::POINTS_PER_CIRCLE = 360;
+const double PhantomElement::SCALE_PELEM_EXTENT=0.005;  // increase pelem limits by 0.5% 
 
 const int Phantom::PHM_INVALID = -1;
 const int Phantom::PHM_HERMAN = 0;
-const int Phantom::PHM_BHERMAN = 1;
-const int Phantom::PHM_ROWLAND = 2;
-const int Phantom::PHM_BROWLAND = 3;
+const int Phantom::PHM_B_HERMAN = 1;
+const int Phantom::PHM_SHEPP_LOGAN = 2;
+const int Phantom::PHM_B_SHEPP_LOGAN = 3;
 const int Phantom::PHM_UNITPULSE = 4;
 
 const char* Phantom::s_aszPhantomName[] = 
 {
   {"herman"},
-  {"bherman"},
-  {"rowland"},
-  {"browland"},
+  {"herman-b"},
+  {"shepp-logan"},
+  {"shepp-logan-b"},
   {"unitpulse"},
 };
 
 const char* Phantom::s_aszPhantomTitle[] = 
 {
   {"Herman Head"},
-  {"Herman Head Bordered"},
-  {"Rowland Head"},
-  {"Rowland Head Bordered"},
+  {"Herman Head (Bordered)"},
+  {"Shepp-Logan"},
+  {"Shepp-Logan (Bordered)"},
   {"Unit Pulse"},
 };
 
@@ -155,14 +157,14 @@ Phantom::createFromPhantom (const int phmid)
     case PHM_HERMAN:
       addStdHerman();
       break;
-    case PHM_BHERMAN:
+    case PHM_B_HERMAN:
       addStdHermanBordered();
       break;
-    case PHM_ROWLAND:
-      addStdRowland();
+    case PHM_SHEPP_LOGAN:
+      addStdSheppLogan();
       break;
-    case PHM_BROWLAND:
-      addStdRowlandBordered ();
+    case PHM_B_SHEPP_LOGAN:
+      addStdSheppLoganBordered();
       break;
     case PHM_UNITPULSE:
       m_composition = P_UNIT_PULSE;
@@ -196,34 +198,32 @@ Phantom::createFromPhantom (const int phmid)
 bool
 Phantom::createFromFile (const char* const fname)
 {
-  bool stoploop = false;
-  bool retval = false;
+  bool bGoodFile = true;
   FILE *fp;
 
   if ((fp = fopen (fname, "r")) == NULL)
     return (false);
 
-  do {
+  m_name = fname;
+
+  while (1) {
     double cx, cy, u, v, rot, dens;
     char pelemtype[80];
-    int n = fscanf (fp, "%79s %lf %lf %lf %lf %lf %lf",
-		pelemtype, &cx, &cy, &u, &v, &rot, &dens);
+
+    int status = fscanf (fp, "%79s %lf %lf %lf %lf %lf %lf", pelemtype, &cx, &cy, &u, &v, &rot, &dens);
     
-    if (n == EOF || n == 0) {	/* end of file */
-      stoploop = true;
-      retval = false;
-    } else if (n != 7) {
-      stoploop = true;
-      retval = false;
-    } else {
-      addPElem (pelemtype, cx, cy, u, v, rot, dens);
-      retval = true;
+    if (status == static_cast<int>(EOF)) 
+      break;
+    else if (status != 7) {
+      sys_error (ERR_WARNING, "Insufficient fields reading phantom file %s [Phantom::createFromFile]", fname);
+      bGoodFile = false;
     }
-  } while (stoploop == false);
+    addPElem (pelemtype, cx, cy, u, v, rot, dens);
+  }
   
   fclose (fp);
 
-  return (retval);
+  return (bGoodFile);
 }
 
 
@@ -306,29 +306,30 @@ Phantom::print () const
 void 
 Phantom::show () const
 {
-  double wsize = m_xmax - m_xmin;
-  double xmin = m_xmin;
-  double ymin = m_ymin;
-  double xmax, ymax;
-  SGP_ID gid;
+  SGPDriver driverSGP ("Phantom Show");
+  SGP sgp (driverSGP);
 
-  if ((m_ymax - m_ymin) > wsize)
-      wsize = m_ymax - m_ymin;
-  wsize *= 1.1;
+  show (sgp);
 
-  xmax = xmin + wsize;
-  ymax = ymin + wsize; 
-  
-  printf("Drawing Phantom:\n\n");
-  printf("    data limits: %9.3g, %9.3g, %9.3g, %9.3g\n", m_xmin, m_ymin, m_xmax, m_ymax);
-  printf("    window size: %9.3g, %9.3g, %9.3g, %9.3g\n", xmin, ymin, xmax, ymax);
+  std::cout << "Press return to continue";
+  cio_kb_getc();
+}
 
-  gid = sgp2_init (0, 0, "Phantom Show");
-  sgp2_window (xmin, ymin, xmax, ymax);
+void 
+Phantom::show (SGP& sgp) const
+{
+  double wsize = m_xmax - m_xmin;
+  if ((m_ymax - m_ymin) > wsize) 
+      wsize = m_ymax - m_ymin;
+  wsize *= 1.01;
+  double halfWindow = wsize / 2;
 
-  draw();
+  double xcent = m_xmin + (m_xmax - m_xmin) / 2;
+  double ycent = m_ymin + (m_ymax - m_ymin) / 2;
 
-  sgp2_close (gid);
+  sgp.setWindow (xcent - halfWindow, ycent - halfWindow, xcent + halfWindow, ycent + halfWindow);
+
+  draw (sgp);
 }
 #endif
 
@@ -342,16 +343,16 @@ Phantom::show () const
 
 #ifdef HAVE_SGP
 void 
-Phantom::draw () const
+Phantom::draw (SGP& sgp) const
 {
   for (PElemIterator i = m_listPElem.begin(); i != m_listPElem.end(); i++)
-    sgp2_polyline_abs ((*i)->xOutline(), (*i)->yOutline(), (*i)->nOutlinePoints());
+    sgp.polylineAbs ((*i)->xOutline(), (*i)->yOutline(), (*i)->nOutlinePoints());
 }
 #endif
 
 
 /* NAME
- *   addStdRowland		Make head phantom of S.W. Rowland
+ *   addStdSheppLogan	Make head phantom of Shepp-Logan
  *
  * REFERENCES
  *   S. W. Rowland, "Computer Implementation of Image Reconstruction
@@ -360,7 +361,7 @@ Phantom::draw () const
  */
 
 void 
-Phantom::addStdRowland ()
+Phantom::addStdSheppLogan ()
 {
   addPElem ("ellipse",  0.0000,  0.0000, 0.6900,  0.9200,   0.0,  1.00);
   addPElem ("ellipse",  0.0000, -0.0184, 0.6624,  0.8740,   0.0, -0.98);
@@ -376,10 +377,10 @@ Phantom::addStdRowland ()
 }
 
 void 
-Phantom::addStdRowlandBordered ()
+Phantom::addStdSheppLoganBordered ()
 {
-  addStdRowland ();
-  addPElem ("ellipse", 0.000, 0.0000, 0.7500, 1.000, 0.0, 0.00);
+  addStdSheppLogan ();
+  addPElem ("rectangle", 0.000, 0.0000, 0.8600, 1.150, 0.0, 0.00);
 }
 
 /* NAME
@@ -414,7 +415,7 @@ void
 Phantom::addStdHermanBordered ()
 {
   addStdHerman();
-  addPElem ("ellipse",  0.000, 0.000, 8.650, 8.650,  0.00, 0.000);
+  addPElem ("rectangle",  0.000, 0.00, 10.780, 8.110,  90.00, 0.000);
 }
 
 
@@ -432,11 +433,11 @@ Phantom::addStdHermanBordered ()
 void
 Phantom::convertToImagefile (ImageFile& im, const int in_nsample, const int trace) const
 {
-  convertToImagefile (im, in_nsample, trace, 0, im.nx());
+  convertToImagefile (im, in_nsample, trace, 0, im.nx(), true);
 }
 
 void 
-Phantom::convertToImagefile (ImageFile& im, const int in_nsample, const int trace, const int colStart, const int colCount) const
+Phantom::convertToImagefile (ImageFile& im, const int in_nsample, const int trace, const int colStart, const int colCount, bool bStoreAtColumnPos) const
 {
   int nx = im.nx();
   int ny = im.ny();
@@ -478,8 +479,12 @@ Phantom::convertToImagefile (ImageFile& im, const int in_nsample, const int trac
   ImageFileArray v = im.getArray();
 
   for (int ix = 0; ix < colCount; ix++)
-      for (int iy = 0; iy < ny; iy++)
-	  v[ix][iy] = 0;
+    for (int iy = 0; iy < ny; iy++) {
+      int iColStore = ix;
+      if (bStoreAtColumnPos)
+		iColStore += colStart;
+      v[iColStore][iy] = 0;
+    }
 
   double x_start = xmin + (colStart * xinc);
   for (PElemConstIterator pelem = m_listPElem.begin(); pelem != m_listPElem.end(); pelem++) {
@@ -487,11 +492,14 @@ Phantom::convertToImagefile (ImageFile& im, const int in_nsample, const int trac
       double x, y, xi, yi;
       int ix, iy, kx, ky;
       for (ix = 0, x = x_start; ix < colCount; ix++, x += xinc) {
-	  for (iy = 0, y = ymin; iy < ny; iy++, y += yinc) {
-	      for (kx = 0, xi = x + kxofs; kx < nsample; kx++, xi += kxinc) {
-		  for (ky = 0, yi = y + kyofs; ky < nsample; ky++, yi += kyinc)
-		      if (rPElem.isPointInside (xi, yi, PHM_COORD) == TRUE)
-			  v[ix][iy] += rPElem.atten();
+	int iColStore = ix;
+	if (bStoreAtColumnPos)
+	  iColStore += colStart;
+	for (iy = 0, y = ymin; iy < ny; iy++, y += yinc) {
+	  for (kx = 0, xi = x + kxofs; kx < nsample; kx++, xi += kxinc) {
+	    for (ky = 0, yi = y + kyofs; ky < nsample; ky++, yi += kyinc)
+	      if (rPElem.isPointInside (xi, yi, PHM_COORD) == TRUE)
+			v[iColStore][iy] += rPElem.atten();
 	      } // for kx
 	  } /* for iy */
       }  /* for ix */
@@ -499,11 +507,16 @@ Phantom::convertToImagefile (ImageFile& im, const int in_nsample, const int trac
   
 
   if (nsample > 1) {
-    double factor = 1.0 / (nsample * nsample);
-
-    for (int ix = 0; ix < colCount; ix++)
-      for (int iy = 0; iy < ny; iy++)
-	v[ix][iy] *= factor;
+    double factor = 1.0 / static_cast<double>(nsample * nsample);
+
+
+    for (int ix = 0; ix < colCount; ix++) {
+		int iColStore = ix;
+		if (bStoreAtColumnPos)
+			iColStore += colStart;
+		for (int iy = 0; iy < ny; iy++)
+			v[iColStore][iy] *= factor;
+	}
   }
 }