vtkImageSobel3D.cxx 9.13 KB
Newer Older
1 2 3 4 5
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkImageSobel3D.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.
13 14 15

=========================================================================*/
#include "vtkImageSobel3D.h"
16

Brad King's avatar
Brad King committed
17
#include "vtkImageData.h"
18 19 20 21
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkObjectFactory.h"
#include "vtkStreamingDemandDrivenPipeline.h"
22

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

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

//----------------------------------------------------------------------------
// Construct an instance of vtkImageSobel3D fitler.
vtkImageSobel3D::vtkImageSobel3D()
{
Charles Law's avatar
Charles Law committed
31 32 33 34 35 36 37
  this->KernelSize[0] = 3;
  this->KernelSize[1] = 3;
  this->KernelSize[2] = 3;
  this->KernelMiddle[0] = 1;
  this->KernelMiddle[1] = 1;
  this->KernelMiddle[2] = 1;
  this->HandleBoundaries = 1;
38 39 40
}

//----------------------------------------------------------------------------
41
void vtkImageSobel3D::PrintSelf(ostream& os, vtkIndent indent)
42
{
Brad King's avatar
Brad King committed
43
  this->Superclass::PrintSelf(os, indent);
44 45 46
}

//----------------------------------------------------------------------------
47
int vtkImageSobel3D::RequestInformation (vtkInformation *request,
48 49
                                         vtkInformationVector **inputVector,
                                         vtkInformationVector *outputVector)
50
{
51 52
  int retval =
    this->Superclass::RequestInformation(request, inputVector, outputVector);
53
  vtkInformation *outInfo = outputVector->GetInformationObject(0);
54
  vtkDataObject::SetPointDataActiveScalarInfo(outInfo, VTK_DOUBLE, 3);
55
  return retval;
56 57 58 59
}

//----------------------------------------------------------------------------
// This execute method handles boundaries.
60
// it handles boundaries. Pixels are just replicated to get values
61 62
// out of extent.
template <class T>
Ken Martin's avatar
Ken Martin committed
63
void vtkImageSobel3DExecute(vtkImageSobel3D *self,
64 65
                            vtkImageData *inData, T *inPtr,
                            vtkImageData *outData, int *outExt,
66
                            double *outPtr, int id, vtkInformation *inInfo)
67
{
Ken Martin's avatar
Ken Martin committed
68
  double r0, r1, r2, *r;
69 70 71
  // For looping though output (and input) pixels.
  int min0, max0, min1, max1, min2, max2;
  int outIdx0, outIdx1, outIdx2;
72
  vtkIdType outInc0, outInc1, outInc2;
Ken Martin's avatar
Ken Martin committed
73
  double *outPtr0, *outPtr1, *outPtr2, *outPtrV;
74
  vtkIdType inInc0, inInc1, inInc2;
75 76
  T *inPtr0, *inPtr1, *inPtr2;
  // For sobel function convolution (Left Right incs for each axis)
77
  vtkIdType inInc0L, inInc0R, inInc1L, inInc1R, inInc2L, inInc2R;
78
  T *inPtrL, *inPtrR;
Ken Martin's avatar
Ken Martin committed
79
  double sum;
80
  // Boundary of input image
81 82
  int inWholeMin0, inWholeMax0, inWholeMin1, inWholeMax1;
  int inWholeMin2, inWholeMax2;
83
  int inWholeExt[6];
84 85
  unsigned long count = 0;
  unsigned long target;
86

87
  // Get boundary information
88 89 90 91 92 93 94
  inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), inWholeExt);
  inWholeMin0 = inWholeExt[0];
  inWholeMax0 = inWholeExt[1];
  inWholeMin1 = inWholeExt[2];
  inWholeMax1 = inWholeExt[3];
  inWholeMin2 = inWholeExt[4];
  inWholeMax2 = inWholeExt[5];
95

96
  // Get information to march through data (skip component)
97 98
  inData->GetIncrements(inInc0, inInc1, inInc2);
  outData->GetIncrements(outInc0, outInc1, outInc2);
Charles Law's avatar
Charles Law committed
99 100 101
  min0 = outExt[0];   max0 = outExt[1];
  min1 = outExt[2];   max1 = outExt[3];
  min2 = outExt[4];   max2 = outExt[5];
102

103
  // We want the input pixel to correspond to output
104
  inPtr = static_cast<T *>(inData->GetScalarPointer(min0,min1,min2));
105

Will Schroeder's avatar
Will Schroeder committed
106
  // The data spacing is important for computing the gradient.
107
  // Scale so it has the same range as gradient.
Charles Law's avatar
Charles Law committed
108 109 110 111
  r = inData->GetSpacing();
  r0 = 0.060445 / r[0];
  r1 = 0.060445 / r[1];
  r2 = 0.060445 / r[2];
112

113
  target = static_cast<unsigned long>((max2-min2+1)*(max1-min1+1)/50.0);
114 115
  target++;

116 117 118 119 120
  // loop through pixels of output
  outPtr2 = outPtr;
  inPtr2 = inPtr;
  for (outIdx2 = min2; outIdx2 <= max2; ++outIdx2)
    {
121 122
    inInc2L = (outIdx2 == inWholeMin2) ? 0 : -inInc2;
    inInc2R = (outIdx2 == inWholeMax2) ? 0 : inInc2;
123 124 125

    outPtr1 = outPtr2;
    inPtr1 = inPtr2;
126
    for (outIdx1 = min1; !self->AbortExecute && outIdx1 <= max1; ++outIdx1)
127
      {
128
      if (!id)
129 130 131 132 133 134 135
        {
        if (!(count%target))
          {
          self->UpdateProgress(count/(50.0*target));
          }
        count++;
        }
136 137
      inInc1L = (outIdx1 == inWholeMin1) ? 0 : -inInc1;
      inInc1R = (outIdx1 == inWholeMax1) ? 0 : inInc1;
138 139 140 141

      outPtr0 = outPtr1;
      inPtr0 = inPtr1;
      for (outIdx0 = min0; outIdx0 <= max0; ++outIdx0)
142 143 144 145 146 147 148 149 150 151
        {
        inInc0L = (outIdx0 == inWholeMin0) ? 0 : -inInc0;
        inInc0R = (outIdx0 == inWholeMax0) ? 0 : inInc0;

        // compute vector.
        outPtrV = outPtr0;
        // 12 Plane
        inPtrL = inPtr0 + inInc0L;
        inPtrR = inPtr0 + inInc0R;
        sum = 2.0 * (*inPtrR - *inPtrL);
152
        sum += static_cast<double>(inPtrR[inInc1L] + inPtrR[inInc1R]
153
                + inPtrR[inInc2L] + inPtrR[inInc2R]);
154
        sum += static_cast<double>(0.586 * (inPtrR[inInc1L+inInc2L] + inPtrR[inInc1L+inInc2R]
155
                         + inPtrR[inInc1R+inInc2L] + inPtrR[inInc1R+inInc2R]));
156
        sum -= static_cast<double>(inPtrL[inInc1L] + inPtrL[inInc1R]
157
                + inPtrL[inInc2L] + inPtrL[inInc2R]);
158
        sum -= static_cast<double>(0.586 * (inPtrL[inInc1L+inInc2L] + inPtrL[inInc1L+inInc2R]
159 160 161 162 163 164 165
                         + inPtrL[inInc1R+inInc2L] + inPtrL[inInc1R+inInc2R]));
        *outPtrV = sum * r0;
        ++outPtrV;
        // 02 Plane
        inPtrL = inPtr0 + inInc1L;
        inPtrR = inPtr0 + inInc1R;
        sum = 2.0 * (*inPtrR - *inPtrL);
166
        sum += static_cast<double>(inPtrR[inInc0L] + inPtrR[inInc0R]
167
                + inPtrR[inInc2L] + inPtrR[inInc2R]);
168
        sum += static_cast<double>(0.586 * (inPtrR[inInc0L+inInc2L] + inPtrR[inInc0L+inInc2R]
169
                         + inPtrR[inInc0R+inInc2L] + inPtrR[inInc0R+inInc2R]));
170
        sum -= static_cast<double>(inPtrL[inInc0L] + inPtrL[inInc0R]
171
                + inPtrL[inInc2L] + inPtrL[inInc2R]);
172
        sum -= static_cast<double>(0.586 * (inPtrL[inInc0L+inInc2L] + inPtrL[inInc0L+inInc2R]
173 174 175 176 177 178 179
                         + inPtrL[inInc0R+inInc2L] + inPtrL[inInc0R+inInc2R]));
        *outPtrV = sum * r1;
        ++outPtrV;
        // 01 Plane
        inPtrL = inPtr0 + inInc2L;
        inPtrR = inPtr0 + inInc2R;
        sum = 2.0 * (*inPtrR - *inPtrL);
180
        sum += static_cast<double>(inPtrR[inInc0L] + inPtrR[inInc0R]
181
                + inPtrR[inInc1L] + inPtrR[inInc1R]);
182
        sum += static_cast<double>(0.586 * (inPtrR[inInc0L+inInc1L] + inPtrR[inInc0L+inInc1R]
183
                         + inPtrR[inInc0R+inInc1L] + inPtrR[inInc0R+inInc1R]));
184
        sum -= static_cast<double>(inPtrL[inInc0L] + inPtrL[inInc0R]
185
                + inPtrL[inInc1L] + inPtrL[inInc1R]);
186
        sum -= static_cast<double>(0.586 * (inPtrL[inInc0L+inInc1L] + inPtrL[inInc0L+inInc1R]
187
                         + inPtrL[inInc0R+inInc1L] + inPtrL[inInc0R+inInc1R]));
Ken Martin's avatar
Ken Martin committed
188
        *outPtrV = static_cast<double>(sum * r2);
189 190 191 192 193
        ++outPtrV;

        outPtr0 += outInc0;
        inPtr0 += inInc0;
        }
194 195 196 197 198 199 200 201 202 203
      outPtr1 += outInc1;
      inPtr1 += inInc1;
      }
    outPtr2 += outInc2;
    inPtr2 += inInc2;
    }
}

//----------------------------------------------------------------------------
// This method contains a switch statement that calls the correct
Charles Law's avatar
Charles Law committed
204
// templated function for the input Data type.  The output Data
Ken Martin's avatar
Ken Martin committed
205
// must be of type double.  This method does handle boundary conditions.
206
// The third axis is the component axis for the output.
207 208 209 210 211 212 213
void vtkImageSobel3D::ThreadedRequestData(
  vtkInformation *vtkNotUsed(request),
  vtkInformationVector **inputVector,
  vtkInformationVector *vtkNotUsed(outputVector),
  vtkImageData ***inData,
  vtkImageData **outData,
  int outExt[6], int id)
214
{
Charles Law's avatar
Charles Law committed
215
  void *inPtr, *outPtr;
216 217 218 219 220 221 222 223
  int inExt[6], wholeExt[6];

  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), wholeExt);
  this->InternalRequestUpdateExtent(inExt, outExt, wholeExt);

  inPtr = inData[0][0]->GetScalarPointerForExtent(inExt);
  outPtr = outData[0]->GetScalarPointerForExtent(outExt);
Charles Law's avatar
Charles Law committed
224 225

  // this filter cannot handle multi component input.
226
  if (inData[0][0]->GetNumberOfScalarComponents() != 1)
Charles Law's avatar
Charles Law committed
227 228 229
    {
    vtkWarningMacro("Expecting input with only one compenent.\n");
    }
230

Ken Martin's avatar
Ken Martin committed
231
  // this filter expects that output is type double.
232
  if (outData[0]->GetScalarType() != VTK_DOUBLE)
233 234
    {
    vtkErrorMacro(<< "Execute: output ScalarType, "
235
                  << vtkImageScalarTypeNameMacro(outData[0]->GetScalarType())
Ken Martin's avatar
Ken Martin committed
236
                  << ", must be double");
237 238
    return;
    }
239

240
  switch (inData[0][0]->GetScalarType())
241
    {
242
    vtkTemplateMacro(
243 244 245
      vtkImageSobel3DExecute(this, inData[0][0],
                             static_cast<VTK_TT *>(inPtr), outData[0], outExt,
                             static_cast<double *>(outPtr),id, inInfo));
246 247 248 249 250
    default:
      vtkErrorMacro(<< "Execute: Unknown ScalarType");
      return;
    }
}