Interim work on creating projections from an imagefile

diff --git a/include/phantom.h b/include/phantom.h
index 0b66b49077e2bfd8dc81b60972c33c78d9c98363..b1cb50fe9c7132584299c34bd896d1e802940946 100644 (file)
--- a/include/phantom.h
+++ b/include/phantom.h
@@ -103,11 +103,8 @@ class PhantomElement
     static const char* const convertTypeToName (PhmElemType iType);
 
     void makeTransformMatrices ();
-
     void makeVectorOutline ();
-
     void calcArcPoints (double x[], double y[], const int pts, const double xcent, const double ycent, const double r, const double start, const double stop);
-
     void calcEllipsePoints (double x[], double y[], const int pts, const double u, const double v);
 
     static int numCirclePoints (double theta);
@@ -150,15 +147,12 @@ class Phantom
         { return m_composition; }
 
     bool createFromPhantom (const char* const phmName);
-
     bool createFromPhantom (const int phmid);
-
-    bool createFromFile (const char* const fname);
-
+    bool createFromPhmFile (const char* const fname);
+    bool createFromImageFile (const char* const fname);
     bool fileWrite (const char* const fname);
 
     void addPElem (const PhantomElement& pelem);
-
     void addPElem (const char* const composition, const double cx, const double cy, const double u, const double v, const double rot, const double atten);
 
     void convertToImagefile (ImageFile& im, double dViewRatio, const int in_nsample, const int trace) const;
@@ -200,6 +194,8 @@ class Phantom
           std::list<PhantomElement*>& listPElem() {return m_listPElem;}
     const std::list<PhantomElement*>& listPElem() const {return m_listPElem;}
     const int nPElem() const {return m_nPElem;}
+    const bool isImagefile(void) const { return m_im != NULL; }
+    const ImageFile* getImagefile() const { return m_im; }
 
     static const int getPhantomCount() {return s_iPhantomCount;}
     static const char** getPhantomNameArray() {return s_aszPhantomName;}
@@ -220,6 +216,7 @@ class Phantom
     static const char* s_aszPhantomName[];
     static const char* s_aszPhantomTitle[];
     static const int s_iPhantomCount;
+    ImageFile* m_im;                  // If defining a phantom from an ImageFile
 
     void init();
 

diff --git a/libctsim/phantom.cpp b/libctsim/phantom.cpp
index bf1165499445375848e475861060d00a8762c443..2b969e87e25128808cc399a2cd23bd95aa5d709a 100644 (file)
--- a/libctsim/phantom.cpp
+++ b/libctsim/phantom.cpp
@@ -78,6 +78,7 @@ Phantom::init ()
   m_composition = P_PELEMS;
   m_fail = false;
   m_id = PHM_INVALID;
+  m_im = NULL;
 }
 
 Phantom::~Phantom ()
@@ -85,9 +86,10 @@ Phantom::~Phantom ()
   for (PElemIterator i = m_listPElem.begin(); i != m_listPElem.end(); i++) {
     delete *i;
   }
+  if (m_im)
+    delete m_im;
 }
 
-
 const char*
 Phantom::convertPhantomIDToName (int phmID)
 {
@@ -171,10 +173,10 @@ Phantom::createFromPhantom (const int phmid)
 
 
 /* METHOD IDENTIFICATION
-*   createFromFile          Add PhantomElements from file
+*   createFromPhmFile          Add PhantomElements from file
 *
 * SYNOPSIS
-*   createFromFile (filename)
+*   createFromPhmFile (filename)
 *
 * RETURNS
 *   true if pelem were added
@@ -182,7 +184,7 @@ Phantom::createFromPhantom (const int phmid)
 */
 
 bool
-Phantom::createFromFile (const char* const fname)
+Phantom::createFromPhmFile (const char* const fname)
 {
   bool bGoodFile = true;
   FILE *fp;
@@ -201,7 +203,7 @@ Phantom::createFromFile (const char* const fname)
     if (status == static_cast<int>(EOF))
       break;
     else if (status != 7) {
-      sys_error (ERR_WARNING, "Insufficient fields reading phantom file %s [Phantom::createFromFile]", fname);
+      sys_error (ERR_WARNING, "Insufficient fields reading phantom file %s [Phantom::createFromPhmFile]", fname);
       bGoodFile = false;
     }
     addPElem (pelemtype, cx, cy, u, v, rot, dens);
@@ -212,6 +214,32 @@ Phantom::createFromFile (const char* const fname)
   return (bGoodFile);
 }
 
+bool
+Phantom::createFromImageFile (const char* const fname)
+{
+  bool bGoodFile = true;
+
+  m_im = new ImageFile ();
+  if (! m_im || ! m_im->fileRead (fname)) {
+    sys_error (ERR_SEVERE, "Unable to read image file %s", fname);
+    bGoodFile = false;
+    delete m_im;
+    m_im = NULL;
+    m_fail = true;
+  } else {
+    m_name = fname;
+    m_id = -1;
+    m_fail = false;
+
+    if (! m_im->getAxisExtent(m_xmin, m_xmax, m_ymin, m_ymax)) {
+      m_xmax = m_im->nx() / 2.;  m_xmin = -m_xmax;
+      m_ymax = m_im->ny() / 2.;  m_ymin = -m_ymax;
+    }
+  }
+
+  return (bGoodFile);
+}
+
 bool
 Phantom::fileWrite (const char* const fname)
 {

diff --git a/libctsim/scanner.cpp b/libctsim/scanner.cpp
index 850b99eef9f8a4c9da0185ce89bd507cbc622913..34864fb018a6ce94a1ba63bd1cb2904777226614 100644 (file)
--- a/libctsim/scanner.cpp
+++ b/libctsim/scanner.cpp
@@ -312,10 +312,10 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS
 #ifdef HAVE_OPENMP
   #pragma omp parallel for
 #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);
@@ -326,12 +326,12 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS
       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 );
 
@@ -370,7 +370,6 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS
 
       m_pSGP->setMarker (SGP::MARKER_BDIAMOND);
 
-      
       m_pSGP->setColor (C_BLACK);
       m_pSGP->setPenWidth (2);
       if (m_idGeometry == GEOMETRY_PARALLEL) {
@@ -605,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);
 }
 

diff --git a/src/docs.cpp b/src/docs.cpp
index a47d8d993fe4f952e93edace295168acef3098aa..3f38883388fd67a5cd37bd2de07c11f51848f783 100644 (file)
--- a/src/docs.cpp
+++ b/src/docs.cpp
@@ -275,7 +275,7 @@ PhantomFileDocument::OnOpenDocument(const wxString& constFilename)
   wxString filename (constFilename);
 
   if (wxFile::Exists (filename)) {
-    m_phantom.createFromFile (filename.mb_str(wxConvUTF8));
+    m_phantom.createFromPhmFile (filename.mb_str(wxConvUTF8));
     if (theApp->getVerboseLogging())
       *theApp->getLog() << _T("Read phantom file ") << filename << _T("\n");
   } else {

diff --git a/tools/phm2helix.cpp b/tools/phm2helix.cpp
index 93a8390f56bc86614e5170d01626280f74f29ffa..5ead089814673a9678415ae5f1ca5648411925a2 100644 (file)
--- a/tools/phm2helix.cpp
+++ b/tools/phm2helix.cpp
@@ -281,7 +281,7 @@ phm2helix_main (int argc, char* const argv[])
           if (stat != 0 )
                         std::cerr << "Error executing external phantom program " << opt_PhmProg << " with command " << extcommand << std::endl;
 
-          phm.createFromFile (opt_PhmFileName.c_str());
+          phm.createFromPhmFile (opt_PhmFileName.c_str());
           remove(opt_PhmFileName.c_str());
 
           Scanner scanner (phm, optGeometryName.c_str(), opt_ndet, opt_nview,
@@ -309,7 +309,7 @@ phm2helix_main (int argc, char* const argv[])
                 if (stat != 0 )
                         std::cerr << "Error executing external phantom program " << opt_PhmProg << " with command " << extcommand << std::endl;
                 Phantom phmtmp;
-                phmtmp.createFromFile (opt_PhmFileName.c_str());
+                phmtmp.createFromPhmFile (opt_PhmFileName.c_str());
 
                 scanner.collectProjections (pjGlobal, phmtmp, iView,
                           1, scanner.offsetView(), true, opt_trace);

diff --git a/tools/phm2if.cpp b/tools/phm2if.cpp
index 040f7e17e6d264bddcb05e57bfd12e18c610ca18..13f5b922705fad8074ea5004e34441169e9e9a04 100644 (file)
--- a/tools/phm2if.cpp
+++ b/tools/phm2if.cpp
@@ -241,33 +241,34 @@ phm2if_main (int argc, char* const argv[])
 
     std::ostringstream oss;
     oss << "phm2if: nx=" << opt_nx << ", ny=" << opt_ny << ", viewRatio=" << optViewRatio << ", nsample=" << opt_nsample << ", ";
-    if (optPhmFilename != "")
+    if (!optPhmFilename.empty()) {
       oss << "phantomFile=" << optPhmFilename;
-    else if (optPhmName != "")
+    } else if (!optPhmName.empty()) {
       oss << "phantom=" << optPhmName;
-    else if (optFilterName != "") {
+    } else if (!optFilterName.empty()) {
       oss << "filter=" << optFilterName << " - " << optDomainName;
     }
-    if (optDesc != "")
+    if (!optDesc.empty())
       oss << ": " << optDesc;
     optDesc = oss.str();
 
-    if (optPhmName != "") {
+    if (! optPhmName.empty()) {
       phm.createFromPhantom (optPhmName.c_str());
       if (phm.fail()) {
         std::cout << phm.failMessage() << std::endl << std::endl;
         phm2if_usage(argv[0]);
         return (1);
       }
-    }
-
-    if (optPhmFilename != "") {
-      phm.createFromFile(optPhmFilename.c_str());
+    } else if (!optPhmFilename.empty()) {
+      phm.createFromPhmFile(optPhmFilename.c_str());
 #ifdef HAVE_MPI
       if (mpiWorld.getRank() == 0)
         std::cerr << "Can't use phantom from file in MPI mode\n";
       return (1);
 #endif
+    } else {
+      std::cerr << "Internal error\n";
+      return (1);
     }
 
     if (optVerbose)

diff --git a/tools/phm2pj.cpp b/tools/phm2pj.cpp
index 7fa9bde79fe8133c34b3472905d70b4e51237533..2bd7e3facb045a8acd4877ac8288b3feffd72832 100644 (file)
--- a/tools/phm2pj.cpp
+++ b/tools/phm2pj.cpp
@@ -27,13 +27,14 @@
 #include "timer.h"
 
 
-enum { O_PHANTOM, O_DESC, O_NRAY, O_ROTANGLE, O_PHMFILE, O_GEOMETRY, O_FOCAL_LENGTH, O_CENTER_DETECTOR_LENGTH,
+enum { O_PHANTOM, O_DESC, O_NRAY, O_ROTANGLE, O_PHMFILE, O_IMAGEFILE, O_GEOMETRY, O_FOCAL_LENGTH, O_CENTER_DETECTOR_LENGTH,
 O_VIEW_RATIO, O_SCAN_RATIO, O_OFFSETVIEW, O_TRACE, O_VERBOSE, O_HELP, O_DEBUG, O_VERSION };
 
 static struct option phm2pj_options[] =
 {
   {"phantom", 1, 0, O_PHANTOM},
   {"phmfile", 1, 0, O_PHMFILE},
+  {"imagefile", 1, 0, O_IMAGEFILE},
   {"desc", 1, 0, O_DESC},
   {"nray", 1, 0, O_NRAY},
   {"rotangle", 1, 0, O_ROTANGLE},
@@ -62,12 +63,14 @@ phm2pj_usage (const char *program)
   std::cout << "     outfile          Name of output file for projections\n";
   std::cout << "     ndet             Number of detectors\n";
   std::cout << "     nview            Number of rotated views\n";
+  std::cout << "     SELECT PHANTOM   Use either --phantom, --phmfile, or --imagefile\n";
   std::cout << "     --phantom        Phantom to use for projection\n";
   std::cout << "        herman        Herman head phantom\n";
   std::cout << "        shepp-logan   Shepp-Logan head phantom\n";
   std::cout << "        unit-pulse     Unit pulse phantom\n";
   std::cout << "     --phmfile        Get Phantom from phantom file\n";
-  std::cout << "     --desc           Description of raysum\n";
+  std::cout << "     --imagefile      Get Phantom from imagefile\n";
+  std::cout << "     --desc           Add description for projection file\n";
   std::cout << "     --nray           Number of rays per detector (default = 1)\n";
   std::cout << "     --rotangle       Angle to rotate view through (fraction of a circle)\n";
   std::cout << "                      (default = select appropriate for geometry)\n";
@@ -106,6 +109,7 @@ phm2pj_main (int argc, char* const argv[])
   std::string opt_desc;
   std::string optPhmName;
   std::string optPhmFileName;
+  std::string optImageFileName;
   int opt_ndet;
   int opt_nview;
   int opt_offsetview = 0;
@@ -121,7 +125,7 @@ phm2pj_main (int argc, char* const argv[])
   char* endptr = NULL;
   char* endstr;
   UNUSED(opt_debug);
-  
+
 #ifdef HAVE_MPI
   MPIWorld mpiWorld (argc, argv);
 #endif
@@ -144,6 +148,9 @@ phm2pj_main (int argc, char* const argv[])
       case O_PHMFILE:
         optPhmFileName = optarg;
         break;
+      case O_IMAGEFILE:
+        optImageFileName = optarg;
+        break;
       case O_VERBOSE:
         opt_verbose = 1;
         break;
@@ -244,11 +251,19 @@ phm2pj_main (int argc, char* const argv[])
       }
     }
 
-    if (optPhmName == "" && optPhmFileName == "") {
+    if (optPhmName.empty() && optPhmFileName.empty() && optImageFileName.empty()) {
       std::cerr << "No phantom defined\n" << std::endl;
       phm2pj_usage(argv[0]);
       return (1);
     }
+    if ((! optPhmName.empty() && !optPhmFileName.empty()) ||
+        (! optPhmName.empty() && !optImageFileName.empty()) ||
+        (! optPhmFileName.empty() && !optImageFileName.empty())) {
+      std::cerr << "Multiple phantoms define. Use only --phantom, --phmfile, or --imagefile\n" << std::endl;
+      phm2pj_usage(argv[0]);
+      return (1);
+    }
+
     if (optind + 3 != argc) {
       phm2pj_usage(argv[0]);
       return (1);
@@ -279,31 +294,38 @@ phm2pj_main (int argc, char* const argv[])
 
     std::ostringstream desc;
     desc << "phm2pj: NDet=" << opt_ndet << ", Nview=" << opt_nview << ", NRay=" << opt_nray << ", RotAngle=" << opt_rotangle << ", OffsetView =" << opt_offsetview << ", Geometry=" << optGeometryName << ", ";
-    if (optPhmFileName.length()) {
+
+    if (!optPhmFileName.empty()) {
       desc << "PhantomFile=" << optPhmFileName;
-    } else if (optPhmName != "") {
+    } else if (!optPhmName.empty()) {
       desc << "Phantom=" << optPhmName;
+    } else if (!optImageFileName.empty()) {
+      desc << "Imagefile=" << optImageFileName;
     }
     if (opt_desc.length()) {
       desc << ": " << opt_desc;
     }
     opt_desc = desc.str();
 
-    if (optPhmName != "") {
+    if (!optPhmName.empty()) {
       phm.createFromPhantom (optPhmName.c_str());
       if (phm.fail()) {
         std::cout << phm.failMessage() << std::endl << std::endl;
         phm2pj_usage(argv[0]);
         return (1);
       }
-    }
-
-    if (optPhmFileName != "") {
+    } else if (!optPhmFileName.empty()) {
 #ifdef HAVE_MPI
       std::cerr << "Can not read phantom from file in MPI mode\n";
       return (1);
 #endif
-      phm.createFromFile (optPhmFileName.c_str());
+      phm.createFromPhmFile (optPhmFileName.c_str());
+    } else if (!optImageFileName.empty()) {
+#ifdef HAVE_MPI
+      std::cerr << "Can not read image from file in MPI mode\n";
+      return (1);
+#endif
+      phm.createFromImageFile (optImageFileName.c_str());
     }
 
 #ifdef HAVE_MPI