vtkImageLaplacian.cxx 7.17 KB
Newer Older
Charles Law's avatar
Charles Law committed
1 2 3 4 5
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkImageLaplacian.cxx

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.
Charles Law's avatar
Charles Law committed
13 14 15

=========================================================================*/
#include "vtkImageLaplacian.h"
16 17

#include "vtkImageData.h"
18 19
#include "vtkInformation.h"
#include "vtkInformationVector.h"
20
#include "vtkObjectFactory.h"
21
#include "vtkStreamingDemandDrivenPipeline.h"
22

Sean McBride's avatar
Sean McBride committed
23
#include <cmath>
24

Brad King's avatar
Brad King committed
25
vtkStandardNewMacro(vtkImageLaplacian);
Charles Law's avatar
Charles Law committed
26 27 28 29 30

//----------------------------------------------------------------------------
// Construct an instance of vtkImageLaplacian fitler.
vtkImageLaplacian::vtkImageLaplacian()
{
31
  this->Dimensionality = 2;
Charles Law's avatar
Charles Law committed
32 33 34
}

//----------------------------------------------------------------------------
35
void vtkImageLaplacian::PrintSelf(ostream& os, vtkIndent indent)
Charles Law's avatar
Charles Law committed
36
{
Brad King's avatar
Brad King committed
37
  this->Superclass::PrintSelf(os, indent);
38
  os << indent << "Dimensionality: " << this->Dimensionality;
Charles Law's avatar
Charles Law committed
39 40 41
}

//----------------------------------------------------------------------------
42
// Just clip the request.  The subclass may need to overwrite this method.
43
int vtkImageLaplacian::RequestUpdateExtent (
44 45 46
  vtkInformation * vtkNotUsed(request),
  vtkInformationVector **inputVector,
  vtkInformationVector *outputVector)
Charles Law's avatar
Charles Law committed
47
{
48 49 50 51
  // get the info objects
  vtkInformation* outInfo = outputVector->GetInformationObject(0);
  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);

52
  int idx;
53
  int wholeExtent[6], inUExt[6];
54

55 56 57
  inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), wholeExtent);
  outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), inUExt);

58 59
  // update and Clip
  for (idx = 0; idx < 3; ++idx)
Charles Law's avatar
Charles Law committed
60
    {
61 62 63
    --inUExt[idx*2];
    ++inUExt[idx*2+1];
    if (inUExt[idx*2] < wholeExtent[idx*2])
Charles Law's avatar
Charles Law committed
64
      {
65
      inUExt[idx*2] = wholeExtent[idx*2];
Charles Law's avatar
Charles Law committed
66
      }
67
    if (inUExt[idx*2] > wholeExtent[idx*2 + 1])
Charles Law's avatar
Charles Law committed
68
      {
69
      inUExt[idx*2] = wholeExtent[idx*2 + 1];
70
      }
71
    if (inUExt[idx*2+1] < wholeExtent[idx*2])
72
      {
73
      inUExt[idx*2+1] = wholeExtent[idx*2];
74
      }
75
    if (inUExt[idx*2 + 1] > wholeExtent[idx*2 + 1])
76
      {
77
      inUExt[idx*2 + 1] = wholeExtent[idx*2 + 1];
Charles Law's avatar
Charles Law committed
78 79
      }
    }
80
  inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), inUExt, 6);
Charles Law's avatar
Charles Law committed
81

82 83
  return 1;
}
Charles Law's avatar
Charles Law committed
84 85 86

//----------------------------------------------------------------------------
// This execute method handles boundaries.
87
// it handles boundaries. Pixels are just replicated to get values
Charles Law's avatar
Charles Law committed
88 89
// out of extent.
template <class T>
Ken Martin's avatar
Ken Martin committed
90 91 92 93
void vtkImageLaplacianExecute(vtkImageLaplacian *self,
                              vtkImageData *inData, T *inPtr,
                              vtkImageData *outData, T *outPtr,
                              int outExt[6], int id)
Charles Law's avatar
Charles Law committed
94
{
95 96
  int idxC, idxX, idxY, idxZ;
  int maxC, maxX, maxY, maxZ;
97 98
  vtkIdType inIncX, inIncY, inIncZ;
  vtkIdType outIncX, outIncY, outIncZ;
99 100 101
  unsigned long count = 0;
  unsigned long target;
  int axesNum;
102 103
  int *wholeExtent;
  vtkIdType *inIncs;
Ken Martin's avatar
Ken Martin committed
104
  double r[3], d, sum;
105
  int useZMin, useZMax, useYMin, useYMax, useXMin, useXMax;
106

107 108 109
  // find the region to loop over
  maxC = inData->GetNumberOfScalarComponents();
  maxX = outExt[1] - outExt[0];
110
  maxY = outExt[3] - outExt[2];
111
  maxZ = outExt[5] - outExt[4];
112
  target = static_cast<unsigned long>((maxZ+1)*(maxY+1)/50.0);
113 114 115 116
  target++;

  // Get the dimensionality of the gradient.
  axesNum = self->GetDimensionality();
117 118

  // Get increments to march through data
119 120 121 122 123 124 125 126 127 128 129
  inData->GetContinuousIncrements(outExt, inIncX, inIncY, inIncZ);
  outData->GetContinuousIncrements(outExt, outIncX, outIncY, outIncZ);

  // The data spacing is important for computing the Laplacian.
  // Divided by dx twice (second derivative).
  inData->GetSpacing(r);
  r[0] = 1.0 / (r[0] * r[0]);
  r[1] = 1.0 / (r[1] * r[1]);
  r[2] = 1.0 / (r[2] * r[2]);

  // get some other info we need
130 131
  inIncs = inData->GetIncrements();
  wholeExtent = inData->GetExtent();
132

Sven Buijssen's avatar
Sven Buijssen committed
133
  // Loop through output pixels
134
  for (idxZ = 0; idxZ <= maxZ; idxZ++)
Charles Law's avatar
Charles Law committed
135
    {
136 137
    useZMin = ((idxZ + outExt[4]) <= wholeExtent[4]) ? 0 : -inIncs[2];
    useZMax = ((idxZ + outExt[4]) >= wholeExtent[5]) ? 0 : inIncs[2];
138
    for (idxY = 0; !self->AbortExecute && idxY <= maxY; idxY++)
Charles Law's avatar
Charles Law committed
139
      {
140
      if (!id)
141 142 143 144 145 146 147
        {
        if (!(count%target))
          {
          self->UpdateProgress(count/(50.0*target));
          }
        count++;
        }
148 149 150
      useYMin = ((idxY + outExt[2]) <= wholeExtent[2]) ? 0 : -inIncs[1];
      useYMax = ((idxY + outExt[2]) >= wholeExtent[3]) ? 0 : inIncs[1];
      for (idxX = 0; idxX <= maxX; idxX++)
151 152 153 154 155 156 157
        {
        useXMin = ((idxX + outExt[0]) <= wholeExtent[0]) ? 0 : -inIncs[0];
        useXMax = ((idxX + outExt[0]) >= wholeExtent[1]) ? 0 : inIncs[0];
        for (idxC = 0; idxC < maxC; idxC++)
          {
          // do X axis
          d = -2.0*(*inPtr);
158 159
          d += static_cast<double>(inPtr[useXMin]);
          d += static_cast<double>(inPtr[useXMax]);
160 161 162 163
          sum = d * r[0];

          // do y axis
          d = -2.0*(*inPtr);
164 165
          d += static_cast<double>(inPtr[useYMin]);
          d += static_cast<double>(inPtr[useYMax]);
166 167 168 169 170
          sum = sum + d * r[1];
          if (axesNum == 3)
            {
            // do z axis
            d = -2.0*(*inPtr);
171 172
            d += static_cast<double>(inPtr[useZMin]);
            d += static_cast<double>(inPtr[useZMax]);
173 174
            sum = sum + d * r[2];
            }
175
          *outPtr = static_cast<T>(sum);
176 177 178 179
          inPtr++;
          outPtr++;
          }
        }
180 181
      outPtr += outIncY;
      inPtr += inIncY;
Charles Law's avatar
Charles Law committed
182
      }
183 184
    outPtr += outIncZ;
    inPtr += inIncZ;
Charles Law's avatar
Charles Law committed
185 186 187 188 189
    }
}

//----------------------------------------------------------------------------
// This method contains a switch statement that calls the correct
190 191
// templated function for the input data type.  The output data
// must match input type.  This method does handle boundary conditions.
192 193
void vtkImageLaplacian::ThreadedRequestData(
  vtkInformation *vtkNotUsed(request),
194 195
  vtkInformationVector **vtkNotUsed(inputVector),
  vtkInformationVector *vtkNotUsed(outputVector),
196
  vtkImageData ***inData,
197 198
  vtkImageData **outData,
  int outExt[6], int id)
Charles Law's avatar
Charles Law committed
199
{
200 201 202
  void *inPtr = inData[0][0]->GetScalarPointerForExtent(outExt);
  void *outPtr = outData[0]->GetScalarPointerForExtent(outExt);

203
  // this filter expects that input is the same type as output.
204 205
  if (inData[0][0]->GetScalarType() !=
      outData[0]->GetScalarType())
Charles Law's avatar
Charles Law committed
206
    {
207
    vtkErrorMacro(<< "Execute: input ScalarType, "
208 209
      << inData[0][0]->GetScalarType() << ", must match out ScalarType "
      << outData[0]->GetScalarType());
Charles Law's avatar
Charles Law committed
210 211
    return;
    }
212

213
  switch (inData[0][0]->GetScalarType())
Charles Law's avatar
Charles Law committed
214
    {
215 216
    vtkTemplateMacro(
      vtkImageLaplacianExecute( this, inData[0][0],
217 218
                                static_cast<VTK_TT *>(inPtr), outData[0],
                                static_cast<VTK_TT *>(outPtr),
219
                                outExt, id));
Charles Law's avatar
Charles Law committed
220 221 222 223 224 225
    default:
      vtkErrorMacro(<< "Execute: Unknown ScalarType");
      return;
    }
}