[vtkusers] Once more -- how do I get the outline contours of regions in a mask/binary image
Bill Lorensen
bill.lorensen at gmail.com
Thu Mar 27 15:57:35 EDT 2008
Kent,
In the past, I did this sort of thing:
extract contours from an image with vtkContourFilter with an isovalue of .5.
convert disconnected line segments into polylines using vtkStripper
Bill
On Thu, Mar 27, 2008 at 12:17 PM, kent williams
<nkwmailinglists at gmail.com> wrote:
> Here's what I want:
>
> Input: Image where pixel == 0 is outside mask pixel != 0 is inside mask
> Output: Polygon/Polyline outline of all regions of nonzero pixels in image.
>
> Reason: My program uses vtkContourWidget to trace anatomical features.
> I save the contours in an ITK xml file. I can also save a binary
> image -- i.e. an image file the same dimensions as the image being
> traced. If a voxel is inside the traced contour, then it has a value
> of 255, otherwise zero.
>
> What I want to do is to read a binary image, and then find all the
> contours around the set pixels. Basically the inverse of the
> operation above -- rather than using outline contours to generate a
> binary image, I want to use a binary image to get a set of outline
> contours.
>
> I've written what I think should work, below, but it looks really,
> really wrong when I load it -- it seems to be some woven agglomeration
> of many too many conour outlines.
>
> #include "MaskToContoursCLP.h"
> #include "itkBrains2MaskImageIO.h"
> #include "itkBrains2MaskImageIOFactory.h"
> #include "itkIO.h"
> #include "vtkKWImage.h"
> #include <vtkAppendPolyData.h>
> #include <vtkContourFilter.h>
> #include <vtkXMLPolyDataWriter.h>
> #include <vtkPolyData.h>
> #include <vtkImageData.h>
> #include <vtkDataObject.h>
> #include <vtkMath.h>
> #include <vtkIdTypeArray.h>
> #include <vtkCellArray.h>
> #include <itkExtractImageFilter.h>
>
> //
> // ConverPointSequenceToPolyData
> // does what the name says -- code borrowed from someon
> // on the VTK Mailing list
> //
> // logically what you end up is a collection of lines; at some point
> // we might want something different. At which point I'll have to dive
> // in and learn what the hell is going on with vtkCells.
> void ConvertPointSequenceToPolyData(vtkPoints* inPts,
> const int& closed,
> vtkPolyData* outPoly)
> {
> if ( !inPts || !outPoly )
> {
> return;
> }
>
> int npts = inPts->GetNumberOfPoints();
>
> if ( npts < 2 )
> {
> return;
> }
>
> double p0[3];
> double p1[3];
> inPts->GetPoint( 0, p0 );
> inPts->GetPoint( npts - 1, p1 );
> if ( p0[0] == p1[0] && p0[1] == p1[1] && p0[2] == p1[2] && closed )
> {
> --npts;
> }
>
> vtkPoints* temp = vtkPoints::New();
> temp->SetNumberOfPoints( npts );
> for ( int i = 0; i < npts; ++i )
> {
> temp->SetPoint( i, inPts->GetPoint( i ) );
> }
>
> vtkCellArray *cells = vtkCellArray::New();
> cells->Allocate( cells->EstimateSize( npts + closed, 2 ) );
>
> cells->InsertNextCell( npts + closed );
>
> for ( int i = 0; i < npts; ++i )
> {
> cells->InsertCellPoint( i );
> }
>
> if ( closed )
> {
> cells->InsertCellPoint( 0 );
> }
>
> outPoly->SetPoints( temp );
> temp->Delete();
> outPoly->SetLines( cells );
> cells->Delete();
> }
> //--------------------------------------------------------------------------
> // assumes a piecwise linear polyline with points at discrete locations
> //
> int
> ReducePolyData2D(vtkPolyData* inPoly,
> vtkPolyData* outPoly, const int& closed )
> {
> if ( !inPoly || !outPoly )
> {
> return 0;
> }
>
> vtkPoints* inPts = inPoly->GetPoints();
> if ( !inPts )
> {
> return 0;
> }
> int n = inPts->GetNumberOfPoints();
> if ( n < 3 )
> {
> return 0;
> }
>
> double p0[3];
> inPts->GetPoint( 0, p0 );
> double p1[3];
> inPts->GetPoint( n - 1, p1 );
> bool minusNth = ( p0[0] == p1[0] && p0[1] == p1[1] && p0[2] == p1[2] );
> if ( minusNth && closed )
> {
> --n;
> }
>
> struct frenet
> {
> double t[3]; // unit tangent vector
> bool state; // state of kappa: zero or non-zero T/F
> };
>
> frenet* f;
> f = new frenet[n];
> double tL;
> // calculate the tangent vector by forward differences
> for ( int i = 0; i < n; ++i )
> {
> inPts->GetPoint( i, p0 );
> inPts->GetPoint( ( i + 1 ) % n, p1 );
> tL = sqrt( vtkMath::Distance2BetweenPoints( p0, p1 ) );
> if ( tL == 0.0 ){ tL = 1.0; }
> for ( int j = 0 ; j < 3; ++j )
> {
> f[i].t[j] = (p1[j] - p0[j]) / tL;
> }
> }
>
> // calculate kappa from tangent vectors by forward differences
> // mark those points that have very low curvature
> double eps = 1.e-10;
>
> for ( int i = 0; i < n; ++i )
> {
> f[i].state = ( fabs( vtkMath::Dot( f[i].t, f[( i + 1 ) % n].t ) - 1.0 )
> < eps );
> }
>
> vtkPoints* tempPts = vtkPoints::New();
>
> // mark keepers
> vtkIdTypeArray* ids = vtkIdTypeArray::New();
>
> // for now, insist on keeping the first point for closure
> ids->InsertNextValue( 0 );
>
> for ( int i = 1; i < n; ++i )
> {
> bool pre = f[( i - 1 + n ) % n].state; // means fik != 1
> bool cur = f[i].state; // means fik = 1
> bool nex = f[( i + 1 ) % n].state;
>
> // possible vertex bend patterns for keep: pre cur nex
> // 0 0 1
> // 0 1 1
> // 0 0 0
> // 0 1 0
>
> // definite delete pattern
> // 1 1 1
>
> bool keep = false;
>
> if ( pre && cur && nex ) { keep = false; }
> else if ( !pre && !cur && nex ) { keep = true; }
> else if ( !pre && cur && nex ) { keep = true; }
> else if ( !pre && !cur && !nex ) { keep = true; }
> else if ( !pre && cur && !nex ) { keep = true; }
>
> if ( keep ) // not a current sure thing but the preceding delete was
> {
> ids->InsertNextValue( i );
> }
> }
>
> for ( int i = 0; i < ids->GetNumberOfTuples(); ++i )
> {
> tempPts->InsertNextPoint( inPts->GetPoint( ids->GetValue( i ) ) );
> }
>
> if ( closed )
> {
> tempPts->InsertNextPoint( inPts->GetPoint( ids->GetValue( 0 ) ) );
> }
>
> ConvertPointSequenceToPolyData( tempPts, closed, outPoly );
>
> ids->Delete();
> tempPts->Delete();
> delete [] f;
> return 1;
> }
>
> /*This program generates a signed distance image for a binary image.
> The intensity values inside the structues will be positive and outside
> the structure will be negative in the signed distance image. */
>
> int main( int argc, char * argv[] )
> {
> PARSE_ARGS;
> itk::Brains2MaskImageIOFactory::RegisterOneFactory();
>
> const unsigned int Dimension = 3;
> typedef float PixelType;
> typedef itk::Image< PixelType, Dimension > ImageType;
> typedef ImageType SliceImageType;
> // typedef itk::Image<PixelType,2> SliceImageType;
> typedef itk::ExtractImageFilter<ImageType,SliceImageType> ExtractFilterType;
>
> ImageType::Pointer inputImage(itkUtil::ReadImage<ImageType>(InputImage));
>
> ImageType::SizeType size = inputImage->GetLargestPossibleRegion().GetSize();
> //
> // accumulate PolyData
> vtkAppendPolyData *accumulator = vtkAppendPolyData::New();
>
> //
> // for each slice
> for(unsigned i = 0; i < size[2]; i++)
> {
> //
> // extract a slice from the image.
> ImageType::RegionType extractRegion;
>
> ImageType::SizeType extractSize(size);
> extractSize[2] = 1;
>
> ImageType::IndexType extractIndex;
> extractIndex[0] = 0;
> extractIndex[1] = 0;
> extractIndex[2] = i;
>
> extractRegion.SetSize(extractSize);
> extractRegion.SetIndex(extractIndex);
>
> ExtractFilterType::Pointer extractFilter = ExtractFilterType::New();
> extractFilter->SetInput(inputImage);
> extractFilter->SetExtractionRegion(extractRegion);
> extractFilter->Update();
>
> SliceImageType::Pointer sliceImage(extractFilter->GetOutput());
>
> vtkKWImage *adapter = vtkKWImage::New();
> adapter->SetITKImageBase(sliceImage);
> vtkContourFilter *contourFilter = vtkContourFilter::New();
> contourFilter->SetValue(0,0.0);
> contourFilter->SetValue(1,1.0);
> contourFilter->SetInput(vtkDataObject::SafeDownCast(adapter->GetVTKImage()));
> contourFilter->Update();
> vtkPolyData *pd = contourFilter->GetOutput();
> vtkPolyData *pd2 = vtkPolyData::New();
> ReducePolyData2D(pd,pd2,true);
> accumulator->AddInput(pd2);
> adapter->Delete();
> contourFilter->Delete();
> }
> accumulator->Update();
> vtkXMLPolyDataWriter *writer = vtkXMLPolyDataWriter::New();
> writer->SetInput(accumulator->GetOutput());
> writer->SetDataModeToAscii();
> writer->SetFileName(OutputContours.c_str());
> writer->Write();
> writer->Delete();
>
> return EXIT_SUCCESS;
>
> }
> _______________________________________________
> This is the private VTK discussion list.
> Please keep messages on-topic. Check the FAQ at: http://www.vtk.org/Wiki/VTK_FAQ
> Follow this link to subscribe/unsubscribe:
> http://www.vtk.org/mailman/listinfo/vtkusers
>
More information about the vtkusers
mailing list