From 747a2ec9e0f3c49723b36da0cc77270fbecc9dfe Mon Sep 17 00:00:00 2001 From: "Kevin M. Rosenberg" Date: Wed, 21 Mar 2018 16:04:44 -0600 Subject: [PATCH] Interim work on creating projections from an imagefile --- include/phantom.h | 13 +++++------- libctsim/phantom.cpp | 38 +++++++++++++++++++++++++++++----- libctsim/scanner.cpp | 49 ++++++++++++++++++++++++++++++++++++-------- src/docs.cpp | 2 +- tools/phm2helix.cpp | 4 ++-- tools/phm2if.cpp | 19 +++++++++-------- tools/phm2pj.cpp | 44 +++++++++++++++++++++++++++++---------- 7 files changed, 124 insertions(+), 45 deletions(-) diff --git a/include/phantom.h b/include/phantom.h index 0b66b49..b1cb50f 100644 --- 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& listPElem() {return m_listPElem;} const std::list& 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 bf11654..2b969e8 100644 --- 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(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 850b99e..34864fb 100644 --- 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 a47d8d9..3f38883 100644 --- 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 93a8390..5ead089 100644 --- 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 040f7e1..13f5b92 100644 --- 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 7fa9bde..2bd7e3f 100644 --- 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 -- 2.34.1