From: Kevin M. Rosenberg Date: Wed, 21 Mar 2001 21:45:31 +0000 (+0000) Subject: r641: no message X-Git-Tag: debian-4.5.3-3~376 X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=commitdiff_plain;h=3ea498d51ce4597e9649cd21f155b51175ea0bea r641: no message --- diff --git a/include/fourier.h b/include/fourier.h index 5fce2f6..b956fcc 100644 --- a/include/fourier.h +++ b/include/fourier.h @@ -9,7 +9,7 @@ ** 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 @@ -34,12 +34,73 @@ public: 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 + 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 + 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* 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* pdVector, const int n); +#endif }; diff --git a/include/imagefile.h b/include/imagefile.h index 140955c..56c255b 100644 --- a/include/imagefile.h +++ b/include/imagefile.h @@ -9,7 +9,7 @@ ** 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 @@ -214,9 +214,8 @@ public: #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); diff --git a/include/interpolator.h b/include/interpolator.h index def913e..cea4961 100644 --- a/include/interpolator.h +++ b/include/interpolator.h @@ -2,7 +2,7 @@ ** 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 @@ -47,3 +47,114 @@ public: double interpolate (double x); }; + +template +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 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; + } +}; + diff --git a/include/projections.h b/include/projections.h index 9c0cce0..d2fe89c 100644 --- a/include/projections.h +++ b/include/projections.h @@ -9,7 +9,7 @@ ** 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 @@ -134,7 +134,7 @@ class Projections bool convertPolar (ImageFile& rIF, int iInterpolation); bool convertFFTPolar (ImageFile& rIF, int iInterpolation, int iZeropad); void calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double** ppdView, double** ppdDet, - int iNumDetWithZeros, double dZeropadRatio); + int iNumDetWithZeros, double dZeropadRatio, double dDetInc); void interpolatePolar (ImageFileArray& v, ImageFileArray& vImag, unsigned int nx, unsigned int ny, std::complex** ppcDetValue, double** ppdDet, double** ppdView, unsigned int nView, unsigned int nDet, unsigned int nDetWithZeros, int iInterpolate); diff --git a/libctgraphics/bresenham.cpp b/libctgraphics/bresenham.cpp index e9936fe..070f44e 100644 --- a/libctgraphics/bresenham.cpp +++ b/libctgraphics/bresenham.cpp @@ -1,278 +1,95 @@ -/* ; 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); +} + + diff --git a/libctsim/fourier.cpp b/libctsim/fourier.cpp index f7357cd..6fec9f5 100644 --- a/libctsim/fourier.cpp +++ b/libctsim/fourier.cpp @@ -9,7 +9,7 @@ ** 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 @@ -98,182 +98,3 @@ Fourier::shuffleNaturalToFourierOrder (ImageFile& im) } 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* pdVector, const int n) -{ - std::complex* pdTemp = new std::complex [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* pdVector, const int n) -{ - std::complex* pdTemp = new std::complex [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; -} - - diff --git a/libctsim/imagefile.cpp b/libctsim/imagefile.cpp index 1c8a063..bab87bd 100644 --- a/libctsim/imagefile.cpp +++ b/libctsim/imagefile.cpp @@ -9,7 +9,7 @@ ** 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 @@ -150,7 +150,8 @@ ImageFile::getCenterCoordinates (unsigned int& iXCenter, unsigned int& iYCenter) 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); @@ -166,60 +167,6 @@ ImageFile::filterResponse (const char* const domainName, double bw, const char* } } -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((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 // @@ -739,40 +686,22 @@ ImageFile::scaleImage (ImageFile& result) const ImageFileArray vResult = result.getArray(); ImageFileArray vResultImag = result.getImaginaryArray(); + BilinearInterpolator realInterp (vReal, nx, ny); + BilinearInterpolator 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 (dXPos); - unsigned int scaleNY = static_cast (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; } @@ -797,7 +726,7 @@ ImageFile::fft (ImageFile& result) const 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()) @@ -806,29 +735,28 @@ ImageFile::fft (ImageFile& result) const 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; } @@ -850,7 +778,7 @@ ImageFile::ifft (ImageFile& result) const 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()) @@ -858,36 +786,38 @@ ImageFile::ifft (ImageFile& result) const 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 @@ -902,18 +832,16 @@ ImageFile::fftRows (ImageFile& result) const 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* pcRow = new std::complex [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* pcRow = new std::complex [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()) @@ -959,7 +887,7 @@ ImageFile::ifftRows (ImageFile& result) const 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* pcRow = new std::complex [m_nx]; unsigned int ix, iy; @@ -1006,18 +934,16 @@ ImageFile::fftCols (ImageFile& result) const 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* pcCol = new std::complex [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* pcCol = new std::complex [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()) @@ -1063,7 +989,7 @@ ImageFile::ifftCols (ImageFile& result) const 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* pcCol = new std::complex [m_ny]; unsigned int ix, iy; @@ -1263,7 +1189,7 @@ ImageFile::magnitude (ImageFile& result) const } if (result.isComplex()) - result.convertComplexToReal(); + result.reallocComplexToReal(); return true; } @@ -1280,18 +1206,67 @@ ImageFile::phase (ImageFile& result) const 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 @@ -1324,8 +1299,7 @@ ImageFile::square (ImageFile& result) const vResult[ix][iy] = vLHS[ix][iy] * vLHS[ix][iy]; } } - - + return true; } @@ -1484,7 +1458,10 @@ ImageFile::readImagePPM (const char* const pszFile) } 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++) { @@ -1511,23 +1488,37 @@ ImageFile::readImagePPM (const char* const pszFile) 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; } } @@ -1676,7 +1667,7 @@ ImageFile::exportImage (const char* const pszFormat, const char* const pszFilena return bSuccess; } #endif - + sys_error (ERR_SEVERE, "Invalid format %s [ImageFile::exportImage]", pszFormat); return false; } @@ -1767,7 +1758,7 @@ ImageFile::writeImageText (const char* const outfile) int ny = m_ny; ImageFileArray v = getArray(); ImageFileArray vImag = getImaginaryArray(); - + if ((fp = fopen (outfile, "w")) == NULL) return false; diff --git a/libctsim/procsignal.cpp b/libctsim/procsignal.cpp index 2d6600d..1c4f25f 100644 --- a/libctsim/procsignal.cpp +++ b/libctsim/procsignal.cpp @@ -9,7 +9,7 @@ ** 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 @@ -418,15 +418,15 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw } 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++) diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp index 1aee958..4d4a02d 100644 --- a/libctsim/projections.cpp +++ b/libctsim/projections.cpp @@ -8,7 +8,7 @@ ** 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 @@ -697,8 +697,6 @@ Projections::convertPolar (ImageFile& rIF, int iInterpolationID) double** ppdView = adView.getArray(); double** ppdDet = adDet.getArray(); - pProj->calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet, m_nDet, 1.); - std::complex** ppcDetValue = new std::complex* [m_nView]; unsigned int iView; for (iView = 0; iView < m_nView; iView++) { @@ -708,7 +706,10 @@ Projections::convertPolar (ImageFile& rIF, int iInterpolationID) ppcDetValue[iView][iDet] = std::complex(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]; @@ -724,7 +725,8 @@ Projections::convertPolar (ImageFile& rIF, int iInterpolationID) 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(); @@ -741,43 +743,50 @@ Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad sys_error (ERR_WARNING, "convertFFTPolar supports Parallel only"); return false; } - - Array2d adView (nx, ny); - Array2d adDet (nx, ny); - double** ppdView = adView.getArray(); - double** ppdDet = adDet.getArray(); - std::complex** ppcDetValue = new std::complex* [m_nView]; + int iInterpDet = nx; + int iNumInterpDetWithZeros = ProcessSignal::addZeropadFactor (iInterpDet, iZeropad); - int iNumDetWithZeros = ProcessSignal::addZeropadFactor (m_nDet, iZeropad); - double dZeropadRatio = iNumDetWithZeros / static_cast(m_nDet); + double dZeropadRatio = static_cast(iNumInterpDetWithZeros) / static_cast(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** ppcDetValue = new std::complex* [m_nView]; + double dInterpScale = (m_nDet-1) / static_cast(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 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 [iNumDetWithZeros]; - for (unsigned int iD = 0; iD < iNumDetWithZeros; iD++) - ppcDetValue[iView][iD] = std::complex (pcIn[iD].re, pcIn[iD].im); -//ppcDetValue[iView][iD] = std::complex (std::abs(std::complex(pcIn[iD].re, pcIn[iD].im)), 0); - Fourier::shuffleFourierToNaturalOrder (ppcDetValue[iView], iNumDetWithZeros); + + ppcDetValue[iView] = new std::complex [iNumInterpDetWithZeros]; + for (unsigned int iD = 0; iD < iNumInterpDetWithZeros; iD++) + ppcDetValue[iView][iD] = std::complex (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 adView (nx, ny); + Array2d 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++) @@ -791,12 +800,14 @@ Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad 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; @@ -830,7 +841,7 @@ Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double } ppdView[ix][iy] = (phi - m_rotStart) / m_rotInc; - ppdDet[ix][iy] = (r / m_detInc) + iDetCenter; + ppdDet[ix][iy] = (r / dDetInc) + iDetCenter; } } } @@ -840,29 +851,36 @@ Projections::interpolatePolar (ImageFileArray& v, ImageFileArray& vImag, unsigned int nx, unsigned int ny, std::complex** ppcDetValue, double** ppdView, double** ppdDet, unsigned int nView, unsigned int nDet, unsigned int nDetWithZeros, int iInterpolationID) { + typedef std::complex complexValue; + BilinearInterpolator 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 (ppdView[ix][iy]); unsigned int iDet = nearest (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(ppdView[ix][iy]); +#if 1 + std::complex 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(ppdDet[ix][iy]); + int iFloorDet = ::floor (ppdDet[ix][iy]); double dFracDet = ppdDet[ix][iy] - iFloorDet; if (iFloorDet >= 0 && iFloorView >= 0) { @@ -882,6 +900,7 @@ Projections::interpolatePolar (ImageFileArray& v, ImageFileArray& vImag, v3 = v2; else v3 = ppcDetValue[0][iFloorDet+1]; + std::complex vInterp = (1 - dFracView) * (1 - dFracDet) * v1 + dFracView * (1 - dFracDet) * v2 + dFracView * dFracDet * v3 + dFracDet * (1 - dFracView) * v4; @@ -889,12 +908,12 @@ Projections::interpolatePolar (ImageFileArray& v, ImageFileArray& vImag, 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) @@ -1013,14 +1032,14 @@ Projections::interpolateToParallel () const // 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 interp (pdThetaValuesForT, pdRaysumsForT, pProjNew->nView()); + detValues[iD] = interp.interpolate (dViewAngle, &iLastFloor); } } delete pdThetaValuesForT; @@ -1042,9 +1061,9 @@ Projections::interpolateToParallel () const 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 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; @@ -1192,36 +1211,3 @@ ParallelRaysums::getDetPositions (double* pdDetPos) 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])); -} - diff --git a/libctsupport/interpolator.cpp b/libctsupport/interpolator.cpp index bbc7f27..d5167ad 100644 --- a/libctsupport/interpolator.cpp +++ b/libctsupport/interpolator.cpp @@ -2,7 +2,7 @@ ** 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 @@ -136,3 +136,6 @@ CubicSplineInterpolator::interpolate (double x) + + + diff --git a/libctsupport/mathfuncs.cpp b/libctsupport/mathfuncs.cpp index 2ff4e0b..6bead31 100644 --- a/libctsupport/mathfuncs.cpp +++ b/libctsupport/mathfuncs.cpp @@ -2,7 +2,7 @@ ** 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 @@ -135,3 +135,5 @@ vectorNumericStatistics (std::vector vec, const int nPoints, double& min else // Even median = (vec[ (nPoints / 2) - 1 ] + vec[ nPoints / 2 ]) / 2; } + + diff --git a/msvc/ctsim/ctsim.plg b/msvc/ctsim/ctsim.plg index 4c73f26..f45aeea 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\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 @@ -48,16 +48,16 @@ Creating temporary file "C:\DOCUME~1\kevin\LOCALS~1\Temp\RSP5E6.tmp" with conten .\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"

Output Window

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

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

Command Lines

-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 @@ -84,7 +84,7 @@ winmm.lib rpcrt4.lib ws2_32.lib ../libctsim/Debug/libctsim.lib libcmtd.lib ..\.. \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"

Output Window

Linking... diff --git a/src/ctsim.cpp b/src/ctsim.cpp index d42bab2..9d7e870 100644 --- a/src/ctsim.cpp +++ b/src/ctsim.cpp @@ -9,7 +9,7 @@ ** 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 @@ -70,7 +70,7 @@ #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[] = { @@ -293,6 +293,12 @@ CTSimApp::openConfig() 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 @@ -304,7 +310,13 @@ CTSimApp::closeConfig() 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; } diff --git a/src/ctsim.h b/src/ctsim.h index 06d4cc3..4b992da 100644 --- a/src/ctsim.h +++ b/src/ctsim.h @@ -9,7 +9,7 @@ ** 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 @@ -360,6 +360,8 @@ enum { IFMENU_FILTER_IFFT_COLS, IFMENU_FILTER_MAGNITUDE, IFMENU_FILTER_PHASE, + IFMENU_FILTER_REAL, + IFMENU_FILTER_IMAGINARY, IFMENU_FILTER_SHUFFLENATURALTOFOURIERORDER, IFMENU_FILTER_SHUFFLEFOURIERTONATURALORDER, @@ -372,7 +374,7 @@ enum { PLOTMENU_VIEW_SCALE_AUTO, PLOTMENU_VIEW_SCALE_FULL, - GRAPH3D_VIEW_SURFACE, + GRAPH3D_VIEW_WIREFRAME, GRAPH3D_VIEW_COLOR, GRAPH3D_VIEW_LIGHTING, GRAPH3D_VIEW_SMOOTH, diff --git a/src/graph3dview.cpp b/src/graph3dview.cpp index 526ab64..15d560f 100644 --- a/src/graph3dview.cpp +++ b/src/graph3dview.cpp @@ -9,7 +9,7 @@ ** 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 @@ -140,7 +140,7 @@ EVT_MENU(IFMENU_FILE_PROPERTIES, Graph3dFileView::OnProperties) 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) @@ -151,16 +151,15 @@ Graph3dFileView::Graph3dFileView () { 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() @@ -194,7 +193,7 @@ Graph3dFileView::OnCreate (wxDocument *doc, long WXUNUSED(flags) ) 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; } @@ -241,15 +240,7 @@ Graph3dFileView::DrawSurface() glTranslated (-static_cast(nx - 1) / 2, 0.0, -static_cast(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); @@ -261,6 +252,13 @@ Graph3dFileView::DrawSurface() 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); } } @@ -366,13 +364,6 @@ Graph3dFileView::CreateDisplayList() } -void -Graph3dFileView::calculateSurfaceBackground () -{ - if (m_bCalculatedSurfaceBackground) - return; -} - void Graph3dFileView::OnProperties (wxCommandEvent& event) { @@ -392,10 +383,10 @@ Graph3dFileView::OnLighting (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(); } @@ -777,7 +768,7 @@ Graph3dFileView::CreateChildFrame (wxDocument *doc, wxView *view) 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); @@ -802,7 +793,7 @@ Graph3dFileView::CreateChildFrame (wxDocument *doc, wxView *view) subframe->Centre(wxBOTH); wxAcceleratorEntry accelEntries[7]; - accelEntries[0].Set (wxACCEL_CTRL, static_cast('R'), GRAPH3D_VIEW_SURFACE); + accelEntries[0].Set (wxACCEL_CTRL, static_cast('R'), GRAPH3D_VIEW_WIREFRAME); accelEntries[1].Set (wxACCEL_CTRL, static_cast('L'), GRAPH3D_VIEW_COLOR); accelEntries[2].Set (wxACCEL_CTRL, static_cast('G'), GRAPH3D_VIEW_LIGHTING); accelEntries[3].Set (wxACCEL_CTRL, static_cast('M'), GRAPH3D_VIEW_SMOOTH); @@ -901,8 +892,8 @@ Graph3dFileCanvas::OnChar(wxKeyEvent& event) 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); diff --git a/src/graph3dview.h b/src/graph3dview.h index 44adb86..7933f17 100644 --- a/src/graph3dview.h +++ b/src/graph3dview.h @@ -9,7 +9,7 @@ ** 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 @@ -56,7 +56,7 @@ private: bool m_bDoubleBuffer; bool m_bSmooth; bool m_bLighting; - bool m_bSurface; + bool m_bWireframe; bool m_bColor; enum { DISPLAYLIST_COLOR = 1, @@ -69,7 +69,6 @@ private: double m_dColorScaleMax; bool m_bColorScaleMinSet; bool m_bColorScaleMaxSet; - bool m_bCalculatedSurfaceBackground; void Draw(); void DrawSurface(); @@ -109,14 +108,13 @@ public: 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 diff --git a/src/views.cpp b/src/views.cpp index 6c59f8e..4de349e 100644 --- a/src/views.cpp +++ b/src/views.cpp @@ -9,7 +9,7 @@ ** 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 @@ -241,6 +241,8 @@ EVT_MENU(IFMENU_FILTER_IFFT_COLS, ImageFileView::OnIFFTCols) #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 @@ -251,7 +253,7 @@ EVT_MENU(IFMENU_PLOT_HISTOGRAM, ImageFileView::OnPlotHistogram) 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; } @@ -817,6 +819,38 @@ ImageFileView::OnPhase (wxCommandEvent& event) 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) @@ -849,7 +883,6 @@ ImageFileView::CreateChildFrame(wxDocument *doc, wxView *view) 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"); @@ -886,30 +919,33 @@ ImageFileView::CreateChildFrame(wxDocument *doc, wxView *view) 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..."); @@ -952,7 +988,7 @@ ImageFileView::CreateChildFrame(wxDocument *doc, wxView *view) 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"); @@ -960,7 +996,7 @@ ImageFileView::CreateChildFrame(wxDocument *doc, wxView *view) subframe->Centre(wxBOTH); - wxAcceleratorEntry accelEntries[8]; + wxAcceleratorEntry accelEntries[10]; accelEntries[0].Set (wxACCEL_CTRL, static_cast('A'), IFMENU_VIEW_SCALE_AUTO); accelEntries[1].Set (wxACCEL_CTRL, static_cast('U'), IFMENU_VIEW_SCALE_FULL); accelEntries[2].Set (wxACCEL_CTRL, static_cast('E'), IFMENU_VIEW_SCALE_MINMAX); @@ -968,12 +1004,15 @@ ImageFileView::CreateChildFrame(wxDocument *doc, wxView *view) accelEntries[4].Set (wxACCEL_CTRL, static_cast('C'), IFMENU_EDIT_COPY); accelEntries[5].Set (wxACCEL_CTRL, static_cast('X'), IFMENU_EDIT_CUT); accelEntries[6].Set (wxACCEL_CTRL, static_cast('V'), IFMENU_EDIT_PASTE); + int iEntry = 7; +#ifdef HAVE_FFT + accelEntries[iEntry++].Set (wxACCEL_CTRL, static_cast('2'), IFMENU_FILTER_FFT); + accelEntries[iEntry++].Set (wxACCEL_ALT, static_cast('2'), IFMENU_FILTER_IFFT); +#endif #ifdef wxUSE_GLCANVAS - accelEntries[7].Set (wxACCEL_CTRL, static_cast('3'), IFMENU_IMAGE_CONVERT3D); - wxAcceleratorTable accelTable (8, accelEntries); -#else - wxAcceleratorTable accelTable (7, accelEntries); + accelEntries[iEntry++].Set (wxACCEL_CTRL, static_cast('3'), IFMENU_IMAGE_CONVERT3D); #endif + wxAcceleratorTable accelTable (iEntry, accelEntries); subframe->SetAcceleratorTable (accelTable); @@ -1026,6 +1065,17 @@ void 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(); @@ -1142,17 +1192,26 @@ ImageFileView::OnEditPaste (wxCommandEvent& event) 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); @@ -1472,7 +1531,7 @@ ImageFileView::OnPlotFFTRow (wxCommandEvent& event) 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); @@ -1482,8 +1541,8 @@ ImageFileView::OnPlotFFTRow (wxCommandEvent& event) 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); @@ -1575,7 +1634,7 @@ ImageFileView::OnPlotFFTCol (wxCommandEvent& event) 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); @@ -1585,8 +1644,8 @@ ImageFileView::OnPlotFFTCol (wxCommandEvent& event) 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); } diff --git a/src/views.h b/src/views.h index 3f902cf..99551f2 100644 --- a/src/views.h +++ b/src/views.h @@ -9,7 +9,7 @@ ** 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 @@ -62,6 +62,7 @@ private: ImageFileCanvas *m_pCanvas; wxMenu* m_pFileMenu; + wxMenu* m_pFilterMenu; bool m_bMinSpecified; bool m_bMaxSpecified; double m_dMinPixel; @@ -129,6 +130,8 @@ public: void OnMagnitude (wxCommandEvent& event); void OnPhase (wxCommandEvent& event); + void OnReal (wxCommandEvent& event); + void OnImaginary (wxCommandEvent& event); void OnScaleAuto (wxCommandEvent& event); void OnScaleMinMax (wxCommandEvent& event);