vtkImageButterworthLowPass.cxx 5.98 KB
Newer Older
1
2
3
/*=========================================================================

  Program:   Visualization Toolkit
Charles Law's avatar
Charles Law committed
4
  Module:    vtkImageButterworthLowPass.cxx
5

6
  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7
8
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
Will Schroeder's avatar
Will Schroeder committed
9

10
11
     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12
     PURPOSE.  See the above copyright notice for more information.
13
14

=========================================================================*/
Charles Law's avatar
Charles Law committed
15
#include "vtkImageButterworthLowPass.h"
16
17

#include "vtkImageData.h"
18
19
#include "vtkObjectFactory.h"

Brad King's avatar
Brad King committed
20
#include <math.h>
21

Ken Martin's avatar
Ken Martin committed
22
vtkCxxRevisionMacro(vtkImageButterworthLowPass, "1.24");
Brad King's avatar
Brad King committed
23
vtkStandardNewMacro(vtkImageButterworthLowPass);
24
25

//----------------------------------------------------------------------------
Charles Law's avatar
Charles Law committed
26
vtkImageButterworthLowPass::vtkImageButterworthLowPass()
27
{
Ken Martin's avatar
Ken Martin committed
28
  this->CutOff[0] = this->CutOff[1] = this->CutOff[2] = VTK_FLOAT_MAX;
Charles Law's avatar
Charles Law committed
29
  this->Order = 1;
30
31
32
33
}


//----------------------------------------------------------------------------
Ken Martin's avatar
Ken Martin committed
34
void vtkImageButterworthLowPass::SetXCutOff(double cutOff)
35
{
36
  if (cutOff == this->CutOff[0])
37
    {
38
    return;
39
    }
40
  this->CutOff[0] = cutOff;
41
42
43
  this->Modified();
}
//----------------------------------------------------------------------------
Ken Martin's avatar
Ken Martin committed
44
void vtkImageButterworthLowPass::SetYCutOff(double cutOff)
45
{
46
  if (cutOff == this->CutOff[1])
47
    {
48
    return;
49
    }
50
  this->CutOff[1] = cutOff;
51
52
53
  this->Modified();
}
//----------------------------------------------------------------------------
Ken Martin's avatar
Ken Martin committed
54
void vtkImageButterworthLowPass::SetZCutOff(double cutOff)
55
{
56
  if (cutOff == this->CutOff[2])
57
    {
58
    return;
59
    }
60
61
  this->CutOff[2] = cutOff;
  this->Modified();
62
63
64
}

//----------------------------------------------------------------------------
65
void vtkImageButterworthLowPass::ThreadedExecute(vtkImageData *inData, 
66
67
                                                 vtkImageData *outData,
                                                 int ext[6], int id)
68
{
69
70
  int idx0, idx1, idx2;
  int min0, max0;
Ken Martin's avatar
Ken Martin committed
71
72
  double *inPtr;
  double *outPtr;
73
  int *wholeExtent;
Ken Martin's avatar
Ken Martin committed
74
  double *spacing;
75
76
  int inInc0, inInc1, inInc2;
  int outInc0, outInc1, outInc2;
Ken Martin's avatar
Ken Martin committed
77
  double temp0, temp1, temp2, mid0, mid1, mid2;
78
  // normalization factors
Ken Martin's avatar
Ken Martin committed
79
80
  double norm0, norm1, norm2;
  double sum1, sum0;
81
82
  unsigned long count = 0;
  unsigned long target;
83
84
85

  // Error checking
  if (inData->GetNumberOfScalarComponents() != 2)
86
    {
87
    vtkErrorMacro("Expecting 2 components not " 
88
                  << inData->GetNumberOfScalarComponents());
89
90
    return;
    }
Ken Martin's avatar
Ken Martin committed
91
92
  if (inData->GetScalarType() != VTK_DOUBLE || 
      outData->GetScalarType() != VTK_DOUBLE)
93
    {
Ken Martin's avatar
Ken Martin committed
94
    vtkErrorMacro("Expecting input and output to be of type double");
95
96
    return;
    }
97
  
98
  wholeExtent = this->GetInput()->GetWholeExtent();
99
  spacing = inData->GetSpacing();
100

Ken Martin's avatar
Ken Martin committed
101
102
  inPtr = (double *)(inData->GetScalarPointerForExtent(ext));
  outPtr = (double *)(outData->GetScalarPointerForExtent(ext));
103
104
105
106
107
108

  inData->GetContinuousIncrements(ext, inInc0, inInc1, inInc2);
  outData->GetContinuousIncrements(ext, outInc0, outInc1, outInc2);  
  
  min0 = ext[0];
  max0 = ext[1];
Ken Martin's avatar
Ken Martin committed
109
110
111
  mid0 = (double)(wholeExtent[0] + wholeExtent[1] + 1) / 2.0;
  mid1 = (double)(wholeExtent[2] + wholeExtent[3] + 1) / 2.0;
  mid2 = (double)(wholeExtent[4] + wholeExtent[5] + 1) / 2.0;
112
113
  if ( this->CutOff[0] == 0.0)
    {
Ken Martin's avatar
Ken Martin committed
114
    norm0 = VTK_FLOAT_MAX;
115
116
117
118
119
120
    }
  else
    {
    norm0 = 1.0 / ((spacing[0] * 2.0 * mid0) * this->CutOff[0]);
    }
  if ( this->CutOff[1] == 0.0)
121
    {
Ken Martin's avatar
Ken Martin committed
122
    norm1 = VTK_FLOAT_MAX;
123
124
125
126
127
128
129
    }
  else
    {
    norm1 = 1.0 / ((spacing[1] * 2.0 * mid1) * this->CutOff[1]);
    }
  if ( this->CutOff[2] == 0.0)
    {
Ken Martin's avatar
Ken Martin committed
130
    norm2 = VTK_FLOAT_MAX;
131
132
133
134
135
136
    }
  else
    {
    norm2 = 1.0 / ((spacing[2] * 2.0 * mid2) * this->CutOff[2]);
    }
  
137
138
139
  target = (unsigned long)((ext[5]-ext[4]+1)*(ext[3]-ext[2]+1)/50.0);
  target++;

140
141
142
143
  // loop over all the pixels (keeping track of normalized distance to origin.
  for (idx2 = ext[4]; idx2 <= ext[5]; ++idx2)
    {
    // distance to min (this axis' contribution)
Ken Martin's avatar
Ken Martin committed
144
    temp2 = (double)idx2;
145
    // Wrap back to 0.
146
    if (temp2 > mid2)
147
      {
148
      temp2 = mid2 + mid2 - temp2;
149
      }
150
151
152
    // Convert location into normalized cycles/world unit
    temp2 = temp2 * norm2;

Ken Martin's avatar
Ken Martin committed
153
    for (idx1 = ext[2]; !this->AbortExecute && idx1 <= ext[3]; ++idx1)
154
      {
155
      if (!id) 
156
157
158
159
160
161
162
        {
        if (!(count%target))
          {
          this->UpdateProgress(count/(50.0*target));
          }
        count++;
        }
163
      // distance to min (this axis' contribution)
Ken Martin's avatar
Ken Martin committed
164
      temp1 = (double)idx1;
165
166
      // Wrap back to 0.
      if (temp1 > mid1)
167
168
169
        {
        temp1 = mid1 + mid1 - temp1;
        }
170
171
172
173
174
      // Convert location into cycles / world unit
      temp1 = temp1 * norm1;
      sum1 = temp2 * temp2 + temp1 * temp1;
      
      for (idx0 = min0; idx0 <= max0; ++idx0)
175
176
        {
        // distance to min (this axis' contribution)
Ken Martin's avatar
Ken Martin committed
177
        temp0 = (double)idx0;
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
        // Wrap back to 0.
        if (temp0 > mid0)
          {
          temp0 = mid0 + mid0 - temp0;
          }
        // Convert location into cycles / world unit
        temp0 = temp0 * norm0;
        sum0 = sum1 + temp0 * temp0;
        
        // compute Butterworth1D function from sum = d^2
        if (this->Order == 1)
          {
          sum0 = 1.0 / (1.0 + sum0);
          }
        else
          {
Ken Martin's avatar
Ken Martin committed
194
          sum0 = 1.0 / (1.0 + pow(sum0, (double)this->Order));
195
196
197
198
199
200
201
202
          }     
        
        // real component
        *outPtr++ = *inPtr++ * sum0;
        // imaginary component  
        *outPtr++ = *inPtr++ * sum0;
        
        }
203
204
      inPtr += inInc1;
      outPtr += outInc1;
205
      }
206
207
    inPtr += inInc2;
    outPtr += outInc2;    
208
209
210
    }
}

211
void vtkImageButterworthLowPass::PrintSelf(ostream& os, vtkIndent indent)
212
{
Brad King's avatar
Brad King committed
213
  this->Superclass::PrintSelf(os,indent);
214

215
216
217
218
219
220
  os << indent << "Order: " << this->Order << "\n";

  os << indent << "CutOff: ( "
     << this->CutOff[0] << ", "
     << this->CutOff[1] << ", "
     << this->CutOff[2] << " )\n";
221
}
222