[Insight-users] Creating a 2D binary image from a polygon

[Luis Ibanez]

Hi AurelieC,


You may want to try the following filters:

http://www.itk.org/Insight/Doxygen/html/classitk_1_1PolylineMask2DImageFilter.html
http://www.itk.org/Insight/Doxygen/html/classitk_1_1PolylineMaskImageFilter.html


or you could also do it by combining the PolygonSpatialObject with the
SpatialObjectToImageFilter:

http://www.itk.org/Insight/Doxygen/html/classitk_1_1PolygonSpatialObject.html
http://www.itk.org/Insight/Doxygen/html/classitk_1_1SpatialObjectToImageFilter.html

as illustrated in the recently added examples:
http://www.itk.org/cgi-bin/viewcvs.cgi/Examples/Filtering/SpatialObjectToImage1.cxx?root=Insight&sortby=date&view=log

   Insight/Examples/Filtering/SpatialObjectToImage1.cxx




   Regards,


       Luis



-----------------------
AurelieC wrote:
> Hello,
> 
> I have a 2D polygonal region of interest defined by a list of points (x and
> y coordinates). I am trying to create a new image of a given size which will
> be white inside the polygon and black outside.
> 
> What I tried to do : 
> I create a mesh, I add in the mesh the points, then LineCell between the
> points and vertexCell at each points.
> Then I create a MeshSpatialObject which take the created mesh.
> Then I try to create a binary image using the MeshSpatialObject : first I
> was using the MeshToImageFilter but it did not work so I tried to scan all
> the image and for each pixel I just check if it is inside or outside the
> polygon, if it is inside it gets the value 1.
> The code compiles, but nothing happens, no point is considered as inside the
> mesh, the outside image is black, so I think I made a mistake while creating
> the mesh. Maybe I should use PolygonCell, or PolygonGroupSpatialObject?
> 
> Please help me.
> 
> This is my code:
> 
> 
> // creating the mesh
> 
> typedef   float   PixelType;  
> const unsigned int Dimension = 2;
> typedef itk::Mesh< PixelType, Dimension>   MeshType;
> MeshType::Pointer  mesh = MeshType::New();
> 
> // adding points and creating the cells
> 
> MeshType::PointType p;
> 
> typedef MeshType::CellType CellType;
> typedef itk::LineCell< CellType > LineType;
> typedef itk::VertexCell< CellType > VertexType;
> 
> // tab contained the coordinates of each points
> for (int i = 0; i < nbPoint; i++)
> {
>     p[0] = tab[i*2];
>     p[1] = tab[i*2+1];
>     mesh->SetPoint(i,p);
> }
> 
> //creating lines between points
> CellType::CellAutoPointer cellpointer;
> 
> for(int cellId=0; cellId < nbPoint-1; cellId++)
> {
>     cellpointer.TakeOwnership( new LineType );
>     cellpointer->SetPointId( 0, cellId ); // first point
>     cellpointer->SetPointId( 1, cellId+1 ); // second point
>     mesh->SetCell( cellId, cellpointer ); // insert the cell
> }
> 
> cellpointer.TakeOwnership( new LineType );
> cellpointer->SetPointId( 0, (nbPoint-1) ); // first point
> cellpointer->SetPointId( 1, 0 ); // second point
> mesh->SetCell( nbPoint-1, cellpointer); // insert the cell
> 
> for(int cellId=0; cellId < nbPoint; cellId++)
> {
>     cellpointer.TakeOwnership( new VertexType );
>     cellpointer->SetPointId( 0, cellId );
>     mesh->SetCell( cellId + nbPoint, cellpointer );
> }
> 
> 
> std::cout << "# Points= " << mesh->GetNumberOfPoints() << std::endl;
> std::cout << "# Cell  = " << mesh->GetNumberOfCells() << std::endl;
> 
> 
> // creating the MeshSpatialObject
> typedef itk::MeshSpatialObject<MeshType>     MeshSpatialObjectType;
> MeshSpatialObjectType::Pointer myMeshSpatialObject = 
> MeshSpatialObjectType::New();
> 
> myMeshSpatialObject->SetMesh(mesh);
> 
> 
> // creating the output image
> 
> typedef itk::Image<float, 2>      ImageType;
> ImageType::Pointer image = ImageType::New();
> 
> // dimX and dimY given dimensions
> size[0] = d->dimX;
> size[1] = d->dimY;
> ImageType::IndexType start;
> start[0] = 0; // first index on X
> start[1] = 0; // first index on Y
> ImageType::RegionType region;
> region.SetSize( size );
> region.SetIndex( start );
> 
> image->SetRegions( region );
> image->Allocate();
> 
> myMeshSpatialObject->Update();
> 
> 
> // region to scan. I reduce the region to scan using the
> mesh->GetBoundingBox()
> 
> typedef itk::FixedArray< float, 4> FixedArraySize;
> FixedArraySize test;
> myMeshSpatialObject->GetBoundingBox()->ComputeBoundingBox();
> test =  myMeshSpatialObject->GetBoundingBox()->GetBounds();
> 
> ImageType::RegionType roi;
> 
> size[0] = int(test[1])+2 -  int(test[0]);
> size[1] =  int(test[3]) +2 - int(test[2]);
> roi.SetSize( size );
> // roi.SetIndex( start );
> 
> 
> // test each pixel of the boundingbox, check if it is inside or outside the
> polygon
> itk::Point<float,2> point;
> typedef itk::ImageRegionIterator<ImageType> Iterator;
> Iterator it( image, roi );
> it.GoToBegin();
>   
>   while( !it.IsAtEnd() )
>     {
> 
>     image->TransformIndexToPhysicalPoint(it.GetIndex(), point );
>     
>     if (myMeshSpatialObject->IsInside(point))
>     {
> 	qDebug() << "isinside "<< endl; // no point is in this case...
> 	image->SetPixel(it.GetIndex(), 1);
>     }
>     else 
> 	image->SetPixel(it.GetIndex(), 0);
> 
>     ++it;
>     }
> 
> 
> 
> // write the result
> 
> typedef itk::ImageFileWriter< ImageType > ImageWriterType;
> ImageWriterType::Pointer imgWriter = ImageWriterType::New();
> imgWriter->SetInput(image);
> imgWriter->SetFileName(filename);
> imgWriter->Update();
> 
> 
> 
> 
> Thank you!