You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
466 lines
15 KiB
466 lines
15 KiB
/*=========================================================================
|
|
|
|
Program: Visualization Toolkit
|
|
Module: $RCSfile: vtkImageGaussianSmooth.cxx,v $
|
|
|
|
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
|
|
All rights reserved.
|
|
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
|
|
|
|
This software is distributed WITHOUT ANY WARRANTY; without even
|
|
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
|
|
PURPOSE. See the above copyright notice for more information.
|
|
|
|
=========================================================================*/
|
|
#include "vtkImageGaussianSmooth.h"
|
|
|
|
#include "vtkImageData.h"
|
|
#include "vtkInformation.h"
|
|
#include "vtkInformationVector.h"
|
|
#include "vtkObjectFactory.h"
|
|
#include "vtkStreamingDemandDrivenPipeline.h"
|
|
|
|
#include <math.h>
|
|
|
|
vtkCxxRevisionMacro(vtkImageGaussianSmooth, "$Revision: 1.46 $");
|
|
vtkStandardNewMacro(vtkImageGaussianSmooth);
|
|
|
|
//----------------------------------------------------------------------------
|
|
vtkImageGaussianSmooth::vtkImageGaussianSmooth()
|
|
{
|
|
this->Dimensionality = 3; // note: this overrides Standard deviation.
|
|
this->StandardDeviations[0] = 2.0;
|
|
this->StandardDeviations[1] = 2.0;
|
|
this->StandardDeviations[2] = 2.0;
|
|
this->RadiusFactors[0] = 1.5;
|
|
this->RadiusFactors[1] = 1.5;
|
|
this->RadiusFactors[2] = 1.5;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
vtkImageGaussianSmooth::~vtkImageGaussianSmooth()
|
|
{
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
void vtkImageGaussianSmooth::PrintSelf(ostream& os, vtkIndent indent)
|
|
{
|
|
this->Superclass::PrintSelf(os, indent);
|
|
|
|
// int idx;
|
|
|
|
//os << indent << "BoundaryRescale: " << this->BoundaryRescale << "\n";
|
|
|
|
os << indent << "Dimensionality: " << this->Dimensionality << "\n";
|
|
|
|
os << indent << "RadiusFactors: ( "
|
|
<< this->RadiusFactors[0] << ", "
|
|
<< this->RadiusFactors[1] << ", "
|
|
<< this->RadiusFactors[2] << " )\n";
|
|
|
|
os << indent << "StandardDeviations: ( "
|
|
<< this->StandardDeviations[0] << ", "
|
|
<< this->StandardDeviations[1] << ", "
|
|
<< this->StandardDeviations[2] << " )\n";
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
void vtkImageGaussianSmooth::ComputeKernel(double *kernel, int min, int max,
|
|
double std)
|
|
{
|
|
int x;
|
|
double sum;
|
|
|
|
// handle special case
|
|
if (std == 0.0)
|
|
{
|
|
kernel[0] = 1.0;
|
|
return;
|
|
}
|
|
|
|
// fill in kernel
|
|
sum = 0.0;
|
|
for (x = min; x <= max; ++x)
|
|
{
|
|
sum += kernel[x-min] =
|
|
exp(- ((double)(x*x)) / (std * std * 2.0));
|
|
}
|
|
|
|
// normalize
|
|
for (x = min; x <= max; ++x)
|
|
{
|
|
kernel[x-min] /= sum;
|
|
}
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
int vtkImageGaussianSmooth::RequestUpdateExtent (
|
|
vtkInformation * vtkNotUsed(request),
|
|
vtkInformationVector **inputVector,
|
|
vtkInformationVector *outputVector)
|
|
{
|
|
// get the info objects
|
|
vtkInformation *outInfo = outputVector->GetInformationObject(0);
|
|
vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
|
|
|
|
int wholeExtent[6], inExt[6];
|
|
|
|
outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), inExt);
|
|
|
|
// Expand filtered axes
|
|
inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), wholeExtent);
|
|
|
|
this->InternalRequestUpdateExtent(inExt, wholeExtent);
|
|
|
|
inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), inExt, 6);
|
|
|
|
return 1;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
void vtkImageGaussianSmooth::InternalRequestUpdateExtent(int *inExt,
|
|
int *wholeExtent)
|
|
{
|
|
int idx, radius;
|
|
|
|
// Expand filtered axes
|
|
for (idx = 0; idx < this->Dimensionality; ++idx)
|
|
{
|
|
radius = (int)(this->StandardDeviations[idx] * this->RadiusFactors[idx]);
|
|
inExt[idx*2] -= radius;
|
|
if (inExt[idx*2] < wholeExtent[idx*2])
|
|
{
|
|
inExt[idx*2] = wholeExtent[idx*2];
|
|
}
|
|
|
|
inExt[idx*2+1] += radius;
|
|
if (inExt[idx*2+1] > wholeExtent[idx*2+1])
|
|
{
|
|
inExt[idx*2+1] = wholeExtent[idx*2+1];
|
|
}
|
|
}
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// For a given position along the convolution axis, this method loops over
|
|
// all other axes, and performs the convolution. Boundary conditions handled
|
|
// previously.
|
|
template <class T>
|
|
void
|
|
vtkImageGaussianSmoothExecute(vtkImageGaussianSmooth *self, int axis,
|
|
double *kernel, int kernelSize,
|
|
vtkImageData *inData, T *inPtrC,
|
|
vtkImageData *outData, int outExt[6],
|
|
T *outPtrC, int *pcycle, int target,
|
|
int *pcount, int total)
|
|
{
|
|
int maxC, max0, max1;
|
|
int idxC, idx0, idx1, idxK;
|
|
vtkIdType *inIncs, *outIncs;
|
|
vtkIdType inInc0, inInc1, inIncK, outInc0, outInc1;
|
|
T *outPtr1, *outPtr0;
|
|
T *inPtr1, *inPtr0, *inPtrK;
|
|
double *ptrK, sum;
|
|
|
|
// avoid warnings
|
|
self = self;
|
|
max0 = max1 = inInc0 = inInc1 = outInc0 = outInc1 = 0;
|
|
|
|
// I am counting on the fact that tight loops (component on outside)
|
|
// is more important than cache misses from shuffled access.
|
|
|
|
// Do the correct shuffling of the axes (increments, extents)
|
|
inIncs = inData->GetIncrements();
|
|
outIncs = outData->GetIncrements();
|
|
inIncK = inIncs[axis];
|
|
maxC = outData->GetNumberOfScalarComponents();
|
|
switch (axis)
|
|
{
|
|
case 0:
|
|
inInc0 = inIncs[1]; inInc1 = inIncs[2];
|
|
outInc0 = outIncs[1]; outInc1 = outIncs[2];
|
|
max0 = outExt[3] - outExt[2] + 1; max1 = outExt[5] - outExt[4] + 1;
|
|
break;
|
|
case 1:
|
|
inInc0 = inIncs[0]; inInc1 = inIncs[2];
|
|
outInc0 = outIncs[0]; outInc1 = outIncs[2];
|
|
max0 = outExt[1] - outExt[0] + 1; max1 = outExt[5] - outExt[4] + 1;
|
|
break;
|
|
case 2:
|
|
inInc0 = inIncs[0]; inInc1 = inIncs[1];
|
|
outInc0 = outIncs[0]; outInc1 = outIncs[1];
|
|
max0 = outExt[1] - outExt[0] + 1; max1 = outExt[3] - outExt[2] + 1;
|
|
break;
|
|
}
|
|
|
|
for (idxC = 0; idxC < maxC; ++idxC)
|
|
{
|
|
inPtr1 = inPtrC;
|
|
outPtr1 = outPtrC;
|
|
for (idx1 = 0; !self->AbortExecute && idx1 < max1; ++idx1)
|
|
{
|
|
inPtr0 = inPtr1;
|
|
outPtr0 = outPtr1;
|
|
for (idx0 = 0; idx0 < max0; ++idx0)
|
|
{
|
|
inPtrK = inPtr0;
|
|
ptrK = kernel;
|
|
sum = 0.0;
|
|
// too bad this short loop has to be the inner most loop
|
|
for (idxK = 0; idxK < kernelSize; ++idxK)
|
|
{
|
|
sum += *ptrK * (double)(*inPtrK);
|
|
++ptrK;
|
|
inPtrK += inIncK;
|
|
}
|
|
*outPtr0 = (T)(sum);
|
|
inPtr0 += inInc0;
|
|
outPtr0 += outInc0;
|
|
}
|
|
inPtr1 += inInc1;
|
|
outPtr1 += outInc1;
|
|
// we finished a row ... do we update ???
|
|
if (total)
|
|
{ // yes this is the main thread
|
|
*pcycle += max0;
|
|
if (*pcycle > target)
|
|
{ // yes
|
|
*pcycle -= target;
|
|
*pcount += target;
|
|
self->UpdateProgress((double)(*pcount) / (double)total);
|
|
//fprintf(stderr, "count: %d, total: %d, progress: %f\n",
|
|
//*pcount, total, (double)(*pcount) / (double)total);
|
|
}
|
|
}
|
|
}
|
|
|
|
++inPtrC;
|
|
++outPtrC;
|
|
}
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
template <class T>
|
|
size_t vtkImageGaussianSmoothGetTypeSize(T*)
|
|
{
|
|
return sizeof(T);
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// This method convolves over one axis. It loops over the convolved axis,
|
|
// and handles boundary conditions.
|
|
void vtkImageGaussianSmooth::ExecuteAxis(int axis,
|
|
vtkImageData *inData, int inExt[6],
|
|
vtkImageData *outData, int outExt[6],
|
|
int *pcycle, int target,
|
|
int *pcount, int total,
|
|
vtkInformation *inInfo)
|
|
{
|
|
int idxA, max;
|
|
int wholeExtent[6], wholeMax, wholeMin;
|
|
double *kernel;
|
|
// previousClip and currentClip rembers that the previous was not clipped
|
|
// keeps from recomputing kernels for center pixels.
|
|
int kernelSize = 0;
|
|
int kernelLeftClip, kernelRightClip;
|
|
int previousClipped, currentClipped;
|
|
int radius, size;
|
|
void *inPtr;
|
|
void *outPtr;
|
|
int coords[3];
|
|
vtkIdType *outIncs, outIncA;
|
|
|
|
// Get the correct starting pointer of the output
|
|
outPtr = outData->GetScalarPointerForExtent(outExt);
|
|
outIncs = outData->GetIncrements();
|
|
outIncA = outIncs[axis];
|
|
|
|
// trick to account for the scalar type of the output(used to be only float)
|
|
switch (outData->GetScalarType())
|
|
{
|
|
vtkTemplateMacro(
|
|
outIncA *= vtkImageGaussianSmoothGetTypeSize(static_cast<VTK_TT*>(0))
|
|
);
|
|
default:
|
|
vtkErrorMacro("Unknown scalar type");
|
|
return;
|
|
}
|
|
|
|
// Determine default starting position of input
|
|
coords[0] = inExt[0];
|
|
coords[1] = inExt[2];
|
|
coords[2] = inExt[4];
|
|
|
|
// get whole extent for boundary checking ...
|
|
inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), wholeExtent);
|
|
wholeMin = wholeExtent[axis*2];
|
|
wholeMax = wholeExtent[axis*2+1];
|
|
|
|
// allocate memory for the kernel
|
|
radius = (int)(this->StandardDeviations[axis] * this->RadiusFactors[axis]);
|
|
size = 2*radius + 1;
|
|
kernel = new double[size];
|
|
|
|
// loop over the convolution axis
|
|
previousClipped = currentClipped = 1;
|
|
max = outExt[axis*2+1];
|
|
for (idxA = outExt[axis*2]; idxA <= max; ++idxA)
|
|
{
|
|
// left boundary condition
|
|
coords[axis] = idxA - radius;
|
|
kernelLeftClip = wholeMin - coords[axis];
|
|
if (kernelLeftClip > 0)
|
|
{ // front of kernel is cut off ("kernelStart" samples)
|
|
coords[axis] += kernelLeftClip;
|
|
}
|
|
else
|
|
{
|
|
kernelLeftClip = 0;
|
|
}
|
|
// Right boundary condition
|
|
kernelRightClip = (idxA + radius) - wholeMax;
|
|
if (kernelRightClip < 0)
|
|
{
|
|
kernelRightClip = 0;
|
|
}
|
|
|
|
// We can only use previous kernel if it is not clipped and new
|
|
// kernel is also not clipped.
|
|
currentClipped = kernelLeftClip + kernelRightClip;
|
|
if (currentClipped || previousClipped)
|
|
{
|
|
this->ComputeKernel(kernel, -radius+kernelLeftClip,
|
|
radius-kernelRightClip,
|
|
(double)(this->StandardDeviations[axis]));
|
|
kernelSize = size - kernelLeftClip - kernelRightClip;
|
|
}
|
|
previousClipped = currentClipped;
|
|
|
|
/* now do the convolution on the rest of the axes */
|
|
inPtr = inData->GetScalarPointer(coords);
|
|
switch (inData->GetScalarType())
|
|
{
|
|
vtkTemplateMacro(
|
|
vtkImageGaussianSmoothExecute(this, axis, kernel, kernelSize,
|
|
inData, (VTK_TT*)(inPtr),
|
|
outData, outExt, (VTK_TT*)(outPtr),
|
|
pcycle, target, pcount, total)
|
|
);
|
|
default:
|
|
vtkErrorMacro("Unknown scalar type");
|
|
return;
|
|
}
|
|
outPtr = (void *)((unsigned char *)outPtr + outIncA);
|
|
}
|
|
|
|
// get rid of temporary kernel
|
|
delete [] kernel;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
// This method decomposes the gaussian and smooths along each axis.
|
|
void vtkImageGaussianSmooth::ThreadedRequestData(
|
|
vtkInformation *vtkNotUsed(request),
|
|
vtkInformationVector **inputVector,
|
|
vtkInformationVector *outputVector,
|
|
vtkImageData ***inData,
|
|
vtkImageData **outData,
|
|
int outExt[6], int id)
|
|
{
|
|
int inExt[6];
|
|
int target, count, total, cycle;
|
|
|
|
// for feed back, determine line target to get 50 progress update
|
|
// update is called every target lines. Progress is computed from
|
|
// the number of pixels processed so far.
|
|
count = 0; target = 0; total = 0; cycle = 0;
|
|
if (id == 0)
|
|
{
|
|
// determine the number of pixels.
|
|
total = this->Dimensionality * (outExt[1] - outExt[0] + 1)
|
|
* (outExt[3] - outExt[2] + 1) * (outExt[5] - outExt[4] + 1)
|
|
* inData[0][0]->GetNumberOfScalarComponents();
|
|
// pixels per update (50 updates)
|
|
target = total / 50;
|
|
}
|
|
|
|
// this filter expects that input is the same type as output.
|
|
if (inData[0][0]->GetScalarType() != outData[0]->GetScalarType())
|
|
{
|
|
vtkErrorMacro("Execute: input ScalarType, "
|
|
<< inData[0][0]->GetScalarType()
|
|
<< ", must match out ScalarType "
|
|
<< outData[0]->GetScalarType());
|
|
return;
|
|
}
|
|
|
|
// Decompose
|
|
vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
|
|
vtkInformation *outInfo = outputVector->GetInformationObject(0);
|
|
int wholeExt[6];
|
|
inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), wholeExt);
|
|
outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), inExt);
|
|
this->InternalRequestUpdateExtent(inExt, wholeExt);
|
|
|
|
switch (this->Dimensionality)
|
|
{
|
|
case 1:
|
|
this->ExecuteAxis(0, inData[0][0], inExt, outData[0], outExt,
|
|
&cycle, target, &count, total, inInfo);
|
|
break;
|
|
case 2:
|
|
int tempExt[6];
|
|
vtkImageData *tempData;
|
|
// compute intermediate extent
|
|
tempExt[0] = inExt[0]; tempExt[1] = inExt[1];
|
|
tempExt[2] = outExt[2]; tempExt[3] = outExt[3];
|
|
tempExt[4] = inExt[4]; tempExt[5] = inExt[5];
|
|
// create a temp data for intermediate results
|
|
tempData = vtkImageData::New();
|
|
tempData->SetExtent(tempExt);
|
|
tempData->SetNumberOfScalarComponents(
|
|
inData[0][0]->GetNumberOfScalarComponents());
|
|
tempData->SetScalarType(inData[0][0]->GetScalarType());
|
|
this->ExecuteAxis(1, inData[0][0], inExt, tempData, tempExt,
|
|
&cycle, target, &count, total, inInfo);
|
|
this->ExecuteAxis(0, tempData, tempExt, outData[0], outExt,
|
|
&cycle, target, &count, total, inInfo);
|
|
// release temporary data
|
|
tempData->Delete();
|
|
break;
|
|
case 3:
|
|
// we do z first because it is most likely smallest
|
|
int temp0Ext[6], temp1Ext[6];
|
|
vtkImageData *temp0Data, *temp1Data;
|
|
// compute intermediate extents
|
|
temp0Ext[0] = inExt[0]; temp0Ext[1] = inExt[1];
|
|
temp0Ext[2] = inExt[2]; temp0Ext[3] = inExt[3];
|
|
temp0Ext[4] = outExt[4]; temp0Ext[5] = outExt[5];
|
|
|
|
temp1Ext[0] = inExt[0]; temp1Ext[1] = inExt[1];
|
|
temp1Ext[2] = outExt[2]; temp1Ext[3] = outExt[3];
|
|
temp1Ext[4] = outExt[4]; temp1Ext[5] = outExt[5];
|
|
|
|
// create a temp data for intermediate results
|
|
temp0Data = vtkImageData::New();
|
|
temp0Data->SetExtent(temp0Ext);
|
|
temp0Data->SetNumberOfScalarComponents(
|
|
inData[0][0]->GetNumberOfScalarComponents());
|
|
temp0Data->SetScalarType(inData[0][0]->GetScalarType());
|
|
|
|
temp1Data = vtkImageData::New();
|
|
temp1Data->SetExtent(temp1Ext);
|
|
temp1Data->SetNumberOfScalarComponents(
|
|
inData[0][0]->GetNumberOfScalarComponents());
|
|
temp1Data->SetScalarType(inData[0][0]->GetScalarType());
|
|
this->ExecuteAxis(2, inData[0][0], inExt, temp0Data, temp0Ext,
|
|
&cycle, target, &count, total, inInfo);
|
|
this->ExecuteAxis(1, temp0Data, temp0Ext, temp1Data, temp1Ext,
|
|
&cycle, target, &count, total, inInfo);
|
|
temp0Data->Delete();
|
|
this->ExecuteAxis(0, temp1Data, temp1Ext, outData[0], outExt,
|
|
&cycle, target, &count, total, inInfo);
|
|
temp1Data->Delete();
|
|
break;
|
|
}
|
|
}
|
|
|