X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fphantom.cpp;h=dba6576c8cf2b78761c38a47fca530bab4b8ffec;hp=1508bea745912e3ffc07661c057f52c2096db4b1;hb=3fba6928127cd65870bdcd96c8114ad5894247ae;hpb=99dd1d6ed10db1f669a5fe6af71225a50fc0ddfb diff --git a/libctsim/phantom.cpp b/libctsim/phantom.cpp index 1508bea..dba6576 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.1 2000/06/19 02:59:34 kevin Exp $ +** $Id: phantom.cpp,v 1.6 2000/07/13 07:03:21 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; } @@ -249,10 +327,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,7 +336,7 @@ 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); @@ -277,17 +352,14 @@ Phantom::std_rowland (void) } 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,7 +367,7 @@ 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); @@ -314,6 +386,96 @@ Phantom::std_herman (void) addPElem("ellipse", 0.000, 0.00, 7.875, 5.7187, 90.00, -0.206); } +void +Phantom::addStdHermanBordered (void) +{ + addPElem("ellipse", 0., 0., 6.6, 5.9, 90., 0.); +} + + +/* NAME + * convertToImagefile Make image array from Phantom + * + * SYNOPSIS + * pic_to_imagefile (pic, im, nsample) + * Phantom& pic Phantom definitions + * ImageFile *im Computed pixel array + * int nsample Number of samples along each axis for each pixel + * (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 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; + + double dx = m_xmax - m_xmin; + double dy = m_ymax - m_ymin; + double xcent = m_xmin + dx / 2; + double ycent = m_ymin + dy / 2; + double phmlen = (dx > dy ? dx : dy); + + double phmradius = phmlen / 2; + + double xmin = xcent - phmradius; + double xmax = xcent + phmradius; + double ymin = ycent - phmradius; + double ymax = ycent + phmradius; + + // Each pixel holds the average of the intensity of the cell with (ix,iy) at the center of the pixel + // Set major increments so that the last cell v[nx-1][ny-1] will start at xmax - xinc, ymax - yinc). + // Set minor increments so that sample points are centered in cell + + double xinc = (xmax - xmin) / nx; + double yinc = (ymax - ymin) / ny; + + double kxinc = xinc / nsample; /* interval between samples */ + double kyinc = yinc / nsample; + double kxofs = kxinc / 2; /* offset of 1st point */ + double kyofs = kyinc / 2; + + im.setAxisExtent (xmin, xmax, ymin, ymax); + im.setAxisIncrement (xinc, yinc); + + ImageFileArray v = im.getArray(); + + double x_start = xmin + (colStart * xinc); + 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) { + 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 ((*pelem)->isPointInside (xi, yi, PHM_COORD) == TRUE) + v[ix][iy] += (*pelem)->atten(); + } // for kx + } /* for iy */ + } /* for ix */ + } /* for pelem */ + + + 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; + } +} //////////////////////////////////////////////////////////////////////////////////////////////////////// // CLASS IDENTIFICATION @@ -330,20 +492,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 @@ -368,6 +517,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) { @@ -408,8 +578,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) { @@ -500,7 +668,6 @@ PhantomElement::makeVectorOutline (void) } - /* NAME * calc_arc Calculate outline of a arc of a circle * @@ -564,10 +731,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); }