[vtkusers] Threading issue when extracting surfaces from a labeled ImageData
Jérôme
jerome.velut at gmail.com
Thu Mar 31 06:34:18 EDT 2011
Hi Sébastien !
I played around with your problem, as I had to use MultiThreader in VTK. It
seems that vtkThreshold is not thread-safe. Attached is your modified file and
a CMakeLists.
In a nutshell, I integrated a MutexLock in you MyThreadHelper class, then
I locked the update of the threshold. Well, yes: it means that Threshold is no
more multithreaded. But Geometry and Writer are!
I am sure there is other workflow that would work. Take a look a
vtkContourFilter: It allows you to select several scalar values on which you
want to extract an isosurface. Maybe it will be of interest for you.
HTH,
Jerome
Le 30 mars 2011 17:59, Sébastien Valette
<sebastien.valette at creatis.insa-lyon.fr> a écrit :
> Hi,
>
> I'm trying to extract the surfaces from labeled volume data. The source code
> is pasted at the end of this message. The code works when the number of
> threads is set to 1 (uncomment the only commented line) but when using
> multithreading, the output is randomly corrupted (but parts of the surfaces
> are still correct).
>
> My workflow is as follows :
> ImageData->vtkThreshold->vtkDataSetSurfaceFilter->vtkPolyDataWriter.
>
> Does anyone knows where the problem could come from? Or is there any better
> workflow for this? I'm using vtk 5.4.2 on linux (Ubuntu 10.10).
>
> Thanks,
>
> Sebastien
>
> #include <sstream>
> #include <vtkImageData.h>
> #include <vtkMetaImageReader.h>
> #include <vtkIdList.h>
> #include <vtkThreshold.h>
> #include <vtkDataSetSurfaceFilter.h>
> #include <vtkPolyDataWriter.h>
> #include <vtkPointData.h>
> #include <vtkCellData.h>
> #include <vtkMultiThreader.h>
>
> class MyThreaderHelperClass
> {
> public:
>
> vtkImageData *Image;
> vtkIdList *Labels;
> };
>
> VTK_THREAD_RETURN_TYPE ThreadedSurfaceExtraction (void *arg)
> {
> vtkMultiThreader::ThreadInfo *Info = (vtkMultiThreader::ThreadInfo*) arg;
> MyThreaderHelperClass *Helper = (MyThreaderHelperClass *) Info->UserData;
> int MyId = Info->ThreadID;
> int NumberOfThreads=Info->NumberOfThreads;
>
> vtkImageData *Image=vtkImageData::New();
> Image->ShallowCopy(Helper->Image);
> vtkIdList *Labels=Helper->Labels;
>
> vtkThreshold *Threshold=vtkThreshold::New();
> Threshold->SetInput(Image);
>
> vtkDataSetSurfaceFilter *Geometry=vtkDataSetSurfaceFilter::New();
> vtkPolyDataWriter *Writer=vtkPolyDataWriter::New();
>
> for (int i=MyId;i<Labels->GetNumberOfIds();i+=NumberOfThreads)
> {
> double Level=Labels->GetId(i);
> Threshold->ThresholdBetween(Level-0.5,Level+0.5);
> Geometry->SetInput((vtkDataObject *)Threshold->GetOutput());
> Geometry->Update();
> std::stringstream Name;
> Name<<Level<<".vtk";
> Writer->SetFileName(Name.str().c_str());
> Geometry->GetOutput()->GetCellData()->Initialize();
> Geometry->GetOutput()->GetPointData()->Initialize();
> Writer->SetInput(Geometry->GetOutput());
> Writer->Write();
> }
> Threshold->Delete();
> Geometry->Delete();
> Writer->Delete();
> return (VTK_THREAD_RETURN_VALUE);
> }
>
> int main( int argc, char *argv[] )
> {
>
> if (argc<2)
> {
> cout<<"Usage : VolumeAnalysis file.mhd"<<endl;
> exit(1);
> }
>
> // Load Volume
> cout <<"load : "<<argv[1]<<endl;
>
> vtkMetaImageReader *Reader=vtkMetaImageReader::New();
> Reader->SetFileName(argv[1]);
> Reader->Update();
> vtkImageData *Image=Reader->GetOutput();
>
> vtkIdList *Labels=vtkIdList::New();
> int Dimensions[3];
>
> Image->GetDimensions(Dimensions);
> std::cout<<"Image dimensions : "<<Dimensions[0]<<" "<<Dimensions[1]<<"
> "<<Dimensions[2]<<endl;
>
> //collect labels
> for (int z=0;z<Dimensions[2];z++)
> {
> for (int y=0;y<Dimensions[1];y++)
> {
> for (int x=0;x<Dimensions[0];x++)
> {
> int Label=Image->GetScalarComponentAsFloat (x,y,z,0);
> if (Labels->IsId(Label)==-1)
> Labels->InsertNextId(Label);
> }
> }
> }
> cout<<"Found "<<Labels->GetNumberOfIds()<<" labels"<<endl;
>
> MyThreaderHelperClass Helper;
> Helper.Image=Image;
> Helper.Labels=Labels;
>
> vtkMultiThreader *Threader=vtkMultiThreader::New();
> // Threader->SetNumberOfThreads(1);
> Threader->SetSingleMethod (ThreadedSurfaceExtraction, (void *) &Helper);
> Threader->SingleMethodExecute ();
> return (0);
> }
> _______________________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Please keep messages on-topic and check the VTK FAQ at:
> http://www.vtk.org/Wiki/VTK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://www.vtk.org/mailman/listinfo/vtkusers
>
-------------- next part --------------
#include <sstream>
#include <vtkImageData.h>
#include <vtkMetaImageReader.h>
#include <vtkIdList.h>
#include <vtkThreshold.h>
#include <vtkDataSetSurfaceFilter.h>
#include <vtkPolyDataWriter.h>
#include <vtkPointData.h>
#include <vtkCellData.h>
#include <vtkMultiThreader.h>
#include <vtkMutexLock.h>
class MyThreaderHelperClass
{
public:
vtkImageData *Image;
vtkIdList *Labels;
vtkMutexLock* Mutex;
};
VTK_THREAD_RETURN_TYPE ThreadedSurfaceExtraction (void *arg)
{
vtkMultiThreader::ThreadInfo *Info = (vtkMultiThreader::ThreadInfo*) arg;
MyThreaderHelperClass *Helper = (MyThreaderHelperClass *) Info->UserData;
int MyId = Info->ThreadID;
int NumberOfThreads=Info->NumberOfThreads;
vtkIdList *Labels=Helper->Labels;
vtkImageData *Image=vtkImageData::New();
Image->ShallowCopy(Helper->Image);
vtkThreshold *Threshold=vtkThreshold::New();
Threshold->SetInput(Image);
vtkDataSetSurfaceFilter *Geometry=vtkDataSetSurfaceFilter::New();
Geometry->SetInputConnection(Threshold->GetOutputPort());
vtkPolyDataWriter *Writer=vtkPolyDataWriter::New();
Writer->SetInputConnection(Geometry->GetOutputPort());
for (int i=MyId; i<Labels->GetNumberOfIds(); i+=NumberOfThreads)
{
double Level=Labels->GetId(i);
Threshold->ThresholdBetween(Level-0.5,Level+0.5);
Helper->Mutex->Lock();
cout << "Thread#" << MyId << " - Label#" << i << " - Level="<<Level << endl;
Threshold->Update( );
Helper->Mutex->Unlock();
std::stringstream Name;
Name<<Level<<".vtk";
Writer->SetFileName(Name.str().c_str());
Writer->Write();
}
Threshold->Delete();
Geometry->Delete();
Writer->Delete();
Image->Delete();
return (VTK_THREAD_RETURN_VALUE);
}
int main( int argc, char *argv[] )
{
if (argc<2)
{
cout<<"Usage : VolumeAnalysis file.mhd"<<endl;
exit(1);
}
// Load Volume
cout <<"load : "<<argv[1]<<endl;
vtkMetaImageReader *Reader=vtkMetaImageReader::New();
Reader->SetFileName(argv[1]);
Reader->Update();
vtkImageData *Image=Reader->GetOutput();
vtkIdList *Labels=vtkIdList::New();
int Dimensions[3];
Image->GetDimensions(Dimensions);
std::cout<<"Image dimensions : "<<Dimensions[0]<<" "<<Dimensions[1]<<" "<<Dimensions[2]<<endl;
//collect labels
for (int z=0; z<Dimensions[2]; z++)
{
for (int y=0; y<Dimensions[1]; y++)
{
for (int x=0; x<Dimensions[0]; x++)
{
int Label=Image->GetScalarComponentAsFloat (x,y,z,0);
if (Labels->IsId(Label)==-1)
Labels->InsertNextId(Label);
}
}
}
cout<<"Found "<<Labels->GetNumberOfIds()<<" labels"<<endl;
MyThreaderHelperClass Helper;
Helper.Image=Image;
Helper.Labels=Labels;
Helper.Mutex=vtkMutexLock::New();
vtkMultiThreader *Threader=vtkMultiThreader::New();
//Threader->SetNumberOfThreads(1);
Threader->SetSingleMethod (ThreadedSurfaceExtraction, (void *) &Helper);
Threader->SingleMethodExecute ();
return (0);
}
-------------- next part --------------
project(VolumeAnalysis)
### The following are needed to avoid warnings. See CMake doc. ######
CMAKE_MINIMUM_REQUIRED(VERSION 2.8)
if(COMMAND cmake_policy)
cmake_policy(SET CMP0003 NEW)
endif(COMMAND cmake_policy)
FIND_PACKAGE(VTK REQUIRED)
INCLUDE(${VTK_USE_FILE})
add_executable(VolumeAnalysis VolumeAnalysis.cxx)
TARGET_LINK_LIBRARIES( VolumeAnalysis vtkCommon vtkGraphics vtkImaging vtkIO)
More information about the vtkusers
mailing list