r641: no message
authorKevin M. Rosenberg <kevin@rosenberg.net>
Wed, 21 Mar 2001 21:45:31 +0000 (21:45 +0000)
committerKevin M. Rosenberg <kevin@rosenberg.net>
Wed, 21 Mar 2001 21:45:31 +0000 (21:45 +0000)
18 files changed:
include/fourier.h
include/imagefile.h
include/interpolator.h
include/projections.h
libctgraphics/bresenham.cpp
libctsim/fourier.cpp
libctsim/imagefile.cpp
libctsim/procsignal.cpp
libctsim/projections.cpp
libctsupport/interpolator.cpp
libctsupport/mathfuncs.cpp
msvc/ctsim/ctsim.plg
src/ctsim.cpp
src/ctsim.h
src/graph3dview.cpp
src/graph3dview.h
src/views.cpp
src/views.h

index 5fce2f6..b956fcc 100644 (file)
@@ -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<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
 };
 
 
index 140955c..56c255b 100644 (file)
@@ -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);
 
index def913e..cea4961 100644 (file)
@@ -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 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;
+  }
+};
+
index 9c0cce0..d2fe89c 100644 (file)
@@ -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<double>** ppcDetValue,
     double** ppdDet, double** ppdView, unsigned int nView, unsigned int nDet, unsigned int nDetWithZeros, 
     int iInterpolate);
index e9936fe..070f44e 100644 (file)
-/* ; 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);
+}
+  
+  
index f7357cd..6fec9f5 100644 (file)
@@ -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<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;
-}
-
-
index 1c8a063..bab87bd 100644 (file)
@@ -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<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
 //
@@ -739,40 +686,22 @@ ImageFile::scaleImage (ImageFile& result) const
   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;
 }
 
@@ -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<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())
@@ -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<double>* pcRow = new std::complex<double> [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<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())
@@ -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<double>* pcCol = new std::complex<double> [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;
   
index 2d6600d..1c4f25f 100644 (file)
@@ -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++) 
index 1aee958..4d4a02d 100644 (file)
@@ -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<double>** ppcDetValue = new std::complex<double>* [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<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];
@@ -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<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++)
@@ -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<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) { 
@@ -882,6 +900,7 @@ Projections::interpolatePolar (ImageFileArray& v, ImageFileArray& vImag,
             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;
@@ -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<double> 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<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;
@@ -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]));
-}
-
index bbc7f27..d5167ad 100644 (file)
@@ -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)
 
 
 
+  
+
+
index 2ff4e0b..6bead31 100644 (file)
@@ -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<double> vec, const int nPoints, double& min
   else        // Even
     median = (vec[ (nPoints / 2) - 1 ] + vec[ nPoints / 2 ]) / 2;
 }
+
+
index 4c73f26..f45aeea 100644 (file)
@@ -6,13 +6,13 @@
 --------------------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
@@ -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"
 <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
@@ -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"
 <h3>Output Window</h3>
 Linking...
 
index d42bab2..9d7e870 100644 (file)
@@ -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;
 }
 
index 06d4cc3..4b992da 100644 (file)
@@ -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,
index 526ab64..15d560f 100644 (file)
@@ -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<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);
@@ -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<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);
@@ -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);
index 44adb86..7933f17 100644 (file)
@@ -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
index 6c59f8e..4de349e 100644 (file)
@@ -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&lter");
+  menu_bar->Append(m_pFilterMenu, "Fi&lter");
   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<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);
@@ -968,12 +1004,15 @@ ImageFileView::CreateChildFrame(wxDocument *doc, wxView *view)
   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);
   
@@ -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);
     }
     
index 3f902cf..99551f2 100644 (file)
@@ -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);