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.
473 lines
15 KiB
473 lines
15 KiB
/*=========================================================================
|
|
|
|
Program: Visualization Toolkit
|
|
Module: $RCSfile: vtkImageDifference.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 "vtkImageDifference.h"
|
|
|
|
#include "vtkImageData.h"
|
|
#include "vtkInformation.h"
|
|
#include "vtkInformationVector.h"
|
|
#include "vtkObjectFactory.h"
|
|
#include "vtkStreamingDemandDrivenPipeline.h"
|
|
|
|
vtkCxxRevisionMacro(vtkImageDifference, "$Revision: 1.39 $");
|
|
vtkStandardNewMacro(vtkImageDifference);
|
|
|
|
// Construct object to extract all of the input data.
|
|
vtkImageDifference::vtkImageDifference()
|
|
{
|
|
int i;
|
|
for ( i = 0; i < this->NumberOfThreads; i++ )
|
|
{
|
|
this->ErrorPerThread[i] = 0;
|
|
this->ThresholdedErrorPerThread[i] = 0.0;
|
|
}
|
|
this->Threshold = 16;
|
|
this->AllowShift = 1;
|
|
this->Averaging = 1;
|
|
this->SetNumberOfInputPorts(2);
|
|
}
|
|
|
|
|
|
|
|
// not so simple macro for calculating error
|
|
#define vtkImageDifferenceComputeError(c1,c2) \
|
|
/* compute the pixel to pixel difference first */ \
|
|
r1 = abs(((int)(c1)[0] - (int)(c2)[0])); \
|
|
g1 = abs(((int)(c1)[1] - (int)(c2)[1])); \
|
|
b1 = abs(((int)(c1)[2] - (int)(c2)[2]));\
|
|
if ((r1+g1+b1) < (tr+tg+tb)) { tr = r1; tg = g1; tb = b1; } \
|
|
/* if averaging is on and we have neighbor info then compute */ \
|
|
/* input1 to avg(input2) */ \
|
|
if (this->Averaging && \
|
|
(idx0 > inMinX + 1) && (idx0 < inMaxX - 1) && \
|
|
(idx1 > inMinY + 1) && (idx1 < inMaxY - 1)) \
|
|
{\
|
|
ar1 = (int)(c1)[0]; \
|
|
ag1 = (int)(c1)[1]; \
|
|
ab1 = (int)(c1)[2]; \
|
|
ar2 = (int)(c2)[0] + (int)(c2 - in2Inc0)[0] + (int)(c2 + in2Inc0)[0] + \
|
|
(int)(c2-in2Inc1)[0] + (int)(c2-in2Inc1-in2Inc0)[0] + (int)(c2-in2Inc1+in2Inc0)[0] + \
|
|
(int)(c2+in2Inc1)[0] + (int)(c2+in2Inc1-in2Inc0)[0] + (int)(c2+in2Inc1+in2Inc0)[0]; \
|
|
ag2 = (int)(c2)[1] + (int)(c2 - in2Inc0)[1] + (int)(c2 + in2Inc0)[1] + \
|
|
(int)(c2-in2Inc1)[1] + (int)(c2-in2Inc1-in2Inc0)[1] + (int)(c2-in2Inc1+in2Inc0)[1] + \
|
|
(int)(c2+in2Inc1)[1] + (int)(c2+in2Inc1-in2Inc0)[1] + (int)(c2+in2Inc1+in2Inc0)[1]; \
|
|
ab2 = (int)(c2)[2] + (int)(c2 - in2Inc0)[2] + (int)(c2 + in2Inc0)[2] + \
|
|
(int)(c2-in2Inc1)[2] + (int)(c2-in2Inc1-in2Inc0)[2] + (int)(c2-in2Inc1+in2Inc0)[2] + \
|
|
(int)(c2+in2Inc1)[2] + (int)(c2+in2Inc1-in2Inc0)[2] + (int)(c2+in2Inc1+in2Inc0)[2]; \
|
|
r1 = abs(ar1 - ar2/9); \
|
|
g1 = abs(ag1 - ag2/9); \
|
|
b1 = abs(ab1 - ab2/9); \
|
|
if ((r1+g1+b1) < (tr+tg+tb)) { tr = r1; tg = g1; tb = b1; } \
|
|
/* Now compute the avg(input1) to avg(input2) comparison */ \
|
|
ar1 = (int)(c1)[0] + (int)(c1 - in1Inc0)[0] + (int)(c1 + in1Inc0)[0] + \
|
|
(int)(c1-in1Inc1)[0] + (int)(c1-in1Inc1-in1Inc0)[0] + (int)(c1-in1Inc1+in1Inc0)[0] + \
|
|
(int)(c1+in1Inc1)[0] + (int)(c1+in1Inc1-in1Inc0)[0] + (int)(c1+in1Inc1+in1Inc0)[0]; \
|
|
ag1 = (int)(int)(c1)[1] + (int)(c1 - in1Inc0)[1] + (int)(c1 + in1Inc0)[1] + \
|
|
(int)(c1-in1Inc1)[1] + (int)(c1-in1Inc1-in1Inc0)[1] + (int)(c1-in1Inc1+in1Inc0)[1] + \
|
|
(int)(c1+in1Inc1)[1] + (int)(c1+in1Inc1-in1Inc0)[1] + (int)(c1+in1Inc1+in1Inc0)[1]; \
|
|
ab1 = (int)(int)(c1)[2] + (int)(c1 - in1Inc0)[2] + (int)(c1 + in1Inc0)[2] + \
|
|
(int)(c1-in1Inc1)[2] + (int)(c1-in1Inc1-in1Inc0)[2] + (int)(c1-in1Inc1+in1Inc0)[2] + \
|
|
(int)(c1+in1Inc1)[2] + (int)(c1+in1Inc1-in1Inc0)[2] + (int)(c1+in1Inc1+in1Inc0)[2]; \
|
|
r1 = abs(ar1/9 - ar2/9); \
|
|
g1 = abs(ag1/9 - ag2/9); \
|
|
b1 = abs(ab1/9 - ab2/9); \
|
|
if ((r1+g1+b1) < (tr+tg+tb)) { tr = r1; tg = g1; tb = b1; } \
|
|
/* finally compute avg(input1) to input2) */ \
|
|
ar2 = (int)(c2)[0]; \
|
|
ag2 = (int)(c2)[1]; \
|
|
ab2 = (int)(c2)[2]; \
|
|
r1 = abs(ar1/9 - ar2); \
|
|
g1 = abs(ag1/9 - ag2); \
|
|
b1 = abs(ab1/9 - ab2); \
|
|
if ((r1+g1+b1) < (tr+tg+tb)) { tr = r1; tg = g1; tb = b1; } \
|
|
}
|
|
|
|
|
|
|
|
|
|
//----------------------------------------------------------------------------
|
|
// This method computes the input extent necessary to generate the output.
|
|
int vtkImageDifference::RequestUpdateExtent(
|
|
vtkInformation * vtkNotUsed(request),
|
|
vtkInformationVector **inputVector,
|
|
vtkInformationVector *outputVector)
|
|
{
|
|
int idx;
|
|
|
|
// get the info objects
|
|
vtkInformation* outInfo = outputVector->GetInformationObject(0);
|
|
vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
|
|
|
|
int *wholeExtent =
|
|
inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT());
|
|
|
|
int uExt[6];
|
|
outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(),uExt);
|
|
|
|
// grow input whole extent.
|
|
for (idx = 0; idx < 2; ++idx)
|
|
{
|
|
uExt[idx*2] -= 2;
|
|
uExt[idx*2+1] += 2;
|
|
|
|
// we must clip extent with whole extent is we hanlde boundaries.
|
|
if (uExt[idx*2] < wholeExtent[idx*2])
|
|
{
|
|
uExt[idx*2] = wholeExtent[idx*2];
|
|
}
|
|
if (uExt[idx*2 + 1] > wholeExtent[idx*2 + 1])
|
|
{
|
|
uExt[idx*2 + 1] = wholeExtent[idx*2 + 1];
|
|
}
|
|
}
|
|
inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(),uExt,6);
|
|
|
|
// now do the second input
|
|
inInfo = inputVector[1]->GetInformationObject(0);
|
|
wholeExtent =
|
|
inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT());
|
|
outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(),uExt);
|
|
|
|
// grow input whole extent.
|
|
for (idx = 0; idx < 2; ++idx)
|
|
{
|
|
uExt[idx*2] -= 2;
|
|
uExt[idx*2+1] += 2;
|
|
|
|
// we must clip extent with whole extent is we hanlde boundaries.
|
|
if (uExt[idx*2] < wholeExtent[idx*2])
|
|
{
|
|
uExt[idx*2] = wholeExtent[idx*2];
|
|
}
|
|
if (uExt[idx*2 + 1] > wholeExtent[idx*2 + 1])
|
|
{
|
|
uExt[idx*2 + 1] = wholeExtent[idx*2 + 1];
|
|
}
|
|
}
|
|
inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(),uExt,6);
|
|
|
|
return 1;
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
void vtkImageDifference::ThreadedRequestData(
|
|
vtkInformation * vtkNotUsed( request ),
|
|
vtkInformationVector ** vtkNotUsed( inputVector ),
|
|
vtkInformationVector * vtkNotUsed( outputVector ),
|
|
vtkImageData ***inData,
|
|
vtkImageData **outData,
|
|
int outExt[6], int id)
|
|
{
|
|
unsigned char *in1Ptr0, *in1Ptr1, *in1Ptr2;
|
|
unsigned char *in2Ptr0, *in2Ptr1, *in2Ptr2;
|
|
unsigned char *outPtr0, *outPtr1, *outPtr2;
|
|
int min0, max0, min1, max1, min2, max2;
|
|
int idx0, idx1, idx2;
|
|
vtkIdType in1Inc0, in1Inc1, in1Inc2;
|
|
vtkIdType in2Inc0, in2Inc1, in2Inc2;
|
|
vtkIdType outInc0, outInc1, outInc2;
|
|
int tr, tg, tb, r1, g1, b1;
|
|
int ar1, ag1, ab1, ar2, ag2, ab2;
|
|
int inMinX, inMaxX, inMinY, inMaxY;
|
|
int *inExt;
|
|
int matched;
|
|
unsigned long count = 0;
|
|
unsigned long target;
|
|
|
|
this->ErrorPerThread[id] = 0;
|
|
this->ThresholdedErrorPerThread[id] = 0;
|
|
|
|
if (inData[0] == NULL || inData[1] == NULL || outData == NULL)
|
|
{
|
|
if (!id)
|
|
{
|
|
vtkErrorMacro(<< "Execute: Missing data");
|
|
}
|
|
this->ErrorPerThread[id] = 1000;
|
|
this->ThresholdedErrorPerThread[id] = 1000;
|
|
return;
|
|
}
|
|
|
|
if (inData[0][0]->GetNumberOfScalarComponents() != 3 ||
|
|
inData[1][0]->GetNumberOfScalarComponents() != 3 ||
|
|
outData[0]->GetNumberOfScalarComponents() != 3)
|
|
{
|
|
if (!id)
|
|
{
|
|
vtkErrorMacro(<< "Execute: Expecting 3 components (RGB)");
|
|
}
|
|
this->ErrorPerThread[id] = 1000;
|
|
this->ThresholdedErrorPerThread[id] = 1000;
|
|
return;
|
|
}
|
|
|
|
// this filter expects that input is the same type as output.
|
|
if (inData[0][0]->GetScalarType() != VTK_UNSIGNED_CHAR ||
|
|
inData[1][0]->GetScalarType() != VTK_UNSIGNED_CHAR ||
|
|
outData[0]->GetScalarType() != VTK_UNSIGNED_CHAR)
|
|
{
|
|
if (!id)
|
|
{
|
|
vtkErrorMacro(<< "Execute: All ScalarTypes must be unsigned char");
|
|
}
|
|
this->ErrorPerThread[id] = 1000;
|
|
this->ThresholdedErrorPerThread[id] = 1000;
|
|
return;
|
|
}
|
|
|
|
in1Ptr2 = (unsigned char *) inData[0][0]->GetScalarPointerForExtent(outExt);
|
|
in2Ptr2 = (unsigned char *) inData[1][0]->GetScalarPointerForExtent(outExt);
|
|
outPtr2 = (unsigned char *) outData[0]->GetScalarPointerForExtent(outExt);
|
|
|
|
inData[0][0]->GetIncrements(in1Inc0, in1Inc1, in1Inc2);
|
|
inData[1][0]->GetIncrements(in2Inc0, in2Inc1, in2Inc2);
|
|
outData[0]->GetIncrements(outInc0, outInc1, outInc2);
|
|
|
|
min0 = outExt[0]; max0 = outExt[1];
|
|
min1 = outExt[2]; max1 = outExt[3];
|
|
min2 = outExt[4]; max2 = outExt[5];
|
|
|
|
inExt = inData[0][0]->GetExtent();
|
|
// we set min and Max to be one pixel in from actual values to support
|
|
// the 3x3 averaging we do
|
|
inMinX = inExt[0]; inMaxX = inExt[1];
|
|
inMinY = inExt[2]; inMaxY = inExt[3];
|
|
|
|
target = (unsigned long)((max2 - min2 +1)*(max1 - min1 + 1)/50.0);
|
|
target++;
|
|
|
|
for (idx2 = min2; idx2 <= max2; ++idx2)
|
|
{
|
|
in1Ptr1 = in1Ptr2;
|
|
in2Ptr1 = in2Ptr2;
|
|
outPtr1 = outPtr2;
|
|
for (idx1 = min1; !this->AbortExecute && idx1 <= max1; ++idx1)
|
|
{
|
|
if (!id)
|
|
{
|
|
if (!(count%target))
|
|
{
|
|
this->UpdateProgress(count/(50.0*target));
|
|
}
|
|
count++;
|
|
}
|
|
in1Ptr0 = in1Ptr1;
|
|
in2Ptr0 = in2Ptr1;
|
|
outPtr0 = outPtr1;
|
|
for (idx0 = min0; idx0 <= max0; ++idx0)
|
|
{
|
|
tr = 1000;
|
|
tg = 1000;
|
|
tb = 1000;
|
|
|
|
/* check the exact match pixel */
|
|
vtkImageDifferenceComputeError(in1Ptr0,in2Ptr0);
|
|
|
|
// do a quick check to see if this match is exact, if so
|
|
// we can save some seious time by skipping the eight
|
|
// connected neighbors
|
|
matched = 0;
|
|
if ((tr <= 0)&&(tg <= 0)&&(tb <= 0))
|
|
{
|
|
matched = 1;
|
|
}
|
|
|
|
/* If AllowShift, then we examine neighboring pixels to
|
|
find the least difference. This feature is used to
|
|
allow images to shift slightly between different graphics
|
|
systems, like between opengl and starbase. */
|
|
if (!matched && this->AllowShift)
|
|
{
|
|
/* lower row */
|
|
if (idx1 > inMinY)
|
|
{
|
|
vtkImageDifferenceComputeError(in1Ptr0-in1Inc1,in2Ptr0);
|
|
if (idx0 > inMinX)
|
|
{
|
|
vtkImageDifferenceComputeError(in1Ptr0-in1Inc0-in1Inc1,in2Ptr0);
|
|
}
|
|
if (idx0 < inMaxX)
|
|
{
|
|
vtkImageDifferenceComputeError(in1Ptr0+in1Inc0-in1Inc1,in2Ptr0);
|
|
}
|
|
}
|
|
/* middle row (center already considered) */
|
|
if (idx0 > inMinX)
|
|
{
|
|
vtkImageDifferenceComputeError(in1Ptr0-in1Inc0,in2Ptr0);
|
|
}
|
|
if (idx0 < inMaxX)
|
|
{
|
|
vtkImageDifferenceComputeError(in1Ptr0+in1Inc0,in2Ptr0);
|
|
}
|
|
/* upper row */
|
|
if (idx1 < inMaxY)
|
|
{
|
|
vtkImageDifferenceComputeError(in1Ptr0+in1Inc1,in2Ptr0);
|
|
if (idx0 > inMinX)
|
|
{
|
|
vtkImageDifferenceComputeError(in1Ptr0-in1Inc0+in1Inc1,in2Ptr0);
|
|
}
|
|
if (idx0 < inMaxX)
|
|
{
|
|
vtkImageDifferenceComputeError(in1Ptr0+in1Inc0+in1Inc1,in2Ptr0);
|
|
}
|
|
}
|
|
}
|
|
|
|
this->ErrorPerThread[id] = this->ErrorPerThread[id] + (tr + tg + tb)/(3.0*255);
|
|
tr -= this->Threshold;
|
|
if (tr < 0)
|
|
{
|
|
tr = 0;
|
|
}
|
|
tg -= this->Threshold;
|
|
if (tg < 0)
|
|
{
|
|
tg = 0;
|
|
}
|
|
tb -= this->Threshold;
|
|
if (tb < 0)
|
|
{
|
|
tb = 0;
|
|
}
|
|
*outPtr0++ = (unsigned char)tr;
|
|
*outPtr0++ = (unsigned char)tg;
|
|
*outPtr0++ = (unsigned char)tb;
|
|
this->ThresholdedErrorPerThread[id] += (tr + tg + tb)/(3.0*255.0);
|
|
|
|
in1Ptr0 += in1Inc0;
|
|
in2Ptr0 += in2Inc0;
|
|
}
|
|
outPtr1 += outInc1;
|
|
in1Ptr1 += in1Inc1;
|
|
in2Ptr1 += in2Inc1;
|
|
}
|
|
outPtr2 += outInc2;
|
|
in1Ptr2 += in1Inc2;
|
|
in2Ptr2 += in2Inc2;
|
|
}
|
|
}
|
|
|
|
//----------------------------------------------------------------------------
|
|
//Make the output the intersection of the inputs, of course the inputs better
|
|
//be the same size
|
|
int vtkImageDifference::RequestInformation (
|
|
vtkInformation * vtkNotUsed(request),
|
|
vtkInformationVector **inputVector,
|
|
vtkInformationVector *outputVector)
|
|
{
|
|
// get the info objects
|
|
vtkInformation* outInfo = outputVector->GetInformationObject(0);
|
|
vtkInformation *inInfo1 = inputVector[0]->GetInformationObject(0);
|
|
vtkInformation *inInfo2 = inputVector[1]->GetInformationObject(0);
|
|
|
|
int *in1Ext = inInfo1->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT());
|
|
int *in2Ext = inInfo2->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT());
|
|
|
|
int i;
|
|
if (in1Ext[0] != in2Ext[0] || in1Ext[1] != in2Ext[1] ||
|
|
in1Ext[2] != in2Ext[2] || in1Ext[3] != in2Ext[3] ||
|
|
in1Ext[4] != in2Ext[4] || in1Ext[5] != in2Ext[5])
|
|
{
|
|
for (i = 0; i < this->NumberOfThreads; i++)
|
|
{
|
|
this->ErrorPerThread[i] = 1000;
|
|
this->ThresholdedErrorPerThread[i] = 1000;
|
|
}
|
|
vtkErrorMacro("ExecuteInformation: Input are not the same size.\n"
|
|
<< " Input1 is: " << in1Ext[0] << "," << in1Ext[1] << ","
|
|
<< in1Ext[2] << "," << in1Ext[3] << ","
|
|
<< in1Ext[4] << "," << in1Ext[5] << "\n"
|
|
<< " Input2 is: " << in2Ext[0] << "," << in2Ext[1] << ","
|
|
<< in2Ext[2] << "," << in2Ext[3] << ","
|
|
<< in2Ext[4] << "," << in2Ext[5] );
|
|
}
|
|
|
|
// We still need to set the whole extent to be the intersection.
|
|
// Otherwise the execute may crash.
|
|
int ext[6];
|
|
for (i = 0; i < 3; ++i)
|
|
{
|
|
ext[i*2] = in1Ext[i*2];
|
|
if (ext[i*2] < in2Ext[i*2])
|
|
{
|
|
ext[i*2] = in2Ext[i*2];
|
|
}
|
|
ext[i*2+1] = in1Ext[i*2+1];
|
|
if (ext[i*2+1] > in2Ext[i*2+1])
|
|
{
|
|
ext[i*2+1] = in2Ext[i*2+1];
|
|
}
|
|
}
|
|
outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(),ext,6);
|
|
|
|
return 1;
|
|
}
|
|
|
|
double vtkImageDifference::GetError()
|
|
{
|
|
double error = 0.0;
|
|
int i;
|
|
|
|
for ( i= 0; i < this->NumberOfThreads; i++ )
|
|
{
|
|
error += this->ErrorPerThread[i];
|
|
}
|
|
|
|
return error;
|
|
}
|
|
|
|
double vtkImageDifference::GetThresholdedError()
|
|
{
|
|
double error = 0.0;
|
|
int i;
|
|
|
|
for ( i= 0; i < this->NumberOfThreads; i++ )
|
|
{
|
|
error += this->ThresholdedErrorPerThread[i];
|
|
}
|
|
|
|
return error;
|
|
}
|
|
|
|
vtkImageData *vtkImageDifference::GetImage()
|
|
{
|
|
if (this->GetNumberOfInputConnections(1) < 1)
|
|
{
|
|
return 0;
|
|
}
|
|
return vtkImageData::SafeDownCast(
|
|
this->GetExecutive()->GetInputData(1, 0));
|
|
}
|
|
|
|
|
|
void vtkImageDifference::PrintSelf(ostream& os, vtkIndent indent)
|
|
{
|
|
this->Superclass::PrintSelf(os,indent);
|
|
|
|
int i;
|
|
|
|
for ( i= 0; i < this->NumberOfThreads; i++ )
|
|
{
|
|
os << indent << "Error for thread " << i << ": " << this->ErrorPerThread[i] << "\n";
|
|
os << indent << "ThresholdedError for thread " << i << ": "
|
|
<< this->ThresholdedErrorPerThread[i] << "\n";
|
|
}
|
|
os << indent << "Threshold: " << this->Threshold << "\n";
|
|
os << indent << "AllowShift: " << this->AllowShift << "\n";
|
|
os << indent << "Averaging: " << this->Averaging << "\n";
|
|
}
|
|
|
|
|
|
|