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);
{ 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;
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;}
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();
m_composition = P_PELEMS;
m_fail = false;
m_id = PHM_INVALID;
+ m_im = NULL;
}
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)
{
/* METHOD IDENTIFICATION
-* createFromFile Add PhantomElements from file
+* createFromPhmFile Add PhantomElements from file
*
* SYNOPSIS
-* createFromFile (filename)
+* createFromPhmFile (filename)
*
* RETURNS
* true if pelem were added
*/
bool
-Phantom::createFromFile (const char* const fname)
+Phantom::createFromPhmFile (const char* const fname)
{
bool bGoodFile = true;
FILE *fp;
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);
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)
{
#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);
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 );
m_pSGP->setMarker (SGP::MARKER_BDIAMOND);
-
m_pSGP->setColor (C_BLACK);
m_pSGP->setPenWidth (2);
if (m_idGeometry == GEOMETRY_PARALLEL) {
*/
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);
}
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 {
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,
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);
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)
#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},
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";
std::string opt_desc;
std::string optPhmName;
std::string optPhmFileName;
+ std::string optImageFileName;
int opt_ndet;
int opt_nview;
int opt_offsetview = 0;
char* endptr = NULL;
char* endstr;
UNUSED(opt_debug);
-
+
#ifdef HAVE_MPI
MPIWorld mpiWorld (argc, argv);
#endif
case O_PHMFILE:
optPhmFileName = optarg;
break;
+ case O_IMAGEFILE:
+ optImageFileName = optarg;
+ break;
case O_VERBOSE:
opt_verbose = 1;
break;
}
}
- 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);
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