[vtkusers] Threading issue when extracting surfaces from a labeled ImageData

[Sébastien Valette]

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);
}