X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fphantom.cpp;h=7acfc3fc9ccba40168f889867ed5b404b04bd299;hp=e28074f5917f5555c0c017f9a62a4e6ef3c02d82;hb=e4c1f7f8eb87558c3abf3bf1d20732361f425351;hpb=fff4beb84fcc84e65e4feb457e2ed25c7774cff4 diff --git a/libctsim/phantom.cpp b/libctsim/phantom.cpp index e28074f..7acfc3f 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.2 2000/06/19 19:54:23 kevin Exp $ +** $Id: phantom.cpp,v 1.7 2000/07/15 08:36:13 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 @@ -33,6 +33,19 @@ // Phantom::Phantom (void) +{ + init (); +} + + +Phantom::Phantom (const char* const phmName) +{ + init (); + createFromPhantom (phmName); +} + +void +Phantom::init (void) { m_nPElem = 0; m_xmin = 1E30; @@ -41,9 +54,10 @@ Phantom::Phantom (void) m_ymax = -1E30; m_diameter = 0; m_composition = P_PELEMS; + m_fail = false; + m_id = PHM_INVALID; } - Phantom::~Phantom (void) { for (PElemIterator i = m_listPElem.begin(); i != m_listPElem.end(); i++) { @@ -52,29 +66,93 @@ Phantom::~Phantom (void) } -void -Phantom::create (const int phmid) +const char* +Phantom::convertPhantomIDToName (const PhantomID phmID) +{ + const char *name = ""; + + if (phmID == PHM_HERMAN) + name = PHM_HERMAN_STR; + else if (phmID == PHM_BHERMAN) + name = PHM_BHERMAN_STR; + else if (phmID == PHM_ROWLAND) + name = PHM_ROWLAND_STR; + else if (phmID == PHM_BROWLAND) + name = PHM_BROWLAND_STR; + else if (phmID == PHM_UNITPULSE) + name = PHM_UNITPULSE_STR; + + return (name); +} + +Phantom::PhantomID +Phantom::convertNameToPhantomID (const char* const phmName) +{ + PhantomID id = PHM_INVALID; + + if (strcasecmp (phmName, PHM_HERMAN_STR) == 0) + id = PHM_HERMAN; + else if (strcasecmp (phmName, PHM_BHERMAN_STR) == 0) + id = PHM_BHERMAN; + else if (strcasecmp (phmName, PHM_ROWLAND_STR) == 0) + id = PHM_ROWLAND; + else if (strcasecmp (phmName, PHM_BROWLAND_STR) == 0) + id = PHM_BROWLAND; + else if (strcasecmp (phmName, PHM_UNITPULSE_STR) == 0) + id = PHM_UNITPULSE; + + return (id); +} + + +bool +Phantom::createFromPhantom (const char* const phmName) +{ + PhantomID phmid = convertNameToPhantomID (phmName); + if (phmid == PHM_INVALID) { + m_fail = true; + m_failMessage = "Invalid phantom name "; + m_failMessage += phmName; + return false; + } + + m_name = phmName; + createFromPhantom (phmid); + return true; +} + +bool +Phantom::createFromPhantom (const PhantomID phmid) { switch (phmid) { - case O_PHM_HERMAN: - std_herman(); + case PHM_HERMAN: + addStdHerman(); + break; + case PHM_BHERMAN: + addStdHermanBordered(); break; - case O_PHM_ROWLAND: - std_rowland(); + case PHM_ROWLAND: + addStdRowland(); break; - case O_PHM_BROWLAND: - std_rowland_bordered (); + case PHM_BROWLAND: + addStdRowlandBordered (); break; - case O_PHM_UNITPULSE: + case PHM_UNITPULSE: m_composition = P_UNIT_PULSE; addPElem ("rectangle", 0., 0., 100., 100., 0., 0.); // outline addPElem ("ellipse", 0., 0., 1., 1., 0., 1.); // pulse break; default: - sys_error (ERR_WARNING, "Illegal phantom id %d\n", phmid); - break; + m_fail = true; + m_failMessage = "Illegal phantom id "; + m_failMessage += phmid; + return false; } + + m_id = phmid; + + return true; } @@ -216,17 +294,15 @@ Phantom::show (void) const 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); + 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); - gid = sgp2_init(0, 0, "Phantom Show"); + gid = sgp2_init (0, 0, "Phantom Show"); sgp2_window (xmin, ymin, xmax, ymax); draw(); - termgrf2(); + sgp2_close (gid); } #endif @@ -249,10 +325,7 @@ Phantom::draw (void) const /* NAME - * std_rowland Make head phantom of S.W. Rowland - * - * SYNOPSIS - * std_rowland () + * addStdRowland Make head phantom of S.W. Rowland * * REFERENCES * S. W. Rowland, "Computer Implementation of Image Reconstruction @@ -261,33 +334,30 @@ Phantom::draw (void) const */ void -Phantom::std_rowland (void) +Phantom::addStdRowland (void) { - 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); - addPElem("ellipse", 0.2200, 0.0000, 0.1100, 0.3100, -18.0, -0.02); - addPElem("ellipse", -0.2200, 0.0000, 0.1600, 0.4100, 18.0, -0.02); - addPElem("ellipse", 0.0000, 0.3500, 0.2100, 0.2500, 0.0, 0.01); - addPElem("ellipse", 0.0000, 0.1000, 0.0460, 0.0460, 0.0, 0.01); - addPElem("ellipse", 0.0000, -0.1000, 0.0460, 0.0460, 0.0, 0.01); - addPElem("ellipse", -0.0800, -0.6050, 0.0460, 0.0230, 0.0, 0.01); - addPElem("ellipse", 0.0000, -0.6050, 0.0230, 0.0230, 0.0, 0.01); - addPElem("ellipse", 0.0600, -0.6050, 0.0230, 0.0230, 0.0, 0.01); - addPElem("ellipse", 0.5538, -0.3858, 0.0330, 0.2060, -18.0, 0.03); + 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); + addPElem ("ellipse", 0.2200, 0.0000, 0.1100, 0.3100, -18.0, -0.02); + addPElem ("ellipse", -0.2200, 0.0000, 0.1600, 0.4100, 18.0, -0.02); + addPElem ("ellipse", 0.0000, 0.3500, 0.2100, 0.2500, 0.0, 0.01); + addPElem ("ellipse", 0.0000, 0.1000, 0.0460, 0.0460, 0.0, 0.01); + addPElem ("ellipse", 0.0000, -0.1000, 0.0460, 0.0460, 0.0, 0.01); + addPElem ("ellipse", -0.0800, -0.6050, 0.0460, 0.0230, 0.0, 0.01); + addPElem ("ellipse", 0.0000, -0.6050, 0.0230, 0.0230, 0.0, 0.01); + addPElem ("ellipse", 0.0600, -0.6050, 0.0230, 0.0230, 0.0, 0.01); + addPElem ("ellipse", 0.5538, -0.3858, 0.0330, 0.2060, -18.0, 0.03); } void -Phantom::std_rowland_bordered (void) +Phantom::addStdRowlandBordered (void) { - std_rowland (); + addStdRowland (); addPElem ("ellipse", 0.000, 0.0000, 0.7500, 1.000, 0.0, 0.00); } /* NAME - * std_herman Standard head phantom of G. T. Herman - * - * SYNOPSIS - * std_herman () + * addStdHerman Standard head phantom of G. T. Herman * * REFERENCES * G. T. Herman, "Image Reconstructions from Projections: The Fundementals @@ -295,23 +365,30 @@ Phantom::std_rowland_bordered (void) */ void -Phantom::std_herman (void) +Phantom::addStdHerman (void) +{ + addPElem ("ellipse", 0.000, 1.50, 0.375, 0.3000, 90.00, -0.003); + addPElem ("ellipse", 0.675, -0.75, 0.225, 0.1500, 140.00, 0.010); + addPElem ("ellipse", 0.750, 1.50, 0.375, 0.2250, 50.00, 0.003); + addPElem ("segment", 1.375, -7.50, 1.100, 0.6250, 19.20, -0.204); + addPElem ("segment", 1.375, -7.50, 1.100, 4.3200, 19.21, 0.204); + addPElem ("segment", 0.000, -2.25, 1.125, 0.3750, 0.00, -0.003); + addPElem ("segment", 0.000, -2.25, 1.125, 3.0000, 0.00, 0.003); + addPElem ("segment", -1.000, 3.75, 1.000, 0.5000, 135.00, -0.003); + addPElem ("segment", -1.000, 3.75, 1.000, 3.0000, 135.00, 0.003); + addPElem ("segment", 1.000, 3.75, 1.000, 0.5000, 225.00, -0.003); + addPElem ("segment", 1.000, 3.75, 1.000, 3.0000, 225.00, 0.003); + addPElem ("triangle", 5.025, 3.75, 1.125, 0.5000, 110.75, 0.206); + addPElem ("triangle",-5.025, 3.75, 1.125, 0.9000,-110.75, 0.206); + addPElem ("ellipse", 0.000, 0.00, 8.625, 6.4687, 90.00, 0.416); + addPElem ("ellipse", 0.000, 0.00, 7.875, 5.7187, 90.00, -0.206); +} + +void +Phantom::addStdHermanBordered (void) { - addPElem("ellipse", 0.000, 1.50, 0.375, 0.3000, 90.00, -0.003); - addPElem("ellipse", 0.675, -0.75, 0.225, 0.1500, 140.00, 0.010); - addPElem("ellipse", 0.750, 1.50, 0.375, 0.2250, 50.00, 0.003); - addPElem("segment", 1.375, -7.50, 1.100, 0.6250, 19.20, -0.204); - addPElem("segment", 1.375, -7.50, 1.100, 4.3200, 19.21, 0.204); - addPElem("segment", 0.000, -2.25, 1.125, 0.3750, 0.00, -0.003); - addPElem("segment", 0.000, -2.25, 1.125, 3.0000, 0.00, 0.003); - addPElem("segment", -1.000, 3.75, 1.000, 0.5000, 135.00, -0.003); - addPElem("segment", -1.000, 3.75, 1.000, 3.0000, 135.00, 0.003); - addPElem("segment", 1.000, 3.75, 1.000, 0.5000, 225.00, -0.003); - addPElem("segment", 1.000, 3.75, 1.000, 3.0000, 225.00, 0.003); - addPElem("triangle", 5.025, 3.75, 1.125, 0.5000, 110.75, 0.206); - addPElem("triangle",-5.025, 3.75, 1.125, 0.9000,-110.75, 0.206); - addPElem("ellipse", 0.000, 0.00, 8.625, 6.4687, 90.00, 0.416); - addPElem("ellipse", 0.000, 0.00, 7.875, 5.7187, 90.00, -0.206); + addStdHerman(); + addPElem ("ellipse", 0.000, 0.000, 8.650, 8.650, 0.00, 0.000); } @@ -326,18 +403,24 @@ Phantom::std_herman (void) * (total samples per pixel = nsample * nsample) */ +void +Phantom::convertToImagefile (ImageFile& im, const int in_nsample, const int trace) const +{ + convertToImagefile (im, in_nsample, trace, 0, im.nx()); +} + void -Phantom::convertToImagefile (ImageFile& im, const int colStart, const int colCount, const int in_nsample, const int trace) const +Phantom::convertToImagefile (ImageFile& im, const int in_nsample, const int trace, const int colStart, const int colCount) const { int nx = im.nx(); int ny = im.ny(); if (nx < 2 || ny < 2) return; + int nsample = in_nsample; if (nsample < 1) - nsample = 1; + nsample = 1; - int nsample = in_nsample; double dx = m_xmax - m_xmin; double dy = m_ymax - m_ymin; double xcent = m_xmin + dx / 2; @@ -369,7 +452,7 @@ Phantom::convertToImagefile (ImageFile& im, const int colStart, const int colCou ImageFileArray v = im.getArray(); double x_start = xmin + (colStart * xinc); - for (PElemConstIterator pelem = m_listPElem.begin; pelem != m_listPElem.end; pelem++) { + for (PElemConstIterator pelem = m_listPElem.begin(); pelem != m_listPElem.end(); pelem++) { double x, y, xi, yi; int ix, iy, kx, ky; for (ix = 0, x = x_start; ix < colCount; ix++, x += xinc) { @@ -408,20 +491,7 @@ PhantomElement::PhantomElement (const char *type, const double cx, const double { m_rot = convertDegreesToRadians (rot); // convert angle to radians - if (strcasecmp (type, "rectangle") == 0) - m_type = PELEM_RECTANGLE; - else if (strcasecmp (type, "triangle") == 0) - m_type = PELEM_TRIANGLE; - else if (strcasecmp (type, "ellipse") == 0) - m_type = PELEM_ELLIPSE; - else if (strcasecmp (type, "sector") == 0) - m_type = PELEM_SECTOR; - else if (strcasecmp (type, "segment") == 0) - m_type = PELEM_SEGMENT; - else { - sys_error (ERR_WARNING, "Unknown PhantomElement type %s [PhantomElement::PhantomElement]", type); - m_type = PELEM_INVALID; - } + m_type = convertNameToType (type); makeTransformMatrices (); // calc transform matrices between phantom and normalized phantomelement makeVectorOutline (); // calculate vector outline of pelem @@ -446,6 +516,27 @@ PhantomElement::~PhantomElement (void) delete m_yOutline; } +PhmElemType +PhantomElement::convertNameToType (const char* const typeName) +{ + PhmElemType type = PELEM_INVALID; + + if (strcasecmp (typeName, "rectangle") == 0) + type = PELEM_RECTANGLE; + else if (strcasecmp (typeName, "triangle") == 0) + type = PELEM_TRIANGLE; + else if (strcasecmp (typeName, "ellipse") == 0) + type = PELEM_ELLIPSE; + else if (strcasecmp (typeName, "sector") == 0) + type = PELEM_SECTOR; + else if (strcasecmp (typeName, "segment") == 0) + type = PELEM_SEGMENT; + else + sys_error (ERR_WARNING, "Unknown PhantomElement type %s [PhantomElement::PhantomElement]", type); + + return (type); +} + void PhantomElement::makeTransformMatrices (void) { @@ -486,8 +577,6 @@ PhantomElement::makeTransformMatrices (void) * Called by phm_add_pelem() */ -static const double SCALE_PELEM_EXTENT=0.005; // increase pelem limits by 0.5% - void PhantomElement::makeVectorOutline (void) { @@ -578,7 +667,6 @@ PhantomElement::makeVectorOutline (void) } - /* NAME * calc_arc Calculate outline of a arc of a circle * @@ -642,10 +730,9 @@ PhantomElement::calcEllipsePoints (double x[], double y[], const int pts, const */ int -PhantomElement::numCirclePoints (double theta) const +PhantomElement::numCirclePoints (double theta) { - if (theta < 0.0 || theta > TWOPI) - sys_error(ERR_WARNING, "illegal values sent to circle_pts"); + theta = clamp (theta, 0., TWOPI); return static_cast (POINTS_PER_CIRCLE * theta / TWOPI + 1.5); }