《itk实用demo》-轮廓提取

轮廓提取

    typedef itk::SimpleContourExtractorImageFilter <ImageType, ImageType>
    SimpleContourExtractorImageFilterType;
    SimpleContourExtractorImageFilterType::Pointer contourFilter
    = SimpleContourExtractorImageFilterType::New();
    contourFilter->SetInput(slicedata);
    contourFilter->Update();
--------------------------------------------------------------------------------
vtkPolyData* getContour(SliceType *slicedata)
{   
    cout<<"start get contour ..."<<endl;

    typedef itk::Image<unsigned char ,2> USliceType ;
    typedef itk::Image<float,2> FSliceType;

    typedef itk::CastImageFilter<SliceType ,USliceType > UCastImageType;
    UCastImageType::Pointer ucastFilter = UCastImageType::New();
    ucastFilter->SetInput(slicedata);
    ucastFilter->Update();
    USliceImageType::Pointer slice2 = ucastFilter->GetOutput();

    //计算每个像素到轮廓的距离
    typedef itk::ApproximateSignedDistanceMapImageFilter<USliceType ,FSliceType> DistanceMapType;
    DistanceMapType::Pointer distanceFilter = DistanceMapType::New();
    distanceFilter->SetInput(ucastFilter->GetOutput());
    distanceFilter->SetInsideValue(255);
    distanceFilter->SetOutsideValue(0);
    try
    {
        distanceFilter->Update();
    }
    catch(itk::ExceptionObject& er)
    {
        cout<<er<<endl;
    }

    //提取轮廓
    typedef itk::ContourExtractor2DImageFilter<FSliceType> ContourExtractorType;
    ContourExtractorType::Pointer contourFilter = ContourExtractorType::New();
    contourFilter->SetInput(distanceFilter->GetOutput());
    contourFilter->SetContourValue(0);
    try
    {
        contourFilter->Update();
    }
    catch(itk::ExceptionObject& er)
    {
        cout<<er<<endl;
    }
    //cout<<"contour num is :"<<contourFilter->GetNumberOfOutputs()<<endl;

    //找出最大轮廓,并依次保存轮廓点
    vtkPolyData* myContour = vtkPolyData::New();

    if (contourFilter->GetNumberOfOutputs()!=0)
    {
        vtkCellArray* cell = vtkCellArray::New();
        vtkPoints* points = vtkPoints::New();
        int pointId = 0 ;   
        int max = contourFilter->GetOutput(0)->GetVertexList()->size();
        int ind = 0;
        for(unsigned int i = 0; i < contourFilter->GetNumberOfOutputs(); ++i)
        {       
            if(contourFilter->GetOutput(i)->GetVertexList()->size()>max)
            {
                max = contourFilter->GetOutput(i)->GetVertexList()->size();
                ind = i;
            }
        }

        //获取轮廓的物理坐标
        //cout<<"contour size is :"<<contourFilter->GetOutput(ind)->GetVertexList()->size()<<endl;
        vtkPolyLine* polyLine = vtkPolyLine::New();
        ContourExtractorType::VertexListType::ConstIterator vertexIterator = contourFilter->GetOutput(ind)->GetVertexList()->Begin();       
        while(vertexIterator != contourFilter->GetOutput(ind)->GetVertexList()->End())
        {
            SliceImageType::PointType pointTemp;
            itk::ContinuousIndex<double,2> continousIndex;
            continousIndex[0] = vertexIterator->Value()[0];
            continousIndex[1] = vertexIterator->Value()[1];
            slicedata->TransformContinuousIndexToPhysicalPoint(continousIndex,pointTemp);
            double point[3];
            point[0] = pointTemp[0]*spacing[0];
            point[1] = pointTemp[1]*spacing[1];
            point[2] =  znum; 
            pointId = points->InsertNextPoint(point);
            polyLine->GetPointIds()->InsertNextId(pointId);
            ++vertexIterator;
        }
        //polyLine->GetPointIds()->InsertNextId(0);
        cell->InsertNextCell(polyLine);
        polyLine->Delete();
        myContour->SetPoints(points);
        myContour->SetLines(cell);
        //myContour->SetPolys(cell);

        cell->Delete();
        myContour->Modified();
        //cout<<"end calc contour ..."<<endl;

        points->Delete();
    }

    return myContour;
}