From 5c6b29ab4885308cc3381af6e0a68f4804956d2e Mon Sep 17 00:00:00 2001 From: "Kevin M. Rosenberg" Date: Fri, 29 Dec 2000 15:45:06 +0000 Subject: [PATCH] r318: *** empty log message *** --- ChangeLog | 15 +- TODO | 8 +- include/ctsupport.h | 4 +- include/ezplot.h | 8 +- include/imagefile.h | 20 +- include/pol.h | 23 +- include/procsignal.h | 14 +- libctgraphics/ezplot.cpp | 247 ++++++------- libctgraphics/ezset.cpp | 4 +- libctgraphics/pol.cpp | 41 ++- libctsim/imagefile.cpp | 242 +++++++++++- libctsim/procsignal.cpp | 748 ++++++++++++++++++++------------------ libctsim/reconstruct.cpp | 12 +- libctsupport/syserror.cpp | 155 ++++---- msvc/ctsim/ctsim.plg | 29 +- src/ctsim.h | 9 +- src/dlgprojections.cpp | 12 +- src/docs.cpp | 19 +- src/views.cpp | 205 ++++++++--- src/views.h | 10 +- tools/if1.cpp | 66 ++-- tools/if2.cpp | 12 +- 22 files changed, 1153 insertions(+), 750 deletions(-) diff --git a/ChangeLog b/ChangeLog index 894573e..91c26bc 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,13 +1,6 @@ TODO Read PlotFile's. - Add help button onto dialog's for reconstruction & projections. - Consider use of wxWindows help file classes. - - Consider changing POL to to read tag tokens. - - Consider being able to write Phantom files as text and view as - text in ctsim. 3.0alpha1 - Released 12/30/00 @@ -21,6 +14,9 @@ TODO * ctsim: Added row and column plot comparisons between two image files. + * ctsim: Added "Process" menu to image file display with math + functions. Added 2-dimensional inverse Fourier to math functions. + * mathfuncs.cpp: Reworked statistics algorithm to share between imagefile and plotfile classes. @@ -29,7 +25,7 @@ TODO if1.cpp and if2.cpp to imagefile.cpp * if1: Updated to handle error conditions, such as sqrt of a - negative number. + negative number. Converted to use new ImageFile math functions. * if2: Updated to output plot files and use new ImageFile class math functions @@ -40,7 +36,8 @@ TODO individual curves. Improved display of labels and ticks. Updated to use POL class member variable. Updated to more C++ conventions. - * pol: converted to C++ class + * pol: converted to C++ class. Extracted HashTable to separate + class. * sgp: Added linestyle settings diff --git a/TODO b/TODO index ad50170..69aa217 100644 --- a/TODO +++ b/TODO @@ -1,4 +1,10 @@ -Convert pol to C++ object. +Add help button onto dialog's for reconstruction & projections. +Consider use of wxWindows help file classes. + +Consider changing POL to to read tag tokens. + +Consider being able to write Phantom files as text and view as +text in ctsim. FFT filtering requires FFTW library. Add a basic FFT routine. diff --git a/include/ctsupport.h b/include/ctsupport.h index 55c25d5..9a9ca83 100644 --- a/include/ctsupport.h +++ b/include/ctsupport.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: ctsupport.h,v 1.17 2000/12/27 20:09:19 kevin Exp $ +** $Id: ctsupport.h,v 1.18 2000/12/29 15:45:06 kevin Exp $ ** ** ** This program is free software; you can redistribute it and/or modify @@ -148,7 +148,7 @@ char *str_upper(char *str); /* syserror.cpp */ void sys_error(int severity, const char *msg, ...); -void sys_verror(int severity, const char *msg, va_list arg); +void sys_verror (std::string& strOutput, int severity, const char *msg, va_list arg); void sys_error_level(int severity); // Math Section diff --git a/include/ezplot.h b/include/ezplot.h index b2512aa..f6dbf77 100644 --- a/include/ezplot.h +++ b/include/ezplot.h @@ -7,7 +7,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: ezplot.h,v 1.22 2000/12/27 03:16:02 kevin Exp $ +** $Id: ezplot.h,v 1.23 2000/12/29 15:45:06 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 @@ -319,7 +319,7 @@ private: void make_numfmt(char *fmtstr, int *fldwid, int *nfrac, double min, double max, int nint); int axis_scale (double min, double max, int nint, double *minp, double *maxp, int *nintp); - SGP& rSGP; + SGP* m_pSGP; POL m_pol; void clearCurves (); @@ -338,7 +338,7 @@ private: { return ygn_min + (y - ygw_min) * m_yWorldScale; } public: - EZPlot (SGP& sgp); + EZPlot (); ~EZPlot (); bool ezset (const std::string& command); @@ -350,7 +350,7 @@ private: void addCurve (const double* y, int n); void addCurve (const float* y, int n); - void plot (); + void plot (SGP* pSGP); }; #endif diff --git a/include/imagefile.h b/include/imagefile.h index a7e69d7..7c7daba 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.23 2000/12/22 04:18:00 kevin Exp $ +** $Id: imagefile.h,v 1.24 2000/12/29 15:45:06 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 @@ -120,25 +120,25 @@ class ImageFile : public ImageFileBase void filterResponse (const char* const domainName, double bw, const char* const filterName, double filt_param); void statistics (double& min, double& max, double& mean, double& mode, double& median, double& stddev) const; - void getMinMax (double& min, double& max) const; - void printStatistics (std::ostream& os) const; - bool comparativeStatistics (const ImageFile& imComp, double& d, double& r, double& e) const; - bool printComparativeStatistics (const ImageFile& imComp, std::ostream& os) const; - bool subtractImages (const ImageFile& rRHS, ImageFile& result) const; - + bool subtractImages (const ImageFile& rRHS, ImageFile& result) const; bool addImages (const ImageFile& rRHS, ImageFile& result) const; - bool multiplyImages (const ImageFile& rRHS, ImageFile& result) const; - bool divideImages (const ImageFile& rRHS, ImageFile& result) const; + bool invertPixelValues (ImageFile& result) const; + bool sqrt (ImageFile& result) const; + bool square (ImageFile& result) const; + bool log (ImageFile& result) const; + bool exp (ImageFile& result) const; + bool FFTMagnitude (ImageFile& result) const; + bool FFTPhase (ImageFile& result) const; + int display (void) const; - int displayScaling (const int scaleFactor, ImageFileValue pmin, ImageFileValue pmax) const; #if HAVE_PNG diff --git a/include/pol.h b/include/pol.h index 262cccb..808b380 100644 --- a/include/pol.h +++ b/include/pol.h @@ -2,7 +2,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: pol.h,v 1.8 2000/12/27 20:09:19 kevin Exp $ +** $Id: pol.h,v 1.9 2000/12/29 15:45:06 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 @@ -21,7 +21,7 @@ #define __H_POL #include "hashtable.h" - +#include class POL { @@ -167,6 +167,17 @@ private: PC_DUMP, }; + enum { + INPUT_STREAM = 1, + INPUT_FILE, + INPUT_STRING, + }; + + int m_iInputType; + std::istream* m_pInputStream; + FILE* m_pInputFile; + char* m_pszInputString; + enum { MAXFILE = 8, }; @@ -178,12 +189,8 @@ private: char inputline[MAXLINE]; /* current input line */ int lineptr; /* current position in inputline */ - enum { - BUFSIZE = 100, - }; - int bp; // pointer to next free position - int buf[BUFSIZE]; // pushed back input characters - + std::stack m_stackPushBackInput; + bool skipSingleToken (char term[]); int tok (struct token_st *token); void dumptok (struct token_st *token); diff --git a/include/procsignal.h b/include/procsignal.h index 7a3cb62..a160154 100644 --- a/include/procsignal.h +++ b/include/procsignal.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: procsignal.h,v 1.8 2000/12/16 06:12:47 kevin Exp $ +** $Id: procsignal.h,v 1.9 2000/12/29 15:45:06 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 @@ -57,6 +57,11 @@ class ProcessSignal { static const int FILTER_GENERATION_INVALID; static const int FILTER_GENERATION_DIRECT; static const int FILTER_GENERATION_INVERSE_FOURIER; + + enum { + FORWARD = -1, + BACKWARD = 1, + }; ProcessSignal (const char* szFilterName, const char* szFilterMethodName,double bw, double signalIncrement, int n, double param, const char* szDomainName, const char* szFilterGenerationName, const int zeropad = 0, const int preinterpolationFactor = 1, const int iTraceLevel = Trace::TRACE_NONE, int iGeometry = Scanner::GEOMETRY_PARALLEL, double dFocalLength = 1., SGP* pSGP = NULL); @@ -99,9 +104,10 @@ class ProcessSignal { static void finiteFourierTransform (const std::complex input[], double output[], const int n, const int direction); - static void shuffleNaturalToFourierOrder (double* pdVector, const int n); - - static void shuffleFourierToNaturalOrder (double* pdVector, const int n); + static void shuffleNaturalToFourierOrder (double* pdVector, const int n); + static void shuffleNaturalToFourierOrder (std::complex* pdVector, const int n); + static void shuffleFourierToNaturalOrder (double* pdVector, const int n); + static void shuffleFourierToNaturalOrder (std::complex* pdVector, const int n); private: std::string m_nameFilterMethod; diff --git a/libctgraphics/ezplot.cpp b/libctgraphics/ezplot.cpp index 72f7e51..e4537f0 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.25 2000/12/27 20:09:19 kevin Exp $ +** $Id: ezplot.cpp,v 1.26 2000/12/29 15:45:06 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 @@ -134,8 +134,7 @@ EZPlot::clearCurves () } -EZPlot::EZPlot (SGP& sgp) -: rSGP (sgp) +EZPlot::EZPlot () { initKeywords(); m_pol.addSkipWord ("please"); @@ -158,10 +157,8 @@ EZPlot::EZPlot (SGP& sgp) void EZPlot::initPlotSettings () { - charheight = rSGP.getCharHeight(); - charwidth = rSGP.getCharWidth(); - m_iCurrentCurve = -1; + m_pSGP = NULL; c_xlabel = ""; c_ylabel = ""; @@ -340,17 +337,19 @@ EZPlot::getLegend (unsigned int iCurve) const */ void -EZPlot::plot () +EZPlot::plot (SGP* pSGP) { if (m_vecCurves.size() <= 0) return; - - rSGP.setWindow (0., 0., 1., 1.); + + m_pSGP = pSGP; + + m_pSGP->setWindow (0., 0., 1., 1.); if (s_textsize == TRUE) - rSGP.setTextPointSize (v_textsize); - charheight = rSGP.getCharHeight(); - charwidth = rSGP.getCharWidth(); + m_pSGP->setTextPointSize (v_textsize); + charheight = m_pSGP->getCharHeight(); + charwidth = m_pSGP->getCharWidth(); const EZPlotCurve& firstCurve = *m_vecCurves[0]; double xmin = firstCurve.x[0]; // extent of curves in world coord @@ -457,6 +456,8 @@ EZPlot::plot () title_row = ya_max;; if (c_title.length() > 0) ya_max -= 2 * charheight; + else + ya_max -= 0.7 * charheight; // allow room for yaxis ticklabel // calculate legend box boundaries int max_leg = 0; // longest legend in characters @@ -493,8 +494,8 @@ EZPlot::plot () yl_min = yl_max - leg_height; - rSGP.setColor (clr_legend); - rSGP.drawRect (xl_min, yl_min, xl_max, yl_max); + m_pSGP->setColor (clr_legend); + m_pSGP->drawRect (xl_min, yl_min, xl_max, yl_max); int iLegend = 0; // current legend position for (iCurve = 0; iCurve < m_vecCurves.size(); iCurve++) { @@ -506,21 +507,21 @@ EZPlot::plot () double xmax = xl_max - 1.0 * charwidth; double y = yl_max - (2.0 + iLegend * 3) * charheight; - rSGP.moveAbs (xmin, y + 0.5 * charheight); - rSGP.drawText (pstrLegend->c_str()); - rSGP.setColor (getColor (iCurve)); + m_pSGP->moveAbs (xmin, y + 0.5 * charheight); + m_pSGP->drawText (pstrLegend->c_str()); + m_pSGP->setColor (getColor (iCurve)); int iLS = getLinestyle (iCurve); if (iLS != SGP::LS_NOLINE) { - rSGP.setLineStyle (iLS); - rSGP.moveAbs (xmin, y); - rSGP.lineAbs (xmax, y); + m_pSGP->setLineStyle (iLS); + m_pSGP->moveAbs (xmin, y); + m_pSGP->lineAbs (xmax, y); } int iSymbol = getSymbol (iCurve); if (iSymbol > 0) { double xinc = (xmax - xmin) / (5 - 1); - rSGP.setLineStyle (SGP::LS_SOLID); + m_pSGP->setLineStyle (SGP::LS_SOLID); for (int j = 0; j < 5; j++) { - rSGP.moveAbs (xmin + j * xinc, y); + m_pSGP->moveAbs (xmin + j * xinc, y); symbol(iSymbol, 0.5 * charwidth, 0.5 * charheight); } } @@ -541,13 +542,13 @@ EZPlot::plot () // Y-Label if (c_ylabel.length() > 0) { - rSGP.setTextAngle (HALFPI); - rSGP.setTextSize (1.5 * charheight); + m_pSGP->setTextAngle (HALFPI); + m_pSGP->setTextSize (1.5 * charheight); double xExtent, yExtent; - rSGP.getTextExtent (c_ylabel.c_str(), &xExtent, &yExtent); - rSGP.setTextSize (charheight); + m_pSGP->getTextExtent (c_ylabel.c_str(), &xExtent, &yExtent); + m_pSGP->setTextSize (charheight); xa_min += xExtent; - rSGP.setTextAngle (0.0); + m_pSGP->setTextAngle (0.0); } ylbl_col = xp_min; @@ -650,7 +651,7 @@ EZPlot::plot () // PLOT CURVES - rSGP.setLineStyle (SGP::LS_SOLID); + m_pSGP->setLineStyle (SGP::LS_SOLID); drawAxes(); // size of symbol in NDC @@ -663,12 +664,12 @@ EZPlot::plot () for (iCurve = 0; iCurve < m_vecCurves.size(); iCurve++) { const EZPlotCurve* const pCurve = m_vecCurves [iCurve]; - rSGP.setColor (getColor (iCurve)); + m_pSGP->setColor (getColor (iCurve)); int iSym = getSymbol (iCurve); int iLS = getLinestyle (iCurve); if (iLS != SGP::LS_NOLINE) { - rSGP.setLineStyle (iLS); + m_pSGP->setLineStyle (iLS); double x1 = convertWorldToNDC_X (pCurve->x[0]); double y1 = convertWorldToNDC_Y (pCurve->y[0]); for (int i = 1; i < pCurve->m_iPointCount; i++) { @@ -677,8 +678,8 @@ EZPlot::plot () double x2Clip = x2; double y2Clip = y2; if (clip_rect (x1, y1, x2Clip, y2Clip, clipRect)) { - rSGP.moveAbs (x1, y1); - rSGP.lineAbs (x2Clip, y2Clip); + m_pSGP->moveAbs (x1, y1); + m_pSGP->lineAbs (x2Clip, y2Clip); } x1 = x2; y1 = y2; @@ -686,17 +687,17 @@ EZPlot::plot () } if (iSym > 0) { int iSymFreq = getSymbolFreq (iCurve); - rSGP.setLineStyle (SGP::LS_SOLID); + m_pSGP->setLineStyle (SGP::LS_SOLID); double x = convertWorldToNDC_X (pCurve->x[0]); double y = convertWorldToNDC_Y (pCurve->y[0]); - rSGP.moveAbs (x, y); + m_pSGP->moveAbs (x, y); symbol (iSym, symwidth, symheight); for (int i = 1; i < pCurve->m_iPointCount; i++) if (i % iSymFreq == 0 || i == pCurve->m_iPointCount - 1) { x = convertWorldToNDC_X (pCurve->x[i]); y = convertWorldToNDC_Y (pCurve->y[i]); if (y >= ygn_min && y <= ygn_max) { - rSGP.moveAbs (x, y); + m_pSGP->moveAbs (x, y); symbol (iSym, symwidth, symheight); } } @@ -726,8 +727,8 @@ EZPlot::drawAxes() char str[256]; char *numstr; - rSGP.setTextSize (charheight); - rSGP.setTextColor (1, -1); + m_pSGP->setTextSize (charheight); + m_pSGP->setTextColor (1, -1); if (o_xticks == ABOVE) xticklen = 0.5 * charheight; @@ -741,21 +742,21 @@ EZPlot::drawAxes() if (c_title.length() > 0) { double wText, hText; - rSGP.setTextSize (charheight * 2.0); - rSGP.getTextExtent (c_title.c_str(), &wText, &hText); - rSGP.moveAbs (xa_min + (xa_max-xa_min)/2 - wText/2, title_row); - rSGP.setTextColor (clr_title, -1); - rSGP.drawText (c_title); - rSGP.setTextSize (charheight); + m_pSGP->setTextSize (charheight * 2.0); + m_pSGP->getTextExtent (c_title.c_str(), &wText, &hText); + m_pSGP->moveAbs (xa_min + (xa_max-xa_min)/2 - wText/2, title_row); + m_pSGP->setTextColor (clr_title, -1); + m_pSGP->drawText (c_title); + m_pSGP->setTextSize (charheight); } if (o_grid == TRUE || o_box == TRUE) { - rSGP.setColor (clr_axis); - rSGP.moveAbs (xa_min, ya_min); - rSGP.lineAbs (xa_max, ya_min); - rSGP.lineAbs (xa_max, ya_max); - rSGP.lineAbs (xa_min, ya_max); - rSGP.lineAbs (xa_min, ya_min); + m_pSGP->setColor (clr_axis); + m_pSGP->moveAbs (xa_min, ya_min); + m_pSGP->lineAbs (xa_max, ya_min); + m_pSGP->lineAbs (xa_max, ya_max); + m_pSGP->lineAbs (xa_min, ya_max); + m_pSGP->lineAbs (xa_min, ya_min); } // calculate position of axes @@ -781,44 +782,44 @@ EZPlot::drawAxes() if (o_xaxis == LINEAR) { // draw axis line - rSGP.setColor (clr_axis); + m_pSGP->setColor (clr_axis); if (o_tag && !o_grid && !o_box && s_xcross) { - rSGP.moveAbs (xa_min, xaxispos - charheight); - rSGP.lineAbs (xa_min, xaxispos + charheight); + m_pSGP->moveAbs (xa_min, xaxispos - charheight); + m_pSGP->lineAbs (xa_min, xaxispos + charheight); } - rSGP.moveAbs (xa_min, xaxispos); - rSGP.lineAbs (xa_max, xaxispos); + m_pSGP->moveAbs (xa_min, xaxispos); + m_pSGP->lineAbs (xa_max, xaxispos); if (o_tag && !o_grid && !o_box) { - rSGP.moveAbs (xa_max, xaxispos - charheight); - rSGP.lineAbs (xa_max, xaxispos + charheight); + m_pSGP->moveAbs (xa_max, xaxispos - charheight); + m_pSGP->lineAbs (xa_max, xaxispos + charheight); } if (o_grid == TRUE) { - rSGP.setColor (clr_grid); + m_pSGP->setColor (clr_grid); for (i = 0; i <= x_nint; i++) { - rSGP.moveAbs (xt_min + xn_tickinc * i, ya_max); - rSGP.lineAbs (xt_min + xn_tickinc * i, ya_min); + m_pSGP->moveAbs (xt_min + xn_tickinc * i, ya_max); + m_pSGP->lineAbs (xt_min + xn_tickinc * i, ya_min); } } - rSGP.setTextSize (charheight * 1.5); - rSGP.setTextColor (clr_label, -1); + m_pSGP->setTextSize (charheight * 1.5); + m_pSGP->setTextColor (clr_label, -1); double wText, hText; - rSGP.getTextExtent (c_xlabel.c_str(), &wText, &hText); - rSGP.moveAbs (xa_min + (xa_max-xa_min)/2 - wText / 2, xlbl_row + charheight * 1.5); - rSGP.drawText (c_xlabel); - rSGP.setTextSize (charheight); + m_pSGP->getTextExtent (c_xlabel.c_str(), &wText, &hText); + m_pSGP->moveAbs (xa_min + (xa_max-xa_min)/2 - wText / 2, xlbl_row + charheight * 1.5); + m_pSGP->drawText (c_xlabel); + m_pSGP->setTextSize (charheight); minorinc = xn_tickinc / (o_xminortick + 1); for (i = 0; i <= x_nint; i++) { x = xt_min + xn_tickinc * i; - rSGP.setColor (clr_axis); - rSGP.moveAbs (x, xaxispos); - rSGP.lineAbs (x, xaxispos + xticklen); + m_pSGP->setColor (clr_axis); + m_pSGP->moveAbs (x, xaxispos); + m_pSGP->lineAbs (x, xaxispos + xticklen); if (i != x_nint) for (j = 1; j <= o_xminortick; j++) { x2 = x + minorinc * j; - rSGP.moveAbs (x2, xaxispos); - rSGP.lineAbs (x2, xaxispos + TICKRATIO * xticklen); + m_pSGP->moveAbs (x2, xaxispos); + m_pSGP->lineAbs (x2, xaxispos + TICKRATIO * xticklen); } axis_near = FALSE; if (xaxispos + xtl_ofs > ya_min && o_yaxis != NOAXIS) { @@ -834,9 +835,9 @@ EZPlot::drawAxes() if (o_xtlabel == TRUE && axis_near == FALSE) { snprintf (str, sizeof(str), x_numfmt, xgw_min + xw_tickinc * i); numstr = str_skip_head (str, " "); - rSGP.moveAbs (x-strlen(numstr)*charwidth/2, xaxispos + xtl_ofs); - rSGP.setTextColor (clr_number, -1); - rSGP.drawText (numstr); + m_pSGP->moveAbs (x-strlen(numstr)*charwidth/2, xaxispos + xtl_ofs); + m_pSGP->setTextColor (clr_number, -1); + m_pSGP->drawText (numstr); } } } // X - Axis @@ -848,47 +849,47 @@ EZPlot::drawAxes() if (o_yaxis == LINEAR) { - rSGP.setColor (clr_axis); + m_pSGP->setColor (clr_axis); if (o_tag && !o_grid && !o_box && s_ycross) { - rSGP.moveAbs (yaxispos - charwidth, ya_min); - rSGP.lineAbs (yaxispos + charwidth, ya_min); + m_pSGP->moveAbs (yaxispos - charwidth, ya_min); + m_pSGP->lineAbs (yaxispos + charwidth, ya_min); } - rSGP.moveAbs (yaxispos, ya_min); - rSGP.lineAbs (yaxispos, ya_max); + m_pSGP->moveAbs (yaxispos, ya_min); + m_pSGP->lineAbs (yaxispos, ya_max); if (o_tag && !o_grid && !o_box) { - rSGP.moveAbs (yaxispos - charwidth, ya_max); - rSGP.lineAbs (yaxispos + charwidth, ya_max); + m_pSGP->moveAbs (yaxispos - charwidth, ya_max); + m_pSGP->lineAbs (yaxispos + charwidth, ya_max); } if (o_grid == TRUE) { - rSGP.setColor (clr_grid); + m_pSGP->setColor (clr_grid); for (i = 0; i <= y_nint; i++) { y = yt_min + yn_tickinc * i; - rSGP.moveAbs (xa_max, y); - rSGP.lineAbs (xa_min, y); + m_pSGP->moveAbs (xa_max, y); + m_pSGP->lineAbs (xa_min, y); } } - rSGP.setTextAngle (HALFPI); - rSGP.setTextSize (1.5 * charheight); - rSGP.setTextColor (clr_label, -1); + m_pSGP->setTextAngle (HALFPI); + m_pSGP->setTextSize (1.5 * charheight); + m_pSGP->setTextColor (clr_label, -1); double xExtent, yExtent; - rSGP.getTextExtent (c_ylabel.c_str(), &xExtent, &yExtent); - rSGP.moveAbs (ylbl_col, ya_min + (ya_max-ya_min)/2 - yExtent / 2); - rSGP.drawText (c_ylabel); - rSGP.setTextAngle (0.0); - rSGP.setTextSize (charheight); + m_pSGP->getTextExtent (c_ylabel.c_str(), &xExtent, &yExtent); + m_pSGP->moveAbs (ylbl_col, ya_min + (ya_max-ya_min)/2 - yExtent / 2); + m_pSGP->drawText (c_ylabel); + m_pSGP->setTextAngle (0.0); + m_pSGP->setTextSize (charheight); minorinc = yn_tickinc / (o_yminortick + 1); for (i = 0; i <= y_nint; i++) { y = yt_min + yn_tickinc * i; - rSGP.setColor (clr_axis); - rSGP.moveAbs (yaxispos, y); - rSGP.lineAbs (yaxispos + yticklen, y); + m_pSGP->setColor (clr_axis); + m_pSGP->moveAbs (yaxispos, y); + m_pSGP->lineAbs (yaxispos + yticklen, y); if (i != y_nint) for (j = 1; j <= o_yminortick; j++) { y2 = y + minorinc * j; - rSGP.moveAbs (yaxispos, y2); - rSGP.lineAbs (yaxispos + TICKRATIO * yticklen, y2); + m_pSGP->moveAbs (yaxispos, y2); + m_pSGP->lineAbs (yaxispos + TICKRATIO * yticklen, y2); } axis_near = FALSE; if (yaxispos + ytl_ofs > xa_min && o_xaxis != NOAXIS) { @@ -902,9 +903,9 @@ EZPlot::drawAxes() } if (o_ytlabel == TRUE && axis_near == FALSE) { snprintf (str, sizeof(str), y_numfmt, ygw_min + yw_tickinc * i); - rSGP.moveAbs (yaxispos + ytl_ofs, y + 0.5 * charheight); - rSGP.setTextColor (clr_number, -1); - rSGP.drawText (str); + m_pSGP->moveAbs (yaxispos + ytl_ofs, y + 0.5 * charheight); + m_pSGP->setTextColor (clr_number, -1); + m_pSGP->drawText (str); } } } // Y - Axis @@ -918,34 +919,34 @@ EZPlot::symbol (int sym, double symwidth, double symheight) return; if (sym == SB_CROSS) { - rSGP.moveRel (-0.5 * symwidth, -0.5 * symheight); - rSGP.lineRel (symwidth, symheight); - rSGP.moveRel (-symwidth, 0.0); - rSGP.lineRel (symwidth, -symheight); - rSGP.moveRel (-0.5 * symwidth, 0.5 * symheight); + m_pSGP->moveRel (-0.5 * symwidth, -0.5 * symheight); + m_pSGP->lineRel (symwidth, symheight); + m_pSGP->moveRel (-symwidth, 0.0); + m_pSGP->lineRel (symwidth, -symheight); + m_pSGP->moveRel (-0.5 * symwidth, 0.5 * symheight); } else if (sym == SB_PLUS) { - rSGP.moveRel (-0.5 * symwidth, 0.0); - rSGP.lineRel (symwidth, 0.0); - rSGP.moveRel (-0.5 * symwidth, -0.5 * symheight); - rSGP.lineRel (0.0, symheight); - rSGP.moveRel (0.0, -0.5 * symheight); + m_pSGP->moveRel (-0.5 * symwidth, 0.0); + m_pSGP->lineRel (symwidth, 0.0); + m_pSGP->moveRel (-0.5 * symwidth, -0.5 * symheight); + m_pSGP->lineRel (0.0, symheight); + m_pSGP->moveRel (0.0, -0.5 * symheight); } else if (sym == SB_BOX) { - rSGP.moveRel (-0.5 * symwidth, -0.5 * symheight); - rSGP.lineRel (symwidth, 0.0); - rSGP.lineRel (0.0, symheight); - rSGP.lineRel (-symwidth, 0.0); - rSGP.lineRel (0.0, -symheight); - rSGP.moveRel (0.5 * symwidth, 0.5 * symheight); + m_pSGP->moveRel (-0.5 * symwidth, -0.5 * symheight); + m_pSGP->lineRel (symwidth, 0.0); + m_pSGP->lineRel (0.0, symheight); + m_pSGP->lineRel (-symwidth, 0.0); + m_pSGP->lineRel (0.0, -symheight); + m_pSGP->moveRel (0.5 * symwidth, 0.5 * symheight); } else if (sym == SB_CIRCLE) { - rSGP.drawCircle (symwidth); + m_pSGP->drawCircle (symwidth); } else if (sym == SB_ERRORBAR) { - rSGP.moveRel (-0.5 * symwidth, 0.5 * symheight); - rSGP.lineRel (symwidth, 0.0); - rSGP.moveRel (-0.5 * symwidth, 0.0); - rSGP.lineRel (0.0, -symheight); - rSGP.moveRel (-0.5 * symwidth, 0.0); - rSGP.lineRel (symwidth, 0.0); - rSGP.moveRel (-0.5 * symwidth, 0.5 * symheight); + m_pSGP->moveRel (-0.5 * symwidth, 0.5 * symheight); + m_pSGP->lineRel (symwidth, 0.0); + m_pSGP->moveRel (-0.5 * symwidth, 0.0); + m_pSGP->lineRel (0.0, -symheight); + m_pSGP->moveRel (-0.5 * symwidth, 0.0); + m_pSGP->lineRel (symwidth, 0.0); + m_pSGP->moveRel (-0.5 * symwidth, 0.5 * symheight); } } diff --git a/libctgraphics/ezset.cpp b/libctgraphics/ezset.cpp index 718b250..222b3cd 100644 --- a/libctgraphics/ezset.cpp +++ b/libctgraphics/ezset.cpp @@ -6,7 +6,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: ezset.cpp,v 1.14 2000/12/27 20:09:19 kevin Exp $ +** $Id: ezset.cpp,v 1.15 2000/12/29 15:45:06 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 @@ -81,7 +81,7 @@ EZPlot::do_cmd (int lx) } break; case S_REPLOT: - plot (); + plot (m_pSGP); break; case S_CLEAR: clearCurves (); diff --git a/libctgraphics/pol.cpp b/libctgraphics/pol.cpp index 1d1a13a..9c0e427 100644 --- a/libctgraphics/pol.cpp +++ b/libctgraphics/pol.cpp @@ -6,7 +6,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: pol.cpp,v 1.7 2000/12/27 20:09:19 kevin Exp $ +** $Id: pol.cpp,v 1.8 2000/12/29 15:45:06 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 @@ -54,7 +54,6 @@ const unsigned int POL::NUMCMD = (sizeof(POL::cmdlist) / sizeof (struct POL::Key POL::POL() { - bp = 0; currentf = -1; m_bTrace = false; init(); @@ -832,9 +831,10 @@ POL::usefile (int source, char *fn) sys_error (ERR_SEVERE, "files nested too deeply"); return; } - - bp = 0; /* clear any pushed back input */ - + + while (! m_stackPushBackInput.empty()) + m_stackPushBackInput.pop(); + if (source == P_USE_STR) { filep[currentf] = NULL; } else if (source == P_USE_FILE) { @@ -898,11 +898,13 @@ POL::inchar() int POL::getch (FILE *fp) { - int c; - - if (bp > 0) - return (buf[--bp]); - + int c; + if (m_stackPushBackInput.size() > 0) { + c = m_stackPushBackInput.top(); + m_stackPushBackInput.pop(); + return c; + } + if (fp == NULL) { if ((c = inputline[lineptr]) == EOS) return (EOF); @@ -916,22 +918,21 @@ POL::getch (FILE *fp) return (c); } -/* push character back on input */ +// push character back on input void POL::ungetch (int c) { - if (bp > BUFSIZE) - sys_error (ERR_SEVERE, "too many characters pushed back [ungetch]"); - else - buf[bp++] = c; + m_stackPushBackInput.push (c); } int POL::get_inputline (FILE *fp) { - lineptr = 0; - bp = 0; + while (! m_stackPushBackInput.empty()) + m_stackPushBackInput.pop(); + + lineptr = 0; if (fgets (inputline, MAXLINE, fp) == NULL) return (EOF); else @@ -941,7 +942,9 @@ POL::get_inputline (FILE *fp) void POL::set_inputline (const char* const line) { - lineptr = 0; - bp = 0; + while (! m_stackPushBackInput.empty()) + m_stackPushBackInput.pop(); + strncpy (inputline, line, MAXLINE); + lineptr = 0; } diff --git a/libctsim/imagefile.cpp b/libctsim/imagefile.cpp index a5e9bec..54f38a1 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.24 2000/12/22 04:18:00 kevin Exp $ +** $Id: imagefile.cpp,v 1.25 2000/12/29 15:45:06 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 @@ -62,7 +62,7 @@ ImageFile::filterResponse (const char* const domainName, double bw, const char* for (int i = -hx; i <= hx; i++) { for (int j = -hy; j <= hy; j++) { - double r = sqrt (i * i + j * j); + double r = ::sqrt (i * i + j * j); v[i+hx][j+hy] = filter.response (r); } @@ -170,7 +170,7 @@ ImageFile::comparativeStatistics (const ImageFile& imComp, double& d, double& r, } } - d = sqrt (sqErrorSum / sqDiffFromMeanSum); + d = ::sqrt (sqErrorSum / sqDiffFromMeanSum); r = absErrorSum / absValueSum; int hx = m_nx / 2; @@ -364,6 +364,242 @@ ImageFile::divideImages (const ImageFile& rRHS, ImageFile& result) const } +bool +ImageFile::invertPixelValues (ImageFile& result) const +{ + if (m_nx != result.nx() || m_ny != result.ny()) { + sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]"); + return false; + } + + ImageFileArrayConst vLHS = getArray(); + ImageFileArray vResult = result.getArray(); + + for (int ix = 0; ix < m_nx; ix++) { + ImageFileColumnConst in = vLHS[ix]; + ImageFileColumn out = vResult[ix]; + for (int iy = 0; iy < m_ny; iy++) + *out++ = - *in++; + } + + return true; +} + +bool +ImageFile::sqrt (ImageFile& result) const +{ + if (m_nx != result.nx() || m_ny != result.ny()) { + sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]"); + return false; + } + + ImageFileArrayConst vLHS = getArray(); + ImageFileArray vResult = result.getArray(); + + for (int ix = 0; ix < m_nx; ix++) { + ImageFileColumnConst in = vLHS[ix]; + ImageFileColumn out = vResult[ix]; + for (int iy = 0; iy < m_ny; iy++) + if (*in < 0) + *out++ = -::sqrt(-*in++); + else + *out++ = ::sqrt(*in++); + } + + return true; +} + +bool +ImageFile::log (ImageFile& result) const +{ + if (m_nx != result.nx() || m_ny != result.ny()) { + sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]"); + return false; + } + + ImageFileArrayConst vLHS = getArray(); + ImageFileArray vResult = result.getArray(); + + for (int ix = 0; ix < m_nx; ix++) { + ImageFileColumnConst in = vLHS[ix]; + ImageFileColumn out = vResult[ix]; + for (int iy = 0; iy < m_ny; iy++) + if (*in <= 0) + *out++ = 0; + else + *out++ = ::log(*in++); + } + + return true; +} + +bool +ImageFile::exp (ImageFile& result) const +{ + if (m_nx != result.nx() || m_ny != result.ny()) { + sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]"); + return false; + } + + ImageFileArrayConst vLHS = getArray(); + ImageFileArray vResult = result.getArray(); + + for (int ix = 0; ix < m_nx; ix++) { + ImageFileColumnConst in = vLHS[ix]; + ImageFileColumn out = vResult[ix]; + for (int iy = 0; iy < m_ny; iy++) + *out++ = ::exp (*in++); + } + + return true; +} + +bool +ImageFile::FFTMagnitude (ImageFile& result) const +{ + if (m_nx != result.nx() || m_ny != result.ny()) { + sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]"); + return false; + } + + ImageFileArrayConst vLHS = getArray(); + ImageFileArray vResult = result.getArray(); + + int ix, iy; + double* pY = new double [m_ny]; + std::complex** complexOut = new std::complex* [m_nx]; + for (ix = 0; ix < m_nx; ix++) + complexOut[ix] = new std::complex [m_ny]; + + for (ix = 0; ix < m_nx; ix++) { + for (iy = 0; iy < m_ny; iy++) + pY[iy] = vLHS[ix][iy]; + ProcessSignal::finiteFourierTransform (pY, complexOut[ix], m_ny, ProcessSignal::FORWARD); + } + delete pY; + + std::complex* pX = new std::complex [m_nx]; + std::complex* complexOutCol = new std::complex [m_nx]; + for (iy = 0; iy < m_ny; iy++) { + for (ix = 0; ix < m_nx; ix++) + pX[ix] = complexOut[ix][iy]; + ProcessSignal::finiteFourierTransform (pX, complexOutCol, m_nx, ProcessSignal::FORWARD); + for (ix = 0; ix < m_nx; ix++) + complexOut[ix][iy] = complexOutCol[ix]; + } + delete [] pX; + + for (ix = 0; ix < m_nx; ix++) + ProcessSignal::shuffleFourierToNaturalOrder (complexOut[ix], m_ny); + + + for (iy = 0; iy < m_ny; iy++) { + for (ix = 0; ix < m_nx; ix++) + complexOutCol[ix] = complexOut[ix][iy]; + ProcessSignal::shuffleFourierToNaturalOrder (complexOutCol, m_nx);; + for (ix = 0; ix < m_nx; ix++) + complexOut[ix][iy] = complexOutCol[ix]; + + } + delete [] complexOutCol; + + for (ix = 0; ix < m_nx; ix++) + for (iy = 0; iy < m_ny; iy++) + vResult[ix][iy] = std::abs (complexOut[ix][iy]); + + for (ix = 0; ix < m_nx; ix++) + delete [] complexOut[ix]; + + delete [] complexOut; + + return true; +} + +bool +ImageFile::FFTPhase (ImageFile& result) const +{ + if (m_nx != result.nx() || m_ny != result.ny()) { + sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]"); + return false; + } + + ImageFileArrayConst vLHS = getArray(); + ImageFileArray vResult = result.getArray(); + + int ix, iy; + double* pY = new double [m_ny]; + std::complex** complexOut = new std::complex* [m_nx]; + for (ix = 0; ix < m_nx; ix++) + complexOut[ix] = new std::complex [m_ny]; + + for (ix = 0; ix < m_nx; ix++) { + for (iy = 0; iy < m_ny; iy++) + pY[iy] = vLHS[ix][iy]; + ProcessSignal::finiteFourierTransform (pY, complexOut[ix], m_ny, ProcessSignal::FORWARD); + } + delete pY; + + std::complex* pX = new std::complex [m_nx]; + std::complex* complexOutCol = new std::complex [m_nx]; + for (iy = 0; iy < m_ny; iy++) { + for (ix = 0; ix < m_nx; ix++) + pX[ix] = complexOut[ix][iy]; + ProcessSignal::finiteFourierTransform (pX, complexOutCol, m_nx, ProcessSignal::FORWARD); + for (ix = 0; ix < m_nx; ix++) + complexOut[ix][iy] = complexOutCol[ix]; + } + delete [] pX; + + for (ix = 0; ix < m_nx; ix++) + ProcessSignal::shuffleFourierToNaturalOrder (complexOut[ix], m_ny); + + + for (iy = 0; iy < m_ny; iy++) { + for (ix = 0; ix < m_nx; ix++) + complexOutCol[ix] = complexOut[ix][iy]; + ProcessSignal::shuffleFourierToNaturalOrder (complexOutCol, m_nx);; + for (ix = 0; ix < m_nx; ix++) + complexOut[ix][iy] = complexOutCol[ix]; + + } + delete [] complexOutCol; + + for (ix = 0; ix < m_nx; ix++) + for (iy = 0; iy < m_ny; iy++) + vResult[ix][iy] = atan (complexOut[ix][iy].imag / complexOut[ix][iy].real); + + for (ix = 0; ix < m_nx; ix++) + delete [] complexOut[ix]; + + delete [] complexOut; + + return true; +} + +bool +ImageFile::square (ImageFile& result) const +{ + if (m_nx != result.nx() || m_ny != result.ny()) { + sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]"); + return false; + } + + ImageFileArrayConst vLHS = getArray(); + ImageFileArray vResult = result.getArray(); + + for (int ix = 0; ix < m_nx; ix++) { + ImageFileColumnConst in = vLHS[ix]; + ImageFileColumn out = vResult[ix]; + for (int iy = 0; iy < m_ny; iy++) { + *out++ = *in * *in; + in++; + } + } + + return true; +} + + void ImageFile::writeImagePGM (const char *outfile, int nxcell, int nycell, double densmin, double densmax) { diff --git a/libctsim/procsignal.cpp b/libctsim/procsignal.cpp index 7bf2711..a7933b2 100644 --- a/libctsim/procsignal.cpp +++ b/libctsim/procsignal.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: procsignal.cpp,v 1.10 2000/12/16 06:12:47 kevin Exp $ +** $Id: procsignal.cpp,v 1.11 2000/12/29 15:45:06 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 @@ -78,7 +78,7 @@ const int ProcessSignal::s_iFilterGenerationCount = sizeof(s_aszFilterGeneration // ProcessSignal // ProcessSignal::ProcessSignal (const char* szFilterName, const char* szFilterMethodName, double dBandwidth, double dSignalIncrement, int nSignalPoints, double dFilterParam, const char* szDomainName, const char* szFilterGenerationName, int iZeropad, int iPreinterpolationFactor, int iTraceLevel, int iGeometry, double dFocalLength, SGP* pSGP) - : m_adFourierCosTable(NULL), m_adFourierSinTable(NULL), m_adFilter(NULL), m_fail(false) +: m_adFourierCosTable(NULL), m_adFourierSinTable(NULL), m_adFilter(NULL), m_fail(false) { m_idFilterMethod = convertFilterMethodNameToID (szFilterMethodName); if (m_idFilterMethod == FILTER_METHOD_INVALID) { @@ -108,7 +108,7 @@ ProcessSignal::ProcessSignal (const char* szFilterName, const char* szFilterMeth m_failMessage += szDomainName; return; } - + init (m_idFilter, m_idFilterMethod, dBandwidth, dSignalIncrement, nSignalPoints, dFilterParam, m_idDomain, m_idFilterGeneration, iZeropad, iPreinterpolationFactor, iTraceLevel, iGeometry, dFocalLength, pSGP); } @@ -123,7 +123,7 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw m_idFilterGeneration = idFilterGeneration; m_idGeometry = iGeometry; m_dFocalLength = dFocalLength; - + if (m_idFilter == SignalFilter::FILTER_INVALID || m_idDomain == SignalFilter::DOMAIN_INVALID || m_idFilterMethod == FILTER_METHOD_INVALID || m_idFilterGeneration == FILTER_GENERATION_INVALID) { m_fail = true; return; @@ -137,14 +137,14 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw m_dFilterParam = dFilterParam; m_iZeropad = iZeropad; m_iPreinterpolationFactor = iPreinterpolationFactor; - + // scale signalInc/BW to signalInc/2 to adjust for imaginary detector // through origin of phantom, see Kak-Slaney Fig 3.22, for Collinear if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) { m_dSignalInc /= 2; m_dBandwidth *= 2; } - + if (m_idFilterMethod == FILTER_METHOD_FFT) { #if HAVE_FFTW m_idFilterMethod = FILTER_METHOD_RFFTW; @@ -154,165 +154,165 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw return; #endif } - + bool m_bFrequencyFiltering = true; if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) m_bFrequencyFiltering = false; - + // Spatial-based filtering if (! m_bFrequencyFiltering) { - + if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) { - m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1; - m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1); - m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1); - m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1); - SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL); - m_adFilter = new double[ m_nFilterPoints ]; - filter.copyFilterData (m_adFilter, 0, m_nFilterPoints); + m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1; + m_dFilterMin = -m_dSignalInc * (m_nSignalPoints - 1); + m_dFilterMax = m_dSignalInc * (m_nSignalPoints - 1); + m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1); + SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL); + m_adFilter = new double[ m_nFilterPoints ]; + filter.copyFilterData (m_adFilter, 0, m_nFilterPoints); } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) { - m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1; - m_dFilterMin = -1. / (2 * m_dSignalInc); - m_dFilterMax = 1. / (2 * m_dSignalInc); - m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1); - SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY); - m_adFilter = new double[ m_nFilterPoints ]; - double* adFrequencyFilter = new double [m_nFilterPoints]; - filter.copyFilterData (adFrequencyFilter, 0, m_nFilterPoints); + m_nFilterPoints = 2 * (m_nSignalPoints - 1) + 1; + m_dFilterMin = -1. / (2 * m_dSignalInc); + m_dFilterMax = 1. / (2 * m_dSignalInc); + m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1); + SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY); + m_adFilter = new double[ m_nFilterPoints ]; + double* adFrequencyFilter = new double [m_nFilterPoints]; + filter.copyFilterData (adFrequencyFilter, 0, m_nFilterPoints); #ifdef HAVE_SGP - EZPlot* pEZPlot = NULL; - if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) { - pEZPlot = new EZPlot (*pSGP); - pEZPlot->ezset ("title Filter Response: Natural Order"); - pEZPlot->ezset ("ylength 0.25"); - pEZPlot->addCurve (adFrequencyFilter, m_nFilterPoints); - pEZPlot->plot(); - } + EZPlot* pEZPlot = NULL; + if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) { + pEZPlot = new EZPlot (); + pEZPlot->ezset ("title Filter Response: Natural Order"); + pEZPlot->ezset ("ylength 0.25"); + pEZPlot->addCurve (adFrequencyFilter, m_nFilterPoints); + pEZPlot->plot (pSGP); + } #endif - shuffleNaturalToFourierOrder (adFrequencyFilter, m_nFilterPoints); + shuffleNaturalToFourierOrder (adFrequencyFilter, m_nFilterPoints); #ifdef HAVE_SGP - if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) { - pEZPlot->ezset ("title Filter Response: Fourier Order"); - pEZPlot->ezset ("ylength 0.25"); - pEZPlot->ezset ("yporigin 0.25"); - pEZPlot->addCurve (adFrequencyFilter, m_nFilterPoints); - pEZPlot->plot(); - } + if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) { + pEZPlot->ezset ("title Filter Response: Fourier Order"); + pEZPlot->ezset ("ylength 0.25"); + pEZPlot->ezset ("yporigin 0.25"); + pEZPlot->addCurve (adFrequencyFilter, m_nFilterPoints); + pEZPlot->plot (pSGP); + } #endif - ProcessSignal::finiteFourierTransform (adFrequencyFilter, m_adFilter, m_nFilterPoints, -1); - delete adFrequencyFilter; + ProcessSignal::finiteFourierTransform (adFrequencyFilter, m_adFilter, m_nFilterPoints, -1); + delete adFrequencyFilter; #ifdef HAVE_SGP - if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) { - pEZPlot->ezset ("title Inverse Fourier Frequency: Fourier Order"); - pEZPlot->ezset ("ylength 0.25"); - pEZPlot->ezset ("yporigin 0.50"); - pEZPlot->addCurve (m_adFilter, m_nFilterPoints); - pEZPlot->plot(); - } + if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) { + pEZPlot->ezset ("title Inverse Fourier Frequency: Fourier Order"); + pEZPlot->ezset ("ylength 0.25"); + pEZPlot->ezset ("yporigin 0.50"); + pEZPlot->addCurve (m_adFilter, m_nFilterPoints); + pEZPlot->plot (pSGP); + } #endif - shuffleFourierToNaturalOrder (m_adFilter, m_nFilterPoints); + shuffleFourierToNaturalOrder (m_adFilter, m_nFilterPoints); #ifdef HAVE_SGP - if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) { - pEZPlot->ezset ("title Inverse Fourier Frequency: Natural Order"); - pEZPlot->ezset ("ylength 0.25"); - pEZPlot->ezset ("yporigin 0.75"); - pEZPlot->addCurve (m_adFilter, m_nFilterPoints); - pEZPlot->plot(); - delete pEZPlot; - } + if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) { + pEZPlot->ezset ("title Inverse Fourier Frequency: Natural Order"); + pEZPlot->ezset ("ylength 0.25"); + pEZPlot->ezset ("yporigin 0.75"); + pEZPlot->addCurve (m_adFilter, m_nFilterPoints); + pEZPlot->plot (pSGP); + delete pEZPlot; + } #endif - for (i = 0; i < m_nFilterPoints; i++) { - m_adFilter[i] /= m_dSignalInc; - } + for (i = 0; i < m_nFilterPoints; i++) { + m_adFilter[i] /= m_dSignalInc; + } } if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) { - for (i = 0; i < m_nFilterPoints; i++) - m_adFilter[i] *= 0.5; + for (i = 0; i < m_nFilterPoints; i++) + m_adFilter[i] *= 0.5; } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) { - for (i = 0; i < m_nFilterPoints; i++) { - int iDetFromZero = i - ((m_nFilterPoints - 1) / 2); - double sinScale = sin (iDetFromZero * m_dSignalInc); - if (fabs(sinScale) < 1E-7) - sinScale = 1; - else - sinScale = (iDetFromZero * m_dSignalInc) / sinScale; - double dScale = 0.5 * sinScale * sinScale; - m_adFilter[i] *= dScale; - } + for (i = 0; i < m_nFilterPoints; i++) { + int iDetFromZero = i - ((m_nFilterPoints - 1) / 2); + double sinScale = sin (iDetFromZero * m_dSignalInc); + if (fabs(sinScale) < 1E-7) + sinScale = 1; + else + sinScale = (iDetFromZero * m_dSignalInc) / sinScale; + double dScale = 0.5 * sinScale * sinScale; + m_adFilter[i] *= dScale; + } } // if (geometry) } // if (spatial filtering) - + else if (m_bFrequencyFiltering) { // Frequency-based filtering - + if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) { // calculate number of filter points with zeropadding m_nFilterPoints = m_nSignalPoints; if (m_iZeropad > 0) { - double logBase2 = log(m_nFilterPoints) / log(2); - int nextPowerOf2 = static_cast(floor(logBase2)); - if (logBase2 != floor(logBase2)) - nextPowerOf2++; - nextPowerOf2 += (m_iZeropad - 1); - m_nFilterPoints = 1 << nextPowerOf2; + double logBase2 = log(m_nFilterPoints) / log(2); + int nextPowerOf2 = static_cast(floor(logBase2)); + if (logBase2 != floor(logBase2)) + nextPowerOf2++; + nextPowerOf2 += (m_iZeropad - 1); + m_nFilterPoints = 1 << nextPowerOf2; #ifdef DEBUG - if (m_traceLevel >= Trace::TRACE_CONSOLE) - std::cout << "nFilterPoints = " << m_nFilterPoints << endl; + if (m_traceLevel >= Trace::TRACE_CONSOLE) + std::cout << "nFilterPoints = " << m_nFilterPoints << endl; #endif } m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor; - + if (m_nFilterPoints % 2) { // Odd - m_dFilterMin = -1. / (2 * m_dSignalInc); - m_dFilterMax = 1. / (2 * m_dSignalInc); - m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1); + m_dFilterMin = -1. / (2 * m_dSignalInc); + m_dFilterMax = 1. / (2 * m_dSignalInc); + m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (m_nFilterPoints - 1); } else { // Even - m_dFilterMin = -1. / (2 * m_dSignalInc); - m_dFilterMax = 1. / (2 * m_dSignalInc); - m_dFilterInc = (m_dFilterMax - m_dFilterMin) / m_nFilterPoints; - m_dFilterMax -= m_dFilterInc; + m_dFilterMin = -1. / (2 * m_dSignalInc); + m_dFilterMax = 1. / (2 * m_dSignalInc); + m_dFilterInc = (m_dFilterMax - m_dFilterMin) / m_nFilterPoints; + m_dFilterMax -= m_dFilterInc; } - + SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY); m_adFilter = new double [m_nFilterPoints]; filter.copyFilterData (m_adFilter, 0, m_nFilterPoints); - + // This doesn't work! // Need to add filtering for divergent geometries & Frequency/Direct filtering if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) { - for (i = 0; i < m_nFilterPoints; i++) - m_adFilter[i] *= 0.5; + for (i = 0; i < m_nFilterPoints; i++) + m_adFilter[i] *= 0.5; } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) { - for (i = 0; i < m_nFilterPoints; i++) { - int iDetFromZero = i - ((m_nFilterPoints - 1) / 2); - double sinScale = sin (iDetFromZero * m_dSignalInc); - if (fabs(sinScale) < 1E-7) - sinScale = 1; - else - sinScale = (iDetFromZero * m_dSignalInc) / sinScale; - double dScale = 0.5 * sinScale * sinScale; - m_adFilter[i] *= dScale; - } + for (i = 0; i < m_nFilterPoints; i++) { + int iDetFromZero = i - ((m_nFilterPoints - 1) / 2); + double sinScale = sin (iDetFromZero * m_dSignalInc); + if (fabs(sinScale) < 1E-7) + sinScale = 1; + else + sinScale = (iDetFromZero * m_dSignalInc) / sinScale; + double dScale = 0.5 * sinScale * sinScale; + m_adFilter[i] *= dScale; + } } #ifdef HAVE_SGP EZPlot* pEZPlot = NULL; if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) { - pEZPlot = new EZPlot (*pSGP); - pEZPlot->ezset ("title Filter Filter: Natural Order"); - pEZPlot->ezset ("ylength 0.50"); - pEZPlot->ezset ("yporigin 0.00"); - pEZPlot->addCurve (m_adFilter, m_nFilterPoints); - pEZPlot->plot(); + pEZPlot = new EZPlot; + pEZPlot->ezset ("title Filter Filter: Natural Order"); + pEZPlot->ezset ("ylength 0.50"); + pEZPlot->ezset ("yporigin 0.00"); + pEZPlot->addCurve (m_adFilter, m_nFilterPoints); + pEZPlot->plot (pSGP); } #endif shuffleNaturalToFourierOrder (m_adFilter, m_nFilterPoints); #ifdef HAVE_SGP if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) { - pEZPlot->ezset ("title Filter Filter: Fourier Order"); - pEZPlot->ezset ("ylength 0.50"); - pEZPlot->ezset ("yporigin 0.50"); - pEZPlot->addCurve (m_adFilter, m_nFilterPoints); - pEZPlot->plot(); - delete pEZPlot; + pEZPlot->ezset ("title Filter Filter: Fourier Order"); + pEZPlot->ezset ("ylength 0.50"); + pEZPlot->ezset ("yporigin 0.50"); + pEZPlot->addCurve (m_adFilter, m_nFilterPoints); + pEZPlot->plot (pSGP); + delete pEZPlot; } #endif } else if (m_idFilterGeneration == FILTER_GENERATION_INVERSE_FOURIER) { @@ -323,17 +323,17 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw m_dFilterInc = (m_dFilterMax - m_dFilterMin) / (nSpatialPoints - 1); m_nFilterPoints = nSpatialPoints; if (m_iZeropad > 0) { - double logBase2 = log(nSpatialPoints) / log(2); - int nextPowerOf2 = static_cast(floor(logBase2)); - if (logBase2 != floor(logBase2)) - nextPowerOf2++; - nextPowerOf2 += (m_iZeropad - 1); - m_nFilterPoints = 1 << nextPowerOf2; + double logBase2 = log(nSpatialPoints) / log(2); + int nextPowerOf2 = static_cast(floor(logBase2)); + if (logBase2 != floor(logBase2)) + nextPowerOf2++; + nextPowerOf2 += (m_iZeropad - 1); + m_nFilterPoints = 1 << nextPowerOf2; } m_nOutputPoints = m_nFilterPoints * m_iPreinterpolationFactor; #ifdef DEBUG if (m_traceLevel >= Trace::TRACE_CONSOLE) - std::cout << "nFilterPoints = " << m_nFilterPoints << endl; + std::cout << "nFilterPoints = " << m_nFilterPoints << endl; #endif double* adSpatialFilter = new double [m_nFilterPoints]; SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, nSpatialPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_SPATIAL); @@ -341,48 +341,48 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw #ifdef HAVE_SGP EZPlot* pEZPlot = NULL; if (pSGP && m_traceLevel >= Trace::TRACE_PLOT) { - pEZPlot = new EZPlot (*pSGP); - pEZPlot->ezset ("title Spatial Filter: Natural Order"); - pEZPlot->ezset ("ylength 0.50"); - pEZPlot->ezset ("yporigin 0.00"); - pEZPlot->addCurve (adSpatialFilter, nSpatialPoints); - pEZPlot->plot(); - delete pEZPlot; + pEZPlot = new EZPlot; + pEZPlot->ezset ("title Spatial Filter: Natural Order"); + pEZPlot->ezset ("ylength 0.50"); + pEZPlot->ezset ("yporigin 0.00"); + pEZPlot->addCurve (adSpatialFilter, nSpatialPoints); + pEZPlot->plot (pSGP); + delete pEZPlot; } #endif if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) { - for (i = 0; i < m_nFilterPoints; i++) - adSpatialFilter[i] *= 0.5; + for (i = 0; i < m_nFilterPoints; i++) + adSpatialFilter[i] *= 0.5; } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) { - for (i = 0; i < m_nFilterPoints; i++) { - int iDetFromZero = i - ((m_nFilterPoints - 1) / 2); - double sinScale = sin (iDetFromZero * m_dSignalInc); - if (fabs(sinScale) < 1E-7) - sinScale = 1; - else - sinScale = (iDetFromZero * m_dSignalInc) / sinScale; - double dScale = 0.5 * sinScale * sinScale; - adSpatialFilter[i] *= dScale; - } + for (i = 0; i < m_nFilterPoints; i++) { + int iDetFromZero = i - ((m_nFilterPoints - 1) / 2); + double sinScale = sin (iDetFromZero * m_dSignalInc); + if (fabs(sinScale) < 1E-7) + sinScale = 1; + else + sinScale = (iDetFromZero * m_dSignalInc) / sinScale; + double dScale = 0.5 * sinScale * sinScale; + adSpatialFilter[i] *= dScale; + } } for (i = nSpatialPoints; i < m_nFilterPoints; i++) - adSpatialFilter[i] = 0; - + adSpatialFilter[i] = 0; + m_adFilter = new double [m_nFilterPoints]; std::complex* acInverseFilter = new std::complex [m_nFilterPoints]; finiteFourierTransform (adSpatialFilter, acInverseFilter, m_nFilterPoints, 1); - delete adSpatialFilter; - for (i = 0; i < m_nFilterPoints; i++) - m_adFilter[i] = std::abs(acInverseFilter[i]) * m_dSignalInc; - delete acInverseFilter; + delete adSpatialFilter; + for (i = 0; i < m_nFilterPoints; i++) + m_adFilter[i] = std::abs(acInverseFilter[i]) * m_dSignalInc; + delete acInverseFilter; #ifdef HAVE_SGP if (pEZPlot && m_traceLevel >= Trace::TRACE_PLOT) { - pEZPlot->ezset ("title Spatial Filter: Inverse"); - pEZPlot->ezset ("ylength 0.50"); - pEZPlot->ezset ("yporigin 0.50"); - pEZPlot->addCurve (m_adFilter, m_nFilterPoints); - pEZPlot->plot(); - delete pEZPlot; + pEZPlot->ezset ("title Spatial Filter: Inverse"); + pEZPlot->ezset ("ylength 0.50"); + pEZPlot->ezset ("yporigin 0.50"); + pEZPlot->addCurve (m_adFilter, m_nFilterPoints); + pEZPlot->plot (pSGP); + delete pEZPlot; } #endif } @@ -401,13 +401,13 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw angle += angleIncrement; } } - + #if HAVE_FFTW if (m_idFilterMethod == FILTER_METHOD_FFTW || m_idFilterMethod == FILTER_METHOD_RFFTW) { for (i = 0; i < m_nFilterPoints; i++) //fftw uses unnormalized fft m_adFilter[i] /= m_nFilterPoints; } - + if (m_idFilterMethod == FILTER_METHOD_RFFTW) { m_realPlanForward = rfftw_create_plan (m_nFilterPoints, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE); m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE); @@ -431,23 +431,23 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw ProcessSignal::~ProcessSignal (void) { - delete [] m_adFourierSinTable; - delete [] m_adFourierCosTable; - delete [] m_adFilter; - + delete [] m_adFourierSinTable; + delete [] m_adFourierCosTable; + delete [] m_adFilter; + #if HAVE_FFTW - if (m_idFilterMethod == FILTER_METHOD_FFTW) { - fftw_destroy_plan(m_complexPlanForward); - fftw_destroy_plan(m_complexPlanBackward); - delete [] m_adComplexFftInput; - delete [] m_adComplexFftSignal; - } - if (m_idFilterMethod == FILTER_METHOD_RFFTW) { - rfftw_destroy_plan(m_realPlanForward); - rfftw_destroy_plan(m_realPlanBackward); - delete [] m_adRealFftInput; - delete [] m_adRealFftSignal; - } + if (m_idFilterMethod == FILTER_METHOD_FFTW) { + fftw_destroy_plan(m_complexPlanForward); + fftw_destroy_plan(m_complexPlanBackward); + delete [] m_adComplexFftInput; + delete [] m_adComplexFftSignal; + } + if (m_idFilterMethod == FILTER_METHOD_RFFTW) { + rfftw_destroy_plan(m_realPlanForward); + rfftw_destroy_plan(m_realPlanBackward); + delete [] m_adRealFftInput; + delete [] m_adRealFftSignal; + } #endif } @@ -455,24 +455,24 @@ int ProcessSignal::convertFilterMethodNameToID (const char* const filterMethodName) { int fmID = FILTER_METHOD_INVALID; - + for (int i = 0; i < s_iFilterMethodCount; i++) - if (strcasecmp (filterMethodName, s_aszFilterMethodName[i]) == 0) { + if (strcasecmp (filterMethodName, s_aszFilterMethodName[i]) == 0) { fmID = i; break; } - - return (fmID); + + return (fmID); } const char * ProcessSignal::convertFilterMethodIDToName (const int fmID) { static const char *name = ""; - + if (fmID >= 0 && fmID < s_iFilterMethodCount) - return (s_aszFilterMethodName [fmID]); - + return (s_aszFilterMethodName [fmID]); + return (name); } @@ -480,10 +480,10 @@ const char * ProcessSignal::convertFilterMethodIDToTitle (const int fmID) { static const char *title = ""; - + if (fmID >= 0 && fmID < s_iFilterMethodCount) - return (s_aszFilterMethodTitle [fmID]); - + return (s_aszFilterMethodTitle [fmID]); + return (title); } @@ -492,24 +492,24 @@ int ProcessSignal::convertFilterGenerationNameToID (const char* const fgName) { int fgID = FILTER_GENERATION_INVALID; - + for (int i = 0; i < s_iFilterGenerationCount; i++) - if (strcasecmp (fgName, s_aszFilterGenerationName[i]) == 0) { + if (strcasecmp (fgName, s_aszFilterGenerationName[i]) == 0) { fgID = i; break; } - - return (fgID); + + return (fgID); } const char * ProcessSignal::convertFilterGenerationIDToName (const int fgID) { static const char *name = ""; - + if (fgID >= 0 && fgID < s_iFilterGenerationCount) - return (s_aszFilterGenerationName [fgID]); - + return (s_aszFilterGenerationName [fgID]); + return (name); } @@ -517,10 +517,10 @@ const char * ProcessSignal::convertFilterGenerationIDToTitle (const int fgID) { static const char *name = ""; - + if (fgID >= 0 && fgID < s_iFilterGenerationCount) - return (s_aszFilterGenerationTitle [fgID]); - + return (s_aszFilterGenerationTitle [fgID]); + return (name); } @@ -531,7 +531,7 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const int i; for (i = 0; i < m_nSignalPoints; i++) input[i] = constInput[i]; - + if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) { for (int i = 0; i < m_nSignalPoints; i++) { int iDetFromCenter = i - (m_nSignalPoints / 2); @@ -544,8 +544,8 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const } } if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) { - for (i = 0; i < m_nSignalPoints; i++) - output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints); + for (i = 0; i < m_nSignalPoints; i++) + output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints); } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) { double* inputSignal = new double [m_nFilterPoints]; for (i = 0; i < m_nSignalPoints; i++) @@ -554,15 +554,15 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const inputSignal[i] = 0; // zeropad std::complex* fftSignal = new std::complex [m_nFilterPoints]; finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, -1); - delete inputSignal; + delete inputSignal; for (i = 0; i < m_nFilterPoints; i++) fftSignal[i] *= m_adFilter[i]; double* inverseFourier = new double [m_nFilterPoints]; finiteFourierTransform (fftSignal, inverseFourier, m_nFilterPoints, 1); - delete fftSignal; + delete fftSignal; for (i = 0; i < m_nSignalPoints; i++) output[i] = inverseFourier[i]; - delete inverseFourier; + delete inverseFourier; } else if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) { double* inputSignal = new double [m_nFilterPoints]; for (i = 0; i < m_nSignalPoints; i++) @@ -571,50 +571,50 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const inputSignal[i] = 0; // zeropad std::complex* fftSignal = new std::complex [m_nFilterPoints]; finiteFourierTransform (inputSignal, fftSignal, -1); - delete inputSignal; + delete inputSignal; for (i = 0; i < m_nFilterPoints; i++) fftSignal[i] *= m_adFilter[i]; double* inverseFourier = new double [m_nFilterPoints]; finiteFourierTransform (fftSignal, inverseFourier, 1); - delete fftSignal; + delete fftSignal; for (i = 0; i < m_nSignalPoints; i++) output[i] = inverseFourier[i]; - delete inverseFourier; + delete inverseFourier; } #if HAVE_FFTW else if (m_idFilterMethod == FILTER_METHOD_RFFTW) { - for (i = 0; i < m_nSignalPoints; i++) - m_adRealFftInput[i] = input[i]; - - fftw_real* fftOutput = new fftw_real [ m_nFilterPoints ]; - rfftw_one (m_realPlanForward, m_adRealFftInput, fftOutput); - for (i = 0; i < m_nFilterPoints; i++) - m_adRealFftSignal[i] = m_adFilter[i] * fftOutput[i]; - delete [] fftOutput; - for (i = m_nFilterPoints; i < m_nOutputPoints; i++) + for (i = 0; i < m_nSignalPoints; i++) + m_adRealFftInput[i] = input[i]; + + fftw_real* fftOutput = new fftw_real [ m_nFilterPoints ]; + rfftw_one (m_realPlanForward, m_adRealFftInput, fftOutput); + for (i = 0; i < m_nFilterPoints; i++) + m_adRealFftSignal[i] = m_adFilter[i] * fftOutput[i]; + delete [] fftOutput; + for (i = m_nFilterPoints; i < m_nOutputPoints; i++) m_adRealFftSignal[i] = 0; - - fftw_real* ifftOutput = new fftw_real [ m_nOutputPoints ]; - rfftw_one (m_realPlanBackward, m_adRealFftSignal, ifftOutput); - for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++) - output[i] = ifftOutput[i]; - delete [] ifftOutput; + + fftw_real* ifftOutput = new fftw_real [ m_nOutputPoints ]; + rfftw_one (m_realPlanBackward, m_adRealFftSignal, ifftOutput); + for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++) + output[i] = ifftOutput[i]; + delete [] ifftOutput; } else if (m_idFilterMethod == FILTER_METHOD_FFTW) { - for (i = 0; i < m_nSignalPoints; i++) - m_adComplexFftInput[i].re = input[i]; - - fftw_complex* fftOutput = new fftw_complex [ m_nFilterPoints ]; - fftw_one (m_complexPlanForward, m_adComplexFftInput, fftOutput); - for (i = 0; i < m_nFilterPoints; i++) { - m_adComplexFftSignal[i].re = m_adFilter[i] * fftOutput[i].re; - m_adComplexFftSignal[i].im = m_adFilter[i] * fftOutput[i].im; - } - delete [] fftOutput; - fftw_complex* ifftOutput = new fftw_complex [ m_nOutputPoints ]; - fftw_one (m_complexPlanBackward, m_adComplexFftSignal, ifftOutput); - for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++) - output[i] = ifftOutput[i].re; - delete [] ifftOutput; + for (i = 0; i < m_nSignalPoints; i++) + m_adComplexFftInput[i].re = input[i]; + + fftw_complex* fftOutput = new fftw_complex [ m_nFilterPoints ]; + fftw_one (m_complexPlanForward, m_adComplexFftInput, fftOutput); + for (i = 0; i < m_nFilterPoints; i++) { + m_adComplexFftSignal[i].re = m_adFilter[i] * fftOutput[i].re; + m_adComplexFftSignal[i].im = m_adFilter[i] * fftOutput[i].im; + } + delete [] fftOutput; + fftw_complex* ifftOutput = new fftw_complex [ m_nOutputPoints ]; + fftw_one (m_complexPlanBackward, m_adComplexFftSignal, ifftOutput); + for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++) + output[i] = ifftOutput[i].re; + delete [] ifftOutput; } #endif delete input; @@ -622,36 +622,36 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const /* NAME - * convolve Discrete convolution of two functions - * - * SYNOPSIS - * r = convolve (f1, f2, dx, n, np, func_type) - * double r Convolved result - * double f1[], f2[] Functions to be convolved - * double dx Difference between successive x values - * int n Array index to center convolution about - * int np Number of points in f1 array - * int func_type EVEN or ODD or EVEN_AND_ODD function f2 - * - * NOTES - * f1 is the projection data, its indices range from 0 to np - 1. - * The index for f2, the filter, ranges from -(np-1) to (np-1). - * There are 3 ways to handle the negative vertices of f2: - * 1. If we know f2 is an EVEN function, then f2[-n] = f2[n]. - * All filters used in reconstruction are even. - * 2. If we know f2 is an ODD function, then f2[-n] = -f2[n] - * 3. If f2 is both ODD AND EVEN, then we must store the value of f2 - * for negative indices. Since f2 must range from -(np-1) to (np-1), - * if we add (np - 1) to f2's array index, then f2's index will - * range from 0 to 2 * (np - 1), and the origin, x = 0, will be - * stored at f2[np-1]. - */ +* convolve Discrete convolution of two functions +* +* SYNOPSIS +* r = convolve (f1, f2, dx, n, np, func_type) +* double r Convolved result +* double f1[], f2[] Functions to be convolved +* double dx Difference between successive x values +* int n Array index to center convolution about +* int np Number of points in f1 array +* int func_type EVEN or ODD or EVEN_AND_ODD function f2 +* +* NOTES +* f1 is the projection data, its indices range from 0 to np - 1. +* The index for f2, the filter, ranges from -(np-1) to (np-1). +* There are 3 ways to handle the negative vertices of f2: +* 1. If we know f2 is an EVEN function, then f2[-n] = f2[n]. +* All filters used in reconstruction are even. +* 2. If we know f2 is an ODD function, then f2[-n] = -f2[n] +* 3. If f2 is both ODD AND EVEN, then we must store the value of f2 +* for negative indices. Since f2 must range from -(np-1) to (np-1), +* if we add (np - 1) to f2's array index, then f2's index will +* range from 0 to 2 * (np - 1), and the origin, x = 0, will be +* stored at f2[np-1]. +*/ double ProcessSignal::convolve (const double func[], const double dx, const int n, const int np) const { double sum = 0.0; - + #if UNOPTIMIZED_CONVOLUTION for (int i = 0; i < np; i++) sum += func[i] * m_adFilter[n - i + (np - 1)]; @@ -660,7 +660,7 @@ ProcessSignal::convolve (const double func[], const double dx, const int n, cons for (int i = 0; i < np; i++) sum += *func++ * *f2--; #endif - + return (sum * dx); } @@ -669,16 +669,16 @@ double ProcessSignal::convolve (const float func[], const double dx, const int n, const int np) const { double sum = 0.0; - + #if UNOPTIMIZED_CONVOLUTION -for (int i = 0; i < np; i++) - sum += func[i] * m_adFilter[n - i + (np - 1)]; + for (int i = 0; i < np; i++) + sum += func[i] * m_adFilter[n - i + (np - 1)]; #else -double* f2 = m_adFilter + n + (np - 1); -for (int i = 0; i < np; i++) - sum += *func++ * *f2--; + double* f2 = m_adFilter + n + (np - 1); + for (int i = 0; i < np; i++) + sum += *func++ * *f2--; #endif - + return (sum * dx); } @@ -686,12 +686,12 @@ for (int i = 0; i < np; i++) void ProcessSignal::finiteFourierTransform (const double input[], double output[], const int n, int direction) { - std::complex* complexOutput = new std::complex [n]; - - finiteFourierTransform (input, complexOutput, n, direction); - for (int i = 0; i < n; i++) - output[i] = complexOutput[i].real(); - delete [] complexOutput; + std::complex* complexOutput = new std::complex [n]; + + finiteFourierTransform (input, complexOutput, n, direction); + for (int i = 0; i < n; i++) + output[i] = complexOutput[i].real(); + delete [] complexOutput; } void @@ -701,7 +701,7 @@ ProcessSignal::finiteFourierTransform (const double input[], std::complex input[], std:: direction = -1; else direction = 1; - + double angleIncrement = direction * 2 * PI / n; for (int i = 0; i < n; i++) { std::complex sum (0,0); @@ -750,10 +750,10 @@ ProcessSignal::finiteFourierTransform (const std::complex input[], doubl direction = -1; else direction = 1; - + double angleIncrement = direction * 2 * PI / n; for (int i = 0; i < n; i++) { - double sumReal = 0; + double sumReal = 0; for (int j = 0; j < n; j++) { double angle = i * j * angleIncrement; sumReal += input[j].real() * cos(angle) - input[j].imag() * sin(angle); @@ -774,17 +774,17 @@ ProcessSignal::finiteFourierTransform (const double input[], std::complex 0) { - sumReal += input[j] * m_adFourierCosTable[tableIndex]; - sumImag += input[j] * m_adFourierSinTable[tableIndex]; + sumReal += input[j] * m_adFourierCosTable[tableIndex]; + sumImag += input[j] * m_adFourierSinTable[tableIndex]; } else { - sumReal += input[j] * m_adFourierCosTable[tableIndex]; - sumImag -= input[j] * m_adFourierSinTable[tableIndex]; + sumReal += input[j] * m_adFourierCosTable[tableIndex]; + sumImag -= input[j] * m_adFourierSinTable[tableIndex]; } } if (direction < 0) { @@ -803,21 +803,21 @@ ProcessSignal::finiteFourierTransform (const std::complex input[], std:: direction = -1; else direction = 1; - + for (int i = 0; i < m_nFilterPoints; i++) { double sumReal = 0, sumImag = 0; for (int j = 0; j < m_nFilterPoints; j++) { int tableIndex = i * j; if (direction > 0) { - sumReal += input[j].real() * m_adFourierCosTable[tableIndex] - - input[j].imag() * m_adFourierSinTable[tableIndex]; - sumImag += input[j].real() * m_adFourierSinTable[tableIndex] - + input[j].imag() * m_adFourierCosTable[tableIndex]; + sumReal += input[j].real() * m_adFourierCosTable[tableIndex] + - input[j].imag() * m_adFourierSinTable[tableIndex]; + sumImag += input[j].real() * m_adFourierSinTable[tableIndex] + + input[j].imag() * m_adFourierCosTable[tableIndex]; } else { - sumReal += input[j].real() * m_adFourierCosTable[tableIndex] - - input[j].imag() * -m_adFourierSinTable[tableIndex]; - sumImag += input[j].real() * -m_adFourierSinTable[tableIndex] - + input[j].imag() * m_adFourierCosTable[tableIndex]; + sumReal += input[j].real() * m_adFourierCosTable[tableIndex] + - input[j].imag() * -m_adFourierSinTable[tableIndex]; + sumImag += input[j].real() * -m_adFourierSinTable[tableIndex] + + input[j].imag() * m_adFourierCosTable[tableIndex]; } } if (direction < 0) { @@ -835,17 +835,17 @@ ProcessSignal::finiteFourierTransform (const std::complex input[], doubl direction = -1; else direction = 1; - + for (int i = 0; i < m_nFilterPoints; i++) { - double sumReal = 0; + double sumReal = 0; for (int j = 0; j < m_nFilterPoints; j++) { int tableIndex = i * j; if (direction > 0) { - sumReal += input[j].real() * m_adFourierCosTable[tableIndex] - - input[j].imag() * m_adFourierSinTable[tableIndex]; + sumReal += input[j].real() * m_adFourierCosTable[tableIndex] + - input[j].imag() * m_adFourierSinTable[tableIndex]; } else { - sumReal += input[j].real() * m_adFourierCosTable[tableIndex] - - input[j].imag() * -m_adFourierSinTable[tableIndex]; + sumReal += input[j].real() * m_adFourierCosTable[tableIndex] + - input[j].imag() * -m_adFourierSinTable[tableIndex]; } } if (direction < 0) { @@ -862,58 +862,112 @@ ProcessSignal::finiteFourierTransform (const std::complex input[], doubl // Natural Frequency Order: -n/2...-1,0,1...((n/2)-1) // Fourier Frequency Order: 0,1...((n/2)-1),-n/2...-1 -void -ProcessSignal::shuffleNaturalToFourierOrder (double* pdVector, const int n) -{ +void +ProcessSignal::shuffleNaturalToFourierOrder (double* pdVector, const int n) +{ double* pdTemp = new double [n]; - int i; - if (n % 2) { // Odd - int iHalfN = (n - 1) / 2; - - pdTemp[0] = pdVector[iHalfN]; - for (i = 0; i < iHalfN; i++) - pdTemp[i + 1] = pdVector[i + 1 + iHalfN]; - for (i = 0; i < iHalfN; i++) - pdTemp[i + iHalfN + 1] = pdVector[i]; - } else { // Even - int iHalfN = n / 2; - pdTemp[0] = pdVector[iHalfN]; - for (i = 0; i < iHalfN; i++) - pdTemp[i + 1] = pdVector[i + iHalfN]; - for (i = 0; i < iHalfN - 1; i++) - pdTemp[i + iHalfN + 1] = pdVector[i]; - } - - for (i = 0; i < n; i++) - pdVector[i] = pdTemp[i]; - delete pdTemp; -} - - -void -ProcessSignal::shuffleFourierToNaturalOrder (double* pdVector, const int n) -{ + int i; + if (n % 2) { // Odd + int iHalfN = (n - 1) / 2; + + pdTemp[0] = pdVector[iHalfN]; + for (i = 0; i < iHalfN; i++) + pdTemp[i + 1] = pdVector[i + 1 + iHalfN]; + for (i = 0; i < iHalfN; i++) + pdTemp[i + iHalfN + 1] = pdVector[i]; + } else { // Even + int iHalfN = n / 2; + pdTemp[0] = pdVector[iHalfN]; + for (i = 0; i < iHalfN; i++) + pdTemp[i + 1] = pdVector[i + iHalfN]; + for (i = 0; i < iHalfN - 1; i++) + pdTemp[i + iHalfN + 1] = pdVector[i]; + } + + for (i = 0; i < n; i++) + pdVector[i] = pdTemp[i]; + delete pdTemp; +} + +void +ProcessSignal::shuffleNaturalToFourierOrder (std::complex* pdVector, const int n) +{ + std::complex* pdTemp = new std::complex [n]; + int i; + if (n % 2) { // Odd + int iHalfN = (n - 1) / 2; + + pdTemp[0] = pdVector[iHalfN]; + for (i = 0; i < iHalfN; i++) + pdTemp[i + 1] = pdVector[i + 1 + iHalfN]; + for (i = 0; i < iHalfN; i++) + pdTemp[i + iHalfN + 1] = pdVector[i]; + } else { // Even + int iHalfN = n / 2; + pdTemp[0] = pdVector[iHalfN]; + for (i = 0; i < iHalfN; i++) + pdTemp[i + 1] = pdVector[i + iHalfN]; + for (i = 0; i < iHalfN - 1; i++) + pdTemp[i + iHalfN + 1] = pdVector[i]; + } + + for (i = 0; i < n; i++) + pdVector[i] = pdTemp[i]; + delete [] pdTemp; +} + + +void +ProcessSignal::shuffleFourierToNaturalOrder (double* pdVector, const int n) +{ double* pdTemp = new double [n]; - int i; - if (n % 2) { // Odd - int iHalfN = (n - 1) / 2; - + int i; + if (n % 2) { // Odd + int iHalfN = (n - 1) / 2; + pdTemp[iHalfN] = pdVector[0]; - for (i = 0; i < iHalfN; i++) - pdTemp[i + 1 + iHalfN] = pdVector[i + 1]; - for (i = 0; i < iHalfN; i++) - pdTemp[i] = pdVector[i + iHalfN + 1]; - } else { // Even - int iHalfN = n / 2; - pdTemp[iHalfN] = pdVector[0]; - for (i = 0; i < iHalfN; i++) - pdTemp[i] = pdVector[i + iHalfN]; - for (i = 0; i < iHalfN - 1; i++) - pdTemp[i + iHalfN + 1] = pdVector[i+1]; - } - - for (i = 0; i < n; i++) - pdVector[i] = pdTemp[i]; - delete pdTemp; -} - + for (i = 0; i < iHalfN; i++) + pdTemp[i + 1 + iHalfN] = pdVector[i + 1]; + for (i = 0; i < iHalfN; i++) + pdTemp[i] = pdVector[i + iHalfN + 1]; + } else { // Even + int iHalfN = n / 2; + pdTemp[iHalfN] = pdVector[0]; + for (i = 0; i < iHalfN; i++) + pdTemp[i] = pdVector[i + iHalfN]; + for (i = 0; i < iHalfN - 1; i++) + pdTemp[i + iHalfN + 1] = pdVector[i+1]; + } + + for (i = 0; i < n; i++) + pdVector[i] = pdTemp[i]; + delete pdTemp; +} + +void +ProcessSignal::shuffleFourierToNaturalOrder (std::complex* pdVector, const int n) +{ + std::complex* pdTemp = new std::complex [n]; + int i; + if (n % 2) { // Odd + int iHalfN = (n - 1) / 2; + + pdTemp[iHalfN] = pdVector[0]; + for (i = 0; i < iHalfN; i++) + pdTemp[i + 1 + iHalfN] = pdVector[i + 1]; + for (i = 0; i < iHalfN; i++) + pdTemp[i] = pdVector[i + iHalfN + 1]; + } else { // Even + int iHalfN = n / 2; + pdTemp[iHalfN] = pdVector[0]; + for (i = 0; i < iHalfN; i++) + pdTemp[i] = pdVector[i + iHalfN]; + for (i = 0; i < iHalfN - 1; i++) + pdTemp[i + iHalfN + 1] = pdVector[i+1]; + } + + for (i = 0; i < n; i++) + pdVector[i] = pdTemp[i]; + delete [] pdTemp; +} + diff --git a/libctsim/reconstruct.cpp b/libctsim/reconstruct.cpp index ded9420..df45621 100644 --- a/libctsim/reconstruct.cpp +++ b/libctsim/reconstruct.cpp @@ -8,7 +8,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: reconstruct.cpp,v 1.5 2000/12/17 22:30:34 kevin Exp $ +** $Id: reconstruct.cpp,v 1.6 2000/12/29 15:45:06 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 @@ -117,11 +117,11 @@ Reconstructor::plotFilter (SGP* pSGP) adPlotXAxis[i] = f; if (m_pProcessSignal->getFilter()) { - EZPlot ezplot (*pSGP); + EZPlot ezplot; ezplot.ezset ("title Filter Response"); ezplot.addCurve (adPlotXAxis, m_pProcessSignal->getFilter(), nVecFilter); - ezplot.plot(); + ezplot.plot (pSGP); } } delete adPlotXAxis; @@ -171,7 +171,7 @@ Reconstructor::reconstructView (int iStartView, int iViewCount, SGP* pSGP, bool #ifdef HAVE_SGP if (m_iTrace >= Trace::TRACE_PLOT && pSGP) { - EZPlot ezplotProj (*pSGP); + EZPlot ezplotProj; ezplotProj.ezset ("clear"); ezplotProj.ezset ("title Raw Projection"); @@ -183,7 +183,7 @@ Reconstructor::reconstructView (int iStartView, int iViewCount, SGP* pSGP, bool ezplotProj.ezset ("box."); ezplotProj.ezset ("grid."); ezplotProj.addCurve (m_adPlotXAxis, detval, m_rProj.nDet()); - ezplotProj.plot(); + ezplotProj.plot (pSGP); ezplotProj.ezset ("clear"); ezplotProj.ezset ("title Filtered Projection"); ezplotProj.ezset ("xticks major 5"); @@ -193,7 +193,7 @@ Reconstructor::reconstructView (int iStartView, int iViewCount, SGP* pSGP, bool ezplotProj.ezset ("box"); ezplotProj.ezset ("grid"); ezplotProj.addCurve (m_adPlotXAxis, adFilteredProj, m_nFilteredProjections); - ezplotProj.plot(); + ezplotProj.plot (pSGP); } #endif //HAVE_SGP } diff --git a/libctsupport/syserror.cpp b/libctsupport/syserror.cpp index 8721de2..8468f82 100644 --- a/libctsupport/syserror.cpp +++ b/libctsupport/syserror.cpp @@ -2,7 +2,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: syserror.cpp,v 1.9 2000/12/27 20:09:19 kevin Exp $ +** $Id: syserror.cpp,v 1.10 2000/12/29 15:45:06 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 @@ -22,76 +22,87 @@ #include #include #include -#include -#include "ctsupport.h" +#include +#include +#include "ct.h" /* NAME - * sys_error System error handler - * - * SYNOPSIS - * sys_error (severity, msg, args . . .) - * int severity Severity of error - * char *msg Error message - * args Argument list, direct transfer to printf stack - * Can take 24 byte transfer - */ - -static int errorlevel = ERR_WARNING; /* Set error reporting level */ - +* sys_error System error handler +* +* SYNOPSIS +* sys_error (severity, msg, args . . .) +* int severity Severity of error +* char *msg Error message +* args Argument list, direct transfer to printf stack +* Can take 24 byte transfer +*/ + +static int s_reportErrorLevel = ERR_WARNING; // Set error reporting level + +bool g_bRunningWXWindows = false; void sys_error (int severity, const char *msg, ...) { va_list arg; - + va_start(arg, msg); - - sys_verror (severity, msg, arg); + + std::string strOutput; + sys_verror (strOutput, severity, msg, arg); + +// if (g_bRunningWXWindows) +// theApp->getLog() << strOutput.c_str(); +// else + std::cout << strOutput; va_end(arg); } -static int nErrorCount = 0; +static int s_nErrorCount = 0; const static int MAX_ERROR_COUNT = 20; -void sys_verror (int severity, const char *msg, va_list arg) +void sys_verror (std::string& strOutput, int severity, const char *msg, va_list arg) { - if (severity < errorlevel) - return; // ignore error if less than max level + if (severity < s_reportErrorLevel) + return; // ignore error if less than reporting level + + std::ostringstream os; - nErrorCount++; + s_nErrorCount++; if (severity != ERR_FATAL) { - if (nErrorCount > MAX_ERROR_COUNT) + if (s_nErrorCount > MAX_ERROR_COUNT) return; - else if (nErrorCount == MAX_ERROR_COUNT) { - std::cout << "*****************************************************************\n"; - std::cout << "*** M A X I M U M E R R O R C O U N T R E A C H E D ***\n"; - std::cout << "*** ***\n"; - std::cout << "*** No further errors will be reported ***\n"; - std::cout << "*****************************************************************\n"; + else if (s_nErrorCount == MAX_ERROR_COUNT) { + os << "*****************************************************************\n"; + os << "*** M A X I M U M E R R O R C O U N T R E A C H E D ***\n"; + os << "*** ***\n"; + os << "*** No further errors will be reported ***\n"; + os << "*****************************************************************\n"; + strOutput = os.str(); return; } } - + switch (severity) { case ERR_FATAL: - std::cout << "FATAL ERROR: "; + os << "FATAL ERROR: "; break; case ERR_SEVERE: - std::cout << "SEVERE ERROR: "; + os << "SEVERE ERROR: "; break; case ERR_WARNING: - std::cout << "WARNING ERROR: "; + os << "WARNING ERROR: "; break; case ERR_TRACE: - std::cout << "Trace: "; + os << "Trace: "; break; default: - std::cout << "Illegal error code #" << severity << ": "; + os << "Illegal error severity #" << severity << ": "; } - char errStr[512]; + char errStr[2000]; #if HAVE_VSNPRINTF vsnprintf (errStr, sizeof(errStr), msg, arg); @@ -100,57 +111,59 @@ void sys_verror (int severity, const char *msg, va_list arg) #else strncpy (errStr, sizeof(errStr), "Error message not available on this platform."); #endif - - std::cout << errStr << std::endl; + + os << errStr << std::endl; + strOutput = os.str(); if (severity == ERR_FATAL) - throw std::runtime_error (errStr); + throw std::runtime_error (errStr); #if INTERACTIVE_ERROR_DISPLAY std::cout << "A - Abort C - Continue W - Turn off warnings? "; // fflush(stderr); do - { - int c = cio_kb_waitc("AaBbCcWw", TRUE); /* get code from keyboard */ - c = tolower (c); - fputc (c, stderr); - fputc (NEWLINE, stderr); - - if (c == 'a') - exit (1); - else if (c == 'c') - return; - else if (c == 'w') - { - sys_error_level (ERR_SEVERE); /* report severe & fatal errors */ - break; - } - } while (TRUE); + { + int c = cio_kb_waitc("AaBbCcWw", TRUE); /* get code from keyboard */ + c = tolower (c); + fputc (c, stderr); + fputc (NEWLINE, stderr); + + if (c == 'a') + exit (1); + else if (c == 'c') + return; + else if (c == 'w') + { + sys_error_level (ERR_SEVERE); /* report severe & fatal errors */ + break; + } + } while (TRUE); #endif } /* NAME - * sys_error_level Set error reporting level - * - * SYNOPSIS - * sys_error_level (severity) - * int severity Report all error as serious as severity and beyond - * - * DESCRIPTION - * Causes the system to ignore all error below the level of severity - * For example, if severity == ERR_SEVERE, then report severe & fatal - * error and ignore warnings - */ +* sys_error_level Set error reporting level +* +* SYNOPSIS +* sys_error_level (severity) +* int severity Report all error as serious as severity and beyond +* +* DESCRIPTION +* Causes the system to ignore all error below the level of severity +* For example, if severity == ERR_SEVERE, then report severe & fatal +* error and ignore warnings +*/ void sys_error_level (int severity) { if (severity == ERR_FATAL || - severity == ERR_SEVERE || - severity == ERR_WARNING) - errorlevel = severity; + severity == ERR_SEVERE || + severity == ERR_WARNING || + severity == ERR_TRACE) + s_reportErrorLevel = severity; else - errorlevel = ERR_WARNING; + s_reportErrorLevel = ERR_WARNING; } diff --git a/msvc/ctsim/ctsim.plg b/msvc/ctsim/ctsim.plg index d695bc3..1487b05 100644 --- a/msvc/ctsim/ctsim.plg +++ b/msvc/ctsim/ctsim.plg @@ -6,13 +6,13 @@ --------------------Configuration: libctsim - Win32 Debug--------------------

Command Lines

-Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP7B.tmp" with contents +Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP34.tmp" with contents [ /nologo /G6 /MTd /W3 /Gm /GR /GX /ZI /Od /I "..\..\..\zlib" /I "..\..\INCLUDE" /I "..\..\getopt" /I "..\..\..\lpng108" /I "..\..\..\fftw-2.1.3\fftw" /I "..\..\..\fftw-2.1.3\rfftw" /I "..\..\..\wx2\include" /D "_DEBUG" /D "HAVE_WXWIN" /D "HAVE_STRING_H" /D "HAVE_GETOPT_H" /D "WIN32" /D "_MBCS" /D "_LIB" /D "MSVC" /D "HAVE_FFTW" /D "HAVE_PNG" /D "HAVE_SGP" /D "HAVE_WXWINDOWS" /D "__WXMSW__" /FR"Debug/" /Fp"Debug/libctsim.pch" /YX /Fo"Debug/" /Fd"Debug/" /FD /GZ /c -"C:\ctsim-2.0.6\libctgraphics\pol.cpp" +"C:\ctsim-2.0.6\libctsim\imagefile.cpp" ] -Creating command line "cl.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP7B.tmp" -Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP7C.tmp" with contents +Creating command line "cl.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP34.tmp" +Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP35.tmp" with contents [ /nologo /out:"Debug\libctsim.lib" ".\Debug\array2dfile.obj" @@ -26,6 +26,7 @@ Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP7C.tmp" with content ".\Debug\fnetorderstream.obj" ".\Debug\getopt.obj" ".\Debug\getopt1.obj" +".\Debug\hashtable.obj" ".\Debug\imagefile.obj" ".\Debug\mathfuncs.obj" ".\Debug\phantom.obj" @@ -41,24 +42,17 @@ Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP7C.tmp" with content ".\Debug\trace.obj" ".\Debug\transformmatrix.obj" ".\Debug\xform.obj" -".\Debug\hashtable.obj" ] -Creating command line "link.exe -lib @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP7C.tmp" +Creating command line "link.exe -lib @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP35.tmp"

Output Window

Compiling... -pol.cpp +imagefile.cpp Creating library...

--------------------Configuration: ctsim - Win32 Debug--------------------

Command Lines

-Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP7D.tmp" with contents -[ -/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=\"2.5.0\" /FR"Debug/" /Fp"Debug/ctsim.pch" /YX /Fo"Debug/" /Fd"Debug/" /FD /GZ /c -"C:\ctsim-2.0.6\src\ctsim.cpp" -] -Creating command line "cl.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP7D.tmp" -Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP7E.tmp" with contents +Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP36.tmp" with contents [ 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" @@ -73,14 +67,9 @@ comctl32.lib winmm.lib rpcrt4.lib ws2_32.lib kernel32.lib user32.lib gdi32.lib w "\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\RSP7E.tmp" +Creating command line "link.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP36.tmp"

Output Window

-Compiling... -ctsim.cpp Linking... -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... -

Output Window

diff --git a/src/ctsim.h b/src/ctsim.h index 24ef0a2..bbd7f83 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.13 2000/12/25 21:54:26 kevin Exp $ +** $Id: ctsim.h,v 1.14 2000/12/29 15:45:06 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 @@ -147,6 +147,13 @@ enum { IFMENU_COMPARE_IMAGES, IFMENU_COMPARE_ROW, IFMENU_COMPARE_COL, + IFMENU_PROCESS_INVERTVALUES, + IFMENU_PROCESS_SQRT, + IFMENU_PROCESS_SQUARE, + IFMENU_PROCESS_LOG, + IFMENU_PROCESS_EXP, + IFMENU_PROCESS_FFT_MAGNITUDE, + IFMENU_PROCESS_FFT_PHASE, PHMMENU_PROCESS_RASTERIZE, PHMMENU_PROCESS_PROJECTIONS, PLOTMENU_VIEW_SCALE_MINMAX, diff --git a/src/dlgprojections.cpp b/src/dlgprojections.cpp index a169de8..489d912 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.13 2000/12/18 12:29:41 kevin Exp $ +** $Id: dlgprojections.cpp,v 1.14 2000/12/29 15:45:06 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,13 +163,13 @@ ProjectionsDialog::showView (int iViewNumber) double* detPos = new double [detArray.nDet()]; for (int i = 0; i < detArray.nDet(); i++) detPos[i] = i; - EZPlot ezplot (*m_pSGP); - ezplot.ezset("grid"); - ezplot.ezset("box"); - ezplot.ezset("yticks left"); + EZPlot ezplot; + ezplot.ezset ("grid"); + ezplot.ezset ("box"); + ezplot.ezset ("yticks left"); ezplot.addCurve (detValues, detPos, detArray.nDet()); m_pSGP->setViewport (0.67, 0.1, 1., 1.); - ezplot.plot(); + ezplot.plot (m_pSGP); delete detPos; } } diff --git a/src/docs.cpp b/src/docs.cpp index ac7e3a6..bab1d13 100644 --- a/src/docs.cpp +++ b/src/docs.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: docs.cpp,v 1.8 2000/12/27 03:16:02 kevin Exp $ +** $Id: docs.cpp,v 1.9 2000/12/29 15:45:06 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 @@ -81,9 +81,7 @@ bool ImageFileDocument::OnOpenDocument(const wxString& filename) } Modify(false); UpdateAllViews(); - ImageFileView* ifView = dynamic_cast(GetFirstView()); - if (ifView) - ifView->OnUpdate(ifView, NULL); + GetFirstView()->OnUpdate (GetFirstView(), NULL); return true; } @@ -133,7 +131,8 @@ bool ProjectionFileDocument::OnOpenDocument(const wxString& filename) SetFilename(filename, true); } Modify(false); - UpdateAllViews(); + UpdateAllViews(); + return true; } @@ -176,7 +175,8 @@ bool PhantomDocument::OnOpenDocument(const wxString& filename) } m_idPhantom = m_phantom.id(); Modify(false); - UpdateAllViews(); + UpdateAllViews(); + return true; } @@ -224,14 +224,11 @@ bool PlotFileDocument::OnOpenDocument(const wxString& filename) return false; } *theApp->getLog() << "Read plot file " << filename << "\n"; - SetFilename(filename, true); + SetFilename (filename, true); m_namePlot = filename.c_str(); } - Modify(false); + Modify (false); UpdateAllViews(); - PlotFileView* ifView = dynamic_cast(GetFirstView()); - if (ifView) - ifView->OnUpdate(ifView, NULL); return true; } diff --git a/src/views.cpp b/src/views.cpp index 4659903..8c53e4c 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.39 2000/12/27 03:16:02 kevin Exp $ +** $Id: views.cpp,v 1.40 2000/12/29 15:45:06 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 @@ -164,6 +164,13 @@ EVT_MENU(IFMENU_VIEW_SCALE_AUTO, ImageFileView::OnScaleAuto) EVT_MENU(IFMENU_COMPARE_IMAGES, ImageFileView::OnCompare) EVT_MENU(IFMENU_COMPARE_ROW, ImageFileView::OnCompareRow) EVT_MENU(IFMENU_COMPARE_COL, ImageFileView::OnCompareCol) +EVT_MENU(IFMENU_PROCESS_INVERTVALUES, ImageFileView::OnInvertValues) +EVT_MENU(IFMENU_PROCESS_SQUARE, ImageFileView::OnSquare) +EVT_MENU(IFMENU_PROCESS_SQRT, ImageFileView::OnSquareRoot) +EVT_MENU(IFMENU_PROCESS_LOG, ImageFileView::OnLog) +EVT_MENU(IFMENU_PROCESS_EXP, ImageFileView::OnExp) +EVT_MENU(IFMENU_PROCESS_FFT_MAGNITUDE, ImageFileView::OnFFTMagnitude) +EVT_MENU(IFMENU_PROCESS_FFT_PHASE, ImageFileView::OnFFTPhase) EVT_MENU(IFMENU_PLOT_ROW, ImageFileView::OnPlotRow) EVT_MENU(IFMENU_PLOT_COL, ImageFileView::OnPlotCol) END_EVENT_TABLE() @@ -245,15 +252,15 @@ ImageFileView::OnCompare (wxCommandEvent& event) { std::vector vecIF; theApp->getCompatibleImages (GetDocument(), vecIF); - + if (vecIF.size() == 0) { wxMessageBox("There are no compatible image files open for comparision", "No comparison images"); } else { DialogGetComparisonImage dialogGetCompare(m_frame, "Get Comparison Image", vecIF, true); - + if (dialogGetCompare.ShowModal() == wxID_OK) { - ImageFileDocument* pCompareDoc = dialogGetCompare.getImageFileDocument(); const ImageFile& rIF = GetDocument()->getImageFile(); + ImageFileDocument* pCompareDoc = dialogGetCompare.getImageFileDocument(); const ImageFile& rCompareIF = pCompareDoc->getImageFile(); std::ostringstream os; double min, max, mean, mode, median, stddev; @@ -273,13 +280,13 @@ ImageFileView::OnCompare (wxCommandEvent& event) return; } ImageFile& differenceImage = pDifferenceDoc->getImageFile(); - + differenceImage.setArraySize (rIF.nx(), rIF.ny()); if (! rIF.subtractImages (rCompareIF, differenceImage)) { pDifferenceDoc->DeleteAllViews(); return; } - + if (theApp->getSetModifyNewDocs()) pDifferenceDoc->Modify(true); pDifferenceDoc->UpdateAllViews(this); @@ -290,6 +297,77 @@ ImageFileView::OnCompare (wxCommandEvent& event) } } +void +ImageFileView::OnInvertValues (wxCommandEvent& event) +{ + ImageFile& rIF = GetDocument()->getImageFile(); + rIF.invertPixelValues (rIF); + if (theApp->getSetModifyNewDocs()) + GetDocument()->Modify(TRUE); + GetDocument()->UpdateAllViews(this); +} + +void +ImageFileView::OnSquare (wxCommandEvent& event) +{ + ImageFile& rIF = GetDocument()->getImageFile(); + rIF.square (rIF); + if (theApp->getSetModifyNewDocs()) + GetDocument()->Modify(TRUE); + GetDocument()->UpdateAllViews(this); +} + +void +ImageFileView::OnSquareRoot (wxCommandEvent& event) +{ + ImageFile& rIF = GetDocument()->getImageFile(); + rIF.sqrt (rIF); + if (theApp->getSetModifyNewDocs()) + GetDocument()->Modify(TRUE); + GetDocument()->UpdateAllViews(this); +} + +void +ImageFileView::OnLog (wxCommandEvent& event) +{ + ImageFile& rIF = GetDocument()->getImageFile(); + rIF.log (rIF); + if (theApp->getSetModifyNewDocs()) + GetDocument()->Modify(TRUE); + GetDocument()->UpdateAllViews(this); +} + +void +ImageFileView::OnExp (wxCommandEvent& event) +{ + ImageFile& rIF = GetDocument()->getImageFile(); + rIF.exp (rIF); + if (theApp->getSetModifyNewDocs()) + GetDocument()->Modify(TRUE); + GetDocument()->UpdateAllViews(this); +} + +void +ImageFileView::OnFFTMagnitude (wxCommandEvent& event) +{ + ImageFile& rIF = GetDocument()->getImageFile(); + rIF.FFTMagnitude (rIF); + if (theApp->getSetModifyNewDocs()) + GetDocument()->Modify(TRUE); + GetDocument()->UpdateAllViews(this); +} + +void +ImageFileView::OnFFTPhase (wxCommandEvent& event) +{ + ImageFile& rIF = GetDocument()->getImageFile(); + rIF.FFTPhase (rIF); + if (theApp->getSetModifyNewDocs()) + GetDocument()->Modify(TRUE); + GetDocument()->UpdateAllViews(this); +} + + ImageFileCanvas* ImageFileView::CreateCanvas (wxView *view, wxFrame *parent) { @@ -331,6 +409,15 @@ ImageFileView::CreateChildFrame(wxDocument *doc, wxView *view) view_menu->Append(IFMENU_VIEW_SCALE_MINMAX, "Display Scale &Set..."); view_menu->Append(IFMENU_VIEW_SCALE_AUTO, "Display Scale &Auto..."); + wxMenu* process_menu = new wxMenu; + process_menu->Append (IFMENU_PROCESS_INVERTVALUES, "&Invert Values"); + process_menu->Append (IFMENU_PROCESS_SQUARE, "&Square"); + process_menu->Append (IFMENU_PROCESS_SQRT, "Square &Root"); + process_menu->Append (IFMENU_PROCESS_LOG, "&Log"); + process_menu->Append (IFMENU_PROCESS_EXP, "&Exp"); + process_menu->Append (IFMENU_PROCESS_FFT_MAGNITUDE, "&FFT Magnitude"); + process_menu->Append (IFMENU_PROCESS_FFT_PHASE, "FFT &Phase"); + wxMenu *plot_menu = new wxMenu; plot_menu->Append (IFMENU_PLOT_ROW, "Plot &Row"); plot_menu->Append (IFMENU_PLOT_COL, "Plot &Column"); @@ -339,7 +426,7 @@ ImageFileView::CreateChildFrame(wxDocument *doc, wxView *view) compare_menu->Append (IFMENU_COMPARE_IMAGES, "Compare &Images..."); compare_menu->Append (IFMENU_COMPARE_ROW, "Compare &Row"); compare_menu->Append (IFMENU_COMPARE_COL, "Compare &Column"); - + wxMenu *help_menu = new wxMenu; help_menu->Append(MAINMENU_HELP_ABOUT, "&About"); @@ -347,7 +434,8 @@ ImageFileView::CreateChildFrame(wxDocument *doc, wxView *view) menu_bar->Append(file_menu, "&File"); menu_bar->Append(view_menu, "&View"); - menu_bar->Append(plot_menu, "&Plot"); + menu_bar->Append(process_menu, "&Process"); + menu_bar->Append(plot_menu, "P&lot"); menu_bar->Append(compare_menu, "&Compare"); menu_bar->Append(help_menu, "&Help"); @@ -507,9 +595,9 @@ ImageFileView::OnPlotRow (wxCommandEvent& event) delete pX; delete pY; if (theApp->getSetModifyNewDocs()) - pPlotDoc->Modify(true); + pPlotDoc->Modify(true); + pPlotDoc->UpdateAllViews(); } - } void @@ -556,6 +644,7 @@ ImageFileView::OnPlotCol (wxCommandEvent& event) delete pY; if (theApp->getSetModifyNewDocs()) pPlotDoc->Modify(true); + pPlotDoc->UpdateAllViews(); } } @@ -625,6 +714,7 @@ ImageFileView::OnCompareCol (wxCommandEvent& event) delete pY2; if (theApp->getSetModifyNewDocs()) pPlotDoc->Modify(true); + pPlotDoc->UpdateAllViews(); } } } @@ -645,7 +735,7 @@ ImageFileView::OnCompareRow (wxCommandEvent& event) wxMessageBox ("No compatible images for Row Comparison", "Error"); return; } - + DialogGetComparisonImage dialogGetCompare (m_frame, "Get Comparison Image", vecIFDoc, false); if (dialogGetCompare.ShowModal() == wxID_OK) { @@ -697,6 +787,7 @@ ImageFileView::OnCompareRow (wxCommandEvent& event) delete pY2; if (theApp->getSetModifyNewDocs()) pPlotDoc->Modify(true); + pPlotDoc->UpdateAllViews(); } } } @@ -1371,12 +1462,14 @@ EVT_MENU(PLOTMENU_VIEW_SCALE_AUTO, PlotFileView::OnScaleAuto) END_EVENT_TABLE() PlotFileView::PlotFileView(void) -: wxView(), m_canvas(NULL), m_frame(NULL) +: wxView(), m_canvas(NULL), m_frame(NULL), m_pEZPlot(NULL) { } PlotFileView::~PlotFileView(void) -{ +{ + if (m_pEZPlot) + delete m_pEZPlot; } void @@ -1502,7 +1595,7 @@ PlotFileView::CreateChildFrame(wxDocument *doc, wxView *view) bool -PlotFileView::OnCreate(wxDocument *doc, long WXUNUSED(flags) ) +PlotFileView::OnCreate (wxDocument *doc, long WXUNUSED(flags) ) { m_frame = CreateChildFrame(doc, this); SetFrame(m_frame); @@ -1524,7 +1617,7 @@ PlotFileView::OnCreate(wxDocument *doc, long WXUNUSED(flags) ) m_frame->Show(true); Activate(true); - + return true; } @@ -1540,48 +1633,58 @@ PlotFileView::OnDraw (wxDC* dc) m_canvas->GetClientSize (&xsize, &ysize); SGPDriver driver (dc, xsize, ysize); SGP sgp (driver); - EZPlot plot (sgp); - - for (int iEzset = 0; iEzset < rPlotFile.getNumEzsetCommands(); iEzset++) - plot.ezset (rPlotFile.getEzsetCommand (iEzset)); - - if (m_bMinSpecified) { - std::ostringstream os; - os << "ymin " << m_dMinPixel; - plot.ezset (os.str()); - } - - if (m_bMaxSpecified) { - std::ostringstream os; - os << "ymax " << m_dMaxPixel; - plot.ezset (os.str()); - } - - plot.ezset("box"); - plot.ezset("grid"); - - double* pdXaxis = new double [iNRecords]; - rPlotFile.getColumn (0, pdXaxis); - - double* pdY = new double [iNRecords]; - for (int iCol = 1; iCol < iNColumns; iCol++) { - rPlotFile.getColumn (iCol, pdY); - plot.addCurve (pdXaxis, pdY, iNRecords); - } - - delete pdXaxis; - delete pdY; - - plot.plot(); + if (m_pEZPlot) + m_pEZPlot->plot (&sgp); } } void -PlotFileView::OnUpdate(wxView *WXUNUSED(sender), wxObject *WXUNUSED(hint) ) +PlotFileView::OnUpdate (wxView *WXUNUSED(sender), wxObject *WXUNUSED(hint) ) { - if (m_canvas) - m_canvas->Refresh(); + if (m_canvas) { + const PlotFile& rPlotFile = GetDocument()->getPlotFile(); + const int iNColumns = rPlotFile.getNumColumns(); + const int iNRecords = rPlotFile.getNumRecords(); + + if (iNColumns > 0 && iNRecords > 0) { + if (m_pEZPlot) + delete m_pEZPlot; + m_pEZPlot = new EZPlot; + + for (int iEzset = 0; iEzset < rPlotFile.getNumEzsetCommands(); iEzset++) + m_pEZPlot->ezset (rPlotFile.getEzsetCommand (iEzset)); + + if (m_bMinSpecified) { + std::ostringstream os; + os << "ymin " << m_dMinPixel; + m_pEZPlot->ezset (os.str()); + } + + if (m_bMaxSpecified) { + std::ostringstream os; + os << "ymax " << m_dMaxPixel; + m_pEZPlot->ezset (os.str()); + } + + m_pEZPlot->ezset("box"); + m_pEZPlot->ezset("grid"); + + double* pdXaxis = new double [iNRecords]; + rPlotFile.getColumn (0, pdXaxis); + + double* pdY = new double [iNRecords]; + for (int iCol = 1; iCol < iNColumns; iCol++) { + rPlotFile.getColumn (iCol, pdY); + m_pEZPlot->addCurve (pdXaxis, pdY, iNRecords); + } + + delete pdXaxis; + delete pdY; + } + + m_canvas->Refresh(); + } } bool diff --git a/src/views.h b/src/views.h index 1457e84..8bb5e6d 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.16 2000/12/22 04:18:00 kevin Exp $ +** $Id: views.h,v 1.17 2000/12/29 15:45:06 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 @@ -66,6 +66,13 @@ public: bool OnClose (bool deleteWindow = true); void OnProperties (wxCommandEvent& event); void OnCompare (wxCommandEvent& event); + void OnInvertValues (wxCommandEvent& event); + void OnSquare (wxCommandEvent& event); + void OnSquareRoot (wxCommandEvent& event); + void OnLog (wxCommandEvent& event); + void OnExp (wxCommandEvent& event); + void OnFFTMagnitude (wxCommandEvent& event); + void OnFFTPhase (wxCommandEvent& event); void OnScaleAuto (wxCommandEvent& event); void OnScaleMinMax (wxCommandEvent& event); void OnPlotRow (wxCommandEvent& event); @@ -217,6 +224,7 @@ private: PlotFileCanvas *CreateCanvas(wxView *view, wxFrame *parent); wxFrame *CreateChildFrame(wxDocument *doc, wxView *view); + EZPlot* m_pEZPlot; PlotFileCanvas *m_canvas; wxFrame *m_frame; bool m_bMinSpecified; diff --git a/tools/if1.cpp b/tools/if1.cpp index dcd1dc1..9e4c169 100644 --- a/tools/if1.cpp +++ b/tools/if1.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: if1.cpp,v 1.2 2000/12/23 18:12:35 kevin Exp $ +** $Id: if1.cpp,v 1.3 2000/12/29 15:45:06 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 @@ -46,7 +46,7 @@ static struct option my_options[] = {0, 0, 0, 0} }; -static const char* g_szIdStr = "$Id: if1.cpp,v 1.2 2000/12/23 18:12:35 kevin Exp $"; +static const char* g_szIdStr = "$Id: if1.cpp,v 1.3 2000/12/29 15:45:06 kevin Exp $"; void if1_usage (const char *program) @@ -67,8 +67,6 @@ if1_usage (const char *program) int if1_main (int argc, char *const argv[]) { - ImageFile *im_in; - ImageFile *im_out; char *in_file; char *out_file; int opt_verbose = 0; @@ -131,61 +129,39 @@ if1_main (int argc, char *const argv[]) in_file = argv[optind]; out_file = argv[optind + 1]; - std::string histString; if (opt_invert || opt_log || opt_exp || opt_sqr || opt_sqrt) { - int ix, iy; - - im_in = new ImageFile (); - im_in->fileRead (in_file); - int nx = im_in->nx(); - int ny = im_in->ny(); - im_out = new ImageFile (nx, ny); - - ImageFileArray vIn = im_in->getArray(); - ImageFileArray vOut = im_out->getArray(); - - if (opt_invert) { - for (ix = 0; ix < nx; ix++) - for (iy = 0; iy < ny; iy++) - vOut[ix][iy] = - vIn[ix][iy]; + ImageFile im_in; + im_in.fileRead (in_file); + int nx = im_in.nx(); + int ny = im_in.ny(); + ImageFile im_out (nx, ny); + + if (opt_invert) { + im_in.invertPixelValues (im_out); histString = "Invert transformation"; } - if (opt_log) { - for (ix = 0; ix < nx; ix++) - for (iy = 0; iy < ny; iy++) - if (vIn[ix][iy] < 0) - vOut[ix][iy] = 0; - else - vOut[ix][iy] = log (vIn[ix][iy]); + if (opt_log) { + im_in.log (im_out); histString = "Logrithmic transformation"; } - if (opt_exp) { - for (ix = 0; ix < nx; ix++) - for (iy = 0; iy < ny; iy++) - vOut[ix][iy] = exp (vIn[ix][iy]); + if (opt_exp) { + im_in.exp (im_out); histString = "Exponential transformation"; } - if (opt_sqr) { - for (ix = 0; ix < nx; ix++) - for (iy = 0; iy < ny; iy++) - vOut[ix][iy] = vIn[ix][iy] * vIn[ix][iy]; + if (opt_sqr) { + im_in.square (im_out); histString = "Square transformation"; } - if (opt_sqrt) { - for (ix = 0; ix < nx; ix++) - for (iy = 0; iy < ny; iy++) - if (vIn[ix][iy] < 0) - vOut[ix][iy] = sqrt (-vIn[ix][iy]); - else - vOut[ix][iy] = sqrt (vIn[ix][iy]); + if (opt_sqrt) { + im_in.sqrt (im_out); histString = "Square root transformation"; } - im_out->labelsCopy (*im_in); - im_out->labelAdd (Array2dFileLabel::L_HISTORY, histString.c_str()); - im_out->fileWrite (out_file); + im_out.labelsCopy (im_in); + im_out.labelAdd (Array2dFileLabel::L_HISTORY, histString.c_str()); + im_out.fileWrite (out_file); } return (0); diff --git a/tools/if2.cpp b/tools/if2.cpp index 571e8d6..0f36b59 100644 --- a/tools/if2.cpp +++ b/tools/if2.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: if2.cpp,v 1.5 2000/12/23 18:12:35 kevin Exp $ +** $Id: if2.cpp,v 1.6 2000/12/29 15:45:06 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 @@ -46,7 +46,7 @@ static struct option my_options[] = {0, 0, 0, 0} }; -static const char* g_szIdStr = "$Id: if2.cpp,v 1.5 2000/12/23 18:12:35 kevin Exp $"; +static const char* g_szIdStr = "$Id: if2.cpp,v 1.6 2000/12/29 15:45:06 kevin Exp $"; void if2_usage (const char *program) @@ -254,7 +254,7 @@ if2_main (int argc, char *const argv[]) #if HAVE_SGP SGPDriver driver ("Column Plot"); SGP sgp (driver); - EZPlot ezplot (sgp); + EZPlot ezplot; ezplot.ezset ("clear."); ezplot.ezset ("xticks major 5."); ezplot.ezset ("xlabel Column"); @@ -264,7 +264,7 @@ if2_main (int argc, char *const argv[]) ezplot.ezset ("grid."); ezplot.addCurve (plot_xaxis, v1[opt_columnPlot], im_in1.ny()); ezplot.addCurve (plot_xaxis, v2[opt_columnPlot], im_in2.ny()); - ezplot.plot(); + ezplot.plot (&sgp); std::cout << "Press enter to continue" << flush; cio_kb_getc(); #endif @@ -309,7 +309,7 @@ if2_main (int argc, char *const argv[]) #if HAVE_SGP SGPDriver driver ("Row Plot"); SGP sgp (driver); - EZPlot ezplot (sgp); + EZPlot ezplot; ezplot.ezset ("clear."); ezplot.ezset ("xticks major 5."); ezplot.ezset ("title Row Plot"); @@ -319,7 +319,7 @@ if2_main (int argc, char *const argv[]) ezplot.ezset ("grid."); ezplot.addCurve (plot_xaxis, v1Row, im_in1.nx()); ezplot.addCurve (plot_xaxis, v2Row, im_in2.nx()); - ezplot.plot(); + ezplot.plot (&sgp); std::cout << "Press enter to continue" << flush; cio_kb_getc(); #endif -- 2.34.1