From 23f5654dacb1952c15bda92c2606fae3a55e48ad Mon Sep 17 00:00:00 2001 From: "Kevin M. Rosenberg" Date: Thu, 4 Jan 2001 21:28:41 +0000 Subject: [PATCH] r354: Added Projection Polar conversions --- ChangeLog | 8 + include/imagefile.h | 4 +- include/projections.h | 26 +- include/sgp.h | 5 +- libctgraphics/ezplot.cpp | 20 +- libctgraphics/sgp.cpp | 42 ++- libctsim/backprojectors.cpp | 410 +++++++++++++++--------------- libctsim/imagefile.cpp | 20 +- libctsim/projections.cpp | 494 ++++++++++++++++++++++++++---------- libctsim/scanner.cpp | 394 ++++++++++++++-------------- msvc/ctsim/ctsim.plg | 47 ++-- src/ctsim.h | 6 +- src/dialogs.cpp | 203 +++++++++++---- src/dialogs.h | 29 ++- src/dlgprojections.cpp | 10 +- src/views.cpp | 108 ++++++-- src/views.h | 14 +- 17 files changed, 1178 insertions(+), 662 deletions(-) diff --git a/ChangeLog b/ChangeLog index a64a063..5e0220b 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,6 +1,14 @@ 3.0.0beta1 - Released * ctsim: Fixed initialization of min/max for PlotFiles + + * ctsim: Added reset to full-intensity scale menu item + + * ctsim: Add convert projections to polar plot + + * ezplot: Cleaned up y-tick label placement + + * sgp: Added better support for projection/reconstruction animation 3.0.0alpha3 - Released 1/02/01 diff --git a/include/imagefile.h b/include/imagefile.h index eaa69af..aabb903 100644 --- a/include/imagefile.h +++ b/include/imagefile.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: imagefile.h,v 1.28 2001/01/02 16:02:12 kevin Exp $ +** $Id: imagefile.h,v 1.29 2001/01/04 21:28:41 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 @@ -151,6 +151,8 @@ public: : ImageFileBase () {} + void getCenterCoordinates (unsigned int& iXCenter, unsigned int& iYCenter); + bool convertRealToComplex (); bool convertComplexToReal (); diff --git a/include/projections.h b/include/projections.h index 338f376..4c93304 100644 --- a/include/projections.h +++ b/include/projections.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: projections.h,v 1.19 2001/01/03 22:00:46 kevin Exp $ +** $Id: projections.h,v 1.20 2001/01/04 21:28:41 kevin Exp $ ** ** ** This program is free software; you can redistribute it and/or modify @@ -32,17 +32,34 @@ class Scanner; class DetectorArray; class Array2dFileLabel; +class fnetorderstream; +#include "array2dfile.h" +#include "imagefile.h" +#include // Projections class Projections { public: + + static const int POLAR_INTERP_INVALID; + static const int POLAR_INTERP_NEAREST; + static const int POLAR_INTERP_BILINEAR; + static const int POLAR_INTERP_BICUBIC; + Projections (const int nView, const int nDet); Projections (const Scanner& scanner); Projections (); ~Projections (); + static const int getInterpCount() {return s_iInterpCount;} + static const char** getInterpNameArray() {return s_aszInterpName;} + static const char** getInterpTitleArray() {return s_aszInterpTitle;} + static int convertInterpNameToID (const char* const interpName); + static const char* convertInterpIDToName (const int interpID); + static const char* convertInterpIDToTitle (const int interpID); + void initFromScanner (const Scanner& scanner); void printProjectionData (int startView, int endView); @@ -58,6 +75,9 @@ class Projections bool convertPolar (ImageFile& rIF, int iInterpolation); bool convertFFTPolar (ImageFile& rIF, int iInterpolation, int iZeropad); + void calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double** ppdView, double** ppdDet); + void interpolatePolar (ImageFileArray& v, ImageFileArray& vImag, unsigned int nx, unsigned int ny, std::complex** ppcDetValue, + double** ppdDet, double** ppdView, unsigned int nView, unsigned int nDet, int iInterpolate); bool reconstruct (ImageFile& im, const char* const filterName, double filt_param, const char* const filterMethodName, const int zeropad, const char* frequencyFilterName, const char* const interpName, int interp_param, const char* const backprojName, const int trace) const; @@ -124,6 +144,10 @@ class Projections const static kuint16 m_signature; + static const char* s_aszInterpName[]; + static const char* s_aszInterpTitle[]; + static const int s_iInterpCount; + bool headerWrite (fnetorderstream& fs); bool headerRead (fnetorderstream& fs); void newProjData (); diff --git a/include/sgp.h b/include/sgp.h index 2d83b2e..611654e 100644 --- a/include/sgp.h +++ b/include/sgp.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: sgp.h,v 1.24 2001/01/02 16:02:13 kevin Exp $ +** $Id: sgp.h,v 1.25 2001/01/04 21:28:41 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 @@ -125,6 +125,7 @@ private: double m_dCurrentWorldX; double m_dCurrentWorldY; double m_dTextAngle; + int m_iTextPointSize; bool m_bRecalcTransform; double m_dPointsPerPixel; // points (72pt/in) per screen pixel; int m_iLinestyle; @@ -145,6 +146,8 @@ private: #if HAVE_WXWINDOWS wxPen m_pen; wxFont* m_pFont; + + void initFromDC (wxDC* pDC); #endif public: diff --git a/libctgraphics/ezplot.cpp b/libctgraphics/ezplot.cpp index 260d5ba..ea29c76 100644 --- a/libctgraphics/ezplot.cpp +++ b/libctgraphics/ezplot.cpp @@ -6,7 +6,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: ezplot.cpp,v 1.27 2001/01/02 09:58:11 kevin Exp $ +** $Id: ezplot.cpp,v 1.28 2001/01/04 21:28:41 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 @@ -723,18 +723,20 @@ EZPlot::plot (SGP* pSGP) if (o_yaxis == NOAXIS || o_ytlabel == FALSE) ytl_ofs = 0.0; else if (o_yticks == LEFT) { - double xExtent, yExtent; + double xExtentMin, xExtentMax, yExtent; char s[1024]; - snprintf (s, sizeof(s), y_numfmt, 0); - m_pSGP->getTextExtent (s, &xExtent, &yExtent); - ytl_ofs = -2.0 * charwidth - xExtent; + snprintf (s, sizeof(s), y_numfmt, ymin); + m_pSGP->getTextExtent (s, &xExtentMin, &yExtent); + snprintf (s, sizeof(s), y_numfmt, ymax); + m_pSGP->getTextExtent (s, &xExtentMax, &yExtent); + if (xExtentMin > xExtentMax) + xExtentMax = xExtentMin; + ytl_ofs = -xExtentMax; } else if (o_yticks == RIGHT) ytl_ofs = 1.5 * charwidth; xa_max -= 0.7 * x_fldwid * charwidth; // make room for last x tick label - - xt_min = xa_min; yt_min = ya_min; xt_max = xa_max; @@ -754,11 +756,11 @@ EZPlot::plot (SGP* pSGP) if (ytl_ofs != 0.0 && s_xcross == FALSE) { if (o_yticks == LEFT) { - xa_min += (2 + y_fldwid) * charwidth; + xa_min += 2*charwidth - ytl_ofs; // (2 + y_fldwid) * charwidth; xt_min = xa_min; } else if (o_yticks == RIGHT) { xa_min += 0.0; - xt_min = xa_min + ytl_ofs + y_fldwid * charwidth; + xt_min = xa_min + ytl_ofs; // + y_fldwid * charwidth; } } else xt_min = xa_min; diff --git a/libctgraphics/sgp.cpp b/libctgraphics/sgp.cpp index ec129d2..263101f 100644 --- a/libctgraphics/sgp.cpp +++ b/libctgraphics/sgp.cpp @@ -7,7 +7,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: sgp.cpp,v 1.25 2001/01/02 16:02:13 kevin Exp $ +** $Id: sgp.cpp,v 1.26 2001/01/04 21:28:41 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 @@ -91,7 +91,26 @@ SGP::SGP (const SGPDriver& driver) ndc_to_mc.setIdentity(); m_ctm.setIdentity(); -#if HAVE_WXWINDOWS +#ifdef HAVE_WXWINDOWS + initFromDC (driver.idWX()); +#endif + + setWindow (0., 0., 1., 1.); + setViewport (0., 0., 1., 1.); + moveAbs (0., 0.); + stylusNDC (0., 0., false); + + setTextAngle (0.); + setTextPointSize (12); + setColor (C_BLACK); + setLineStyle (LS_SOLID); +} + + +#ifdef HAVE_WXWINDOWS +void +SGP::initFromDC (wxDC* pDC) +{ m_pen.SetWidth (1); if (m_driver.isWX()) { @@ -112,18 +131,9 @@ SGP::SGP (const SGPDriver& driver) m_dPointsPerPixel = iTestPointSize / dTestCharHeight; m_driver.idWX()->SetBackground (*wxWHITE_BRUSH); } +} #endif - setWindow (0., 0., 1., 1.); - setViewport (0., 0., 1., 1.); - moveAbs (0., 0.); - stylusNDC (0., 0., false); - - setTextAngle (0.); - setTextPointSize (12); - setColor (C_BLACK); - setLineStyle (LS_SOLID); -} SGP::~SGP() { @@ -519,7 +529,8 @@ SGP::setTextPointSize (double height) #endif #if HAVE_WXWINDOWS if (m_driver.isWX()) { - m_pFont->SetPointSize (static_cast(height+0.5)); + m_iTextPointSize = static_cast(height+0.5); + m_pFont->SetPointSize (m_iTextPointSize); m_driver.idWX()->SetFont (*m_pFont); } #endif @@ -941,7 +952,10 @@ const unsigned char SGP::MARKER_BITMAP[MARK_COUNT][5] = void SGP::setDC (wxDC* pDC) { - if (m_driver.isWX()) + if (m_driver.isWX()) { m_driver.setDC(pDC); + initFromDC (pDC); + setTextPointSize (m_iTextPointSize); + } } #endif diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp index 1f8a1d2..8343314 100644 --- a/libctsim/backprojectors.cpp +++ b/libctsim/backprojectors.cpp @@ -8,7 +8,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: backprojectors.cpp,v 1.22 2001/01/02 16:02:13 kevin Exp $ +** $Id: backprojectors.cpp,v 1.23 2001/01/04 21:28:41 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 @@ -101,7 +101,7 @@ Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char { m_fail = false; m_pBackprojectImplem = NULL; - + initBackprojector (proj, im, backprojName, interpName, interpFactor); } @@ -142,35 +142,35 @@ Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const m_failMessage = "Invalid interpolation name "; m_failMessage += interpName; } - + if (m_fail || m_idBackproject == BPROJ_INVALID || m_idInterpolation == INTERP_INVALID) { m_fail = true; return false; } - + if (proj.geometry() == Scanner::GEOMETRY_EQUILINEAR) - m_pBackprojectImplem = static_cast(new BackprojectEquilinear(proj, im, m_idInterpolation, interpFactor)); + m_pBackprojectImplem = static_cast(new BackprojectEquilinear(proj, im, m_idInterpolation, interpFactor)); else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR) - m_pBackprojectImplem = static_cast(new BackprojectEquiangular(proj, im, m_idInterpolation, interpFactor)); + m_pBackprojectImplem = static_cast(new BackprojectEquiangular(proj, im, m_idInterpolation, interpFactor)); else if (proj.geometry() == Scanner::GEOMETRY_PARALLEL) { - if (m_idBackproject == BPROJ_TRIG) - m_pBackprojectImplem = static_cast(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor)); - else if (m_idBackproject == BPROJ_TABLE) - m_pBackprojectImplem = static_cast(new BackprojectTable (proj, im, m_idInterpolation, interpFactor)); - else if (m_idBackproject == BPROJ_DIFF) - m_pBackprojectImplem = static_cast(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor)); - else if (m_idBackproject == BPROJ_DIFF2) - m_pBackprojectImplem = static_cast(new BackprojectDiff2 (proj, im, m_idInterpolation, interpFactor)); - else if (m_idBackproject == BPROJ_IDIFF2) - m_pBackprojectImplem = static_cast(new BackprojectIntDiff2 (proj, im, m_idInterpolation, interpFactor)); - else if (m_idBackproject == BPROJ_IDIFF3) - m_pBackprojectImplem = static_cast(new BackprojectIntDiff3 (proj, im, m_idInterpolation, interpFactor)); + if (m_idBackproject == BPROJ_TRIG) + m_pBackprojectImplem = static_cast(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor)); + else if (m_idBackproject == BPROJ_TABLE) + m_pBackprojectImplem = static_cast(new BackprojectTable (proj, im, m_idInterpolation, interpFactor)); + else if (m_idBackproject == BPROJ_DIFF) + m_pBackprojectImplem = static_cast(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor)); + else if (m_idBackproject == BPROJ_DIFF2) + m_pBackprojectImplem = static_cast(new BackprojectDiff2 (proj, im, m_idInterpolation, interpFactor)); + else if (m_idBackproject == BPROJ_IDIFF2) + m_pBackprojectImplem = static_cast(new BackprojectIntDiff2 (proj, im, m_idInterpolation, interpFactor)); + else if (m_idBackproject == BPROJ_IDIFF3) + m_pBackprojectImplem = static_cast(new BackprojectIntDiff3 (proj, im, m_idInterpolation, interpFactor)); } else { - m_fail = true; - m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]"; - return false; + m_fail = true; + m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]"; + return false; } - + return true; } @@ -179,24 +179,24 @@ int Backprojector::convertBackprojectNameToID (const char* const backprojName) { int backprojID = BPROJ_INVALID; - + for (int i = 0; i < s_iBackprojectCount; i++) - if (strcasecmp (backprojName, s_aszBackprojectName[i]) == 0) { - backprojID = i; - break; - } - - return (backprojID); + if (strcasecmp (backprojName, s_aszBackprojectName[i]) == 0) { + backprojID = i; + break; + } + + return (backprojID); } const char* Backprojector::convertBackprojectIDToName (int bprojID) { static const char *bprojName = ""; - + if (bprojID >= 0 && bprojID < s_iBackprojectCount) - return (s_aszBackprojectName[bprojID]); - + return (s_aszBackprojectName[bprojID]); + return (bprojName); } @@ -204,10 +204,10 @@ const char* Backprojector::convertBackprojectIDToTitle (const int bprojID) { static const char *bprojTitle = ""; - + if (bprojID >= 0 && bprojID < s_iBackprojectCount) - return (s_aszBackprojectTitle[bprojID]); - + return (s_aszBackprojectTitle[bprojID]); + return (bprojTitle); } @@ -216,24 +216,24 @@ int Backprojector::convertInterpNameToID (const char* const interpName) { int interpID = INTERP_INVALID; - + for (int i = 0; i < s_iInterpCount; i++) - if (strcasecmp (interpName, s_aszInterpName[i]) == 0) { - interpID = i; - break; - } - - return (interpID); + if (strcasecmp (interpName, s_aszInterpName[i]) == 0) { + interpID = i; + break; + } + + return (interpID); } const char* Backprojector::convertInterpIDToName (const int interpID) { static const char *interpName = ""; - + if (interpID >= 0 && interpID < s_iInterpCount) - return (s_aszInterpName[interpID]); - + return (s_aszInterpName[interpID]); + return (interpName); } @@ -241,10 +241,10 @@ const char* Backprojector::convertInterpIDToTitle (const int interpID) { static const char *interpTitle = ""; - + if (interpID >= 0 && interpID < s_iInterpCount) - return (s_aszInterpTitle[interpID]); - + return (s_aszInterpTitle[interpID]); + return (interpTitle); } @@ -257,33 +257,33 @@ Backprojector::convertInterpIDToTitle (const int interpID) // Pure virtual base class for all backprojectors. Backproject::Backproject (const Projections& proj, ImageFile& im, int interpType, const int interpFactor) - : proj(proj), im(im), interpType(interpType), m_interpFactor(interpFactor) +: proj(proj), im(im), interpType(interpType), m_interpFactor(interpFactor) { detInc = proj.detInc(); nDet = proj.nDet(); iDetCenter = (nDet - 1) / 2; // index refering to L=0 projection rotScale = proj.rotInc(); - + if (proj.geometry() == Scanner::GEOMETRY_PARALLEL) - rotScale /= (proj.nView() * proj.rotInc() / PI); // scale by number of PI rotations + rotScale /= (proj.nView() * proj.rotInc() / PI); // scale by number of PI rotations else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR || proj.geometry() == Scanner::GEOMETRY_EQUILINEAR) - rotScale /= (proj.nView() * proj.rotInc() / (2 * PI)); // scale by number of 2PI rotations + rotScale /= (proj.nView() * proj.rotInc() / (2 * PI)); // scale by number of 2PI rotations else - sys_error (ERR_SEVERE, "Invalid geometry type %d [Backproject::Backproject]", proj.geometry()); - + sys_error (ERR_SEVERE, "Invalid geometry type %d [Backproject::Backproject]", proj.geometry()); + v = im.getArray(); nx = im.nx(); ny = im.ny(); im.arrayDataClear(); - + xMin = -proj.phmLen() / 2; // Retangular coords of phantom xMax = xMin + proj.phmLen(); yMin = -proj.phmLen() / 2; yMax = yMin + proj.phmLen(); - + xInc = (xMax - xMin) / nx; // size of cells yInc = (yMax - yMin) / ny; - + m_dFocalLength = proj.focalLength(); } @@ -300,8 +300,8 @@ Backproject::ScaleImageByRotIncrement () void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double r, double phi, double L, int iDetPos) { - sys_error (ERR_WARNING, "r=%f, phi=%f", r, phi); - errorIndexOutsideDetector (ix, iy, theta, L, iDetPos); + sys_error (ERR_WARNING, "r=%f, phi=%f", r, phi); + errorIndexOutsideDetector (ix, iy, theta, L, iDetPos); } void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double L, int iDetPos) @@ -313,7 +313,7 @@ void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, doubl os << "xMin=" << xMin << ", xMax=" << xMax << ", xInc=" << xInc << "\n"; os << "yMin=" << yMin << ", yMax=" << yMax << ", yInc=" << yInc << "\n"; os << "iDetPos index outside bounds: " << iDetPos << " [backprojector]";; - + sys_error (ERR_WARNING, os.str().c_str()); #endif } @@ -329,7 +329,7 @@ void BackprojectTrig::BackprojectView (const double* const filteredProj, const double view_angle) { double theta = view_angle; - + double x = xMin + xInc / 2; // Rectang coords of center of pixel for (int ix = 0; ix < nx; x += xInc, ix++) { double y = yMin + yInc / 2; @@ -337,23 +337,23 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double double r = sqrt (x * x + y * y); // distance of cell from center double phi = atan2 (y, x); // angle of cell from center double L = r * cos (theta - phi); // position on detector - + if (interpType == Backprojector::INTERP_NEAREST) { - int iDetPos = iDetCenter + nearest (L / detInc); // calc'd index in the filter raysum array - - if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos - errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos); - else - v[ix][iy] += rotScale * filteredProj[iDetPos]; + int iDetPos = iDetCenter + nearest (L / detInc); // calc'd index in the filter raysum array + + if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos + errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos); + else + v[ix][iy] += rotScale * filteredProj[iDetPos]; } else if (interpType == Backprojector::INTERP_LINEAR) { - double p = L / detInc; // position along detector - double pFloor = floor (p); - int iDetPos = iDetCenter + static_cast(pFloor); - double frac = p - pFloor; // fraction distance from det - if (iDetPos < 0 || iDetPos >= nDet - 1) // check for impossible: index outside of raysum pos - errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos); - else - v[ix][iy] += rotScale * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); + double p = L / detInc; // position along detector + double pFloor = floor (p); + int iDetPos = iDetCenter + static_cast(pFloor); + double frac = p - pFloor; // fraction distance from det + if (iDetPos < 0 || iDetPos >= nDet - 1) // check for impossible: index outside of raysum pos + errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos); + else + v[ix][iy] += rotScale * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); } } } @@ -367,13 +367,13 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double // Precalculates trigometric function value for each point in image for backprojection. BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, int interpType, const int interpFactor) - : Backproject (proj, im, interpType, interpFactor) +: Backproject (proj, im, interpType, interpFactor) { arrayR.initSetSize (im.nx(), im.ny()); arrayPhi.initSetSize (im.nx(), im.ny()); r = arrayR.getArray(); phi = arrayPhi.getArray(); - + double x, y; // Rectang coords of center of pixel int ix, iy; for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++) @@ -392,29 +392,29 @@ void BackprojectTable::BackprojectView (const double* const filteredProj, const double view_angle) { double theta = view_angle; - + for (int ix = 0; ix < nx; ix++) { ImageFileColumn pImCol = v[ix]; - + for (int iy = 0; iy < ny; iy++) { double L = r[ix][iy] * cos (theta - phi[ix][iy]); - + if (interpType == Backprojector::INTERP_NEAREST) { - int iDetPos = iDetCenter + nearest(L / detInc); // calc index in the filtered raysum vector - - if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos - errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos); - else - pImCol[iy] += filteredProj[iDetPos]; + int iDetPos = iDetCenter + nearest(L / detInc); // calc index in the filtered raysum vector + + if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos + errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos); + else + pImCol[iy] += filteredProj[iDetPos]; } else if (interpType == Backprojector::INTERP_LINEAR) { - double dPos = L / detInc; // position along detector - double dPosFloor = floor (dPos); - int iDetPos = iDetCenter + static_cast(dPosFloor); - double frac = dPos - dPosFloor; // fraction distance from det - if (iDetPos < 0 || iDetPos >= nDet - 1) - errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos); - else - pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); + double dPos = L / detInc; // position along detector + double dPosFloor = floor (dPos); + int iDetPos = iDetCenter + static_cast(dPosFloor); + double frac = dPos - dPosFloor; // fraction distance from det + if (iDetPos < 0 || iDetPos >= nDet - 1) + errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos); + else + pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); } } // end for y } // end for x @@ -429,14 +429,14 @@ BackprojectTable::BackprojectView (const double* const filteredProj, const doubl // Iterates in x & y direction by adding difference in L position BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, int interpType, const int interpFactor) - : Backproject (proj, im, interpType, interpFactor) +: Backproject (proj, im, interpType, interpFactor) { // calculate center of first pixel v[0][0] double x = xMin + xInc / 2; double y = yMin + yInc / 2; start_r = sqrt (x * x + y * y); start_phi = atan2 (y, x); - + im.arrayDataClear(); } @@ -452,31 +452,31 @@ BackprojectDiff::BackprojectView (const double* const filteredProj, const double double det_dx = xInc * cos (theta); double det_dy = yInc * sin (theta); double lColStart = start_r * cos (theta - start_phi); // calculate L for first point in image - + for (int ix = 0; ix < nx; ix++, lColStart += det_dx) { double curDetPos = lColStart; ImageFileColumn pImCol = v[ix]; - + for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { #ifdef DEBUG printf ("[%2d,%2d]: %8.5f ", ix, iy, curDetPos); #endif if (interpType == Backprojector::INTERP_NEAREST) { - int iDetPos = iDetCenter + nearest(curDetPos / detInc); // calc index in the filtered raysum vector - - if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos - errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); - else - pImCol[iy] += filteredProj[iDetPos]; + int iDetPos = iDetCenter + nearest(curDetPos / detInc); // calc index in the filtered raysum vector + + if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos + errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); + else + pImCol[iy] += filteredProj[iDetPos]; } else if (interpType == Backprojector::INTERP_LINEAR) { - double detPos = curDetPos / detInc; // position along detector - double detPosFloor = floor (detPos); - int iDetPos = iDetCenter + static_cast(detPosFloor); - double frac = detPos - detPosFloor; // fraction distance from det - if (iDetPos < 0 || iDetPos >= nDet - 1) - errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); - else - pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); + double detPos = curDetPos / detInc; // position along detector + double detPosFloor = floor (detPos); + int iDetPos = iDetCenter + static_cast(detPosFloor); + double frac = detPos - detPosFloor; // fraction distance from det + if (iDetPos < 0 || iDetPos >= nDet - 1) + errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); + else + pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); } } // end for y } // end for x @@ -493,40 +493,40 @@ void BackprojectDiff2::BackprojectView (const double* const filteredProj, const double view_angle) { double theta = view_angle; - + // Distance betw. detectors for an angle given in units of detectors double det_dx = xInc * cos (theta) / detInc; double det_dy = yInc * sin (theta) / detInc; - + // calculate detPosition for first point in image (ix=0, iy=0) double detPosColStart = start_r * cos (theta - start_phi) / detInc; - + #ifdef DEBUG printf ("start_r=%8.5f, start_phi=%8.5f, rotScale=%8.5f\n", start_r, start_phi, rotScale); #endif for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) { double curDetPos = detPosColStart; ImageFileColumn pImCol = v[ix]; - + for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { #ifdef DEBUG printf ("[%2d,%2d]: %8.5f %8.5f\n", ix, iy, curDetPos, filteredProj[iDetCenter + nearest(curDetPos)]); #endif if (interpType == Backprojector::INTERP_NEAREST) { - int iDetPos = iDetCenter + nearest (curDetPos); // calc index in the filtered raysum vector - - if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos - errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); - else - *pImCol++ += filteredProj[iDetPos]; + int iDetPos = iDetCenter + nearest (curDetPos); // calc index in the filtered raysum vector + + if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos + errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); + else + *pImCol++ += filteredProj[iDetPos]; } else if (interpType == Backprojector::INTERP_LINEAR) { - double detPosFloor = floor (curDetPos); - int iDetPos = iDetCenter + static_cast(detPosFloor); - double frac = curDetPos - detPosFloor; // fraction distance from det - if (iDetPos < 0 || iDetPos >= nDet - 1) - errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); - else - *pImCol++ += filteredProj[iDetPos] + (frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])); + double detPosFloor = floor (curDetPos); + int iDetPos = iDetCenter + static_cast(detPosFloor); + double frac = curDetPos - detPosFloor; // fraction distance from det + if (iDetPos < 0 || iDetPos >= nDet - 1) + errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); + else + *pImCol++ += filteredProj[iDetPos] + (frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])); } } // end for y } // end for x @@ -542,43 +542,43 @@ void BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const double view_angle) { double theta = view_angle; - + static const kint32 scale = 1 << 16; static const double dScale = scale; static const kint32 halfScale = scale / 2; - + const kint32 det_dx = nearest (xInc * cos (theta) / detInc * scale); const kint32 det_dy = nearest (yInc * sin (theta) / detInc * scale); - + // calculate L for first point in image (0, 0) kint32 detPosColStart = nearest (start_r * cos (theta - start_phi) / detInc * scale); - + for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) { kint32 curDetPos = detPosColStart; ImageFileColumn pImCol = v[ix]; - + for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { if (interpType == Backprojector::INTERP_NEAREST) { - int detPosNearest = (curDetPos >= 0 ? ((curDetPos + halfScale) / scale) : ((curDetPos - halfScale) / scale)); - int iDetPos = iDetCenter + detPosNearest; // calc index in the filtered raysum vector - - if (iDetPos < 0 || iDetPos >= nDet) // check for index outside of raysum pos - errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); - else - *pImCol++ += filteredProj[iDetPos]; + int detPosNearest = (curDetPos >= 0 ? ((curDetPos + halfScale) / scale) : ((curDetPos - halfScale) / scale)); + int iDetPos = iDetCenter + detPosNearest; // calc index in the filtered raysum vector + + if (iDetPos < 0 || iDetPos >= nDet) // check for index outside of raysum pos + errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); + else + *pImCol++ += filteredProj[iDetPos]; } else if (interpType == Backprojector::INTERP_LINEAR) { - kint32 detPosFloor = curDetPos / scale; - kint32 detPosRemainder = curDetPos % scale; - if (detPosRemainder < 0) { - detPosFloor--; - detPosRemainder += scale; - } - int iDetPos = iDetCenter + detPosFloor; - double frac = detPosRemainder / dScale; - if (iDetPos < 0 || iDetPos >= nDet - 1) - errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); - else - *pImCol++ += ((1.-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); + kint32 detPosFloor = curDetPos / scale; + kint32 detPosRemainder = curDetPos % scale; + if (detPosRemainder < 0) { + detPosFloor--; + detPosRemainder += scale; + } + int iDetPos = iDetCenter + detPosFloor; + double frac = detPosRemainder / dScale; + if (iDetPos < 0 || iDetPos >= nDet - 1) + errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); + else + *pImCol++ += ((1.-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); } } // end for y } // end for x @@ -599,13 +599,13 @@ BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const do static const kint32 scaleBitmask = scale - 1; static const kint32 halfScale = scale / 2; static const double dInvScale = 1. / scale; - + const kint32 det_dx = nearest (xInc * cos (theta) / detInc * scale); const kint32 det_dy = nearest (yInc * sin (theta) / detInc * scale); - + // calculate L for first point in image (0, 0) kint32 detPosColStart = nearest ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale); - + // precalculate scaled difference for linear interpolation double* deltaFilteredProj = new double [nDet]; if (interpType == Backprojector::INTERP_LINEAR) { @@ -613,34 +613,34 @@ BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const do deltaFilteredProj[i] = (filteredProj[i+1] - filteredProj[i]) * dInvScale; } deltaFilteredProj[nDet - 1] = 0; // last detector - + int iLastDet = nDet - 1; for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) { kint32 curDetPos = detPosColStart; ImageFileColumn pImCol = v[ix]; - + if (interpType == Backprojector::INTERP_NEAREST) { for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { - const int iDetPos = (curDetPos + halfScale) >> 16; - if (iDetPos >= 0 && iDetPos <= iLastDet) - *pImCol++ += filteredProj[iDetPos]; + const int iDetPos = (curDetPos + halfScale) >> 16; + if (iDetPos >= 0 && iDetPos <= iLastDet) + *pImCol++ += filteredProj[iDetPos]; } // end for iy } else if (interpType == Backprojector::INTERP_FREQ_PREINTERPOLATION) { for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { - const int iDetPos = ((curDetPos + halfScale) >> 16) * m_interpFactor; - if (iDetPos >= 0 && iDetPos <= iLastDet) - *pImCol++ += filteredProj[iDetPos]; + const int iDetPos = ((curDetPos + halfScale) >> 16) * m_interpFactor; + if (iDetPos >= 0 && iDetPos <= iLastDet) + *pImCol++ += filteredProj[iDetPos]; } // end for iy } else if (interpType == Backprojector::INTERP_LINEAR) { for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { - const kint32 iDetPos = curDetPos >> scaleShift; - const kint32 detRemainder = curDetPos & scaleBitmask; - if (iDetPos >= 0 && iDetPos <= iLastDet) - *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]); + const kint32 iDetPos = curDetPos >> scaleShift; + const kint32 detRemainder = curDetPos & scaleBitmask; + if (iDetPos >= 0 && iDetPos <= iLastDet) + *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]); } // end for iy } //end linear } // end for ix - + delete deltaFilteredProj; } @@ -649,10 +649,10 @@ void BackprojectEquiangular::BackprojectView (const double* const filteredProj, const double view_angle) { double beta = view_angle; - + for (int ix = 0; ix < nx; ix++) { ImageFileColumn pImCol = v[ix]; - + for (int iy = 0; iy < ny; iy++) { double dAngleDiff = beta - phi[ix][iy]; double rcos_t = r[ix][iy] * cos (dAngleDiff); @@ -660,23 +660,23 @@ BackprojectEquiangular::BackprojectView (const double* const filteredProj, const double dFLPlusSin = m_dFocalLength + rsin_t; double gamma = atan (rcos_t / dFLPlusSin); double dL2 = dFLPlusSin * dFLPlusSin + (rcos_t * rcos_t); - + if (interpType == Backprojector::INTERP_NEAREST) { - int iDetPos =iDetCenter + nearest(gamma / detInc); // calc index in the filtered raysum vector - - if (iDetPos < 0 || iDetPos >= nDet) { // check for impossible: index outside of raysum pos - ; // errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], gamma, iDetPos); - } else - pImCol[iy] += filteredProj[iDetPos] / dL2; + int iDetPos =iDetCenter + nearest(gamma / detInc); // calc index in the filtered raysum vector + + if (iDetPos < 0 || iDetPos >= nDet) { // check for impossible: index outside of raysum pos + ; // errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], gamma, iDetPos); + } else + pImCol[iy] += filteredProj[iDetPos] / dL2; } else if (interpType == Backprojector::INTERP_LINEAR) { - double dPos = gamma / detInc; // position along detector - double dPosFloor = floor (dPos); - int iDetPos = iDetCenter + static_cast(dPosFloor); - double frac = dPos - dPosFloor; // fraction distance from det - if (iDetPos < 0 || iDetPos >= nDet - 1) { - ; // errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], gamma, iDetPos); - } else - pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / dL2; + double dPos = gamma / detInc; // position along detector + double dPosFloor = floor (dPos); + int iDetPos = iDetCenter + static_cast(dPosFloor); + double frac = dPos - dPosFloor; // fraction distance from det + if (iDetPos < 0 || iDetPos >= nDet - 1) { + ; // errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], gamma, iDetPos); + } else + pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / dL2; } } // end for y } // end for x @@ -686,37 +686,37 @@ void BackprojectEquilinear::BackprojectView (const double* const filteredProj, const double view_angle) { double beta = view_angle; - + for (int ix = 0; ix < nx; ix++) { ImageFileColumn pImCol = v[ix]; - + for (int iy = 0; iy < ny; iy++) { double dAngleDiff = beta - phi[ix][iy]; double rcos_t = r[ix][iy] * cos (dAngleDiff); double rsin_t = r[ix][iy] * sin (dAngleDiff); - + double dU = (m_dFocalLength + rsin_t) / m_dFocalLength; double dDetPos = rcos_t / dU; // double to scale for imaginary detector that passes through origin // of phantom, see Kak-Slaney Figure 3.22 dDetPos *= 2; - + if (interpType == Backprojector::INTERP_NEAREST) { - int iDetPos = iDetCenter + nearest(dDetPos / detInc); // calc index in the filtered raysum vector - - if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos - ; /// errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], dDetPos, iDetPos); - else - pImCol[iy] += (filteredProj[iDetPos] / (dU * dU)); + int iDetPos = iDetCenter + nearest(dDetPos / detInc); // calc index in the filtered raysum vector + + if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos + ; /// errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], dDetPos, iDetPos); + else + pImCol[iy] += (filteredProj[iDetPos] / (dU * dU)); } else if (interpType == Backprojector::INTERP_LINEAR) { - double dPos = dDetPos / detInc; // position along detector - double dPosFloor = floor (dPos); - int iDetPos = iDetCenter + static_cast(dPosFloor); - double frac = dPos - dPosFloor; // fraction distance from det - if (iDetPos < 0 || iDetPos >= nDet - 1) - ; // errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], dDetPos, iDetPos); - else - pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / (dU * dU); + double dPos = dDetPos / detInc; // position along detector + double dPosFloor = floor (dPos); + int iDetPos = iDetCenter + static_cast(dPosFloor); + double frac = dPos - dPosFloor; // fraction distance from det + if (iDetPos < 0 || iDetPos >= nDet - 1) + ; // errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], dDetPos, iDetPos); + else + pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / (dU * dU); } } // end for y } // end for x diff --git a/libctsim/imagefile.cpp b/libctsim/imagefile.cpp index df0e51e..946e4db 100644 --- a/libctsim/imagefile.cpp +++ b/libctsim/imagefile.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: imagefile.cpp,v 1.33 2001/01/02 16:02:13 kevin Exp $ +** $Id: imagefile.cpp,v 1.34 2001/01/04 21:28:41 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 @@ -84,21 +84,29 @@ F64Image::F64Image (void) } void -ImageFile::filterResponse (const char* const domainName, double bw, const char* const filterName, double filt_param, double dInputScale, double dOutputScale) +ImageFile::getCenterCoordinates (unsigned int& iXCenter, unsigned int& iYCenter) { - ImageFileArray v = getArray(); - SignalFilter filter (filterName, domainName, bw, filt_param); - - int iXCenter, iYCenter; if (isEven (m_nx)) iXCenter = m_nx / 2; else iXCenter = (m_nx - 1) / 2; + if (isEven (m_ny)) iYCenter = m_ny / 2; else iYCenter = (m_ny - 1) / 2; +} + + +void +ImageFile::filterResponse (const char* const domainName, double bw, const char* const filterName, double filt_param, double dInputScale, double dOutputScale) +{ + ImageFileArray v = getArray(); + SignalFilter filter (filterName, domainName, bw, filt_param); + unsigned int iXCenter, iYCenter; + getCenterCoordinates (iXCenter, iYCenter); + for (unsigned int ix = 0; ix < m_nx; ix++) for (unsigned int iy = 0; iy < m_ny; iy++) { long lD2 = ((ix - iXCenter) * (ix - iXCenter)) + ((iy - iYCenter) * (iy - iYCenter)); diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp index bdcf35e..bdf6073 100644 --- a/libctsim/projections.cpp +++ b/libctsim/projections.cpp @@ -8,7 +8,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: projections.cpp,v 1.36 2001/01/03 22:00:46 kevin Exp $ +** $Id: projections.cpp,v 1.37 2001/01/04 21:28:41 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 @@ -28,32 +28,54 @@ const kuint16 Projections::m_signature = ('P'*256 + 'J'); +const int Projections::POLAR_INTERP_INVALID = -1; +const int Projections::POLAR_INTERP_NEAREST = 0; +const int Projections::POLAR_INTERP_BILINEAR = 1; +const int Projections::POLAR_INTERP_BICUBIC = 2; + +const char* Projections::s_aszInterpName[] = +{ + {"nearest"}, + {"bilinear"}, + {"bicubic"}, +}; + +const char* Projections::s_aszInterpTitle[] = +{ + {"Nearest"}, + {"Bilinear"}, + {"Bicubic"}, +}; + +const int Projections::s_iInterpCount = sizeof(s_aszInterpName) / sizeof(char*); + + /* NAME - * Projections Constructor for projections matrix storage - * - * SYNOPSIS - * proj = projections_create (filename, nView, nDet) - * Projections& proj Allocated projections structure & matrix - * int nView Number of rotated view - * int nDet Number of detectors - * - */ +* Projections Constructor for projections matrix storage +* +* SYNOPSIS +* proj = projections_create (filename, nView, nDet) +* Projections& proj Allocated projections structure & matrix +* int nView Number of rotated view +* int nDet Number of detectors +* +*/ Projections::Projections (const Scanner& scanner) - : m_projData(0) +: m_projData(0) { initFromScanner (scanner); } Projections::Projections (const int nView, const int nDet) - : m_projData(0) +: m_projData(0) { init (nView, nDet); } Projections::Projections (void) - : m_projData(0) +: m_projData(0) { init (0, 0); } @@ -63,6 +85,43 @@ Projections::~Projections (void) deleteProjData(); } +int +Projections::convertInterpNameToID (const char* const interpName) +{ + int interpID = POLAR_INTERP_INVALID; + + for (int i = 0; i < s_iInterpCount; i++) + if (strcasecmp (interpName, s_aszInterpName[i]) == 0) { + interpID = i; + break; + } + + return (interpID); +} + +const char* +Projections::convertInterpIDToName (const int interpID) +{ + static const char *interpName = ""; + + if (interpID >= 0 && interpID < s_iInterpCount) + return (s_aszInterpName[interpID]); + + return (interpName); +} + +const char* +Projections::convertInterpIDToTitle (const int interpID) +{ + static const char *interpTitle = ""; + + if (interpID >= 0 && interpID < s_iInterpCount) + return (s_aszInterpTitle[interpID]); + + return (interpTitle); +} + + void Projections::init (const int nView, const int nDet) @@ -71,7 +130,7 @@ Projections::init (const int nView, const int nDet) m_nView = nView; m_nDet = nDet; newProjData (); - + time_t t = time (NULL); tm* lt = localtime (&t); m_year = lt->tm_year; @@ -88,7 +147,7 @@ Projections::initFromScanner (const Scanner& scanner) m_label.setLabelType (Array2dFileLabel::L_HISTORY); deleteProjData(); init (scanner.nView(), scanner.nDet()); - + m_phmLen = scanner.phmLen(); m_rotInc = scanner.rotInc(); m_detInc = scanner.detInc(); @@ -97,12 +156,6 @@ Projections::initFromScanner (const Scanner& scanner) m_fieldOfView = scanner.fieldOfView(); m_rotStart = 0; m_detStart = -(scanner.detLen() / 2); -#if 0 - if (m_geometry == Scanner::GEOMETRY_EQUILINEAR) { - m_detInc /= 2; - std::cout << "Kludge: detInc /= 2 in Projections::initFromScanner" << endl; - } -#endif } void @@ -120,10 +173,10 @@ Projections::newProjData (void) { if (m_projData) sys_error(ERR_WARNING, "m_projData != NULL [newProjData]"); - + if (m_nView > 0 && m_nDet) { m_projData = new DetectorArray* [m_nView]; - + for (int i = 0; i < m_nView; i++) m_projData[i] = new DetectorArray (m_nDet); } @@ -131,12 +184,12 @@ Projections::newProjData (void) /* NAME - * projections_free Free memory allocated to projections - * - * SYNOPSIS - * projections_free(proj) - * Projections& proj Projectionss to be deallocated - */ +* projections_free Free memory allocated to projections +* +* SYNOPSIS +* projections_free(proj) +* Projections& proj Projectionss to be deallocated +*/ void Projections::deleteProjData (void) @@ -144,7 +197,7 @@ Projections::deleteProjData (void) if (m_projData != NULL) { for (int i = 0; i < m_nView; i++) delete m_projData[i]; - + delete m_projData; m_projData = NULL; } @@ -152,9 +205,9 @@ Projections::deleteProjData (void) /* NAME - * Projections::headerWwrite Write data header for projections file - * - */ +* Projections::headerWwrite Write data header for projections file +* +*/ bool Projections::headerWrite (fnetorderstream& fs) @@ -171,7 +224,7 @@ Projections::headerWrite (fnetorderstream& fs) kuint16 _hour = m_hour; kuint16 _minute = m_minute; kuint16 _second = m_second; - + kfloat64 _calcTime = m_calcTime; kfloat64 _rotStart = m_rotStart; kfloat64 _rotInc = m_rotInc; @@ -180,11 +233,11 @@ Projections::headerWrite (fnetorderstream& fs) kfloat64 _phmLen = m_phmLen; kfloat64 _fieldOfView = m_fieldOfView; kfloat64 _focalLength = m_focalLength; - + fs.seekp(0); if (! fs) return false; - + fs.writeInt16 (_hsize); fs.writeInt16 (_signature); fs.writeInt32 (_nView); @@ -206,21 +259,21 @@ Projections::headerWrite (fnetorderstream& fs) fs.writeInt16 (_second); fs.writeInt16 (_remarksize); fs.write (m_remark.c_str(), _remarksize); - + m_headerSize = fs.tellp(); _hsize = m_headerSize; fs.seekp(0); fs.writeInt16 (_hsize); if (! fs) - return false; + return false; return true; } /* NAME - * projections_read_header Read data header for projections file - * - */ +* projections_read_header Read data header for projections file +* +*/ bool Projections::headerRead (fnetorderstream& fs) { @@ -230,8 +283,8 @@ Projections::headerRead (fnetorderstream& fs) fs.seekg(0); if (! fs) - return false; - + return false; + fs.readInt16 (_hsize); fs.readInt16 (_signature); fs.readInt32 (_nView); @@ -252,17 +305,17 @@ Projections::headerRead (fnetorderstream& fs) fs.readInt16 (_minute); fs.readInt16 (_second); fs.readInt16 (_remarksize); - + if (! fs) { - sys_error (ERR_SEVERE, "Error reading header information , _remarksize=%d [projections_read_header]", _remarksize); - return false; + sys_error (ERR_SEVERE, "Error reading header information , _remarksize=%d [projections_read_header]", _remarksize); + return false; } - + if (_signature != m_signature) { sys_error (ERR_SEVERE, "File %s does not have a valid projection file signature", m_filename.c_str()); return false; } - + char* pszRemarkStorage = new char [_remarksize+1]; fs.read (pszRemarkStorage, _remarksize); if (! fs) { @@ -272,11 +325,11 @@ Projections::headerRead (fnetorderstream& fs) pszRemarkStorage[_remarksize] = 0; m_remark = pszRemarkStorage; delete pszRemarkStorage; - + off_t _hsizeread = fs.tellg(); if (!fs || _hsizeread != _hsize) { - sys_error (ERR_WARNING, "File header size read %ld != file header size stored %ld [read_projections_header]\n_remarksize=%ld", (long int) _hsizeread, _hsize, _remarksize); - return false; + sys_error (ERR_WARNING, "File header size read %ld != file header size stored %ld [read_projections_header]\n_remarksize=%ld", (long int) _hsizeread, _hsize, _remarksize); + return false; } m_headerSize = _hsize; @@ -297,12 +350,12 @@ Projections::headerRead (fnetorderstream& fs) m_hour = _hour; m_minute = _minute; m_second = _second; - + m_label.setLabelType (Array2dFileLabel::L_HISTORY); m_label.setLabelString (m_remark); m_label.setCalcTime (m_calcTime); m_label.setDateTime (m_year, m_month, m_day, m_hour, m_minute, m_second); - + return true; } @@ -322,21 +375,21 @@ Projections::read (const char* filename) #else frnetorderstream fileRead (m_filename.c_str(), std::ios::in | std::ios::binary | std::ios::nocreate); #endif - + if (fileRead.fail()) return false; - + if (! headerRead (fileRead)) return false; - + deleteProjData (); newProjData(); - + for (int i = 0; i < m_nView; i++) { if (! detarrayRead (fileRead, *m_projData[i], i)) break; } - + fileRead.close(); return true; } @@ -345,7 +398,7 @@ Projections::read (const char* filename) bool Projections::copyViewData (const std::string& filename, std::ostream& os, int startView, int endView) { - return copyViewData (filename.c_str(), os, startView, endView); + return copyViewData (filename.c_str(), os, startView, endView); } bool @@ -354,56 +407,56 @@ Projections::copyViewData (const char* const filename, std::ostream& os, int sta frnetorderstream is (filename, ios::in | ios::binary); kuint16 sizeHeader, signature; kuint32 _nView, _nDet; - + is.readInt16 (sizeHeader); is.readInt16 (signature); is.readInt32 (_nView); is.readInt32 (_nDet); int nView = _nView; int nDet = _nDet; - + if (signature != m_signature) { sys_error (ERR_FATAL, "Illegal signature in projection file %s", filename); return false; } - + if (startView < 0) - startView = 0; + startView = 0; if (startView > nView - 1) - startView = nView; + startView = nView; if (endView < 0 || endView > nView - 1) - endView = nView - 1; - + endView = nView - 1; + if (startView > endView) { // swap if start > end - int tempView = endView; - endView = startView; - startView = tempView; + int tempView = endView; + endView = startView; + startView = tempView; } - + int sizeView = 8 /* view_angle */ + 4 /* nDet */ + (4 * nDet); unsigned char* pViewData = new unsigned char [sizeView]; - + for (int i = startView; i <= endView; i++) { - is.seekg (sizeHeader + i * sizeView); - is.read (reinterpret_cast(pViewData), sizeView); - os.write (reinterpret_cast(pViewData), sizeView); - if (is.fail() || os.fail()) - break; + is.seekg (sizeHeader + i * sizeView); + is.read (reinterpret_cast(pViewData), sizeView); + os.write (reinterpret_cast(pViewData), sizeView); + if (is.fail() || os.fail()) + break; } - + delete pViewData; if (is.fail()) sys_error (ERR_FATAL, "Error reading projection file"); if (os.fail()) sys_error (ERR_FATAL, "Error writing projection file"); - + return (! (is.fail() | os.fail())); } bool Projections::copyHeader (const std::string& filename, std::ostream& os) { - return copyHeader (filename.c_str(), os); + return copyHeader (filename.c_str(), os); } bool @@ -418,20 +471,20 @@ Projections::copyHeader (const char* const filename, std::ostream& os) sys_error (ERR_FATAL, "Illegal signature in projection file %s", filename); return false; } - + unsigned char* pHdrData = new unsigned char [sizeHeader]; is.read (reinterpret_cast(pHdrData), sizeHeader); if (is.fail()) { - sys_error (ERR_FATAL, "Error reading header"); - return false; + sys_error (ERR_FATAL, "Error reading header"); + return false; } - + os.write (reinterpret_cast(pHdrData), sizeHeader); if (os.fail()) { - sys_error (ERR_FATAL, "Error writing header"); - return false; + sys_error (ERR_FATAL, "Error writing header"); + return false; } - + return true; } @@ -450,32 +503,32 @@ Projections::write (const char* filename) sys_error (ERR_SEVERE, "Error opening file %s for output [projections_create]", filename); return false; } - + if (! headerWrite (fs)) - return false; - + return false; + if (m_projData != NULL) { for (int i = 0; i < m_nView; i++) { if (! detarrayWrite (fs, *m_projData[i], i)) - break; + break; } } if (! fs) return false; - - fs.close(); - + + fs.close(); + return true; } /* NAME - * detarrayRead Read a Detector Array structure from the disk - * - * SYNOPSIS - * detarrayRead (proj, darray, view_num) - * DETARRAY *darray Detector array storage location to be filled - * int view_num View number to read - */ +* detarrayRead Read a Detector Array structure from the disk +* +* SYNOPSIS +* detarrayRead (proj, darray, view_num) +* DETARRAY *darray Detector array storage location to be filled +* int view_num View number to read +*/ bool Projections::detarrayRead (fnetorderstream& fs, DetectorArray& darray, const int iview) @@ -487,39 +540,39 @@ Projections::detarrayRead (fnetorderstream& fs, DetectorArray& darray, const int DetectorValue* detval_ptr = darray.detValues(); kfloat64 view_angle; kuint32 nDet; - + fs.seekg (start_data); - + fs.readFloat64 (view_angle); fs.readInt32 (nDet); darray.setViewAngle (view_angle); // darray.setNDet ( nDet); for (unsigned int i = 0; i < nDet; i++) { - kfloat32 detval; - fs.readFloat32 (detval); - detval_ptr[i] = detval; + kfloat32 detval; + fs.readFloat32 (detval); + detval_ptr[i] = detval; } if (! fs) return false; - + return true; } /* NAME - * detarrayWrite Write detector array data to the disk - * - * SYNOPSIS - * detarrayWrite (darray, view_num) - * DETARRAY *darray Detector array data to be written - * int view_num View number to write - * - * DESCRIPTION - * This routine writes the detarray data from the disk sequentially to - * the file that was opened with open_projections(). Data is written in - * binary format. - */ +* detarrayWrite Write detector array data to the disk +* +* SYNOPSIS +* detarrayWrite (darray, view_num) +* DETARRAY *darray Detector array data to be written +* int view_num View number to write +* +* DESCRIPTION +* This routine writes the detarray data from the disk sequentially to +* the file that was opened with open_projections(). Data is written in +* binary format. +*/ bool Projections::detarrayWrite (fnetorderstream& fs, const DetectorArray& darray, const int iview) @@ -537,27 +590,27 @@ Projections::detarrayWrite (fnetorderstream& fs, const DetectorArray& darray, co sys_error (ERR_SEVERE, "Error seeking detectory array [detarrayWrite]"); return false; } - + fs.writeFloat64 (view_angle); fs.writeInt32 (nDet); - + for (unsigned int i = 0; i < nDet; i++) { kfloat32 detval = detval_ptr[i]; fs.writeFloat32 (detval); } - + if (! fs) return (false); - + return true; } /* NAME - * printProjectionData Print projections data - * - * SYNOPSIS - * printProjectionData () - */ +* printProjectionData Print projections data +* +* SYNOPSIS +* printProjectionData () +*/ void Projections::printProjectionData () @@ -586,7 +639,7 @@ Projections::printProjectionData (int startView, int endView) printf("View %d: angle %f\n", ir, m_projData[ir]->viewAngle()); DetectorValue* detval = m_projData[ir]->detValues(); for (int id = 0; id < m_projData[ir]->nDet(); id++) - printf("%8.4f ", detval[id]); + printf("%8.4f ", detval[id]); printf("\n"); } } @@ -609,14 +662,191 @@ Projections::printScanInfo (std::ostringstream& os) const } -bool Projections::convertPolar (ImageFile& rIF, int iInterpolation) +bool +Projections::convertPolar (ImageFile& rIF, int iInterpolationID) +{ + unsigned int nx = rIF.nx(); + unsigned int ny = rIF.ny(); + ImageFileArray v = rIF.getArray(); + ImageFileArray vImag = rIF.getImaginaryArray(); + + if (! v || nx == 0 || ny == 0) + return false; + + Array2d adView (nx, ny); + Array2d adDet (nx, ny); + double** ppdView = adView.getArray(); + double** ppdDet = adDet.getArray(); + + calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet); + + std::complex** ppcDetValue = new std::complex* [m_nView]; + for (unsigned int iView = 0; iView < m_nView; iView++) { + ppcDetValue[iView] = new std::complex [m_nDet]; + for (unsigned int iDet = 0; iDet < m_nDet; iDet++) + ppcDetValue[iView][iDet] = std::complex(getDetectorArray (iView).detValues()[iDet], 0); + } + + interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, m_nView, m_nDet, iInterpolationID); + + for (iView = 0; iView < m_nView; iView++) + delete [] ppcDetValue[iView]; + delete [] ppcDetValue; + + return true; +} + +void +Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double** ppdView, double** ppdDet) +{ + double xMin = -m_phmLen / 2; + double xMax = xMin + m_phmLen; + double yMin = -m_phmLen / 2; + double yMax = yMin + m_phmLen; + + double xInc = (xMax - xMin) / nx; // size of cells + double yInc = (yMax - yMin) / ny; + + int iDetCenter = (m_nDet - 1) / 2; // index refering to L=0 projection + + if (m_geometry != Scanner::GEOMETRY_PARALLEL) { + sys_error (ERR_WARNING, "convertPolar supports Parallel only"); + } + + double x = xMin + xInc / 2; // Rectang coords of center of pixel + for (unsigned int ix = 0; ix < nx; x += xInc, ix++) { + double y = yMin + yInc / 2; + for (unsigned int iy = 0; iy < ny; y += yInc, iy++) { + double r = ::sqrt (x * x + y * y); + double phi = atan2 (y, x); + + if (phi >= PI) { + phi -= PI; + } else if (phi < 0) { + phi += PI; + } else + r = -r; + + ppdView[ix][iy] = (phi - m_rotStart) / m_rotInc; + ppdDet[ix][iy] = (r / m_detInc) + iDetCenter; + } + } +} + +void +Projections::interpolatePolar (ImageFileArray& v, ImageFileArray& vImag, + unsigned int nx, unsigned int ny, std::complex** ppcDetValue, + double** ppdView, double** ppdDet, unsigned int nView, unsigned int nDet, int iInterpolationID) { - return false; + for (unsigned int ix = 0; ix < ny; ix++) { + for (unsigned int iy = 0; iy < ny; iy++) { + if (iInterpolationID == POLAR_INTERP_NEAREST) { + int iView = nearest (ppdView[ix][iy]); + int iDet = nearest (ppdDet[ix][iy]); + if (iView == nView) + iView = nView - 1; + if (iDet >= 0 && iDet < nDet && iView >= 0 && iView < nView) { + v[ix][iy] = ppcDetValue[iView][iDet].real(); + if (vImag) + vImag[ix][iy] = ppcDetValue[iView][iDet].imag(); + } else { + sys_error (ERR_SEVERE, "Can't find projection data for ix=%d,iy=%d with radView=%f and radDet=%f", + ix, iy, ppdView[ix][iy], ppdDet[ix][iy]); + v[ix][iy] = 0; + } + } else if (iInterpolationID == POLAR_INTERP_BILINEAR) { + int iFloorView = ppdView[ix][iy]; + double dFracView = ppdView[ix][iy] - iFloorView; + int iFloorDet = ppdDet[ix][iy]; + double dFracDet = ppdDet[ix][iy] - iFloorDet; + + if (iFloorDet >= 0 && iFloorView >= 0) { + std::complex v1 = ppcDetValue[iFloorView][iFloorDet]; + std::complex v2, v3, v4; + if (iFloorView < nView - 1) + v2 = ppcDetValue[iFloorView + 1][iFloorDet]; + else + v2 = v1; + if (iFloorDet < nDet - 1) + v4 = ppcDetValue[iFloorView][iFloorDet+1]; + else + v4 = v1; + if (iFloorView < nView - 1 && iFloorDet < nDet - 1) + v3 = ppcDetValue [iFloorView+1][iFloorDet+1]; + else if (iFloorView < nView - 1) + v4 = v2; + else + v4 = v1; + std::complex vInterp = v1 + dFracView * (v2 - v1) + dFracDet * (v4 - v1); + v[ix][iy] = vInterp.real(); + if (vImag) + vImag[ix][iy] = vInterp.imag(); + } else { + sys_error (ERR_SEVERE, "Can't find projection data for ix=%d,iy=%d with radView=%f and radDet=%f", + ix, iy, ppdView[ix][iy], ppdDet[ix][iy]); + v[ix][iy] = 0; + if (vImag) + vImag[ix][iy] = 0; + } + } else if (iInterpolationID == POLAR_INTERP_BICUBIC) { + v[ix][iy] =0; + if (vImag) + vImag[ix][iy] = 0; + } + } + } } -bool Projections::convertFFTPolar (ImageFile& rIF, int iInterpolation, int iZeropad) +bool +Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad) { - return false; + unsigned int nx = rIF.nx(); + unsigned int ny = rIF.ny(); + ImageFileArray v = rIF.getArray(); + if (! rIF.isComplex()) + rIF.convertRealToComplex(); + ImageFileArray vImag = rIF.getImaginaryArray(); + + if (! v || nx == 0 || ny == 0) + return false; + + Array2d adView (nx, ny); + Array2d adDet (nx, ny); + double** ppdView = adView.getArray(); + double** ppdDet = adDet.getArray(); + + std::complex** ppcDetValue = new std::complex* [m_nView]; + unsigned int iView; + double* pdDet = new double [m_nDet]; + fftw_complex* pcIn = new fftw_complex [m_nDet]; + fftw_plan plan = fftw_create_plan (m_nDet, FFTW_FORWARD, FFTW_IN_PLACE); + + for (iView = 0; iView < m_nView; iView++) { + ppcDetValue[iView] = new std::complex [m_nDet]; + unsigned int iDet; + for (iDet = 0; iDet < m_nDet; iDet++) { + pcIn[iDet].re = getDetectorArray(iView).detValues()[iDet]; + pcIn[iDet].im = 0; + } + fftw_one (plan, pcIn, NULL); + for (iDet = 0; iDet < m_nDet; iDet++) + ppcDetValue[iView][iDet] = std::complex (pcIn[iDet].re, pcIn[iDet].im); + Fourier::shuffleFourierToNaturalOrder (ppcDetValue[iView], m_nDet); + } + + fftw_destroy_plan (plan); + + delete [] pcIn; + + calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet); + + interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, m_nView, m_nDet, iInterpolationID); + + for (iView = 0; iView < m_nView; iView++) + delete [] ppcDetValue[iView]; + delete [] ppcDetValue; + + return true; } diff --git a/libctsim/scanner.cpp b/libctsim/scanner.cpp index c593b09..e10446f 100644 --- a/libctsim/scanner.cpp +++ b/libctsim/scanner.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: scanner.cpp,v 1.22 2000/12/18 06:32:13 kevin Exp $ +** $Id: scanner.cpp,v 1.23 2001/01/04 21:28:41 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 @@ -71,21 +71,21 @@ DetectorArray::~DetectorArray (void) /* NAME - * Scanner::Scanner Construct a user specified detector structure - * - * SYNOPSIS - * Scanner (phm, nDet, nView, nSample) - * Phantom& phm PHANTOM that we are making detector for - * int geomety Geometry of detector - * int nDet Number of detector along detector array - * int nView Number of rotated views - * int nSample Number of rays per detector - */ +* Scanner::Scanner Construct a user specified detector structure +* +* SYNOPSIS +* Scanner (phm, nDet, nView, nSample) +* Phantom& phm PHANTOM that we are making detector for +* int geomety Geometry of detector +* int nDet Number of detector along detector array +* int nView Number of rotated views +* int nSample Number of rays per detector +*/ Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet, int nView, int nSample, const double rot_anglen, const double dFocalLengthRatio, const double dFieldOfViewRatio) { m_phmLen = phm.maxAxisLength(); // maximal length along an axis - + m_fail = false; m_idGeometry = convertGeometryNameToID (geometryName); if (m_idGeometry == GEOMETRY_INVALID) { @@ -94,7 +94,7 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet, m_failMessage += geometryName; return; } - + if (nView < 1 || nDet < 1) { m_fail = true; m_failMessage = "nView & nDet must be greater than 0"; @@ -102,7 +102,7 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet, } if (nSample < 1) m_nSample = 1; - + m_nDet = nDet; m_nView = nView; m_nSample = nSample; @@ -110,7 +110,7 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet, m_dFieldOfViewRatio = dFieldOfViewRatio; m_dFocalLength = (m_phmLen * SQRT2 / 2) * dFocalLengthRatio; m_dFieldOfView = m_phmLen * SQRT2 * dFieldOfViewRatio; - + m_dXCenter = phm.xmin() + (phm.xmax() - phm.xmin()) / 2; m_dYCenter = phm.ymin() + (phm.ymax() - phm.ymin()) / 2; m_rotLen = rot_anglen; @@ -145,7 +145,7 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet, double dAngle = atan( dHalfSquare / dFocalPastPhm ); #endif double dHalfDetLen = 2 * m_dFocalLength * tan (dAngle); - + m_detLen = dHalfDetLen * 2; m_detInc = m_detLen / m_nDet; if (m_nDet % 2 == 0) // Adjust for Even number of detectors @@ -181,14 +181,14 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet, m_dAngularDetIncrement = m_detInc * 2; // Angular Position 2x gamma angle m_dAngularDetLen = m_detLen * 2; m_initPos.dAngularDet = -m_dAngularDetLen / 2; - + m_initPos.angle = 0; m_initPos.xs1 = m_dXCenter; m_initPos.ys1 = m_dYCenter + m_dFocalLength;; m_initPos.xs2 = m_dXCenter; m_initPos.ys2 = m_dYCenter + m_dFocalLength; } - + // Calculate incrementatal rotation matrix GRFMTX_2D temp; xlat_mtx2 (m_rotmtxIncrement, -m_dXCenter, -m_dYCenter); @@ -208,10 +208,10 @@ const char* Scanner::convertGeometryIDToName (const int geomID) { const char *name = ""; - + if (geomID >= 0 && geomID < s_iGeometryCount) - return (s_aszGeometryName[geomID]); - + return (s_aszGeometryName[geomID]); + return (name); } @@ -219,37 +219,37 @@ const char* Scanner::convertGeometryIDToTitle (const int geomID) { const char *title = ""; - + if (geomID >= 0 && geomID < s_iGeometryCount) - return (s_aszGeometryName[geomID]); - + return (s_aszGeometryName[geomID]); + return (title); } - + int Scanner::convertGeometryNameToID (const char* const geomName) { int id = GEOMETRY_INVALID; - + for (int i = 0; i < s_iGeometryCount; i++) - if (strcasecmp (geomName, s_aszGeometryName[i]) == 0) { - id = i; - break; - } - - return (id); + if (strcasecmp (geomName, s_aszGeometryName[i]) == 0) { + id = i; + break; + } + + return (id); } - + /* NAME - * collectProjections Calculate projections for a Phantom - * - * SYNOPSIS - * collectProjections (proj, phm, start_view, nView, bStoreViewPos, trace) - * Projectrions& proj Projection storage - * Phantom& phm Phantom for which we collect projections - * bool bStoreViewPos TRUE then storage proj at normal view position - * int trace Trace level +* collectProjections Calculate projections for a Phantom +* +* SYNOPSIS +* collectProjections (proj, phm, start_view, nView, bStoreViewPos, trace) +* Projectrions& proj Projection storage +* Phantom& phm Phantom for which we collect projections +* bool bStoreViewPos TRUE then storage proj at normal view position +* int trace Trace level */ @@ -264,7 +264,7 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS { m_trace = trace; double start_angle = iStartView * proj.rotInc(); - + // Calculate initial rotation matrix GRFMTX_2D rotmtx_initial, temp; xlat_mtx2 (rotmtx_initial, -m_dXCenter, -m_dYCenter); @@ -272,15 +272,15 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS mult_mtx2 (rotmtx_initial, temp, rotmtx_initial); xlat_mtx2 (temp, m_dXCenter, m_dYCenter); mult_mtx2 (rotmtx_initial, temp, rotmtx_initial); - + double xd1=0, yd1=0, xd2=0, yd2=0; if (m_idGeometry != GEOMETRY_EQUIANGULAR) { - xd1 = m_initPos.xd1; - yd1 = m_initPos.yd1; - xd2 = m_initPos.xd2; - yd2 = m_initPos.yd2; - xform_mtx2 (rotmtx_initial, xd1, yd1); // rotate detector endpoints - xform_mtx2 (rotmtx_initial, xd2, yd2); // to initial view_angle + xd1 = m_initPos.xd1; + yd1 = m_initPos.yd1; + xd2 = m_initPos.xd2; + yd2 = m_initPos.yd2; + xform_mtx2 (rotmtx_initial, xd1, yd1); // rotate detector endpoints + xform_mtx2 (rotmtx_initial, xd2, yd2); // to initial view_angle } double xs1 = m_initPos.xs1; @@ -289,18 +289,18 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS double ys2 = m_initPos.ys2; xform_mtx2 (rotmtx_initial, xs1, ys1); // rotate source endpoints to xform_mtx2 (rotmtx_initial, xs2, ys2); // initial view angle - + int iView; double viewAngle; for (iView = 0, viewAngle = start_angle; iView < iNumViews; iView++, viewAngle += proj.rotInc()) { int iStoragePosition = iView; if (bStoreAtViewPosition) iStoragePosition += iStartView; - + DetectorArray& detArray = proj.getDetectorArray( iStoragePosition ); - + #ifdef HAVE_SGP - if (pSGP && m_trace >= Trace::TRACE_PHANTOM) { + if (pSGP && m_trace >= Trace::TRACE_PHANTOM) { m_pSGP = pSGP; double dWindowSize = dmax (m_detLen, m_dFocalLength * 2) * SQRT2; double dHalfWindowSize = dWindowSize / 2; @@ -309,53 +309,53 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS m_dYMinWin = m_dYCenter - dHalfWindowSize; m_dYMaxWin = m_dYCenter + dHalfWindowSize; double dHalfPhmLen = m_phmLen / 2; - - m_pSGP->setWindow (m_dXMinWin, m_dYMinWin, m_dXMaxWin, m_dYMaxWin); - m_pSGP->setRasterOp (RO_COPY); - m_pSGP->setColor (C_RED); - m_pSGP->moveAbs (0., 0.); - m_pSGP->drawRect (m_dXCenter - dHalfPhmLen, m_dYCenter - dHalfPhmLen, m_dXCenter + dHalfPhmLen, m_dYCenter + dHalfPhmLen); - m_pSGP->moveAbs (0., 0.); - m_pSGP->drawCircle (m_dFocalLength); - m_pSGP->setColor (C_BLUE); - phm.draw (*m_pSGP); - m_dTextHeight = m_pSGP->getCharHeight (); - - traceShowParam ("Phantom:", "%s", PROJECTION_TRACE_ROW_PHANT_ID, C_BLACK, phm.name().c_str()); - traceShowParam ("Geometry:", "%s", PROJECTION_TRACE_ROW_GEOMETRY, C_BLUE, convertGeometryIDToName(m_idGeometry)); - traceShowParam ("Focal Length Ratio:", "%.2f", PROJECTION_TRACE_ROW_FOCAL_LENGTH, C_BLUE, m_dFocalLengthRatio); - traceShowParam ("Field Of View Ratio:", "%.2f", PROJECTION_TRACE_ROW_FIELD_OF_VIEW, C_BLUE, m_dFieldOfViewRatio); - traceShowParam ("Num Detectors:", "%d", PROJECTION_TRACE_ROW_NDET, C_BLUE, proj.nDet()); - traceShowParam ("Num Views:", "%d", PROJECTION_TRACE_ROW_NVIEW, C_BLUE, proj.nView()); - traceShowParam ("Samples / Ray:", "%d", PROJECTION_TRACE_ROW_SAMPLES, C_BLUE, m_nSample); - - m_pSGP->setMarker (SGP::MARK_BDIAMOND, C_LTGREEN); - } + + m_pSGP->setWindow (m_dXMinWin, m_dYMinWin, m_dXMaxWin, m_dYMaxWin); + m_pSGP->setRasterOp (RO_COPY); + m_pSGP->setColor (C_RED); + m_pSGP->moveAbs (0., 0.); + m_pSGP->drawRect (m_dXCenter - dHalfPhmLen, m_dYCenter - dHalfPhmLen, m_dXCenter + dHalfPhmLen, m_dYCenter + dHalfPhmLen); + m_pSGP->moveAbs (0., 0.); + m_pSGP->drawCircle (m_dFocalLength); + m_pSGP->setColor (C_BLUE); + phm.draw (*m_pSGP); + m_dTextHeight = m_pSGP->getCharHeight (); + + traceShowParam ("Phantom:", "%s", PROJECTION_TRACE_ROW_PHANT_ID, C_BLACK, phm.name().c_str()); + traceShowParam ("Geometry:", "%s", PROJECTION_TRACE_ROW_GEOMETRY, C_BLUE, convertGeometryIDToName(m_idGeometry)); + traceShowParam ("Focal Length Ratio:", "%.2f", PROJECTION_TRACE_ROW_FOCAL_LENGTH, C_BLUE, m_dFocalLengthRatio); + traceShowParam ("Field Of View Ratio:", "%.2f", PROJECTION_TRACE_ROW_FIELD_OF_VIEW, C_BLUE, m_dFieldOfViewRatio); + traceShowParam ("Num Detectors:", "%d", PROJECTION_TRACE_ROW_NDET, C_BLUE, proj.nDet()); + traceShowParam ("Num Views:", "%d", PROJECTION_TRACE_ROW_NVIEW, C_BLUE, proj.nView()); + traceShowParam ("Samples / Ray:", "%d", PROJECTION_TRACE_ROW_SAMPLES, C_BLUE, m_nSample); + + m_pSGP->setMarker (SGP::MARK_BDIAMOND, C_LTGREEN); + } #endif - + #ifdef HAVE_SGP if (m_pSGP && m_trace >= Trace::TRACE_PHANTOM) { m_pSGP->setColor (C_BLACK); m_pSGP->setPenWidth (2); if (m_idGeometry == GEOMETRY_PARALLEL) { - m_pSGP->moveAbs (xs1, ys1); - m_pSGP->lineAbs (xs2, ys2); - m_pSGP->moveAbs (xd1, yd1); - m_pSGP->lineAbs (xd2, yd2); + m_pSGP->moveAbs (xs1, ys1); + m_pSGP->lineAbs (xs2, ys2); + m_pSGP->moveAbs (xd1, yd1); + m_pSGP->lineAbs (xd2, yd2); } else if (m_idGeometry == GEOMETRY_EQUILINEAR) { - m_pSGP->setPenWidth (4); - m_pSGP->moveAbs (xs1, ys1); - m_pSGP->lineAbs (xs2, ys2); - m_pSGP->setPenWidth (2); - m_pSGP->moveAbs (xd1, yd1); - m_pSGP->lineAbs (xd2, yd2); + m_pSGP->setPenWidth (4); + m_pSGP->moveAbs (xs1, ys1); + m_pSGP->lineAbs (xs2, ys2); + m_pSGP->setPenWidth (2); + m_pSGP->moveAbs (xd1, yd1); + m_pSGP->lineAbs (xd2, yd2); } else if (m_idGeometry == GEOMETRY_EQUIANGULAR) { - m_pSGP->setPenWidth (4); - m_pSGP->moveAbs (xs1, ys1); - m_pSGP->lineAbs (xs2, ys2); - m_pSGP->setPenWidth (2); - m_pSGP->moveAbs (0., 0.); - m_pSGP->drawArc (m_dFocalLength, viewAngle + 3 * HALFPI - (m_dAngularDetLen/2), viewAngle + 3 * HALFPI + (m_dAngularDetLen/2)); + m_pSGP->setPenWidth (4); + m_pSGP->moveAbs (xs1, ys1); + m_pSGP->lineAbs (xs2, ys2); + m_pSGP->setPenWidth (2); + m_pSGP->moveAbs (0., 0.); + m_pSGP->drawArc (m_dFocalLength, viewAngle + 3 * HALFPI - (m_dAngularDetLen/2), viewAngle + 3 * HALFPI + (m_dAngularDetLen/2)); } m_pSGP->setPenWidth (1); } @@ -363,11 +363,11 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS traceShowParam ("Current View:", "%d (%.0f%%)", PROJECTION_TRACE_ROW_CURR_VIEW, C_RED, iView + iStartView, (iView + iStartView) / static_cast(m_nView) * 100.); #endif if (m_trace == Trace::TRACE_CONSOLE) - std::cout << "Current View: " << iView+iStartView << std::endl; - + std::cout << "Current View: " << iView+iStartView << std::endl; + projectSingleView (phm, detArray, xd1, yd1, xd2, yd2, xs1, ys1, xs2, ys2, viewAngle + 3 * HALFPI); detArray.setViewAngle (viewAngle); - + #ifdef HAVE_SGP if (m_pSGP && m_trace >= Trace::TRACE_PHANTOM) { // rs_plot (detArray, xd1, yd1, dXCenter, dYCenter, theta); @@ -376,124 +376,124 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS xform_mtx2 (m_rotmtxIncrement, xs1, ys1); xform_mtx2 (m_rotmtxIncrement, xs2, ys2); if (m_idGeometry != GEOMETRY_EQUIANGULAR) { - xform_mtx2 (m_rotmtxIncrement, xd1, yd1); // rotate detector endpoints - xform_mtx2 (m_rotmtxIncrement, xd2, yd2); + xform_mtx2 (m_rotmtxIncrement, xd1, yd1); // rotate detector endpoints + xform_mtx2 (m_rotmtxIncrement, xd2, yd2); } } /* for each iView */ } /* NAME - * rayview Calculate raysums for a view at any angle - * - * SYNOPSIS - * rayview (phm, detArray, xd1, nSample, yd1, xd2, yd2, xs1, ys1, xs2, ys2) - * Phantom& phm Phantom to scan - * DETARRAY *detArray Storage of values for detector array - * Scanner& det Scanner parameters - * double xd1, yd1, xd2, yd2 Beginning & ending detector positions - * double xs1, ys1, xs2, ys2 Beginning & ending source positions - * - * RAY POSITIONING - * For each detector, have there are a variable number of rays traced. - * The source of each ray is the center of the source x-ray cell. The - * detector positions are equally spaced within the cell - * - * The increments between rays are calculated so that the cells start - * at the beginning of a detector cell and they end on the endpoint - * of the cell. Thus, the last cell starts at (xd2-ddx),(yd2-ddy). - * The exception to this is if there is only one ray per detector. - * In that case, the detector position is the center of the detector cell. - */ +* rayview Calculate raysums for a view at any angle +* +* SYNOPSIS +* rayview (phm, detArray, xd1, nSample, yd1, xd2, yd2, xs1, ys1, xs2, ys2) +* Phantom& phm Phantom to scan +* DETARRAY *detArray Storage of values for detector array +* Scanner& det Scanner parameters +* double xd1, yd1, xd2, yd2 Beginning & ending detector positions +* double xs1, ys1, xs2, ys2 Beginning & ending source positions +* +* RAY POSITIONING +* For each detector, have there are a variable number of rays traced. +* The source of each ray is the center of the source x-ray cell. The +* detector positions are equally spaced within the cell +* +* The increments between rays are calculated so that the cells start +* at the beginning of a detector cell and they end on the endpoint +* of the cell. Thus, the last cell starts at (xd2-ddx),(yd2-ddy). +* The exception to this is if there is only one ray per detector. +* In that case, the detector position is the center of the detector cell. +*/ void Scanner::projectSingleView (const Phantom& phm, DetectorArray& detArray, const double xd1, const double yd1, const double xd2, const double yd2, const double xs1, const double ys1, const double xs2, const double ys2, const double dDetAngle) { - + double sdx = (xs2 - xs1) / detArray.nDet(); // change in coords double sdy = (ys2 - ys1) / detArray.nDet(); // between source double xs_maj = xs1 + (sdx / 2); // put ray source in center of cell double ys_maj = ys1 + (sdy / 2); - + double ddx=0, ddy=0, ddx2=0, ddy2=0, ddx2_ofs=0, ddy2_ofs=0, xd_maj=0, yd_maj=0; double dAngleInc=0, dAngleSampleInc=0, dAngleSampleOffset=0, dAngleMajor=0; if (m_idGeometry == GEOMETRY_EQUIANGULAR) { - dAngleInc = m_dAngularDetIncrement; - dAngleSampleInc = dAngleInc / m_nSample; - dAngleSampleOffset = dAngleSampleInc / 2; - dAngleMajor = dDetAngle - (m_dAngularDetLen/2) + dAngleSampleOffset; + dAngleInc = m_dAngularDetIncrement; + dAngleSampleInc = dAngleInc / m_nSample; + dAngleSampleOffset = dAngleSampleInc / 2; + dAngleMajor = dDetAngle - (m_dAngularDetLen/2) + dAngleSampleOffset; } else { - ddx = (xd2 - xd1) / detArray.nDet(); // change in coords - ddy = (yd2 - yd1) / detArray.nDet(); // between detectors - ddx2 = ddx / m_nSample; // Incr. between rays with detector cell - ddy2 = ddy / m_nSample; // Doesn't include detector endpoints - ddx2_ofs = ddx2 / 2; // offset of 1st ray from start of detector cell - ddy2_ofs = ddy2 / 2; - - xd_maj = xd1 + ddx2_ofs; // Incr. between detector cells - yd_maj = yd1 + ddy2_ofs; + ddx = (xd2 - xd1) / detArray.nDet(); // change in coords + ddy = (yd2 - yd1) / detArray.nDet(); // between detectors + ddx2 = ddx / m_nSample; // Incr. between rays with detector cell + ddy2 = ddy / m_nSample; // Doesn't include detector endpoints + ddx2_ofs = ddx2 / 2; // offset of 1st ray from start of detector cell + ddy2_ofs = ddy2 / 2; + + xd_maj = xd1 + ddx2_ofs; // Incr. between detector cells + yd_maj = yd1 + ddy2_ofs; } - + DetectorValue* detval = detArray.detValues(); - + if (phm.getComposition() == P_UNIT_PULSE) { // put unit pulse in center of view for (int d = 0; d < detArray.nDet(); d++) if (detArray.nDet() / 2 == d && (d % 2) == 1) - detval[d] = 1; + detval[d] = 1; else - detval[d] = 0; + detval[d] = 0; } else { - for (int d = 0; d < detArray.nDet(); d++) { + for (int d = 0; d < detArray.nDet(); d++) { double xs = xs_maj; double ys = ys_maj; double xd=0, yd=0, dAngle=0; if (m_idGeometry == GEOMETRY_EQUIANGULAR) { - dAngle = dAngleMajor; + dAngle = dAngleMajor; } else { - xd = xd_maj; - yd = yd_maj; + xd = xd_maj; + yd = yd_maj; } double sum = 0.0; for (unsigned int i = 0; i < m_nSample; i++) { - if (m_idGeometry == GEOMETRY_EQUIANGULAR) { - xd = m_dFocalLength * cos (dAngle); - yd = m_dFocalLength * sin (dAngle); - } - + if (m_idGeometry == GEOMETRY_EQUIANGULAR) { + xd = m_dFocalLength * cos (dAngle); + yd = m_dFocalLength * sin (dAngle); + } + #ifdef HAVE_SGP - if (m_pSGP && m_trace >= Trace::TRACE_PROJECTIONS) { - m_pSGP->setColor (C_YELLOW); - m_pSGP->setRasterOp (RO_AND); - m_pSGP->moveAbs (xs, ys); - m_pSGP->lineAbs (xd, yd); - } + if (m_pSGP && m_trace >= Trace::TRACE_PROJECTIONS) { + m_pSGP->setColor (C_YELLOW); + m_pSGP->setRasterOp (RO_AND); + m_pSGP->moveAbs (xs, ys); + m_pSGP->lineAbs (xd, yd); + } #endif - - sum += projectSingleLine (phm, xd, yd, xs, ys); - + + sum += projectSingleLine (phm, xd, yd, xs, ys); + #ifdef HAVE_SGP - // if (m_trace >= Trace::TRACE_CLIPPING) { - // traceShowParam ("Attenuation:", "%s", PROJECTION_TRACE_ROW_ATTEN, C_LTMAGENTA, " "); - // traceShowParam ("Attenuation:", "%.3f", PROJECTION_TRACE_ROW_ATTEN, C_LTMAGENTA, sum); - // } + // if (m_trace >= Trace::TRACE_CLIPPING) { + // traceShowParam ("Attenuation:", "%s", PROJECTION_TRACE_ROW_ATTEN, C_LTMAGENTA, " "); + // traceShowParam ("Attenuation:", "%.3f", PROJECTION_TRACE_ROW_ATTEN, C_LTMAGENTA, sum); + // } #endif - if (m_idGeometry == GEOMETRY_EQUIANGULAR) - dAngle += dAngleSampleInc; - else { - xd += ddx2; - yd += ddy2; - } + if (m_idGeometry == GEOMETRY_EQUIANGULAR) + dAngle += dAngleSampleInc; + else { + xd += ddx2; + yd += ddy2; + } } // for each sample in detector - + detval[d] = sum / m_nSample; xs_maj += sdx; ys_maj += sdy; if (m_idGeometry == GEOMETRY_EQUIANGULAR) - dAngleMajor += dAngleInc; - else { - xd_maj += ddx; - yd_maj += ddy; - } + dAngleMajor += dAngleInc; + else { + xd_maj += ddx; + yd_maj += ddy; + } } /* for each detector */ } /* if not unit pulse */ } @@ -529,9 +529,9 @@ void Scanner::traceShowParamRasterOp (int iRasterOp, const char *szLabel, const char *fmt, int row, int color, va_list args) { char szValue[256]; - + vsnprintf (szValue, sizeof(szValue), fmt, args); - + #ifdef HAVE_SGP if (m_pSGP) { m_pSGP->setRasterOp (iRasterOp); @@ -544,24 +544,24 @@ Scanner::traceShowParamRasterOp (int iRasterOp, const char *szLabel, const char m_pSGP->drawText (szValue); } else #endif - { - cio_put_str (szLabel); - cio_put_str (szValue); - cio_put_str ("\n"); - } + { + cio_put_str (szLabel); + cio_put_str (szValue); + cio_put_str ("\n"); + } } /* NAME - * projectSingleLine INTERNAL: Calculates raysum along a line for a Phantom - * - * SYNOPSIS - * rsum = phm_ray_attenuation (phm, x1, y1, x2, y2) - * double rsum Ray sum of Phantom along given line - * Phantom& phm; Phantom from which to calculate raysum - * double *x1, *y1, *x2, y2 Endpoints of ray path (in Phantom coords) - */ +* projectSingleLine INTERNAL: Calculates raysum along a line for a Phantom +* +* SYNOPSIS +* rsum = phm_ray_attenuation (phm, x1, y1, x2, y2) +* double rsum Ray sum of Phantom along given line +* Phantom& phm; Phantom from which to calculate raysum +* double *x1, *y1, *x2, y2 Endpoints of ray path (in Phantom coords) +*/ double Scanner::projectSingleLine (const Phantom& phm, const double x1, const double y1, const double x2, const double y2) @@ -570,20 +570,20 @@ Scanner::projectSingleLine (const Phantom& phm, const double x1, const double y1 double rsum = 0.0; for (PElemConstIterator i = phm.listPElem().begin(); i != phm.listPElem().end(); i++) rsum += projectLineAgainstPElem (**i, x1, y1, x2, y2); - + return (rsum); } /* NAME - * pelem_ray_attenuation Calculate raysum of an pelem along one line - * - * SYNOPSIS - * rsum = pelem_ray_attenuation (pelem, x1, y1, x2, y2) - * double rsum Computed raysum - * PhantomElement& pelem Pelem to scan - * double x1, y1, x2, y2 Endpoints of raysum line - */ +* pelem_ray_attenuation Calculate raysum of an pelem along one line +* +* SYNOPSIS +* rsum = pelem_ray_attenuation (pelem, x1, y1, x2, y2) +* double rsum Computed raysum +* PhantomElement& pelem Pelem to scan +* double x1, y1, x2, y2 Endpoints of raysum line +*/ double Scanner::projectLineAgainstPElem (const PhantomElement& pelem, double x1, double y1, double x2, double y2) @@ -593,7 +593,7 @@ Scanner::projectLineAgainstPElem (const PhantomElement& pelem, double x1, double cio_tone (1000., 0.05); return (0.0); } - + #ifdef HAVE_SGP if (m_pSGP && m_trace == Trace::TRACE_CLIPPING) { m_pSGP->setRasterOp (RO_XOR); @@ -605,7 +605,7 @@ Scanner::projectLineAgainstPElem (const PhantomElement& pelem, double x1, double m_pSGP->setRasterOp (RO_SET); } #endif - + double len = lineLength (x1, y1, x2, y2); return (len * pelem.atten()); } diff --git a/msvc/ctsim/ctsim.plg b/msvc/ctsim/ctsim.plg index 27d403b..5ff72f1 100644 --- a/msvc/ctsim/ctsim.plg +++ b/msvc/ctsim/ctsim.plg @@ -3,41 +3,44 @@ 
 
Build Log

 

---------------------Configuration: ctsim - Win32 Release--------------------
+--------------------Configuration: ctsim - Win32 Debug--------------------
 

 
Command Lines

-Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP2E.tmp" with contents
+Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP1A0.tmp" with contents
 [
-/nologo /G6 /MT /W3 /GR /GX /O2 /I "." /I "..\..\include" /I "..\..\getopt" /I "..\..\..\lpng108" /I "..\..\..\zlib" /I "..\..\..\fftw-2.1.3\fftw" /I "..\..\..\fftw-2.1.3\rfftw" /I "..\..\..\wx2\include" /D "NDEBUG" /D "__WXWIN__" /D "__WXMSW__" /D "HAVE_SGP" /D "HAVE_PNG" /D "HAVE_WXWINDOWS" /D "WIN32" /D "_WINDOWS" /D "_MBCS" /D "HAVE_STRING_H" /D "HAVE_FFTW" /D "HAVE_RFFTW" /D "HAVE_GETOPT_H" /D "MSVC" /D "__WIN95__" /D "__WIN32__" /D WINVER=0x0400 /D "STRICT" /D CTSIMVERSION=\"3.0.0beta1\" /FR"Release/" /Fp"Release/ctsim.pch" /YX /Fo"Release/" /Fd"Release/" /FD /c 
-"C:\ctsim\src\views.cpp"
+/nologo /G6 /MTd /W3 /Gm /GR /GX /ZI /Od /I "\wx2\include" /I "." /I "..\..\include" /I "..\..\getopt" /I "..\..\..\lpng108" /I "..\..\..\zlib" /I "..\..\..\fftw-2.1.3\fftw" /I "..\..\..\fftw-2.1.3\rfftw" /D VERSION=\"2.5.0\" /D "_DEBUG" /D "__WXMSW__" /D "HAVE_SGP" /D "HAVE_PNG" /D "HAVE_WXWINDOWS" /D "WIN32" /D "_WINDOWS" /D "_MBCS" /D "HAVE_STRING_H" /D "HAVE_FFTW" /D "HAVE_RFFTW" /D "HAVE_GETOPT_H" /D "MSVC" /D "__WIN95__" /D "__WIN32__" /D WINVER=0x0400 /D "STRICT" /D CTSIMVERSION=\"3.0.0beta1\" /FR"Debug/" /Fp"Debug/ctsim.pch" /YX /Fo"Debug/" /Fd"Debug/" /FD /GZ /c 
+"C:\ctsim\src\dialogs.cpp"
 ]
-Creating command line "cl.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP2E.tmp" 
-Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP2F.tmp" with contents
+Creating command line "cl.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP1A0.tmp" 
+Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP1A1.tmp" with contents
 [
-kernel32.lib user32.lib wsock32.lib comctl32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib libpng.lib zlib.lib /nologo /subsystem:windows /incremental:no /pdb:"Release/ctsim.pdb" /machine:I386 /out:"Release/ctsim.exe" /libpath:"..\..\..\lpng108\msvc\win32\libpng\lib" /libpath:"..\..\..\lpng108\msvc\win32\zlib\lib" 
-.\Release\ctsim.obj
-.\Release\dialogs.obj
-.\Release\dlgprojections.obj
-.\Release\dlgreconstruct.obj
-.\Release\docs.obj
-.\Release\views.obj
-\ctsim\msvc\libctsim\Release\libctsim.lib
-"\fftw-2.1.3\Win32\FFTW2st\Release\FFTW2st.lib"
-"\fftw-2.1.3\Win32\RFFTW2st\Release\RFFTW2st.lib"
-\wx2\lib\wx.lib
+comctl32.lib winmm.lib rpcrt4.lib ws2_32.lib kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib ../libctsim/Debug/libctsim.lib kernel32.lib user32.lib gdi32.lib winspool.lib comdlg32.lib advapi32.lib shell32.lib ole32.lib oleaut32.lib uuid.lib odbc32.lib odbccp32.lib ..\..\..\lpng108\msvc\win32\libpng\lib_dbg\libpng.lib ..\..\..\lpng108\msvc\win32\zlib\lib_dbg\zlib.lib libcmtd.lib ..\..\..\fftw-2.1.3\Win32\FFTW2st\Debug\FFTW2st.lib ..\..\..\fftw-2.1.3\Win32\RFFTW2st\Debug\RFFTW2st.lib ../../../wx2/lib/wxd.lib /nologo /subsystem:windows /incremental:yes /pdb:"Debug/ctsim.pdb" /debug /machine:I386 /nodefaultlib:"libcd.lib" /nodefaultlib:"libcid.lib" /nodefaultlib:"msvcrtd.lib" /out:"Debug/ctsim.exe" /pdbtype:sept /libpath:"..\..\..\lpng108\msvc\win32\libpng\lib" /libpath:"..\..\..\lpng108\msvc\win32\zlib\lib" 
+.\Debug\ctsim.obj
+.\Debug\dialogs.obj
+.\Debug\dlgprojections.obj
+.\Debug\dlgreconstruct.obj
+.\Debug\docs.obj
+.\Debug\views.obj
+.\Debug\wx.res
+\ctsim\msvc\libctsim\Debug\libctsim.lib
+"\fftw-2.1.3\Win32\FFTW2st\Debug\FFTW2st.lib"
+"\fftw-2.1.3\Win32\RFFTW2st\Debug\RFFTW2st.lib"
+\wx2\lib\wxd.lib
 ]
-Creating command line "link.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP2F.tmp"
+Creating command line "link.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP1A1.tmp"
 
Output Window

 Compiling...
-views.cpp
+dialogs.cpp
 Linking...
-LINK : warning LNK4089: all references to "ADVAPI32.dll" discarded by /OPT:REF
-LINK : warning LNK4089: all references to "WSOCK32.dll" discarded by /OPT:REF
+Creating command line "bscmake.exe /nologo /o"Debug/ctsim.bsc"  .\Debug\ctsim.sbr .\Debug\dialogs.sbr .\Debug\dlgprojections.sbr .\Debug\dlgreconstruct.sbr .\Debug\docs.sbr .\Debug\views.sbr"
+Creating browse info file...
+BSCMAKE: warning BK4503 : minor error in .SBR file '.\Debug\dialogs.sbr' ignored
+
Output Window

 
 
 
 
Results

-ctsim.exe - 0 error(s), 2 warning(s)
+ctsim.exe - 0 error(s), 1 warning(s)
diff --git a/src/ctsim.h b/src/ctsim.h index 31d02a4..e1055ba 100644 --- a/src/ctsim.h +++ b/src/ctsim.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: ctsim.h,v 1.20 2001/01/03 22:00:46 kevin Exp $ +** $Id: ctsim.h,v 1.21 2001/01/04 21:28:41 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 @@ -163,7 +163,6 @@ enum { MAINMENU_FILE_EXIT, IFMENU_FILE_PROPERTIES, - PHMMENU_FILE_PROPERTIES, PJMENU_FILE_PROPERTIES, PJMENU_RECONSTRUCT_FBP, @@ -179,6 +178,7 @@ enum { IFMENU_VIEW_SCALE_AUTO, IFMENU_VIEW_SCALE_MINMAX, + IFMENU_VIEW_SCALE_FULL, IFMENU_COMPARE_IMAGES, IFMENU_COMPARE_ROW, @@ -203,11 +203,13 @@ enum { IFMENU_FILTER_SHUFFLENATURALTOFOURIERORDER, IFMENU_FILTER_SHUFFLEFOURIERTONATURALORDER, + PHMMENU_FILE_PROPERTIES, PHMMENU_PROCESS_RASTERIZE, PHMMENU_PROCESS_PROJECTIONS, PLOTMENU_VIEW_SCALE_MINMAX, PLOTMENU_VIEW_SCALE_AUTO, + PLOTMENU_VIEW_SCALE_FULL, MAINMENU_WINDOW_BASE, }; diff --git a/src/dialogs.cpp b/src/dialogs.cpp index f0ece54..b5b32a5 100644 --- a/src/dialogs.cpp +++ b/src/dialogs.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: dialogs.cpp,v 1.24 2001/01/02 16:02:13 kevin Exp $ +** $Id: dialogs.cpp,v 1.25 2001/01/04 21:28:41 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 @@ -84,7 +84,7 @@ StringValueAndTitleListBox::StringValueAndTitleListBox (wxDialog* pParent, int n }; const char* -StringValueAndTitleListBox::getSelectionStringValue (void) const +StringValueAndTitleListBox::getSelectionStringValue () const { return m_ppszValues[GetSelection()]; } @@ -125,7 +125,7 @@ DialogGetPhantom::DialogGetPhantom (wxFrame* pParent, int iDefaultPhantom) } const char* -DialogGetPhantom::getPhantom(void) +DialogGetPhantom::getPhantom() { return m_pListBoxPhantom->getSelectionStringValue(); } @@ -183,7 +183,7 @@ DialogGetComparisonImage::DialogGetComparisonImage (wxFrame* pParent, const char } ImageFileDocument* -DialogGetComparisonImage::getImageFileDocument(void) +DialogGetComparisonImage::getImageFileDocument() { return m_rVecIF[ m_pListBoxImageChoices->GetSelection() ]; } @@ -241,12 +241,12 @@ DialogGetMinMax::DialogGetMinMax (wxFrame* pParent, const char* const pszTitle, pTopSizer->SetSizeHints (this); } -DialogGetMinMax::~DialogGetMinMax (void) +DialogGetMinMax::~DialogGetMinMax () { } double -DialogGetMinMax::getMinimum (void) +DialogGetMinMax::getMinimum () { wxString strCtrl = m_pTextCtrlMin->GetValue(); double dValue; @@ -257,7 +257,7 @@ DialogGetMinMax::getMinimum (void) } double -DialogGetMinMax::getMaximum (void) +DialogGetMinMax::getMaximum () { wxString strCtrl = m_pTextCtrlMax->GetValue(); double dValue; @@ -401,13 +401,13 @@ DialogGetRasterParameters::DialogGetRasterParameters (wxFrame* pParent, int iDef pTopSizer->SetSizeHints (this); } -DialogGetRasterParameters::~DialogGetRasterParameters (void) +DialogGetRasterParameters::~DialogGetRasterParameters () { } unsigned int -DialogGetRasterParameters::getXSize (void) +DialogGetRasterParameters::getXSize () { wxString strCtrl = m_pTextCtrlXSize->GetValue(); unsigned long lValue; @@ -418,7 +418,7 @@ DialogGetRasterParameters::getXSize (void) } unsigned int -DialogGetRasterParameters::getYSize (void) +DialogGetRasterParameters::getYSize () { wxString strCtrl = m_pTextCtrlYSize->GetValue(); unsigned long lValue; @@ -430,7 +430,7 @@ DialogGetRasterParameters::getYSize (void) unsigned int -DialogGetRasterParameters::getNSamples (void) +DialogGetRasterParameters::getNSamples () { wxString strCtrl = m_pTextCtrlNSamples->GetValue(); unsigned long lValue; @@ -529,13 +529,13 @@ DialogGetProjectionParameters::DialogGetProjectionParameters (wxFrame* pParent, pTopSizer->SetSizeHints (this); } -DialogGetProjectionParameters::~DialogGetProjectionParameters (void) +DialogGetProjectionParameters::~DialogGetProjectionParameters () { } unsigned int -DialogGetProjectionParameters::getNDet (void) +DialogGetProjectionParameters::getNDet () { wxString strCtrl = m_pTextCtrlNDet->GetValue(); unsigned long lValue; @@ -546,7 +546,7 @@ DialogGetProjectionParameters::getNDet (void) } unsigned int -DialogGetProjectionParameters::getNView (void) +DialogGetProjectionParameters::getNView () { wxString strCtrl = m_pTextCtrlNView->GetValue(); unsigned long lValue; @@ -558,7 +558,7 @@ DialogGetProjectionParameters::getNView (void) unsigned int -DialogGetProjectionParameters::getNSamples (void) +DialogGetProjectionParameters::getNSamples () { wxString strCtrl = m_pTextCtrlNSamples->GetValue(); unsigned long lValue; @@ -569,7 +569,7 @@ DialogGetProjectionParameters::getNSamples (void) } double -DialogGetProjectionParameters::getRotAngle (void) +DialogGetProjectionParameters::getRotAngle () { wxString strCtrl = m_pTextCtrlRotAngle->GetValue(); double dValue; @@ -580,7 +580,7 @@ DialogGetProjectionParameters::getRotAngle (void) } double -DialogGetProjectionParameters::getFocalLengthRatio (void) +DialogGetProjectionParameters::getFocalLengthRatio () { wxString strCtrl = m_pTextCtrlFocalLength->GetValue(); double dValue; @@ -591,7 +591,7 @@ DialogGetProjectionParameters::getFocalLengthRatio (void) } double -DialogGetProjectionParameters::getFieldOfViewRatio (void) +DialogGetProjectionParameters::getFieldOfViewRatio () { wxString strCtrl = m_pTextCtrlFieldOfView->GetValue(); double dValue; @@ -602,13 +602,13 @@ DialogGetProjectionParameters::getFieldOfViewRatio (void) } const char* -DialogGetProjectionParameters::getGeometry (void) +DialogGetProjectionParameters::getGeometry () { return m_pListBoxGeometry->getSelectionStringValue(); } int -DialogGetProjectionParameters::getTrace (void) +DialogGetProjectionParameters::getTrace () { return Trace::convertTraceNameToID(m_pListBoxTrace->getSelectionStringValue()); } @@ -713,13 +713,13 @@ DialogGetReconstructionParameters::DialogGetReconstructionParameters (wxFrame* p pTopSizer->SetSizeHints (this); } -DialogGetReconstructionParameters::~DialogGetReconstructionParameters (void) +DialogGetReconstructionParameters::~DialogGetReconstructionParameters () { } unsigned int -DialogGetReconstructionParameters::getXSize (void) +DialogGetReconstructionParameters::getXSize () { wxString strCtrl = m_pTextCtrlXSize->GetValue(); unsigned long lValue; @@ -730,7 +730,7 @@ DialogGetReconstructionParameters::getXSize (void) } unsigned int -DialogGetReconstructionParameters::getYSize (void) +DialogGetReconstructionParameters::getYSize () { wxString strCtrl = m_pTextCtrlYSize->GetValue(); unsigned long lValue; @@ -741,7 +741,7 @@ DialogGetReconstructionParameters::getYSize (void) } unsigned int -DialogGetReconstructionParameters::getZeropad (void) +DialogGetReconstructionParameters::getZeropad () { wxString strCtrl = m_pTextCtrlZeropad->GetValue(); unsigned long lValue; @@ -753,7 +753,7 @@ DialogGetReconstructionParameters::getZeropad (void) unsigned int -DialogGetReconstructionParameters::getInterpParam (void) +DialogGetReconstructionParameters::getInterpParam () { wxString strCtrl = m_pTextCtrlInterpParam->GetValue(); unsigned long lValue; @@ -764,7 +764,7 @@ DialogGetReconstructionParameters::getInterpParam (void) } double -DialogGetReconstructionParameters::getFilterParam (void) +DialogGetReconstructionParameters::getFilterParam () { wxString strCtrl = m_pTextCtrlFilterParam->GetValue(); double dValue; @@ -775,25 +775,25 @@ DialogGetReconstructionParameters::getFilterParam (void) } const char* -DialogGetReconstructionParameters::getFilterName (void) +DialogGetReconstructionParameters::getFilterName () { return m_pListBoxFilter->getSelectionStringValue(); } const char* -DialogGetReconstructionParameters::getFilterMethodName (void) +DialogGetReconstructionParameters::getFilterMethodName () { return m_pListBoxFilterMethod->getSelectionStringValue(); } const char* -DialogGetReconstructionParameters::getInterpName (void) +DialogGetReconstructionParameters::getInterpName () { return m_pListBoxInterp->getSelectionStringValue(); } int -DialogGetReconstructionParameters::getTrace (void) +DialogGetReconstructionParameters::getTrace () { int iTrace = 0; if (strcmp("full", m_pListBoxTrace->getSelectionStringValue()) == 0) @@ -802,13 +802,13 @@ DialogGetReconstructionParameters::getTrace (void) } const char* -DialogGetReconstructionParameters::getBackprojectName (void) +DialogGetReconstructionParameters::getBackprojectName () { return m_pListBoxBackproject->getSelectionStringValue(); } const char* -DialogGetReconstructionParameters::getFilterGenerationName (void) +DialogGetReconstructionParameters::getFilterGenerationName () { return m_pListBoxFilterGeneration->getSelectionStringValue(); } @@ -821,6 +821,7 @@ DialogGetReconstructionParameters::getFilterGenerationName (void) ///////////////////////////////////////////////////////////////////// + DialogGetFilterParameters::DialogGetFilterParameters (wxFrame* pParent, int iDefaultXSize, int iDefaultYSize, int iDefaultFilterID, double dDefaultFilterParam, double dDefaultBandwidth, int iDefaultDomainID, double dDefaultInputScale, double dDefaultOutputScale) : wxDialog (pParent, -1, "Set Filter Parameters", wxDefaultPosition, wxDefaultSize, wxDEFAULT_DIALOG_STYLE | wxCAPTION) { @@ -892,13 +893,13 @@ DialogGetFilterParameters::DialogGetFilterParameters (wxFrame* pParent, int iDef pTopSizer->SetSizeHints (this); } -DialogGetFilterParameters::~DialogGetFilterParameters (void) +DialogGetFilterParameters::~DialogGetFilterParameters () { } unsigned int -DialogGetFilterParameters::getXSize (void) +DialogGetFilterParameters::getXSize () { wxString strCtrl = m_pTextCtrlXSize->GetValue(); unsigned long lValue; @@ -909,7 +910,7 @@ DialogGetFilterParameters::getXSize (void) } unsigned int -DialogGetFilterParameters::getYSize (void) +DialogGetFilterParameters::getYSize () { wxString strCtrl = m_pTextCtrlYSize->GetValue(); unsigned long lValue; @@ -920,7 +921,7 @@ DialogGetFilterParameters::getYSize (void) } double -DialogGetFilterParameters::getBandwidth (void) +DialogGetFilterParameters::getBandwidth () { wxString strCtrl = m_pTextCtrlBandwidth->GetValue(); double dValue; @@ -931,7 +932,7 @@ DialogGetFilterParameters::getBandwidth (void) } double -DialogGetFilterParameters::getFilterParam (void) +DialogGetFilterParameters::getFilterParam () { wxString strCtrl = m_pTextCtrlFilterParam->GetValue(); double dValue; @@ -942,7 +943,7 @@ DialogGetFilterParameters::getFilterParam (void) } double -DialogGetFilterParameters::getInputScale (void) +DialogGetFilterParameters::getInputScale () { wxString strCtrl = m_pTextCtrlInputScale->GetValue(); double dValue; @@ -953,7 +954,7 @@ DialogGetFilterParameters::getInputScale (void) } double -DialogGetFilterParameters::getOutputScale (void) +DialogGetFilterParameters::getOutputScale () { wxString strCtrl = m_pTextCtrlOutputScale->GetValue(); double dValue; @@ -964,13 +965,13 @@ DialogGetFilterParameters::getOutputScale (void) } const char* -DialogGetFilterParameters::getFilterName (void) +DialogGetFilterParameters::getFilterName () { return m_pListBoxFilter->getSelectionStringValue(); } const char* -DialogGetFilterParameters::getDomainName (void) +DialogGetFilterParameters::getDomainName () { return m_pListBoxDomain->getSelectionStringValue(); } @@ -1011,7 +1012,7 @@ DialogExportParameters::DialogExportParameters (wxFrame* pParent, int iDefaultFo } const char* -DialogExportParameters::getFormatName(void) +DialogExportParameters::getFormatName() { return m_pListBoxFormat->getSelectionStringValue(); } @@ -1063,12 +1064,12 @@ DialogGetXYSize::DialogGetXYSize (wxFrame* pParent, const char* const pszTitle, pTopSizer->SetSizeHints (this); } -DialogGetXYSize::~DialogGetXYSize (void) +DialogGetXYSize::~DialogGetXYSize () { } unsigned int -DialogGetXYSize::getXSize (void) +DialogGetXYSize::getXSize () { wxString strCtrl = m_pTextCtrlXSize->GetValue(); long lValue; @@ -1079,7 +1080,7 @@ DialogGetXYSize::getXSize (void) } unsigned int -DialogGetXYSize::getYSize (void) +DialogGetXYSize::getYSize () { wxString strCtrl = m_pTextCtrlYSize->GetValue(); long lValue; @@ -1089,3 +1090,115 @@ DialogGetXYSize::getYSize (void) return (m_iDefaultYSize); } + + +///////////////////////////////////////////////////////////////////// +// CLASS IDENTIFICATION +// +// DialogGetConvertPolarParameters +///////////////////////////////////////////////////////////////////// + +DialogGetConvertPolarParameters::DialogGetConvertPolarParameters (wxFrame* pParent, const char* const pszTitle, + int iDefaultXSize, int iDefaultYSize, int iDefaultInterpolationID, int iDefaultZeropad) +: wxDialog (pParent, -1, pszTitle, wxDefaultPosition, wxDefaultSize, wxDEFAULT_DIALOG_STYLE | wxCAPTION) +{ + m_iDefaultXSize = iDefaultXSize; + m_iDefaultYSize = iDefaultYSize; + m_iDefaultZeropad = iDefaultZeropad; + + wxBoxSizer* pTopSizer = new wxBoxSizer (wxVERTICAL); + + pTopSizer->Add (new wxStaticText (this, -1, pszTitle), 0, wxALIGN_CENTER | wxTOP | wxLEFT | wxRIGHT, 5); + + pTopSizer->Add (new wxStaticLine (this, -1, wxDefaultPosition, wxSize(3,3), wxHORIZONTAL), 0, wxEXPAND | wxALL, 5); + + std::ostringstream os; + os << iDefaultXSize; + m_pTextCtrlXSize = new wxTextCtrl (this, -1, os.str().c_str(), wxDefaultPosition, wxSize(100, 25), 0); + std::ostringstream osYSize; + osYSize << iDefaultYSize; + m_pTextCtrlYSize = new wxTextCtrl (this, -1, osYSize.str().c_str(), wxDefaultPosition, wxSize(100, 25), 0); + + wxFlexGridSizer* pGridSizer = new wxFlexGridSizer (2); + + m_pListBoxInterpolation = new StringValueAndTitleListBox (this, Projections::getInterpCount(), Projections::getInterpTitleArray(), Projections::getInterpNameArray()); + m_pListBoxInterpolation->SetSelection (iDefaultInterpolationID); + pGridSizer->Add (new wxStaticText (this, -1, "Interpolation"), 0, wxALIGN_CENTER_VERTICAL | wxALIGN_RIGHT | wxALL, 5); + pGridSizer->Add (m_pListBoxInterpolation, 0, wxALL | wxALIGN_LEFT | wxEXPAND); + + pGridSizer->Add (new wxStaticText (this, -1, "X Size"), 0, wxALIGN_RIGHT | wxALIGN_CENTER_VERTICAL); + pGridSizer->Add (m_pTextCtrlXSize, 0, wxALIGN_LEFT | wxALIGN_CENTER_VERTICAL); + pGridSizer->Add (new wxStaticText (this, -1, "Y Size"), 0, wxALIGN_RIGHT | wxALIGN_CENTER_VERTICAL); + pGridSizer->Add (m_pTextCtrlYSize, 0, wxALIGN_LEFT | wxALIGN_CENTER_VERTICAL); + if (iDefaultZeropad >= 0) { + std::ostringstream osZeropad; + osZeropad << iDefaultZeropad; + m_pTextCtrlZeropad = new wxTextCtrl (this, -1, osZeropad.str().c_str(), wxDefaultPosition, wxSize(100, 25), 0); + pGridSizer->Add (new wxStaticText (this, -1, "Zeropad"), 0, wxALIGN_RIGHT | wxALIGN_CENTER_VERTICAL); + pGridSizer->Add (m_pTextCtrlZeropad, 0, wxALIGN_LEFT | wxALIGN_CENTER_VERTICAL); + } + + pTopSizer->Add (pGridSizer, 1, wxALL, 3); + + pTopSizer->Add (new wxStaticLine (this, -1, wxDefaultPosition, wxSize(3,3), wxHORIZONTAL), 0, wxEXPAND | wxALL, 5); + + wxBoxSizer* pButtonSizer = new wxBoxSizer (wxHORIZONTAL); + wxButton* pButtonOk = new wxButton (this, wxID_OK, "Okay"); + wxButton* pButtonCancel = new wxButton (this, wxID_CANCEL, "Cancel"); + pButtonSizer->Add (pButtonOk, 0, wxEXPAND | wxALL, 10); + pButtonSizer->Add (pButtonCancel, 0, wxEXPAND | wxALL, 10); + + pTopSizer->Add (pButtonSizer, 0, wxALIGN_CENTER); + + SetAutoLayout (true); + SetSizer (pTopSizer); + pTopSizer->Layout(); + pTopSizer->Fit (this); + pTopSizer->SetSizeHints (this); +} + + +DialogGetConvertPolarParameters::~DialogGetConvertPolarParameters () +{ +} + + +unsigned int +DialogGetConvertPolarParameters::getXSize () +{ + wxString strCtrl = m_pTextCtrlXSize->GetValue(); + unsigned long lValue; + if (strCtrl.ToULong (&lValue)) + return lValue; + else + return (m_iDefaultXSize); +} + +unsigned int +DialogGetConvertPolarParameters::getYSize () +{ + wxString strCtrl = m_pTextCtrlYSize->GetValue(); + unsigned long lValue; + if (strCtrl.ToULong (&lValue)) + return lValue; + else + return (m_iDefaultYSize); +} + +unsigned int +DialogGetConvertPolarParameters::getZeropad () +{ + wxString strCtrl = m_pTextCtrlZeropad->GetValue(); + unsigned long lValue; + if (strCtrl.ToULong (&lValue)) + return lValue; + else + return (m_iDefaultZeropad); +} + +const char* +DialogGetConvertPolarParameters::getInterpolationName () +{ + return m_pListBoxInterpolation->getSelectionStringValue(); +} + diff --git a/src/dialogs.h b/src/dialogs.h index 0553b1a..d16981c 100644 --- a/src/dialogs.h +++ b/src/dialogs.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: dialogs.h,v 1.19 2001/01/03 22:01:50 kevin Exp $ +** $Id: dialogs.h,v 1.20 2001/01/04 21:28:41 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 @@ -36,6 +36,7 @@ #include "phantom.h" #include "procsignal.h" #include "filter.h" +#include "projections.h" // CLASS StringValueAndTitleListBox @@ -291,5 +292,31 @@ class DialogGetXYSize : public wxDialog }; +class DialogGetConvertPolarParameters : public wxDialog +{ + public: + DialogGetConvertPolarParameters (wxFrame* pParent, const char* const pszTitle, int iDefaultXSize = 0, + int iDefaultYSize = 0, int iDefaultInterpolationID = Projections::POLAR_INTERP_BILINEAR, + int iDefaultZeropad = 1); + virtual ~DialogGetConvertPolarParameters (); + + unsigned int getXSize(); + unsigned int getYSize(); + const char* getInterpolationName(); + unsigned int getZeropad(); + + private: + wxTextCtrl* m_pTextCtrlXSize; + wxTextCtrl* m_pTextCtrlYSize; + wxTextCtrl* m_pTextCtrlZeropad; + + StringValueAndTitleListBox* m_pListBoxInterpolation; + + int m_iDefaultXSize; + int m_iDefaultYSize; + int m_iDefaultZeropad; +}; + + #endif diff --git a/src/dlgprojections.cpp b/src/dlgprojections.cpp index fb2f8cd..1d8c5bf 100644 --- a/src/dlgprojections.cpp +++ b/src/dlgprojections.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: dlgprojections.cpp,v 1.15 2001/01/02 16:02:13 kevin Exp $ +** $Id: dlgprojections.cpp,v 1.16 2001/01/04 21:28:41 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 @@ -68,7 +68,8 @@ IMPLEMENT_CLASS(ProjectionsDialog, wxDialog) ProjectionsDialog::ProjectionsDialog (Scanner& rScanner, Projections& rProj, const Phantom& rPhantom, const int iTrace, wxWindow *parent) -: wxDialog(parent, -1, "Collect Projections", wxDefaultPosition), m_rScanner(rScanner), m_rProjections(rProj), m_rPhantom(rPhantom), m_pSGPDriver(NULL), m_pSGP(NULL), m_iTrace(iTrace), m_pDC(NULL), m_btnAbort(0), m_btnPause(0), m_btnStep(0) +: wxDialog(parent, -1, "Collect Projections", wxDefaultPosition), m_rScanner(rScanner), m_rProjections(rProj), m_rPhantom(rPhantom), + m_pSGPDriver(NULL), m_pSGP(NULL), m_iTrace(iTrace), m_pDC(NULL), m_btnAbort(0), m_btnPause(0), m_btnStep(0) { m_state = Continue; m_iLastView = -1; @@ -135,7 +136,7 @@ ProjectionsDialog::ProjectionsDialog (Scanner& rScanner, Projections& rProj, con wxYield(); // Update the display - m_pSGPDriver->idWX()->SetFont(*wxSWISS_FONT); + m_pSGP->setTextPointSize(10); #ifdef __WXMAC__ MacUpdateImmediately(); #endif @@ -228,7 +229,6 @@ ProjectionsDialog::OnPause (wxCommandEvent& event) } else if (m_state == Continue) { m_memoryDC.SelectObject (m_bitmap); // in memoryDC m_pSGP->setDC (&m_memoryDC); - m_memoryDC.SetFont (*wxSWISS_FONT); showView (m_iLastView); m_state = Paused; m_btnPause->SetLabel (wxString("Resume")); @@ -248,7 +248,6 @@ ProjectionsDialog::OnStep (wxCommandEvent& event) } else if (m_state == Continue) { m_memoryDC.SelectObject (m_bitmap); // in memoryDC m_pSGP->setDC (&m_memoryDC); - m_memoryDC.SetFont (*wxSWISS_FONT); showView (m_iLastView); // m_rScanner.collectProjections (m_rProjections, m_rPhantom, m_iLastView, 1, true, m_iTrace, m_pSGP); m_state = Paused; @@ -259,7 +258,6 @@ ProjectionsDialog::OnStep (wxCommandEvent& event) } else if (m_state == Paused) { m_memoryDC.SelectObject (m_bitmap); // in memoryDC m_pSGP->setDC (&m_memoryDC); - m_memoryDC.SetFont (*wxSWISS_FONT); projectView (m_iLastView + 1); m_pSGP->setDC (m_pDC); m_memoryDC.SelectObject(wxNullBitmap); diff --git a/src/views.cpp b/src/views.cpp index f7d8a1f..aba89cc 100644 --- a/src/views.cpp +++ b/src/views.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: views.cpp,v 1.50 2001/01/03 22:00:46 kevin Exp $ +** $Id: views.cpp,v 1.51 2001/01/04 21:28:41 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 @@ -170,6 +170,7 @@ EVT_MENU(IFMENU_FILE_EXPORT, ImageFileView::OnExport) EVT_MENU(IFMENU_FILE_PROPERTIES, ImageFileView::OnProperties) EVT_MENU(IFMENU_VIEW_SCALE_MINMAX, ImageFileView::OnScaleMinMax) EVT_MENU(IFMENU_VIEW_SCALE_AUTO, ImageFileView::OnScaleAuto) +EVT_MENU(IFMENU_VIEW_SCALE_FULL, ImageFileView::OnScaleFull) EVT_MENU(IFMENU_COMPARE_IMAGES, ImageFileView::OnCompare) EVT_MENU(IFMENU_COMPARE_ROW, ImageFileView::OnCompareRow) EVT_MENU(IFMENU_COMPARE_COL, ImageFileView::OnCompareCol) @@ -286,6 +287,16 @@ ImageFileView::OnScaleMinMax (wxCommandEvent& event) } } +void +ImageFileView::OnScaleFull (wxCommandEvent& event) +{ + if (m_bMinSpecified || m_bMaxSpecified) { + m_bMinSpecified = false; + m_bMaxSpecified = false; + OnUpdate (this, NULL); + } +} + void ImageFileView::OnCompare (wxCommandEvent& event) { @@ -304,7 +315,7 @@ ImageFileView::OnCompare (wxCommandEvent& event) std::ostringstream os; double min, max, mean, mode, median, stddev; rIF.statistics (min, max, mean, mode, median, stddev); - os << rIF.getFilename() << ": minimum=" << min << ", maximum=" << max << ", mean=" << mean << ", mode=" << mode << ", median=" << median << ", stddev=" << stddev << "\n"; + os << GetFrame()->GetTitle().c_str() << ": minimum=" << min << ", maximum=" << max << ", mean=" << mean << ", mode=" << mode << ", median=" << median << ", stddev=" << stddev << "\n"; rCompareIF.statistics (min, max, mean, mode, median, stddev); os << pCompareDoc->GetFirstView()->GetFrame()->GetTitle().c_str() << ": minimum=" << min << ", maximum=" << max << ", mean=" << mean << ", mode=" << mode << ", median=" << median << ", stddev=" << stddev << "\n"; os << "\n"; @@ -326,7 +337,7 @@ ImageFileView::OnCompare (wxCommandEvent& event) return; } - wxString s = GetDocument()->GetFirstView()->GetFrame()->GetTitle() + ": "; + wxString s = GetFrame()->GetTitle() + ": "; differenceImage.labelsCopy (rIF, s.c_str()); s = pCompareDoc->GetFirstView()->GetFrame()->GetTitle() + ": "; differenceImage.labelsCopy (rCompareIF, s.c_str()); @@ -561,7 +572,7 @@ ImageFileView::OnDivide (wxCommandEvent& event) } -#ifdef HAVE_FFTW +#ifdef HAVE_FFT void ImageFileView::OnFFT (wxCommandEvent& event) { @@ -604,6 +615,7 @@ ImageFileView::OnFourier (wxCommandEvent& event) GetDocument()->Modify(TRUE); GetDocument()->UpdateAllViews(this); } + void ImageFileView::OnInverseFourier (wxCommandEvent& event) { @@ -717,6 +729,7 @@ ImageFileView::CreateChildFrame(wxDocument *doc, wxView *view) wxMenu *view_menu = new wxMenu; view_menu->Append(IFMENU_VIEW_SCALE_MINMAX, "Display Scale &Set..."); view_menu->Append(IFMENU_VIEW_SCALE_AUTO, "Display Scale &Auto..."); + view_menu->Append(IFMENU_VIEW_SCALE_FULL, "Display &Full Scale"); wxMenu* filter_menu = new wxMenu; filter_menu->Append (IFMENU_FILTER_INVERTVALUES, "&Invert Values"); @@ -725,7 +738,7 @@ ImageFileView::CreateChildFrame(wxDocument *doc, wxView *view) filter_menu->Append (IFMENU_FILTER_LOG, "&Log"); filter_menu->Append (IFMENU_FILTER_EXP, "&Exp"); filter_menu->AppendSeparator(); -#ifdef HAVE_FFTW +#ifdef HAVE_FFT filter_menu->Append (IFMENU_FILTER_FFT, "&FFT"); filter_menu->Append (IFMENU_FILTER_IFFT, "&IFFT"); filter_menu->Append (IFMENU_FILTER_FOURIER, "F&ourier"); @@ -1283,7 +1296,7 @@ PhantomView::PhantomView(void) m_iDefaultNDet = 367; m_iDefaultNView = 320; m_iDefaultNSample = 2; - m_dDefaultRotation = 2; + m_dDefaultRotation = 1; m_dDefaultFocalLength = 2; m_dDefaultFieldOfView = 1; m_iDefaultGeometry = Scanner::GEOMETRY_PARALLEL; @@ -1617,6 +1630,11 @@ ProjectionFileView::ProjectionFileView(void) m_iDefaultInterpolation = Backprojector::INTERP_LINEAR; m_iDefaultInterpParam = 1; m_iDefaultTrace = Trace::TRACE_NONE; + + m_iDefaultPolarNX = 256; + m_iDefaultPolarNY = 256; + m_iDefaultPolarInterpolation = Projections::POLAR_INTERP_BILINEAR; + m_iDefaultPolarZeropad = 1; } ProjectionFileView::~ProjectionFileView(void) @@ -1639,23 +1657,67 @@ void ProjectionFileView::OnConvertPolar (wxCommandEvent& event) { Projections& rProj = GetDocument()->getProjections(); + DialogGetConvertPolarParameters dialogPolar (m_frame, "Convert Polar", m_iDefaultPolarNX, m_iDefaultPolarNY, + m_iDefaultPolarInterpolation, -1); + if (dialogPolar.ShowModal() == wxID_OK) { + wxString strInterpolation (dialogPolar.getInterpolationName()); + m_iDefaultPolarNX = dialogPolar.getXSize(); + m_iDefaultPolarNY = dialogPolar.getYSize(); + ImageFileDocument* pPolarDoc = dynamic_cast(theApp->getDocManager()->CreateDocument("untitled.if", wxDOC_SILENT)); + ImageFile& rIF = pPolarDoc->getImageFile(); + if (! pPolarDoc) { + sys_error (ERR_SEVERE, "Unable to create image file"); + return; + } + rIF.setArraySize (m_iDefaultPolarNX, m_iDefaultPolarNY); + m_iDefaultPolarInterpolation = Projections::convertInterpNameToID (strInterpolation.c_str()); + rProj.convertPolar (rIF, m_iDefaultPolarInterpolation); + rIF.labelAdd (rProj.getLabel().getLabelString().c_str(), rProj.calcTime()); + std::ostringstream os; + os << "Convert projection file " << GetFrame()->GetTitle().c_str() << " to polar image: xSize=" + << m_iDefaultPolarNX << ", ySize=" << m_iDefaultPolarNY << ", interpolation=" + << strInterpolation.c_str(); + *theApp->getLog() << os.str().c_str() << "\n"; + rIF.labelAdd (os.str().c_str()); + if (theApp->getSetModifyNewDocs()) + pPolarDoc->Modify(true); + pPolarDoc->UpdateAllViews(); + pPolarDoc->GetFirstView()->OnUpdate (this, NULL); + } } void ProjectionFileView::OnConvertFFTPolar (wxCommandEvent& event) { Projections& rProj = GetDocument()->getProjections(); - wxMessageBox ("Polar conversion not yet implemented", "Unimplemented function"); -#if 0 - rProj.convertPolar (); - if (theApp->getSetModifyNewDocs()) - GetDocument()->Modify(true); - GetDocument()->UpdateAllViews(); -#ifndef HAVE_FFT - wxMessageBox ("FFT support has not been compiled into this version of CTSim", "Error"); -#endif -#endif -} + DialogGetConvertPolarParameters dialogPolar (m_frame, "Convert to FFT Polar", m_iDefaultPolarNX, m_iDefaultPolarNY, + m_iDefaultPolarInterpolation, m_iDefaultPolarZeropad); + if (dialogPolar.ShowModal() == wxID_OK) { + wxString strInterpolation (dialogPolar.getInterpolationName()); + m_iDefaultPolarNX = dialogPolar.getXSize(); + m_iDefaultPolarNY = dialogPolar.getYSize(); + m_iDefaultPolarZeropad = dialogPolar.getZeropad(); + ImageFileDocument* pPolarDoc = dynamic_cast(theApp->getDocManager()->CreateDocument("untitled.if", wxDOC_SILENT)); + ImageFile& rIF = pPolarDoc->getImageFile(); + if (! pPolarDoc) { + sys_error (ERR_SEVERE, "Unable to create image file"); + return; + } + rIF.setArraySize (m_iDefaultPolarNX, m_iDefaultPolarNY); + m_iDefaultPolarInterpolation = Projections::convertInterpNameToID (strInterpolation.c_str()); + rProj.convertFFTPolar (rIF, m_iDefaultPolarInterpolation, m_iDefaultPolarZeropad); + rIF.labelAdd (rProj.getLabel().getLabelString().c_str(), rProj.calcTime()); + std::ostringstream os; + os << "Convert projection file " << GetFrame()->GetTitle().c_str() << " to FFT polar image: xSize=" + << m_iDefaultPolarNX << ", ySize=" << m_iDefaultPolarNY << ", interpolation=" + << strInterpolation.c_str() << ", zeropad=" << m_iDefaultPolarZeropad; + *theApp->getLog() << os.str().c_str() << "\n"; + rIF.labelAdd (os.str().c_str()); + if (theApp->getSetModifyNewDocs()) + pPolarDoc->Modify(true); + pPolarDoc->UpdateAllViews(); + pPolarDoc->GetFirstView()->OnUpdate (this, NULL); + }} void ProjectionFileView::OnReconstructFourier (wxCommandEvent& event) @@ -1947,6 +2009,7 @@ BEGIN_EVENT_TABLE(PlotFileView, wxView) EVT_MENU(PJMENU_FILE_PROPERTIES, PlotFileView::OnProperties) EVT_MENU(PLOTMENU_VIEW_SCALE_MINMAX, PlotFileView::OnScaleMinMax) EVT_MENU(PLOTMENU_VIEW_SCALE_AUTO, PlotFileView::OnScaleAuto) +EVT_MENU(PLOTMENU_VIEW_SCALE_FULL, PlotFileView::OnScaleFull) END_EVENT_TABLE() PlotFileView::PlotFileView(void) @@ -2026,6 +2089,16 @@ PlotFileView::OnScaleMinMax (wxCommandEvent& event) } } +void +PlotFileView::OnScaleFull (wxCommandEvent& event) +{ + if (m_bMinSpecified || m_bMaxSpecified) { + m_bMinSpecified = false; + m_bMaxSpecified = false; + OnUpdate (this, NULL); + } +} + PlotFileCanvas* PlotFileView::CreateCanvas (wxView *view, wxFrame *parent) @@ -2067,6 +2140,7 @@ PlotFileView::CreateChildFrame(wxDocument *doc, wxView *view) wxMenu *view_menu = new wxMenu; view_menu->Append(PLOTMENU_VIEW_SCALE_MINMAX, "Display Scale &Set..."); view_menu->Append(PLOTMENU_VIEW_SCALE_AUTO, "Display Scale &Auto..."); + view_menu->Append(PLOTMENU_VIEW_SCALE_FULL, "Display &Full Scale"); wxMenu *help_menu = new wxMenu; help_menu->Append(MAINMENU_HELP_CONTENTS, "&Contents"); diff --git a/src/views.h b/src/views.h index 4463689..7ca2e62 100644 --- a/src/views.h +++ b/src/views.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: views.h,v 1.24 2001/01/03 22:00:46 kevin Exp $ +** $Id: views.h,v 1.25 2001/01/04 21:28:41 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 @@ -87,7 +87,7 @@ public: void OnShuffleNaturalToFourierOrder (wxCommandEvent& event); void OnShuffleFourierToNaturalOrder (wxCommandEvent& event); -#ifdef HAVE_FFTW +#ifdef HAVE_FFT void OnFFT (wxCommandEvent& event); void OnIFFT (wxCommandEvent& event); #endif @@ -97,6 +97,7 @@ public: void OnScaleAuto (wxCommandEvent& event); void OnScaleMinMax (wxCommandEvent& event); + void OnScaleFull (wxCommandEvent& event); void OnPlotRow (wxCommandEvent& event); void OnPlotCol (wxCommandEvent& event); void OnPlotHistogram (wxCommandEvent& event); @@ -159,6 +160,11 @@ private: int m_iDefaultBackprojector; int m_iDefaultTrace; + int m_iDefaultPolarNX; + int m_iDefaultPolarNY; + int m_iDefaultPolarInterpolation; + int m_iDefaultPolarZeropad; + public: ProjectionFileView(void); virtual ~ProjectionFileView(void); @@ -266,9 +272,11 @@ public: void OnDraw(wxDC* dc); void OnUpdate(wxView *sender, wxObject *hint = NULL); bool OnClose (bool deleteWindow = true); + void OnProperties (wxCommandEvent& event); - void OnScaleAuto (wxCommandEvent& event); void OnScaleMinMax (wxCommandEvent& event); + void OnScaleAuto (wxCommandEvent& event); + void OnScaleFull (wxCommandEvent& event); wxFrame* getFrame () { return m_frame; } -- 2.34.1