[vtkusers] trigger the execution of a vtkThreadedImageAlgorithm derived filter by calling Update()
Gomez, Alberto
alberto.gomez at kcl.ac.uk
Thu Dec 17 05:50:12 EST 2015
Hi all,
I just fixed it. The problem was that I needed to run Update() on the
image source filters. So my previous code does execute.
There was a bug though: I was trying to calculate a number (the SSD
value) by combining the value obtained from each thread. To achieve
this, My previous code has to be changed to the code below. The main
issue was when to define the return value, how to add return values from
multiple threads and dealing with mutual exclusion.
I would appreciate if someone has a link to a user guide for
implementing these sort of filters in a more efficient way.
Regards,
Alberto
------------------------------------------- CORRECTED CODE
----------------------------------------------
======================================= vtkMeanSSD.h
=======================================
#ifndef vtkMeanSSD_h
#define vtkMeanSSD_h
#include "vtkImagingGeneralModule.h" // For export macro
#include "vtkThreadedImageAlgorithm.h"
#include <vtkMutexLock.h>
class VTKIMAGINGGENERAL_EXPORT vtkMeanSSD : public vtkThreadedImageAlgorithm
{
public:
static vtkMeanSSD *New();
vtkTypeMacro(vtkMeanSSD,vtkThreadedImageAlgorithm);
// Description:
// Set the two inputs to this filter
virtual void SetInput1Data(vtkDataObject *in) {
this->SetInputData(0,in); }
virtual void SetInput2Data(vtkDataObject *in) {
this->SetInputData(1,in); }
// Description:
// Get the metric value from the images. This is
// computed when Update() is called.
double GetValue() {return this->m_value;}
protected:
vtkMeanSSD();
~vtkMeanSSD() { this->m_outputLock->Delete();}
virtual void ThreadedRequestData(vtkInformation *request,
vtkInformationVector **inputVector,
vtkInformationVector *outputVector,
vtkImageData ***inData,
vtkImageData **outData,
int extent[6], int threadId);
//virtual int RequestData(vtkInformation*, vtkInformationVector**,
vtkInformationVector*);
double m_value;
private:
vtkMeanSSD(const vtkMeanSSD&); // Not implemented.
void operator=(const vtkMeanSSD&); // Not implemented.
vtkMutexLock* m_outputLock;
};
#endif
======================================= vtkMeanSSD.cxx
=======================================
#include "vtkMeanSSD.h"
#include "vtkImageData.h"
#include "vtkImageProgressIterator.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkObjectFactory.h"
#include "vtkStreamingDemandDrivenPipeline.h"
vtkStandardNewMacro(vtkMeanSSD);
//----------------------------------------------------------------------------
vtkMeanSSD::vtkMeanSSD()
{
this->SetNumberOfInputPorts(2);
this->m_value = 0;
this->m_outputLock = vtkMutexLock::New();
}
//----------------------------------------------------------------------------
// This templated function executes the filter for any type of data.
// Handles the two input operations
template <class T>
void vtkMeanSSDExecute(vtkMeanSSD *self,
vtkImageData *in1Data,
vtkImageData *in2Data,
double *value,
int outExt[6], int id, T *)
{
vtkImageIterator<T> inIt1(in1Data, outExt);
vtkImageIterator<T> inIt2(in2Data, outExt);
// find the region to loop over
int N = in1Data->GetNumberOfPoints();
float val = 0;
// sum (a -b)*(a-b)
/// Calculate the average of each image
while (!inIt1.IsAtEnd())
{
T* inSI1 = inIt1.BeginSpan();
T* inSI2 = inIt2.BeginSpan();
T* inSI1End = inIt1.EndSpan();
while (inSI1 != inSI1End){
T v1 = static_cast<float>(*inSI1);
T v2 = static_cast<float>(*inSI2);
val += (v1-v2)*(v1-v2);
++inSI1;
++inSI2;
}
inIt1.NextSpan();
inIt2.NextSpan();
}
val/=N;
*value = val;
}
//----------------------------------------------------------------------------
// This method is passed a input and output regions, and executes the filter
// algorithm to fill the output from the inputs.
// It just executes a switch statement to call the correct function for
// the regions data types.
void vtkMeanSSD::ThreadedRequestData(
vtkInformation * vtkNotUsed( request ),
vtkInformationVector ** vtkNotUsed( inputVector ),
vtkInformationVector * vtkNotUsed( outputVector ),
vtkImageData ***inData,
vtkImageData **outData,
int outExt[6], int id)
{
int N = this->GetNumberOfThreads();
// this filter expects that inputs that have the same number of
components
if (inData[0][0]->GetNumberOfScalarComponents() !=
inData[1][0]->GetNumberOfScalarComponents())
{
vtkErrorMacro(<< "Execute: input1 NumberOfScalarComponents, "
<< inData[0][0]->GetNumberOfScalarComponents()
<< ", must match out input2 NumberOfScalarComponents "
<< inData[1][0]->GetNumberOfScalarComponents());
return;
}
double d;
switch (inData[0][0]->GetScalarType())
{
vtkTemplateMacro(vtkMeanSSDExecute(this, inData[0][0],
inData[1][0], &d, outExt, id,
static_cast<VTK_TT *>(0)));
default:
vtkErrorMacro(<< "Execute: Unknown ScalarType");
return;
}
this->m_outputLock->Lock();
this->m_value +=d/N;
this->m_outputLock->Unlock();
}
/*
//----------------------------------------------------------------------------
int vtkMeanSSD::RequestData(
vtkInformation* request, vtkInformationVector** inputVector,
vtkInformationVector* outputVector)
{
// vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
vtkInformation* outInfo = outputVector->GetInformationObject(0);
//vtkImageData *inData =
vtkImageData::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT()));
vtkImageData *outData =
vtkImageData::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
// force reallocation of the scalars
outData->GetPointData()->SetScalars(NULL);
// Allocate the output data
outData->SetExtent(outInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT()));
outData->AllocateScalars(outInfo);
if (outData->GetScalarType() != VTK_FLOAT &&
outData->GetScalarType() != VTK_DOUBLE)
{
vtkErrorMacro(<< "Execute: output data must be be type float or
double.");
return 0;
}
// We are about to call superclass' RequestData which allocates output
// based on the update extent. However, we want the output to be the
// whole extent. So we temprarily override the update extent to be
// the whole extent.
int extentcache[6];
memcpy(extentcache,
outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT()),6*sizeof(int));
outInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(),outInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT()),
6);
// execute over the three directions
for (int i = 0; i < 3; i++)
{
this->Iteration = i;
if (ie[2*i+1] > ie[2*i])
{
if (!this->vtkThreadedImageAlgorithm::RequestData(request,
&outputVector, outputVector))
{
return 0;
}
}
}
// Restore update extent
outInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(),extentcache,6);
return 1;
}
*/
On 16/12/15 11:34, Gomez, Alberto wrote:
> Hi all,
>
> I am trying to write a multithreaded filter (vtkMeanSSD), derived from
> vtkThreadedImageAlgorithm, which will calculate the mean value of the
> difference between the intensities of two images (mean SSD).
>
> What I want to achieve is that when I call Update() on my filter, the
> calculations occur and then I can retrieve the value by calling
> another method GetValue(). This is the snippet from my main file:
>
>
> vtkSmartPointer<vtkMeanSSD> myfilter =
> vtkSmartPointer<vtkMeanSSD>::New();
> myfilter->SetInput1Data(im1);
> myfilter->SetInput2Data(im2);
> myfilter->Update();
> double val = myfilter->GetValue(); // value is 0, which is the default
> by construction
>
> The problem I have is that when I call myfilter->Update() the filter
> is not executed. For some reason, it does execute if I plug it to a
> writer, for example
>
>
> vtkSmartPointer<vtkMeanSSD> myfilter =
> vtkSmartPointer<vtkMeanSSD>::New();
> myfilter->SetInput1Data(im1);
> myfilter->SetInput2Data(im2);
>
> vtkSmartPointer<vtkNIFTIImageWriter> writer
> =vtkSmartPointer<vtkNIFTIImageWriter>::New();
> writer->SetFileName("image_SSD.nii");
> writer->SetInputConnection(myfilter->GetOutputPort());
> writer->Write(); // This triggers the execution
> double val = myfilter->GetValue(); // good value
>
> Any advice how I can trigger the ThreadedRequestData method by calling
> update? Is there any other better way to achieve the functionality I
> want?
>
> Thanks for your help. The full filter code is below.
>
> Alberto
>
>
>
>
> ======================================= vtkMeanSSD.h
> =======================================
>
>
> // .NAME vtkMeanSSD - Simple SSD image-image filter.
> // .SECTION Description
> // This calculates the mean sum of squared differences
> // 1/N*sum (a -b)^2
>
>
> #ifndef vtkMeanSSD_h
> #define vtkMeanSSD_h
>
> #include "vtkImagingGeneralModule.h" // For export macro
> #include "vtkThreadedImageAlgorithm.h"
>
> class VTKIMAGINGGENERAL_EXPORT vtkMeanSSD : public
> vtkThreadedImageAlgorithm
> {
> public:
>
>
> static vtkMeanSSD *New();
> vtkTypeMacro(vtkMeanSSD,vtkThreadedImageAlgorithm);
>
>
> // Description:
> // Set the two inputs to this filter
> virtual void SetInput1Data(vtkDataObject *in) {
> this->SetInputData(0,in); }
> virtual void SetInput2Data(vtkDataObject *in) {
> this->SetInputData(1,in); }
>
> // Description:
> // Get the metric value from the images. This is
> // computed when Update() is called.
> double GetValue() {return this->m_value;}
>
>
> protected:
>
> vtkMeanSSD();
> ~vtkMeanSSD() {}
>
> virtual void ThreadedRequestData(vtkInformation *request,
> vtkInformationVector **inputVector,
> vtkInformationVector *outputVector,
> vtkImageData ***inData,
> vtkImageData **outData,
> int extent[6], int threadId);
>
> double m_value;
>
> private:
>
> vtkMeanSSD(const vtkMeanSSD&); // Not implemented.
> void operator=(const vtkMeanSSD&); // Not implemented.
>
>
> };
>
> #endif
>
>
>
>
>
>
>
> // VTK-HeaderTest-Exclude: vtkMeanSSD.h
>
>
> ======================================= vtkMeanSSD.cxx
> =======================================
>
>
> #include "vtkMeanSSD.h"
>
> #include "vtkImageData.h"
> #include "vtkImageProgressIterator.h"
> #include "vtkInformation.h"
> #include "vtkInformationVector.h"
> #include "vtkObjectFactory.h"
> #include "vtkStreamingDemandDrivenPipeline.h"
>
> vtkStandardNewMacro(vtkMeanSSD);
>
> //----------------------------------------------------------------------------
>
> vtkMeanSSD::vtkMeanSSD()
> {
> this->SetNumberOfInputPorts(2);
> this->m_value = 0;
> }
>
>
> //----------------------------------------------------------------------------
>
> // This templated function executes the filter for any type of data.
> // Handles the two input operations
> template <class T>
> void vtkMeanSSDExecute(vtkMeanSSD *self,
> vtkImageData *in1Data,
> vtkImageData *in2Data,
> double *value,
> int outExt[6], int id, T *)
> {
> vtkImageIterator<T> inIt1(in1Data, outExt);
> vtkImageIterator<T> inIt2(in2Data, outExt);
>
> // find the region to loop over
> int N = in1Data->GetNumberOfPoints();
> float val = 0;
>
> // sum (a -b)*(a-b)
> /// Calculate the average of each image
>
> while (!inIt1.IsAtEnd())
> {
> T* inSI1 = inIt1.BeginSpan();
> T* inSI2 = inIt2.BeginSpan();
> T* inSI1End = inIt1.EndSpan();
>
> while (inSI1 != inSI1End){
>
> T v1 = static_cast<float>(*inSI1);
> T v2 = static_cast<float>(*inSI2);
> val += (v1-v2)*(v1-v2);
> ++inSI1;
> ++inSI2;
> }
> inIt1.NextSpan();
> inIt2.NextSpan();
> }
> val/=N;
>
> *value = val;
>
> }
>
>
>
> //----------------------------------------------------------------------------
>
> // This method is passed a input and output regions, and executes the
> filter
> // algorithm to fill the output from the inputs.
> // It just executes a switch statement to call the correct function for
> // the regions data types.
> void vtkMeanSSD::ThreadedRequestData(
> vtkInformation * vtkNotUsed( request ),
> vtkInformationVector ** vtkNotUsed( inputVector ),
> vtkInformationVector * vtkNotUsed( outputVector ),
> vtkImageData ***inData,
> vtkImageData **outData,
> int outExt[6], int id)
> {
> // this filter expects that input is the same type as output.
> if (inData[0][0]->GetScalarType() != outData[0]->GetScalarType())
> {
> vtkErrorMacro(<< "Execute: input1 ScalarType, "
> << inData[0][0]->GetScalarType()
> << ", must match output ScalarType "
> << outData[0]->GetScalarType());
> return;
> }
>
> if (inData[1][0]->GetScalarType() != outData[0]->GetScalarType())
> {
> vtkErrorMacro(<< "Execute: input2 ScalarType, "
> << inData[1][0]->GetScalarType()
> << ", must match output ScalarType "
> << outData[0]->GetScalarType());
> return;
> }
>
> // this filter expects that inputs that have the same number of
> components
> if (inData[0][0]->GetNumberOfScalarComponents() !=
> inData[1][0]->GetNumberOfScalarComponents())
> {
> vtkErrorMacro(<< "Execute: input1 NumberOfScalarComponents, "
> << inData[0][0]->GetNumberOfScalarComponents()
> << ", must match out input2 NumberOfScalarComponents "
> << inData[1][0]->GetNumberOfScalarComponents());
> return;
> }
>
> switch (inData[0][0]->GetScalarType())
> {
> double d;
> vtkTemplateMacro(
> vtkMeanSSDExecute(this, inData[0][0],
> inData[1][0], &this->m_value, outExt, id,
> static_cast<VTK_TT *>(0)));
> default:
> vtkErrorMacro(<< "Execute: Unknown ScalarType");
> return;
> }
> }
>
--
Dr Alberto Gomez
Research Associate
Department of Biomedical Engineering
King's College London
020 7188 7188 ext 50871
4th North Wing
St Thomas' Hospital
SE1 7EH London, UK
More information about the vtkusers
mailing list