[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