** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: fourier.h,v 1.6 2001/03/18 18:08:25 kevin Exp $
+** $Id: fourier.h,v 1.7 2001/03/21 21:45:31 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
static void shuffleFourierToNaturalOrder (ImageFile& im);
static void shuffleNaturalToFourierOrder (ImageFile& im);
+// Odd Number of Points
+// Natural Frequency Order: -(n-1)/2...-1,0,1...(n-1)/2
+// Fourier Frequency Order: 0, 1..(n-1)/2,-(n-1)/2...-1
+// Even Number of Points
+// Natural Frequency Order: -n/2...-1,0,1...((n/2)-1)
+// Fourier Frequency Order: 0,1...((n/2)-1),-n/2...-1
+ template<class T>
+ static void shuffleNaturalToFourierOrder (T* pVector, const int n)
+ {
+ T* pTemp = new T [n];
+ int i;
+ if (isOdd(n)) { // Odd
+ int iHalfN = (n - 1) / 2;
+
+ pTemp[0] = pVector[iHalfN];
+ for (i = 0; i < iHalfN; i++)
+ pTemp[i + 1] = pVector[i + 1 + iHalfN];
+ for (i = 0; i < iHalfN; i++)
+ pTemp[i + iHalfN + 1] = pVector[i];
+ } else { // Even
+ int iHalfN = n / 2;
+ pTemp[0] = pVector[iHalfN];
+ for (i = 0; i < iHalfN - 1; i++)
+ pTemp[i + 1] = pVector[i + iHalfN + 1];
+ for (i = 0; i < iHalfN; i++)
+ pTemp[i + iHalfN] = pVector[i];
+ }
+
+ for (i = 0; i < n; i++)
+ pVector[i] = pTemp[i];
+ delete pTemp;
+ }
+
+ template<class T>
+ static void shuffleFourierToNaturalOrder (T* pVector, const int n)
+ {
+ T* pTemp = new T [n];
+ int i;
+ if (isOdd(n)) { // Odd
+ int iHalfN = (n - 1) / 2;
+
+ pTemp[iHalfN] = pVector[0];
+ for (i = 0; i < iHalfN; i++)
+ pTemp[i + 1 + iHalfN] = pVector[i + 1];
+ for (i = 0; i < iHalfN; i++)
+ pTemp[i] = pVector[i + iHalfN + 1];
+ } else { // Even
+ int iHalfN = n / 2;
+ pTemp[iHalfN] = pVector[0];
+ for (i = 0; i < iHalfN; i++)
+ pTemp[i] = pVector[i + iHalfN];
+ for (i = 0; i < iHalfN - 1; i++)
+ pTemp[i + iHalfN + 1] = pVector[i+1];
+ }
+
+ for (i = 0; i < n; i++)
+ pVector[i] = pTemp[i];
+ delete pTemp;
+ }
+#if 0
static void shuffleNaturalToFourierOrder (float* pdVector, const int n);
static void shuffleNaturalToFourierOrder (double* pdVector, const int n);
static void shuffleNaturalToFourierOrder (std::complex<double>* pdVector, const int n);
static void shuffleFourierToNaturalOrder (float* pdVector, const int n);
static void shuffleFourierToNaturalOrder (double* pdVector, const int n);
static void shuffleFourierToNaturalOrder (std::complex<double>* pdVector, const int n);
+#endif
};
** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: imagefile.h,v 1.34 2001/03/18 18:08:25 kevin Exp $
+** $Id: imagefile.h,v 1.35 2001/03/21 21:45:31 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
#endif
bool magnitude (ImageFile& result) const;
bool phase (ImageFile& result) const;
-
- int display (void) const;
- int displayScaling (const int scaleFactor, ImageFileValue pmin, ImageFileValue pmax) const;
+ bool real (ImageFile& result) const;
+ bool imaginary (ImageFile& result) const;
bool exportImage (const char* const pszFormat, const char* const pszFilename, int nxcell, int nycell, double densmin, double densmax);
** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: interpolator.h,v 1.1 2001/02/11 21:57:08 kevin Exp $
+** $Id: interpolator.h,v 1.2 2001/03/21 21:45:31 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
double interpolate (double x);
};
+
+template<class T>
+class BilinearInterpolator {
+private:
+ T** const m_ppMatrix;
+ const unsigned int m_nx;
+ const unsigned int m_ny;
+
+public:
+ BilinearInterpolator (T** ppMatrix, unsigned int nx, unsigned int ny)
+ : m_ppMatrix(ppMatrix), m_nx(nx), m_ny(ny)
+ {}
+
+ T interpolate (double dXPos, double dYPos)
+ {
+ int iFloorX = floor (dXPos);
+ int iFloorY = floor (dYPos);
+ double dXFrac = dXPos - iFloorX;
+ double dYFrac = dYPos - iFloorY;
+
+ T result = 0;
+
+ if (iFloorX < 0 || iFloorY < 0 || iFloorX > m_nx-1 || iFloorY > m_ny-1)
+ result = 0;
+ else if (iFloorX == m_nx - 1 && iFloorY == m_ny - 1)
+ result = m_ppMatrix[m_nx-1][m_ny-1];
+ else if (iFloorX == m_nx - 1)
+ result = m_ppMatrix[iFloorX][iFloorY] + dYFrac * (m_ppMatrix[iFloorX][iFloorY+1] - m_ppMatrix[iFloorX][iFloorY]);
+ else if (iFloorY == m_ny - 1)
+ result = m_ppMatrix[iFloorX][iFloorY] + dXFrac * (m_ppMatrix[iFloorX+1][iFloorY] - m_ppMatrix[iFloorX][iFloorY]);
+ else
+ result = (1 - dXFrac) * (1 - dYFrac) * m_ppMatrix[iFloorX][iFloorY] +
+ dXFrac * (1 - dYFrac) * m_ppMatrix[iFloorX+1][iFloorY] +
+ dYFrac * (1 - dXFrac) * m_ppMatrix[iFloorX][iFloorY+1] +
+ dXFrac * dYFrac * m_ppMatrix[iFloorX+1][iFloorY+1];
+
+ return result;
+ }
+};
+
+
+template<class T>
+class LinearInterpolator {
+private:
+ T* const m_pY;
+ T* const m_pX;
+ const unsigned int m_n;
+
+public:
+ LinearInterpolator (T* pY, unsigned int n)
+ : m_pY(pY), m_pX(0), m_n(n)
+ {}
+
+ LinearInterpolator (T* pY, T* pX, unsigned int n)
+ : m_pY(pY), m_pX(pX), m_n(n)
+ {}
+
+ double interpolate (double dX, int* piLastFloor = NULL)
+ {
+ double result = 0;
+
+ if (! m_pX) {
+ if (dX == 0)
+ result = m_pY[0];
+ else if (dX < 0)
+ result = 0;
+ else if (dX == m_n - 1)
+ result = m_pY[m_n-1];
+ else if (dX > m_n - 1)
+ result = 0;
+ else {
+ int iFloor = floor(dX);
+ result = m_pY[iFloor] + (m_pY[iFloor+1] - m_pY[iFloor]) * (dX - iFloor);
+ }
+ } else {
+ int iLower = -1;
+ int iUpper = m_n;
+ if (piLastFloor && *piLastFloor >= 0 && m_pX[*piLastFloor] < dX)
+ iLower = *piLastFloor;
+
+ while (iUpper - iLower > 1) {
+ int iMiddle = (iUpper + iLower) >> 1;
+ if (dX >= m_pX[iMiddle])
+ iLower = iMiddle;
+ else
+ iUpper = iMiddle;
+ }
+ if (dX == m_pX[0])
+ result = m_pY[0];
+ else if (dX < m_pX[0])
+ result = 0;
+ else if (dX == m_pX[m_n-1])
+ result = m_pY[m_n-1];
+ else if (dX > m_pX[m_n-1])
+ result = 0;
+ else {
+ if (iLower < 0 || iLower >= m_n) {
+ sys_error (ERR_SEVERE, "Coordinate out of range [linearInterpolate]");
+ return 0;
+ }
+
+ if (piLastFloor)
+ *piLastFloor = iLower;
+ result = m_pY[iLower] + (m_pY[iUpper] - m_pY[iLower]) * ((dX - m_pX[iLower]) / (m_pX[iUpper] - m_pX[iLower]));
+ }
+ }
+
+ return result;
+ }
+};
+
** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: projections.h,v 1.34 2001/03/13 14:53:44 kevin Exp $
+** $Id: projections.h,v 1.35 2001/03/21 21:45:31 kevin Exp $
**
**
** This program is free software; you can redistribute it and/or modify
bool convertPolar (ImageFile& rIF, int iInterpolation);
bool convertFFTPolar (ImageFile& rIF, int iInterpolation, int iZeropad);
void calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double** ppdView, double** ppdDet,
- int iNumDetWithZeros, double dZeropadRatio);
+ int iNumDetWithZeros, double dZeropadRatio, double dDetInc);
void interpolatePolar (ImageFileArray& v, ImageFileArray& vImag, unsigned int nx, unsigned int ny, std::complex<double>** ppcDetValue,
double** ppdDet, double** ppdView, unsigned int nView, unsigned int nDet, unsigned int nDetWithZeros,
int iInterpolate);
-/* ; FILE IDENTIFICATION
-;
-; Name: CRT_LINE.ASM
-; Purpose: Draw lines on screen using Bresenham's integer algorithm
-; Programmer: Kevin Rosenberg
-; Date Started: 1 Jun 84
-; Last Modification: 15 Jul 85, Version 1.0
-;
-; CORRESPONDENCE
-; Kevin Rosenberg
-; 2515 E. 16th St.
-; Newport Beach, CA 92663
-; (714) 642-0119 (voice)
-; (714) 842-6348 (E-Mail on Consultant's Exchange RBBS)
-; Compuserve 73775,1232
-;
-; LANGAUGE
-; IBM macro v1.000000000
-;
-; LANGAUGE INTERFACE
-; Public routines callable from Lattice C v2.0+
-; Works with S & L memory models. Can easily add support for the
-; other models by changing the STK structure below.
-;
-; NEEDED FILES
-; egadrive.inc Definitions, macros, & procedures for driving EGA
-; hidden.inc Defines OPTIONAL public procedure, crt_hdn_line().
-; dos.mac Lattice's file containing memory model macros
-; See inclusion statement for more documentation
-;
-; GRAPHIC MODES SUPPORTED
-; - Standard CGA modes -- BIOS codes 4, 5, & 6 -- are support through
-; the ROM BIOS pixel setting routine. Things could be speeded up
-; by a factor of 7-15 by writing to the video buffer directly.
-; - You may remove support for the CGA completely by changing the
-; EQU of STD_GRF to 0, a slight speed increase in EGA operations
-; is the gain.
-; - All EGA graphic modes are supported with optimized direct pixel
-; accessing. Note: mode 13 (320x200) is bugged, see LOG date 6-20-85
-;
-; LANGUAGE INTERFACE
-; Lattice C v2.0+
-;
-; HISTORY LOG
-; 7-15-85
-; Released in PD as version 1.0
-; 6-27-85
-; Put EGA routines in INCLUDE file EGADRIVE.INC
-; Put Hidden line routines in OPTIONAL INCLUDE file HIDDEN.INC
-; 6-20-85
-; Tested with video modes 13 (320x200) & 14 (640x200). Mode 14 worked
-; ok, mode 13 was missing every other vertical line: I'm probably
-; not setting one of the myriad of registers correctly. Since I'll
-; never use mode 13, I won't investigate it further.
-; 6-05-85
-; Added & updated documentation
-; 2-28-85
-; Added hidden line routines
-; PUBLIC crt_hdn_line plus internal routines
-; (see function summary for more info)
-; 1-20-85
-; Added global variable STD_GRF
-; Will support standard graphic modes if it is set to non-zero value;
-; however, the slow ROM BIOS is used as well as a slight overhead
-; for the EG code. IT IS MUCH FASTER TO USE THE EGA ENHANCED
-; (350 line) GRAPHIC MODES.
-; 1-09-85
-; Convert to support IBM's Enhanced graphics adapter.
-; Remove all optimized code for color graphics adapter (CGA)
-; 7-16-84
-; Changed calls to ROM that set pixels to speedy routines which use
-; incremental address calculation and direct writing the to graphics
-; buffer.
-; 1-12-84
-; File creation:
-; Converted C function, bresline.c, to assembly langauge.
-; Retained C source of Bresenham's integer rastering algorithm
-; as comments.
-;
-; FUNCTION SUMMARY
-;
-; crt_line_page (page)
-; Set graphics page to plot lines on [default = 0]
-;
-; crt_line_style (style)
-; Set line style [Default = solid line]
-;
-; crt_line (x1, y1, x2, y1, color)
-; Plots a line on the screen using Bresenham's integer algorithm
-;
-; crt_hdn_line (x1, y1, x2, y1, color, col_min, col_max)
-; OPTIONAL routine to plots lines using hidden areas
-; INCLUDE file HIDDEN.INC if you wish to use this routine
-;
-; crt_init_hdn (col_min, col_max)
-; Initializes vectors for use with crt_hdn_line()
-;
-;
-; int x1, y1 - starting point
-; int x2, y2 - end point
-; int color - color to draw line, as defined by IBM bios setting
-; int linestyle - 16 bit wide dot setting pattern
-; int col_min[640] Minimum values for each X column (initialize to 349)
-; int col_max[640] Maximum values for each X column (initialize to 0)
-; (Do not display any points between the min & max,
-; and update these values for each new min & max)
-; If a vector == NULL, then do not clip against that
-; boundary.
-;
-; REFERENCES FOR BRESENHAM'S ALGORITHM
-; Newman & Sproll, "Principals of Interactive Computer Graphics",
-; page 25.
-; Foley & van Dam, "Fundementals of Interactive Computer Graphics",
-; page 433
-*/
-
-;------ Bresenham variables
-
- deltx dw ?
- delty dw ?
- absdx dw ?
- absdy dw ?
- min_inc dw ? ; (-1, 1) indicates direction of minor axis
- dinc1 dw ? ; d increment if (d >= 0)
- dinc2 dw ? ; d increment if (d < 0)
-
-
-
-
-;------------------------------------------------------------------------------
-; PROCEDURE
-; crt_line [PUBLIC] Plot a line on the screen using Bresenham's alg
-;
-; SYNOPSIS
-; crt_line (x1, y1, x2, y2, color)
-
-void bresenham (int x1, int y1, int x2, int y2)
-{
- delta_x = x2 - x1;
- dx_abs = (delta_x >= 0 ? delta_x : -delta_x);
-
- delta_y = y2 - y1;
- dy_abs = (delta_y >= 0 ? delta_y : -delta_y);
-
-if (dx_abs > dy_abs)
- bresx();
-else
- bresy();
-}
-
-lastx= x2;
-lasty = y2;
-
-
-;------------------------------------------------------------------------------
-; MACRO NAME
-; BRES
-;
-; FUNCTION
-; Do inner most loop to draw pictures
-;
-; SYNOPSIS
-; BRES P1,P2
-; P1 Major axis & direction
-; P2 Minor axis
-;
-; EXAMPLES
-; BRES INC_COL,ROW
-;
-; INPUT
-; (DI) Address of local routine to call when setting pixel
-; (BX) Decision variable
-; (CX) Number of points to plot on major axis (loop counter)
-; [bp].min_inc Direction of minor axis increments
-; [bp].dinc1 Bresenham decision increments (see references)
-; [bp].dinc2
-; MACROS ADDR_INC_ROW, _DEC_ROW, _INC_COL, _DEC_COL
-;
-; REGISTERS USED BY PIXEL SETTING ROUTINES
-; (SI) Offset in EGA buffer
-; (AH) Bitmask
-;
-; REGISTERS AVAILABLE FOR USE
-; (AL)
-; (DX)
-;
-; CALLS
-; LS_PSET, this routine uses the addr in (DI) to set a point if
-; the linestyle allows.
-;------------------------------------------------------------------------------
-
-BRES MACRO P1,P2
- LOCAL b0,b1
- cmp [bp].min_inc,0 ;if (min_inc >= 0)
- jl b0
- _BRES P1,INC_&P2 ; inc_minor_axis
- jmp b1 ;else
-b0: _BRES P1,DEC_&P2 ; dec_minor_axis
-b1:
- ENDM
-
-
-_BRES MACRO P1,P2
- LOCAL bres0,bres1,bres2
-bres0: ; do {
- LS_PSET ; ls_pset() /* set pix using linestyle */
-
- ADDR_&P1 ; P1 (change major axis)
-
- cmp bx,0 ;
- jge bres1 ; if (d < 0)
- add bx,[bp].dinc2 ; d += dinc2;
- jmp bres2 ; else {
-bres1: ;
- add bx,[bp].dinc1 ; d += dinc1
- ADDR_&P2 ; P2 (change minor axis)
- ; }
-bres2:
- loop bres0 ; } while (--count > 0)
- ENDM
-
-bresx()
+/*****************************************************************************
+** This is part of the CTSim program
+** Copyright (c) 1983-2001 Kevin Rosenberg
+**
+** $Id: bresenham.cpp,v 1.2 2001/03/21 21:45:31 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
+** published by the Free Software Foundation.
+**
+** This program is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+** GNU General Public License for more details.
+**
+** You should have received a copy of the GNU General Public License
+** along with this program; if not, write to the Free Software
+** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+******************************************************************************/
+
+// REFERENCES FOR BRESENHAM'S ALGORITHM
+// Newman & Sproll, "Principals of Interactive Computer Graphics", page 25.
+// Foley & van Dam, "Fundementals of Interactive Computer Graphics", page 433
+
+
+
+void
+bresenham (int x1, int y1, int x2, int y2)
{
- /* draws a line then abs(dx) >= abs(dy) */
-
- int npix = dx_abs + 1;
-
- /* determine direction of minor axis */
-
- int min_inc = 1;
- if (delta_y < 0)
- min_inc = -1;
-
- /* Put decision variable in d */
-
- int d = dy_abs * 2 - dx_abs;
+ int delta_x = x2 - x1;
+ int dx_abs = (delta_x >= 0 ? delta_x : -delta_x);
+
+ int delta_y = y2 - y1;
+ int dy_abs = (delta_y >= 0 ? delta_y : -delta_y);
+
+ // draws a line when abs(dx) >= abs(dy)
+ if (dx_abs > dy_abs) {
+ int count = dx_abs + 1;
+ int major_inc = (x1 <= x2 ? 1 : -1);
+ int min_inc = (delta_y >= 0 ? 1 : -1); // determine direction of minor axis
+
+ int d = dy_abs * 2 - dx_abs; // Put decision variable in d
int dinc1 = (dy_abs - dx_abs) * 2;
int dinc2 = 2 * dy_abs;
-
- if (x1 <= x2)
- bres(INC_COL, ROW);
- else
- breas (DEC_COL, ROW)l
-
- }
-
-
-void bresy()
-{
- /*; bresy [LOCAL] For plotting lines with abs(dy) > abs(sx)
-
- NOTES
- Logic identical to bresx's, substitution x for y and
- vice versa. I separate the into two routines rather
- than having the main loop spending time determining
- which are the major & minor axis.
- */
-
-
+
+ bresx (x1, y1, major_inc, d, d1, d2);
+ } else { //For plotting lines with abs(dy) > abs(sx)
int count = dy_abs + 1;
-
- /* check direction of minor axis */
-
- min_inc = 1;
- if (delta_x < 0)
- min_inc = -1;
-
-
+
+ int major_inc = (y1 <= y2 ? 1 : -1);
+ int min_inc = (delta_x >= 0 ? 1 : -1); // check direction of minor axis
+
int d = dx_abs * 2 - dy_abs;
dinc1 = (dx_abs - dy_abs) * 2;
dinc2 = 2 * dx_abs;
-
- if (y1 <= y2)
- bres(INC_ROW,COL);
- else
- bres(DEL_ROW,COL);
+
+ bresy (x1, y1, major_inc, min_inc, count, d, d1, d2);
+ }
}
+
+static void
+bresx (int x, int y, int majorinc, int minorinc, int count, int d, int d1, int d2)
+{
+ do {
+ setpixel();
+ x += majorinc;
+
+ if (d < 0)
+ d += dinc2;
+ else {
+ d += dinc1;
+ y += min_inc;
+ }
+ } while (--count > 0);
+}
+
+
+static void
+bresy (int x, int y, int majorinc, int minorinc, int count, int d, int d1, int d2)
+{
+ do {
+ setpixel();
+ y += majorinc;
+
+ if (d < 0)
+ d += dinc2;
+ else {
+ d += dinc1;
+ x += min_inc;
+ }
+ } while (--count > 0);
+}
+
+
** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: fourier.cpp,v 1.5 2001/03/18 18:08:25 kevin Exp $
+** $Id: fourier.cpp,v 1.6 2001/03/21 21:45:31 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
}
delete [] pRow;
}
-
-
-// Odd Number of Points
-// Natural Frequency Order: -(n-1)/2...-1,0,1...(n-1)/2
-// Fourier Frequency Order: 0, 1..(n-1)/2,-(n-1)/2...-1
-// Even Number of Points
-// Natural Frequency Order: -n/2...-1,0,1...((n/2)-1)
-// Fourier Frequency Order: 0,1...((n/2)-1),-n/2...-1
-
-void
-Fourier::shuffleNaturalToFourierOrder (double* pdVector, const int n)
-{
- double* pdTemp = new double [n];
- int i;
- if (isOdd(n)) { // 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 - 1; i++)
- pdTemp[i + 1] = pdVector[i + iHalfN + 1];
- for (i = 0; i < iHalfN; i++)
- pdTemp[i + iHalfN] = pdVector[i];
- }
-
- for (i = 0; i < n; i++)
- pdVector[i] = pdTemp[i];
- delete pdTemp;
-}
-
-void
-Fourier::shuffleNaturalToFourierOrder (std::complex<double>* pdVector, const int n)
-{
- std::complex<double>* pdTemp = new std::complex<double> [n];
- int i;
- if (isOdd(n)) { // 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 - 1; i++)
- pdTemp[i + 1] = pdVector[i + iHalfN + 1];
- for (i = 0; i < iHalfN; i++)
- pdTemp[i + iHalfN] = pdVector[i];
- }
-
- for (i = 0; i < n; i++)
- pdVector[i] = pdTemp[i];
- delete [] pdTemp;
-}
-
-
-void
-Fourier::shuffleNaturalToFourierOrder (float* pdVector, const int n)
-{
- float* pdTemp = new float [n];
- int i;
- if (isOdd (n)) { // 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 - 1; i++)
- pdTemp[i + 1] = pdVector[i + iHalfN + 1];
- for (i = 0; i < iHalfN; i++)
- pdTemp[i + iHalfN] = pdVector[i];
- }
-
- for (i = 0; i < n; i++)
- pdVector[i] = pdTemp[i];
- delete pdTemp;
-}
-
-
-
-void
-Fourier::shuffleFourierToNaturalOrder (double* pdVector, const int n)
-{
- double* pdTemp = new double [n];
- int i;
- if (isOdd(n)) { // 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;
-}
-
-
-void
-Fourier::shuffleFourierToNaturalOrder (std::complex<double>* pdVector, const int n)
-{
- std::complex<double>* pdTemp = new std::complex<double> [n];
- int i;
- if (isOdd(n)) { // 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;
-}
-
-
-
-
-void
-Fourier::shuffleFourierToNaturalOrder (float* pVector, const int n)
-{
- float* pTemp = new float [n];
- int i;
- if (isOdd (n)) { // Odd
- int iHalfN = (n - 1) / 2;
-
- pTemp[iHalfN] = pVector[0];
- for (i = 0; i < iHalfN; i++)
- pTemp[i + 1 + iHalfN] = pVector[i + 1];
- for (i = 0; i < iHalfN; i++)
- pTemp[i] = pVector[i + iHalfN + 1];
- } else { // Even
- int iHalfN = n / 2;
- pTemp[iHalfN] = pVector[0];
- for (i = 0; i < iHalfN; i++)
- pTemp[i] = pVector[i + iHalfN];
- for (i = 0; i < iHalfN - 1; i++)
- pTemp[i + iHalfN + 1] = pVector[i+1];
- }
-
- for (i = 0; i < n; i++)
- pVector[i] = pTemp[i];
- delete [] pTemp;
-}
-
-
** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: imagefile.cpp,v 1.41 2001/03/18 18:08:25 kevin Exp $
+** $Id: imagefile.cpp,v 1.42 2001/03/21 21:45:31 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
void
-ImageFile::filterResponse (const char* const domainName, double bw, const char* const filterName, double filt_param, double dInputScale, double dOutputScale)
+ImageFile::filterResponse (const char* const domainName, double bw, const char* const filterName,
+ double filt_param, double dInputScale, double dOutputScale)
{
ImageFileArray v = getArray();
SignalFilter filter (filterName, domainName, bw, filt_param);
}
}
-int
-ImageFile::display (void) const
-{
- double pmin, pmax;
-
- getMinMax (pmin, pmax);
-
- return (displayScaling (1, pmin, pmax));
-}
-
-int
-ImageFile::displayScaling (const int scale, const ImageFileValue pmin, const ImageFileValue pmax) const
-{
- int nx = m_nx;
- int ny = m_ny;
- ImageFileArrayConst v = getArray();
- if (v == NULL || nx == 0 || ny == 0)
- return 0;
-
-#if HAVE_G2_H
- int* pPens = new int [nx * ny * scale * scale ];
-
- double view_scale = 255 / (pmax - pmin);
- int id_X11 = g2_open_X11 (nx * scale, ny * scale);
- int grayscale[256];
- for (int i = 0; i < 256; i++) {
- double cval = i / 255.;
- grayscale[i] = g2_ink (id_X11, cval, cval, cval);
- }
-
- for (int iy = ny - 1; iy >= 0; iy--) {
- int iRowPos = ((ny - 1 - iy) * scale) * (nx * scale);
- for (int ix = 0; ix < nx; ix++) {
- int cval = static_cast<int>((v[ix][iy] - pmin) * view_scale);
- if (cval < 0)
- cval = 0;
- else if (cval > 255)
- cval = 255;
- for (int sy = 0; sy < scale; sy++)
- for (int sx = 0; sx < scale; sx++)
- pPens[iRowPos+(sy * nx * scale)+(sx + (ix * scale))] = grayscale[cval];
- }
- }
-
- g2_image (id_X11, 0., 0., nx * scale, ny * scale, pPens);
-
- delete pPens;
- return (id_X11);
-#else
- return 0;
-#endif
-}
-
-
// ImageFile::comparativeStatistics Calculate comparative stats
//
ImageFileArray vResult = result.getArray();
ImageFileArray vResultImag = result.getImaginaryArray();
+ BilinearInterpolator<ImageFileValue> realInterp (vReal, nx, ny);
+ BilinearInterpolator<ImageFileValue> imagInterp (vImag, nx, ny);
+
for (unsigned int ix = 0; ix < newNX; ix++) {
for (unsigned int iy = 0; iy < newNY; iy++) {
double dXPos = ix / dXScale;
double dYPos = iy / dYScale;
- unsigned int scaleNX = static_cast<unsigned int> (dXPos);
- unsigned int scaleNY = static_cast<unsigned int> (dYPos);
- double dXFrac = dXPos - scaleNX;
- double dYFrac = dYPos - scaleNY;
- if (scaleNX >= nx - 1 || scaleNY >= ny - 1) {
- vResult[ix][iy] = vReal[scaleNX][scaleNY];
- if (result.isComplex()) {
- if (isComplex())
- vResultImag[ix][iy] = vImag[scaleNX][scaleNY];
- else
- vResultImag[ix][iy] = 0;
- }
- } else {
- vResult[ix][iy] = (1 - dXFrac) * (1 - dYFrac) * vReal[scaleNX][scaleNY] +
- dXFrac * (1 - dYFrac) * vReal[scaleNX+1][scaleNY] +
- dYFrac * (1 - dXFrac) * vReal[scaleNX][scaleNY+1] +
- dXFrac * dYFrac * vReal[scaleNX+1][scaleNY+1];
- if (result.isComplex()) {
- if (isComplex())
- vResultImag[ix][iy] = (1 - dXFrac) * (1 - dYFrac) * vImag[scaleNX][scaleNY] +
- dXFrac * (1 - dYFrac) * vImag[scaleNX+1][scaleNY] +
- dYFrac * (1 - dXFrac) * vImag[scaleNX][scaleNY+1] +
- dXFrac * dYFrac * vImag[scaleNX+1][scaleNY+1];
- else
- vResultImag[ix][iy] = 0;
- }
- }
+ vResult[ix][iy] = realInterp.interpolate (dXPos, dYPos);
+ if (result.isComplex())
+ if (isComplex())
+ vResultImag[ix][iy] = imagInterp.interpolate (dXPos, dYPos);
+ else
+ vResultImag[ix][iy] = 0;
}
}
-
+
return true;
}
unsigned int ix, iy;
unsigned int iArray = 0;
- for (ix = 0; ix < m_nx; ix++)
+ for (ix = 0; ix < m_nx; ix++) {
for (iy = 0; iy < m_ny; iy++) {
in[iArray].re = vReal[ix][iy];
if (isComplex())
in[iArray].im = 0;
iArray++;
}
-
- fftwnd_plan plan = fftw2d_create_plan (m_nx, m_ny, FFTW_FORWARD, FFTW_IN_PLACE);
-
- fftwnd_one (plan, in, NULL);
-
- ImageFileArray vRealResult = result.getArray();
- ImageFileArray vImagResult = result.getImaginaryArray();
- iArray = 0;
- unsigned int iScale = m_nx * m_ny;
- for (ix = 0; ix < m_nx; ix++)
- for (iy = 0; iy < m_ny; iy++) {
- vRealResult[ix][iy] = in[iArray].re / iScale;
- vImagResult[ix][iy] = in[iArray].im / iScale;
- iArray++;
- }
-
- fftwnd_destroy_plan (plan);
- delete in;
-
-
- Fourier::shuffleFourierToNaturalOrder (result);
-
- return true;
+ }
+
+ fftwnd_plan plan = fftw2d_create_plan (m_nx, m_ny, FFTW_FORWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
+ fftwnd_one (plan, in, NULL);
+
+ ImageFileArray vRealResult = result.getArray();
+ ImageFileArray vImagResult = result.getImaginaryArray();
+ iArray = 0;
+ unsigned int iScale = m_nx * m_ny;
+ for (ix = 0; ix < m_nx; ix++) {
+ for (iy = 0; iy < m_ny; iy++) {
+ vRealResult[ix][iy] = in[iArray].re / iScale;
+ vImagResult[ix][iy] = in[iArray].im / iScale;
+ iArray++;
+ }
+ }
+ delete in;
+ fftwnd_destroy_plan (plan);
+
+ Fourier::shuffleFourierToNaturalOrder (result);
+
+ return true;
}
ImageFileArray vRealResult = result.getArray();
ImageFileArray vImagResult = result.getImaginaryArray();
unsigned int ix, iy;
- for (ix = 0; ix < m_nx; ix++)
+ for (ix = 0; ix < m_nx; ix++) {
for (iy = 0; iy < m_ny; iy++) {
vRealResult[ix][iy] = vReal[ix][iy];
if (isComplex())
else
vImagResult[ix][iy] = 0;
}
-
- Fourier::shuffleNaturalToFourierOrder (result);
-
- fftw_complex* in = new fftw_complex [m_nx * m_ny];
-
- unsigned int iArray = 0;
- for (ix = 0; ix < m_nx; ix++)
- for (iy = 0; iy < m_ny; iy++) {
- in[iArray].re = vRealResult[ix][iy];
- in[iArray].im = vImagResult[ix][iy];
- iArray++;
- }
-
- fftwnd_plan plan = fftw2d_create_plan (m_nx, m_ny, FFTW_BACKWARD, FFTW_IN_PLACE);
-
- fftwnd_one (plan, in, NULL);
-
- iArray = 0;
- for (ix = 0; ix < m_nx; ix++)
- for (iy = 0; iy < m_ny; iy++) {
- vRealResult[ix][iy] = in[iArray].re;
- vImagResult[ix][iy] = in[iArray].im;
- iArray++;
- }
-
- fftwnd_destroy_plan (plan);
-
- delete in;
-
- return true;
+ }
+
+ Fourier::shuffleNaturalToFourierOrder (result);
+
+ fftw_complex* in = new fftw_complex [m_nx * m_ny];
+
+ unsigned int iArray = 0;
+ for (ix = 0; ix < m_nx; ix++) {
+ for (iy = 0; iy < m_ny; iy++) {
+ in[iArray].re = vRealResult[ix][iy];
+ in[iArray].im = vImagResult[ix][iy];
+ iArray++;
+ }
+ }
+
+ fftwnd_plan plan = fftw2d_create_plan (m_nx, m_ny, FFTW_BACKWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
+
+ fftwnd_one (plan, in, NULL);
+
+ iArray = 0;
+ for (ix = 0; ix < m_nx; ix++) {
+ for (iy = 0; iy < m_ny; iy++) {
+ vRealResult[ix][iy] = in[iArray].re;
+ vImagResult[ix][iy] = in[iArray].im;
+ iArray++;
+ }
+ }
+ fftwnd_destroy_plan (plan);
+
+ delete in;
+
+ return true;
}
bool
if (! result.convertRealToComplex ())
return false;
}
-
- fftw_complex* in = new fftw_complex [m_nx];
-
+
ImageFileArrayConst vReal = getArray();
ImageFileArrayConst vImag = getImaginaryArray();
- fftw_plan plan = fftw_create_plan (m_nx, FFTW_FORWARD, FFTW_IN_PLACE);
- std::complex<double>* pcRow = new std::complex<double> [m_nx];
-
- unsigned int ix, iy;
- unsigned int iArray = 0;
- for (iy = 0; iy < m_ny; iy++) {
+ fftw_plan plan = fftw_create_plan (m_nx, FFTW_FORWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
+
+ fftw_complex* in = new fftw_complex [m_nx];
+ std::complex<double>* pcRow = new std::complex<double> [m_nx];
+ for (int iy = 0; iy < m_ny; iy++) {
+ unsigned int ix;
for (ix = 0; ix < m_nx; ix++) {
in[ix].re = vReal[ix][iy];
if (isComplex())
ImageFileArrayConst vReal = getArray();
ImageFileArrayConst vImag = getImaginaryArray();
- fftw_plan plan = fftw_create_plan (m_nx, FFTW_BACKWARD, FFTW_IN_PLACE);
+ fftw_plan plan = fftw_create_plan (m_nx, FFTW_BACKWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
std::complex<double>* pcRow = new std::complex<double> [m_nx];
unsigned int ix, iy;
if (! result.convertRealToComplex ())
return false;
}
-
- fftw_complex* in = new fftw_complex [m_ny];
-
+
ImageFileArrayConst vReal = getArray();
ImageFileArrayConst vImag = getImaginaryArray();
- fftw_plan plan = fftw_create_plan (m_ny, FFTW_FORWARD, FFTW_IN_PLACE);
- std::complex<double>* pcCol = new std::complex<double> [m_ny];
-
- unsigned int ix, iy;
- unsigned int iArray = 0;
- for (ix = 0; ix < m_nx; ix++) {
+ fftw_plan plan = fftw_create_plan (m_ny, FFTW_FORWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
+
+ std::complex<double>* pcCol = new std::complex<double> [m_ny];
+ fftw_complex* in = new fftw_complex [m_ny];
+ for (int ix = 0; ix < m_nx; ix++) {
+ unsigned int iy;
for (iy = 0; iy < m_ny; iy++) {
in[iy].re = vReal[ix][iy];
if (isComplex())
ImageFileArrayConst vReal = getArray();
ImageFileArrayConst vImag = getImaginaryArray();
- fftw_plan plan = fftw_create_plan (m_ny, FFTW_BACKWARD, FFTW_IN_PLACE);
+ fftw_plan plan = fftw_create_plan (m_ny, FFTW_BACKWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
std::complex<double>* pcCol = new std::complex<double> [m_ny];
unsigned int ix, iy;
}
if (result.isComplex())
- result.convertComplexToReal();
+ result.reallocComplexToReal();
return true;
}
ImageFileArray vImag = getImaginaryArray();
ImageFileArray vRealResult = result.getArray();
- for (unsigned int ix = 0; ix < m_nx; ix++)
+ for (unsigned int ix = 0; ix < m_nx; ix++) {
for (unsigned int iy = 0; iy < m_ny; iy++) {
if (isComplex())
vRealResult[ix][iy] = ::atan2 (vImag[ix][iy], vReal[ix][iy]);
else
vRealResult[ix][iy] = 0;
}
+ }
+ if (result.isComplex())
+ result.reallocComplexToReal();
- if (result.isComplex())
- result.convertComplexToReal();
+ return true;
+}
+
+bool
+ImageFile::real (ImageFile& result) const
+{
+ if (m_nx != result.nx() || m_ny != result.ny()) {
+ sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
+ return false;
+ }
+
+ ImageFileArray vReal = getArray();
+ ImageFileArray vRealResult = result.getArray();
+
+ for (unsigned int ix = 0; ix < m_nx; ix++) {
+ for (unsigned int iy = 0; iy < m_ny; iy++) {
+ vRealResult[ix][iy] = vReal[ix][iy];
+ }
+ }
+
+ if (result.isComplex())
+ result.reallocComplexToReal();
- return true;
+ return true;
+}
+
+bool
+ImageFile::imaginary (ImageFile& result) const
+{
+ if (m_nx != result.nx() || m_ny != result.ny()) {
+ sys_error (ERR_WARNING, "Difference sizes of images [ImageFile::invertPixelValues]");
+ return false;
+ }
+
+ ImageFileArray vImag = getArray();
+ ImageFileArray vRealResult = result.getArray();
+
+ for (unsigned int ix = 0; ix < m_nx; ix++) {
+ for (unsigned int iy = 0; iy < m_ny; iy++) {
+ if (isComplex())
+ vRealResult[ix][iy] = vImag[ix][iy];
+ else
+ vRealResult[ix][iy] = 0;
+ }
+ }
+
+ if (result.isComplex())
+ result.reallocComplexToReal();
+
+ return true;
}
bool
vResult[ix][iy] = vLHS[ix][iy] * vLHS[ix][iy];
}
}
-
-
+
return true;
}
} else
skipSpacePPM (fp); // ascii may have comments
+ bool bMonochromeImage = false;
double dMaxValue = iMaxValue;
+ double dMaxValue3 = iMaxValue * 3;
+
ImageFileArray v = getArray();
for (int iy = nRows - 1; iy >= 0; iy--) {
for (int ix = 0; ix < nCols; ix++) {
fclose(fp);
return false;
}
- dR = iR / dMaxValue;
- dG = iG / dMaxValue;
- dB = iB / dMaxValue;
- v[ix][iy] = colorToGrayscale (dR, dG, dB);
+ if (ix == 0 && iy == 0 && (iR == iG && iG == iB))
+ bMonochromeImage = true;
+ if (bMonochromeImage)
+ v[ix][iy] = (iR + iG + iB) / dMaxValue3;
+ else {
+ dR = iR / dMaxValue;
+ dG = iG / dMaxValue;
+ dB = iB / dMaxValue;
+ v[ix][iy] = colorToGrayscale (dR, dG, dB);
+ }
break;
case '6':
iR = fgetc(fp);
iG = fgetc(fp);
iB = fgetc(fp);
+
if (iB == EOF) {
fclose(fp);
return false;
}
- dR = iR / dMaxValue;
- dG = iG / dMaxValue;
- dB = iB / dMaxValue;
- v[ix][iy] = colorToGrayscale (dR, dG, dB);
+ if (ix == 0 && iy == 0 && (iR == iG && iG == iB))
+ bMonochromeImage = true;
+
+ if (bMonochromeImage)
+ v[ix][iy] = (iR + iG + iB) / dMaxValue3;
+ else {
+ dR = iR / dMaxValue;
+ dG = iG / dMaxValue;
+ dB = iB / dMaxValue;
+ v[ix][iy] = colorToGrayscale (dR, dG, dB);
+ }
break;
}
}
return bSuccess;
}
#endif
-
+
sys_error (ERR_SEVERE, "Invalid format %s [ImageFile::exportImage]", pszFormat);
return false;
}
int ny = m_ny;
ImageFileArray v = getArray();
ImageFileArray vImag = getImaginaryArray();
-
+
if ((fp = fopen (outfile, "w")) == NULL)
return false;
** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: procsignal.cpp,v 1.29 2001/03/18 18:08:25 kevin Exp $
+** $Id: procsignal.cpp,v 1.30 2001/03/21 21:45:31 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
}
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);
+ m_realPlanForward = rfftw_create_plan (m_nFilterPoints, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE | FFTW_USE_WISDOM);
+ m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE | FFTW_USE_WISDOM);
m_adRealFftInput = new fftw_real [ m_nFilterPoints ];
m_adRealFftSignal = new fftw_real [ m_nOutputPoints ];
for (i = 0; i < m_nFilterPoints; i++)
m_adRealFftInput[i] = 0;
} else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
- m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE);
- m_complexPlanBackward = fftw_create_plan (m_nOutputPoints, FFTW_BACKWARD, FFTW_ESTIMATE);
+ m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE | FFTW_USE_WISDOM);
+ m_complexPlanBackward = fftw_create_plan (m_nOutputPoints, FFTW_BACKWARD, FFTW_ESTIMATE | FFTW_USE_WISDOM);
m_adComplexFftInput = new fftw_complex [ m_nFilterPoints ];
m_adComplexFftSignal = new fftw_complex [ m_nOutputPoints ];
for (i = 0; i < m_nFilterPoints; i++)
** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: projections.cpp,v 1.67 2001/03/18 18:08:25 kevin Exp $
+** $Id: projections.cpp,v 1.68 2001/03/21 21:45:31 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
double** ppdView = adView.getArray();
double** ppdDet = adDet.getArray();
- pProj->calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet, m_nDet, 1.);
-
std::complex<double>** ppcDetValue = new std::complex<double>* [m_nView];
unsigned int iView;
for (iView = 0; iView < m_nView; iView++) {
ppcDetValue[iView][iDet] = std::complex<double>(detval[iDet], 0);
}
- pProj->interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, pProj->m_nView, pProj->m_nDet, pProj->m_nDet, iInterpolationID);
+ pProj->calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet, m_nDet, 1., m_detInc);
+
+ pProj->interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, pProj->m_nView, pProj->m_nDet,
+ pProj->m_nDet, iInterpolationID);
for (iView = 0; iView < m_nView; iView++)
delete [] ppcDetValue[iView];
bool
Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad)
{
-#ifndef HAVE_FFT
+#ifndef HAVE_FFTW
+ rIF.arrayDataClear();
return false;
#else
unsigned int nx = rIF.nx();
sys_error (ERR_WARNING, "convertFFTPolar supports Parallel only");
return false;
}
-
- Array2d<double> adView (nx, ny);
- Array2d<double> adDet (nx, ny);
- double** ppdView = adView.getArray();
- double** ppdDet = adDet.getArray();
- std::complex<double>** ppcDetValue = new std::complex<double>* [m_nView];
+ int iInterpDet = nx;
+ int iNumInterpDetWithZeros = ProcessSignal::addZeropadFactor (iInterpDet, iZeropad);
- int iNumDetWithZeros = ProcessSignal::addZeropadFactor (m_nDet, iZeropad);
- double dZeropadRatio = iNumDetWithZeros / static_cast<double>(m_nDet);
+ double dZeropadRatio = static_cast<double>(iNumInterpDetWithZeros) / static_cast<double>(iInterpDet);
- double* pdDet = new double [iNumDetWithZeros];
- fftw_complex* pcIn = new fftw_complex [iNumDetWithZeros];
- fftw_plan plan = fftw_create_plan (iNumDetWithZeros, FFTW_FORWARD, FFTW_IN_PLACE);
+ fftw_plan plan = fftw_create_plan (iNumInterpDetWithZeros, FFTW_FORWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
+ fftw_complex* pcIn = new fftw_complex [iNumInterpDetWithZeros];
+ std::complex<double>** ppcDetValue = new std::complex<double>* [m_nView];
+ double dInterpScale = (m_nDet-1) / static_cast<double>(iInterpDet-1) / SQRT2;
+
for (unsigned int iView = 0; iView < m_nView; iView++) {
DetectorValue* detval = getDetectorArray(iView).detValues();
- for (unsigned int iDet = 0; iDet < m_nDet; iDet++) {
- pcIn[iDet].re = detval[iDet];
+ LinearInterpolator<DetectorValue> projInterp (detval, m_nDet);
+ for (unsigned int iDet = 0; iDet < iInterpDet; iDet++) {
+// double dInterpPos = iInterpDet * dInterpScale;
+ double dInterpPos = (m_nDet / 2.) + (iDet - iInterpDet/2.) * dInterpScale;
+ pcIn[iDet].re = projInterp.interpolate (dInterpPos);
pcIn[iDet].im = 0;
}
- for (unsigned int iDet2 = m_nDet; iDet2 < iNumDetWithZeros; iDet2++)
+ for (unsigned int iDet2 = iInterpDet; iDet2 < iNumInterpDetWithZeros; iDet2++)
pcIn[iDet2].re = pcIn[iDet2].im = 0;
fftw_one (plan, pcIn, NULL);
- ppcDetValue[iView] = new std::complex<double> [iNumDetWithZeros];
- for (unsigned int iD = 0; iD < iNumDetWithZeros; iD++)
- ppcDetValue[iView][iD] = std::complex<double> (pcIn[iD].re, pcIn[iD].im);
-//ppcDetValue[iView][iD] = std::complex<double> (std::abs(std::complex<double>(pcIn[iD].re, pcIn[iD].im)), 0);
- Fourier::shuffleFourierToNaturalOrder (ppcDetValue[iView], iNumDetWithZeros);
+
+ ppcDetValue[iView] = new std::complex<double> [iNumInterpDetWithZeros];
+ for (unsigned int iD = 0; iD < iNumInterpDetWithZeros; iD++)
+ ppcDetValue[iView][iD] = std::complex<double> (pcIn[iD].re / iInterpDet / (iInterpDet/2), pcIn[iD].im / iInterpDet / (iInterpDet/2));
+
+ Fourier::shuffleFourierToNaturalOrder (ppcDetValue[iView], iNumInterpDetWithZeros);
}
+ delete [] pcIn;
fftw_destroy_plan (plan);
- delete [] pcIn;
- calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet, iNumDetWithZeros, dZeropadRatio);
- interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, m_nView, m_nDet, iNumDetWithZeros,
+ Array2d<double> adView (nx, ny);
+ Array2d<double> adDet (nx, ny);
+ double** ppdView = adView.getArray();
+ double** ppdDet = adDet.getArray();
+ calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet, iNumInterpDetWithZeros, dZeropadRatio,
+ m_detInc * dInterpScale);
+
+ interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, m_nView, m_nDet, iNumInterpDetWithZeros,
iInterpolationID);
for (int i = 0; i < m_nView; i++)
void
Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double** ppdView, double** ppdDet,
- int iNumDetWithZeros, double dZeropadRatio)
+ int iNumDetWithZeros, double dZeropadRatio, double dDetInc)
{
- double xMin = -phmLen() / 2;
- double xMax = xMin + phmLen();
- double yMin = -phmLen() / 2;
- double yMax = yMin + phmLen();
+// double dLength = viewDiameter();
+ double dLength = phmLen();
+ double xMin = -dLength / 2;
+ double xMax = xMin + dLength;
+ double yMin = -dLength / 2;
+ double yMax = yMin + dLength;
double xCent = (xMin + xMax) / 2;
double yCent = (yMin + yMax) / 2;
}
ppdView[ix][iy] = (phi - m_rotStart) / m_rotInc;
- ppdDet[ix][iy] = (r / m_detInc) + iDetCenter;
+ ppdDet[ix][iy] = (r / dDetInc) + iDetCenter;
}
}
}
unsigned int nx, unsigned int ny, std::complex<double>** ppcDetValue, double** ppdView,
double** ppdDet, unsigned int nView, unsigned int nDet, unsigned int nDetWithZeros, int iInterpolationID)
{
+ typedef std::complex<double> complexValue;
+ BilinearInterpolator<complexValue> bilinear (ppcDetValue, nView, nDetWithZeros);
for (unsigned int ix = 0; ix < ny; ix++) {
for (unsigned int iy = 0; iy < ny; iy++) {
+
if (iInterpolationID == POLAR_INTERP_NEAREST) {
unsigned int iView = nearest<int> (ppdView[ix][iy]);
unsigned int iDet = nearest<int> (ppdDet[ix][iy]);
if (iView == nView) {
iView = 0;
- // iDet = m_nDet - iDet;
+ iDet = m_nDet - iDet;
}
if (iDet >= 0 && iDet < nDetWithZeros && iView >= 0 && iView < nView) {
v[ix][iy] = ppcDetValue[iView][iDet].real();
if (vImag)
vImag[ix][iy] = ppcDetValue[iView][iDet].imag();
- } else {
- sys_error (ERR_SEVERE, "Can't find projection data for ix=%d,iy=%d with radView=%f and radDet=%f",
- ix, iy, ppdView[ix][iy], ppdDet[ix][iy]);
+ } else
v[ix][iy] = 0;
- }
+
} else if (iInterpolationID == POLAR_INTERP_BILINEAR) {
- unsigned int iFloorView = static_cast<int>(ppdView[ix][iy]);
+#if 1
+ std::complex<double> vInterp = bilinear.interpolate (ppdView[ix][iy], ppdDet[ix][iy]);
+ v[ix][iy] = vInterp.real();
+ if (vImag)
+ vImag[ix][iy] = vInterp.imag();
+#else
+ int iFloorView = ::floor (ppdView[ix][iy]);
double dFracView = ppdView[ix][iy] - iFloorView;
- unsigned int iFloorDet = static_cast<int>(ppdDet[ix][iy]);
+ int iFloorDet = ::floor (ppdDet[ix][iy]);
double dFracDet = ppdDet[ix][iy] - iFloorDet;
if (iFloorDet >= 0 && iFloorView >= 0) {
v3 = v2;
else
v3 = ppcDetValue[0][iFloorDet+1];
+
std::complex<double> vInterp = (1 - dFracView) * (1 - dFracDet) * v1 +
dFracView * (1 - dFracDet) * v2 + dFracView * dFracDet * v3 +
dFracDet * (1 - dFracView) * v4;
if (vImag)
vImag[ix][iy] = vInterp.imag();
} else {
- sys_error (ERR_SEVERE, "Can't find projection data for ix=%d,iy=%d with radView=%f and radDet=%f",
- ix, iy, ppdView[ix][iy], ppdDet[ix][iy]);
+ // sys_error (ERR_SEVERE, "Can't find projection data for ix=%d,iy=%d with radView=%f and radDet=%f", ix, iy, ppdView[ix][iy], ppdDet[ix][iy]);
v[ix][iy] = 0;
if (vImag)
vImag[ix][iy] = 0;
}
+#endif
} else if (iInterpolationID == POLAR_INTERP_BICUBIC) {
v[ix][iy] =0;
if (vImag)
// interpolate to evenly spaced theta (views)
double dDetPos = pProjNew->m_detStart;
for (int iD = 0; iD < pProjNew->nDet(); iD++, dDetPos += pProjNew->m_detInc) {
- parallel.getThetaAndRaysumsForT (iD, pdThetaValuesForT, pdRaysumsForT);
+ parallel.getThetaAndRaysumsForT (iD, pdThetaValuesForT, pdRaysumsForT);
double dViewAngle = m_rotStart;
int iLastFloor = -1;
for (int iV = 0; iV < pProjNew->nView(); iV++, dViewAngle += pProjNew->m_rotInc) {
DetectorValue* detValues = pProjNew->getDetectorArray (iV).detValues();
-
- detValues[iD] = parallel.interpolate (pdThetaValuesForT, pdRaysumsForT, pProjNew->nView(), dViewAngle, &iLastFloor);
+ LinearInterpolator<double> interp (pdThetaValuesForT, pdRaysumsForT, pProjNew->nView());
+ detValues[iD] = interp.interpolate (dViewAngle, &iLastFloor);
}
}
delete pdThetaValuesForT;
double dDetPos = pProjNew->m_detStart;
int iLastFloor = -1;
- for (int iD = 0; iD < pProjNew->nDet(); iD++, dDetPos += pProjNew->m_detInc) {
- detValues[iD] = parallel.interpolate (pdOriginalDetPositions, pdDetValueCopy, pProjNew->nDet(), dDetPos, &iLastFloor);
- }
+ LinearInterpolator<double> interp (pdOriginalDetPositions, pdDetValueCopy, pProjNew->nDet());
+ for (int iD = 0; iD < pProjNew->nDet(); iD++, dDetPos += pProjNew->m_detInc)
+ detValues[iD] = interp.interpolate (dDetPos, &iLastFloor);
}
delete pdDetValueCopy;
delete pdOriginalDetPositions;
iPos += m_iNumView;
}
}
-
-// locate by bisection, O(log2(n))
-// iLastFloor is used when sequential calls to interpolate with monotonically increasing dX
-double
-ParallelRaysums::interpolate (double* pdX, double* pdY, int n, double dX, int* iLastFloor)
-{
- int iLower = -1;
- int iUpper = n;
- if (iLastFloor && *iLastFloor >= 0 && pdX[*iLastFloor] < dX)
- iLower = *iLastFloor;
-
- while (iUpper - iLower > 1) {
- int iMiddle = (iUpper + iLower) >> 1;
- if (dX >= pdX[iMiddle])
- iLower = iMiddle;
- else
- iUpper = iMiddle;
- }
- if (dX <= pdX[0])
- return pdY[0];
- else if (dX >= pdX[n-1])
- return pdY[1];
-
- if (iLower < 0 || iLower >= n) {
- sys_error (ERR_SEVERE, "Coordinate out of range [locateThetaBase]");
- return 0;
- }
-
- if (iLastFloor)
- *iLastFloor = iLower;
- return pdY[iLower] + (pdY[iUpper] - pdY[iLower]) * ((dX - pdX[iLower]) / (pdX[iUpper] - pdX[iLower]));
-}
-
** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: interpolator.cpp,v 1.2 2001/02/22 11:05:38 kevin Exp $
+** $Id: interpolator.cpp,v 1.3 2001/03/21 21:45:31 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
+
+
+
** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: mathfuncs.cpp,v 1.8 2001/01/28 19:10:18 kevin Exp $
+** $Id: mathfuncs.cpp,v 1.9 2001/03/21 21:45:31 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
else // Even
median = (vec[ (nPoints / 2) - 1 ] + vec[ nPoints / 2 ]) / 2;
}
+
+
--------------------Configuration: libctsim - Win32 Debug--------------------
</h3>
<h3>Command Lines</h3>
-Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP5E5.tmp" with contents
+Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP16.tmp" with contents
[
/nologo /G6 /MTd /W3 /Gm /Gi /GR /GX /Zi /Od /Gy /I "..\..\..\wx2.2.5\src\png" /I "..\..\..\wx2.2.5\src\zlib" /I "..\..\INCLUDE" /I "..\..\getopt" /I "..\..\..\fftw-2.1.3\fftw" /I "..\..\..\fftw-2.1.3\rfftw" /I "..\..\..\wx2.2.5\include" /I "\dicom\ctn\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__" /D "__WIN95__" /D "__WIN32__" /D WINVER=0x0400 /D "HAVE_CTN_DICOM" /D VERSION=\"3.1.0\" /FR"Debug/" /Fp"Debug/libctsim.pch" /YX /Fo"Debug/" /Fd"Debug/" /FD /GZ /c
-"C:\ctsim\libctsim\imagefile.cpp"
+"C:\ctsim\libctsim\projections.cpp"
]
-Creating command line "cl.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP5E5.tmp"
-Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP5E6.tmp" with contents
+Creating command line "cl.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP16.tmp"
+Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP17.tmp" with contents
[
/nologo /out:"Debug\libctsim.lib"
.\Debug\array2dfile.obj
.\Debug\transformmatrix.obj
.\Debug\xform.obj
]
-Creating command line "link.exe -lib @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP5E6.tmp"
+Creating command line "link.exe -lib @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP17.tmp"
<h3>Output Window</h3>
Compiling...
-imagefile.cpp
+projections.cpp
Creating library...
<h3>
--------------------Configuration: ctsim - Win32 Debug--------------------
</h3>
<h3>Command Lines</h3>
-Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP5E7.tmp" with contents
+Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP18.tmp" with contents
[
winmm.lib rpcrt4.lib ws2_32.lib ../libctsim/Debug/libctsim.lib libcmtd.lib ..\..\..\fftw-2.1.3\Win32\FFTW2st\Debug\FFTW2st.lib ..\..\..\fftw-2.1.3\Win32\RFFTW2st\Debug\RFFTW2st.lib wxd.lib xpmd.lib tiffd.lib zlibd.lib pngd.lib comctl32.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 opengl32.lib glu32.lib htmlhelp.lib ctn_lib.lib /nologo /subsystem:windows /incremental:yes /pdb:"Debug/ctsim.pdb" /debug /machine:I386 /out:"Debug/ctsim.exe" /pdbtype:sept /libpath:"\wx2.2.5\lib" /libpath:"\dicom\ctn\winctn\ctn_lib\Debug"
.\Debug\backgroundmgr.obj
\wx2.2.5\lib\zlibd.lib
\wx2.2.5\lib\tiffd.lib
]
-Creating command line "link.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP5E7.tmp"
+Creating command line "link.exe @C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP18.tmp"
<h3>Output Window</h3>
Linking...
** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: ctsim.cpp,v 1.98 2001/03/07 16:34:47 kevin Exp $
+** $Id: ctsim.cpp,v 1.99 2001/03/21 21:45:31 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
#endif
#endif
-static const char* rcsindent = "$Id: ctsim.cpp,v 1.98 2001/03/07 16:34:47 kevin Exp $";
+static const char* rcsindent = "$Id: ctsim.cpp,v 1.99 2001/03/21 21:45:31 kevin Exp $";
struct option CTSimApp::ctsimOptions[] =
{
m_pConfig->Read ("StartupTips", &m_bShowStartupTips);
m_pConfig->Read ("CurrentTip", &m_iCurrentTip);
m_pConfig->Read ("UseBackgroundTasks", &m_bUseBackgroundTasks);
+#ifdef HAVE_FFTW
+ wxString strFftwWisdom;
+ m_pConfig->Read ("FftwWisdom", strFftwWisdom);
+ if (strFftwWisdom.size() > 0)
+ fftw_import_wisdom_from_string (strFftwWisdom.c_str());
+#endif
}
void
m_pConfig->Write ("StartupTips", m_bShowStartupTips);
m_pConfig->Write ("CurrentTip", m_iCurrentTip);
m_pConfig->Write ("UseBackgroundTasks", m_bUseBackgroundTasks);
-
+#ifdef HAVE_FFTW
+ const char* const pszWisdom = fftw_export_wisdom_to_string();
+ wxString strFftwWisdom (pszWisdom);
+ fftw_free ((void*) pszWisdom);
+ m_pConfig->Write ("FftwWisdom", strFftwWisdom);
+#endif
+
delete m_pConfig;
}
** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: ctsim.h,v 1.60 2001/03/18 18:08:26 kevin Exp $
+** $Id: ctsim.h,v 1.61 2001/03/21 21:45:31 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
IFMENU_FILTER_IFFT_COLS,
IFMENU_FILTER_MAGNITUDE,
IFMENU_FILTER_PHASE,
+ IFMENU_FILTER_REAL,
+ IFMENU_FILTER_IMAGINARY,
IFMENU_FILTER_SHUFFLENATURALTOFOURIERORDER,
IFMENU_FILTER_SHUFFLEFOURIERTONATURALORDER,
PLOTMENU_VIEW_SCALE_AUTO,
PLOTMENU_VIEW_SCALE_FULL,
- GRAPH3D_VIEW_SURFACE,
+ GRAPH3D_VIEW_WIREFRAME,
GRAPH3D_VIEW_COLOR,
GRAPH3D_VIEW_LIGHTING,
GRAPH3D_VIEW_SMOOTH,
** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: graph3dview.cpp,v 1.20 2001/03/18 18:08:26 kevin Exp $
+** $Id: graph3dview.cpp,v 1.21 2001/03/21 21:45:31 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
EVT_MENU(GRAPH3D_VIEW_LIGHTING, Graph3dFileView::OnLighting)
EVT_MENU(GRAPH3D_VIEW_COLOR, Graph3dFileView::OnColor)
EVT_MENU(GRAPH3D_VIEW_SMOOTH, Graph3dFileView::OnSmooth)
-EVT_MENU(GRAPH3D_VIEW_SURFACE, Graph3dFileView::OnSurface)
+EVT_MENU(GRAPH3D_VIEW_WIREFRAME, Graph3dFileView::OnWireframe)
EVT_MENU(GRAPH3D_VIEW_SCALE_MINMAX, Graph3dFileView::OnScaleSet)
EVT_MENU(GRAPH3D_VIEW_SCALE_AUTO, Graph3dFileView::OnScaleAuto)
EVT_MENU(GRAPH3D_VIEW_SCALE_FULL, Graph3dFileView::OnScaleFull)
{
m_bDoubleBuffer = true;
m_bSmooth = true;
- m_bLighting = false;
- m_bSurface = false;
- m_bLighting = false;
+ m_bLighting = true;
+ m_bWireframe = false;
+ m_bLighting = true;
m_bColor = true;
m_dXRotate = -180;
m_dYRotate = -210;
m_dZRotate = -195;
m_bColorScaleMinSet = false;
m_bColorScaleMaxSet = false;
- m_bCalculatedSurfaceBackground = false;
}
Graph3dFileView::~Graph3dFileView()
m_pViewMenu->Check (GRAPH3D_VIEW_COLOR, m_bColor);
m_pViewMenu->Check (GRAPH3D_VIEW_LIGHTING, m_bLighting);
m_pViewMenu->Check (GRAPH3D_VIEW_SMOOTH, m_bSmooth);
- m_pViewMenu->Check (GRAPH3D_VIEW_SURFACE, m_bSurface);
+ m_pViewMenu->Check (GRAPH3D_VIEW_WIREFRAME, m_bWireframe);
return true;
}
glTranslated (-static_cast<double>(nx - 1) / 2, 0.0, -static_cast<double>(ny - 1) / 2);
InitMaterials();
- if (m_bSurface) {
- glPolygonMode (GL_FRONT_AND_BACK, GL_FILL);
- if (! m_bColor) {
- glColor3f (1.0f, 1.0f, 1.0f);
- glCallList (DISPLAYLIST_NO_COLOR);
- } else
- glCallList (DISPLAYLIST_COLOR);
- }
- else {
+ if (m_bWireframe) {
if (! m_bColor)
glColor3f (1.0f, 1.0f, 1.0f);
glPolygonOffset (1.0f, 1.0f);
glPolygonOffset (0.0f, 0.0f);
glPolygonMode (GL_FRONT_AND_BACK, GL_LINE);
glCallList (DISPLAYLIST_NO_COLOR);
+ } else {
+ glPolygonMode (GL_FRONT_AND_BACK, GL_FILL);
+ if (! m_bColor) {
+ glColor3f (1.0f, 1.0f, 1.0f);
+ glCallList (DISPLAYLIST_NO_COLOR);
+ } else
+ glCallList (DISPLAYLIST_COLOR);
}
}
}
-void
-Graph3dFileView::calculateSurfaceBackground ()
-{
- if (m_bCalculatedSurfaceBackground)
- return;
-}
-
void
Graph3dFileView::OnProperties (wxCommandEvent& event)
{
}
void
-Graph3dFileView::OnSurface (wxCommandEvent& event)
+Graph3dFileView::OnWireframe (wxCommandEvent& event)
{
- m_bSurface = ! m_bSurface;
- m_pViewMenu->Check (GRAPH3D_VIEW_SURFACE, m_bSurface);
+ m_bWireframe = ! m_bWireframe;
+ m_pViewMenu->Check (GRAPH3D_VIEW_WIREFRAME, m_bWireframe);
m_pCanvas->Refresh();
}
GetDocumentManager()->FileHistoryUseMenu(m_pFileMenu);
m_pViewMenu = new wxMenu;
- m_pViewMenu->Append(GRAPH3D_VIEW_SURFACE, "Su&rface\tCtrl-R", "", true);
+ m_pViewMenu->Append(GRAPH3D_VIEW_WIREFRAME, "Wi&reframe\tCtrl-R", "", true);
m_pViewMenu->Append(GRAPH3D_VIEW_SMOOTH, "S&mooth\tCtrl-M", "", true);
m_pViewMenu->Append(GRAPH3D_VIEW_COLOR, "Co&lor\tCtrl-L", "", true);
m_pViewMenu->Append(GRAPH3D_VIEW_LIGHTING, "Li&ghting\tCtrl-G", "", true);
subframe->Centre(wxBOTH);
wxAcceleratorEntry accelEntries[7];
- accelEntries[0].Set (wxACCEL_CTRL, static_cast<int>('R'), GRAPH3D_VIEW_SURFACE);
+ accelEntries[0].Set (wxACCEL_CTRL, static_cast<int>('R'), GRAPH3D_VIEW_WIREFRAME);
accelEntries[1].Set (wxACCEL_CTRL, static_cast<int>('L'), GRAPH3D_VIEW_COLOR);
accelEntries[2].Set (wxACCEL_CTRL, static_cast<int>('G'), GRAPH3D_VIEW_LIGHTING);
accelEntries[3].Set (wxACCEL_CTRL, static_cast<int>('M'), GRAPH3D_VIEW_SMOOTH);
m_pView->m_dYRotate -= 15.0;
Refresh (false);
break;
- case 'r': case 'R':
- m_pView->OnSurface (dummyEvent);
+ case 'w': case 'W':
+ m_pView->OnWireframe (dummyEvent);
break;
case 's': case 'S':
m_pView->OnSmooth (dummyEvent);
** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: graph3dview.h,v 1.7 2001/03/18 18:08:26 kevin Exp $
+** $Id: graph3dview.h,v 1.8 2001/03/21 21:45:31 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
bool m_bDoubleBuffer;
bool m_bSmooth;
bool m_bLighting;
- bool m_bSurface;
+ bool m_bWireframe;
bool m_bColor;
enum {
DISPLAYLIST_COLOR = 1,
double m_dColorScaleMax;
bool m_bColorScaleMinSet;
bool m_bColorScaleMaxSet;
- bool m_bCalculatedSurfaceBackground;
void Draw();
void DrawSurface();
bool OnClose (bool deleteWindow = true);
void OnProperties (wxCommandEvent& event);
void OnLighting (wxCommandEvent& event);
- void OnSurface (wxCommandEvent& event);
+ void OnWireframe (wxCommandEvent& event);
void OnColor (wxCommandEvent& event);
void OnSmooth (wxCommandEvent& event);
void OnScaleSet (wxCommandEvent& event);
void OnScaleAuto (wxCommandEvent& event);
void OnScaleFull (wxCommandEvent& event);
- void calculateSurfaceBackground();
-
+
#if CTSIM_MDI
wxDocMDIChildFrame* getFrame() { return m_pFrame; }
#else
** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: views.cpp,v 1.139 2001/03/18 18:08:26 kevin Exp $
+** $Id: views.cpp,v 1.140 2001/03/21 21:45:31 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
#endif
EVT_MENU(IFMENU_FILTER_MAGNITUDE, ImageFileView::OnMagnitude)
EVT_MENU(IFMENU_FILTER_PHASE, ImageFileView::OnPhase)
+EVT_MENU(IFMENU_FILTER_REAL, ImageFileView::OnReal)
+EVT_MENU(IFMENU_FILTER_IMAGINARY, ImageFileView::OnImaginary)
EVT_MENU(IFMENU_PLOT_ROW, ImageFileView::OnPlotRow)
EVT_MENU(IFMENU_PLOT_COL, ImageFileView::OnPlotCol)
#ifdef HAVE_FFT
END_EVENT_TABLE()
ImageFileView::ImageFileView()
-: wxView(), m_pFrame(NULL), m_pCanvas(NULL), m_pFileMenu(0), m_bMinSpecified(false), m_bMaxSpecified(false)
+: wxView(), m_pFrame(NULL), m_pCanvas(NULL), m_pFileMenu(0), m_pFilterMenu(0), m_bMinSpecified(false), m_bMaxSpecified(false)
{
m_iDefaultExportFormatID = ImageFile::EXPORT_FORMAT_PNG;
}
GetDocument()->Activate();
}
+void
+ImageFileView::OnReal (wxCommandEvent& event)
+{
+ ImageFile& rIF = GetDocument()->getImageFile();
+ if (rIF.isComplex()) {
+ rIF.real (rIF);
+ rIF.labelAdd ("Real component of complex-image");
+ m_bMinSpecified = false;
+ m_bMaxSpecified = false;
+ if (theApp->getAskDeleteNewDocs())
+ GetDocument()->Modify (true);
+ GetDocument()->UpdateAllViews (this);
+ }
+ GetDocument()->Activate();
+}
+
+void
+ImageFileView::OnImaginary (wxCommandEvent& event)
+{
+ ImageFile& rIF = GetDocument()->getImageFile();
+ if (rIF.isComplex()) {
+ rIF.imaginary (rIF);
+ rIF.labelAdd ("Imaginary component of complex-image");
+ m_bMinSpecified = false;
+ m_bMaxSpecified = false;
+ if (theApp->getAskDeleteNewDocs())
+ GetDocument()->Modify (true);
+ GetDocument()->UpdateAllViews (this);
+ }
+ GetDocument()->Activate();
+}
+
ImageFileCanvas*
ImageFileView::CreateCanvas (wxFrame* parent)
theApp->setIconForFrame (subframe);
m_pFileMenu = new wxMenu;
-
m_pFileMenu->Append(MAINMENU_FILE_CREATE_PHANTOM, "Cr&eate Phantom...\tCtrl-P");
m_pFileMenu->Append(MAINMENU_FILE_CREATE_FILTER, "Create &Filter...\tCtrl-F");
m_pFileMenu->Append(wxID_OPEN, "&Open...\tCtrl-O");
view_menu->Append(IFMENU_VIEW_SCALE_AUTO, "Display Scale &Auto...\tCtrl-A");
view_menu->Append(IFMENU_VIEW_SCALE_FULL, "Display F&ull Scale\tCtrl-U");
- wxMenu* filter_menu = new wxMenu;
- filter_menu->Append (IFMENU_FILTER_INVERTVALUES, "In&vert Values");
- filter_menu->Append (IFMENU_FILTER_SQUARE, "&Square");
- filter_menu->Append (IFMENU_FILTER_SQRT, "Square &Root");
- filter_menu->Append (IFMENU_FILTER_LOG, "&Log");
- filter_menu->Append (IFMENU_FILTER_EXP, "&Exp");
- filter_menu->AppendSeparator();
+ m_pFilterMenu = new wxMenu;
+ m_pFilterMenu->Append (IFMENU_FILTER_INVERTVALUES, "In&vert Values");
+ m_pFilterMenu->Append (IFMENU_FILTER_SQUARE, "&Square");
+ m_pFilterMenu->Append (IFMENU_FILTER_SQRT, "Square &Root");
+ m_pFilterMenu->Append (IFMENU_FILTER_LOG, "&Log");
+ m_pFilterMenu->Append (IFMENU_FILTER_EXP, "E&xp");
+ m_pFilterMenu->AppendSeparator();
#ifdef HAVE_FFT
- filter_menu->Append (IFMENU_FILTER_FFT, "2-D &FFT");
- filter_menu->Append (IFMENU_FILTER_IFFT, "2-D &IFFT");
- filter_menu->Append (IFMENU_FILTER_FFT_ROWS, "FFT Rows");
- filter_menu->Append (IFMENU_FILTER_IFFT_ROWS, "IFFT Rows");
- filter_menu->Append (IFMENU_FILTER_FFT_COLS, "FFT Columns");
- filter_menu->Append (IFMENU_FILTER_IFFT_COLS, "IFFT Columns");
- filter_menu->Append (IFMENU_FILTER_FOURIER, "2-D F&ourier");
- filter_menu->Append (IFMENU_FILTER_INVERSE_FOURIER, "2-D Inverse Fo&urier");
+ m_pFilterMenu->Append (IFMENU_FILTER_FFT, "2-D &FFT\tCtrl-2");
+ m_pFilterMenu->Append (IFMENU_FILTER_IFFT, "2-D &IFFT\tAlt-2");
+ m_pFilterMenu->Append (IFMENU_FILTER_FFT_ROWS, "FFT Rows");
+ m_pFilterMenu->Append (IFMENU_FILTER_IFFT_ROWS, "IFFT Rows");
+ m_pFilterMenu->Append (IFMENU_FILTER_FFT_COLS, "FFT Columns");
+ m_pFilterMenu->Append (IFMENU_FILTER_IFFT_COLS, "IFFT Columns");
+ m_pFilterMenu->Append (IFMENU_FILTER_FOURIER, "2-D F&ourier");
+ m_pFilterMenu->Append (IFMENU_FILTER_INVERSE_FOURIER, "2-D Inverse Fo&urier");
#else
- filter_menu->Append (IFMENU_FILTER_FOURIER, "&Fourier");
- filter_menu->Append (IFMENU_FILTER_INVERSE_FOURIER, "&Inverse Fourier");
+ m_pFilterMenu->Append (IFMENU_FILTER_FOURIER, "&Fourier");
+ m_pFilterMenu->Append (IFMENU_FILTER_INVERSE_FOURIER, "&Inverse Fourier");
#endif
- filter_menu->Append (IFMENU_FILTER_SHUFFLEFOURIERTONATURALORDER, "S&huffle Fourier to Natural Order");
- filter_menu->Append (IFMENU_FILTER_SHUFFLENATURALTOFOURIERORDER, "Shu&ffle Natural to Fourier Order");
- filter_menu->Append (IFMENU_FILTER_MAGNITUDE, "&Magnitude");
- filter_menu->Append (IFMENU_FILTER_PHASE, "&Phase");
+ m_pFilterMenu->Append (IFMENU_FILTER_SHUFFLEFOURIERTONATURALORDER, "Shuffl&e Fourier to Natural Order");
+ m_pFilterMenu->Append (IFMENU_FILTER_SHUFFLENATURALTOFOURIERORDER, "Shuffle &Natural to Fourier Order");
+ m_pFilterMenu->AppendSeparator();
+ m_pFilterMenu->Append (IFMENU_FILTER_MAGNITUDE, "&Magnitude");
+ m_pFilterMenu->Append (IFMENU_FILTER_PHASE, "&Phase");
+ m_pFilterMenu->Append (IFMENU_FILTER_REAL, "Re&al");
+ m_pFilterMenu->Append (IFMENU_FILTER_IMAGINARY, "Ima&ginary");
wxMenu* image_menu = new wxMenu;
image_menu->Append (IFMENU_IMAGE_ADD, "&Add...");
menu_bar->Append(edit_menu, "&Edit");
menu_bar->Append(view_menu, "&View");
menu_bar->Append(image_menu, "&Image");
- menu_bar->Append(filter_menu, "Fi<er");
+ menu_bar->Append(m_pFilterMenu, "Fi<er");
menu_bar->Append(m_pMenuAnalyze, "&Analyze");
menu_bar->Append(help_menu, "&Help");
subframe->Centre(wxBOTH);
- wxAcceleratorEntry accelEntries[8];
+ wxAcceleratorEntry accelEntries[10];
accelEntries[0].Set (wxACCEL_CTRL, static_cast<int>('A'), IFMENU_VIEW_SCALE_AUTO);
accelEntries[1].Set (wxACCEL_CTRL, static_cast<int>('U'), IFMENU_VIEW_SCALE_FULL);
accelEntries[2].Set (wxACCEL_CTRL, static_cast<int>('E'), IFMENU_VIEW_SCALE_MINMAX);
accelEntries[4].Set (wxACCEL_CTRL, static_cast<int>('C'), IFMENU_EDIT_COPY);
accelEntries[5].Set (wxACCEL_CTRL, static_cast<int>('X'), IFMENU_EDIT_CUT);
accelEntries[6].Set (wxACCEL_CTRL, static_cast<int>('V'), IFMENU_EDIT_PASTE);
+ int iEntry = 7;
+#ifdef HAVE_FFT
+ accelEntries[iEntry++].Set (wxACCEL_CTRL, static_cast<int>('2'), IFMENU_FILTER_FFT);
+ accelEntries[iEntry++].Set (wxACCEL_ALT, static_cast<int>('2'), IFMENU_FILTER_IFFT);
+#endif
#ifdef wxUSE_GLCANVAS
- accelEntries[7].Set (wxACCEL_CTRL, static_cast<int>('3'), IFMENU_IMAGE_CONVERT3D);
- wxAcceleratorTable accelTable (8, accelEntries);
-#else
- wxAcceleratorTable accelTable (7, accelEntries);
+ accelEntries[iEntry++].Set (wxACCEL_CTRL, static_cast<int>('3'), IFMENU_IMAGE_CONVERT3D);
#endif
+ wxAcceleratorTable accelTable (iEntry, accelEntries);
subframe->SetAcceleratorTable (accelTable);
ImageFileView::OnUpdate (wxView *WXUNUSED(sender), wxObject *WXUNUSED(hint) )
{
const ImageFile& rIF = GetDocument()->getImageFile();
+ if (m_pFilterMenu && rIF.isComplex()) {
+ m_pFilterMenu->Enable(IFMENU_FILTER_REAL, true);
+ m_pFilterMenu->Enable(IFMENU_FILTER_IMAGINARY, true);
+ m_pFilterMenu->Enable(IFMENU_FILTER_PHASE, true);
+ m_pFilterMenu->Enable(IFMENU_FILTER_MAGNITUDE, true);
+ } else {
+ m_pFilterMenu->Enable(IFMENU_FILTER_REAL, false);
+ m_pFilterMenu->Enable(IFMENU_FILTER_IMAGINARY, false);
+ m_pFilterMenu->Enable(IFMENU_FILTER_PHASE, false);
+ m_pFilterMenu->Enable(IFMENU_FILTER_MAGNITUDE, false);
+ }
ImageFileArrayConst v = rIF.getArray();
int nx = rIF.nx();
int ny = rIF.ny();
int nx = rIF.nx();
int ny = rIF.ny();
+ bool bMonochrome = false;
+
if (bitmap.Ok() == true && bitmap.GetWidth() == nx && bitmap.GetHeight() == ny) {
wxImage image (bitmap);
+ double dScale3 = 3 * 255;
unsigned char* pixels = image.GetData();
ImageFileArray v = rIF.getArray();
for (int ix = 0; ix < rIF.nx(); ix++) {
for (int iy = 0; iy < rIF.ny(); iy++) {
unsigned int iBase = 3 * (iy * nx + ix);
- double dR = pixels[iBase] / 255.;
- double dG = pixels[iBase+1] / 255.;
- double dB = pixels[iBase+2] / 255.;
- v[ix][ny - 1 - iy] = ImageFile::colorToGrayscale (dR, dG, dB);
+ if (ix == 0 && iy == 0 && (pixels[iBase] == pixels[iBase+1] && pixels[iBase+1] == pixels[iBase+2]))
+ bMonochrome = true;
+ if (bMonochrome) {
+ v[ix][ny - 1 - iy] = (pixels[iBase]+pixels[iBase+1]+pixels[iBase+2]) / dScale3;
+ } else {
+ double dR = pixels[iBase] / 255.;
+ double dG = pixels[iBase+1] / 255.;
+ double dB = pixels[iBase+2] / 255.;
+ v[ix][ny - 1 - iy] = ImageFile::colorToGrayscale (dR, dG, dB);
+ }
}
}
OnUpdate(this, NULL);
pcIn[i].im = 0;
}
- fftw_plan plan = fftw_create_plan (nx, FFTW_FORWARD, FFTW_IN_PLACE);
+ fftw_plan plan = fftw_create_plan (nx, FFTW_FORWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
fftw_one (plan, pcIn, NULL);
fftw_destroy_plan (plan);
double* pYMag = new double [nx];
for (i = 0; i < nx; i++) {
pX[i] = i;
- pYReal[i] = pcIn[i].re;
- pYImag[i] = pcIn[i].im;
+ pYReal[i] = pcIn[i].re / nx;
+ pYImag[i] = pcIn[i].im / nx;
pYMag[i] = ::sqrt (pcIn[i].re * pcIn[i].re + pcIn[i].im * pcIn[i].im);
}
Fourier::shuffleFourierToNaturalOrder (pYReal, nx);
for (i = 0; i < ny; i++)
pcIn[i].im = pdTemp[i];
- fftw_plan plan = fftw_create_plan (ny, FFTW_BACKWARD, FFTW_IN_PLACE);
+ fftw_plan plan = fftw_create_plan (ny, FFTW_BACKWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
fftw_one (plan, pcIn, NULL);
fftw_destroy_plan (plan);
double* pYMag = new double [ny];
for (i = 0; i < ny; i++) {
pX[i] = i;
- pYReal[i] = pcIn[i].re;
- pYImag[i] = pcIn[i].im;
+ pYReal[i] = pcIn[i].re / ny;
+ pYImag[i] = pcIn[i].im / ny;
pYMag[i] = ::sqrt (pcIn[i].re * pcIn[i].re + pcIn[i].im * pcIn[i].im);
}
** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: views.h,v 1.52 2001/03/18 18:08:26 kevin Exp $
+** $Id: views.h,v 1.53 2001/03/21 21:45:31 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
ImageFileCanvas *m_pCanvas;
wxMenu* m_pFileMenu;
+ wxMenu* m_pFilterMenu;
bool m_bMinSpecified;
bool m_bMaxSpecified;
double m_dMinPixel;
void OnMagnitude (wxCommandEvent& event);
void OnPhase (wxCommandEvent& event);
+ void OnReal (wxCommandEvent& event);
+ void OnImaginary (wxCommandEvent& event);
void OnScaleAuto (wxCommandEvent& event);
void OnScaleMinMax (wxCommandEvent& event);
projects / ctsim.git / commitdiff