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

[Luis Ibanez]

Hi Aurelie,


             Congratulations !


     You just found a bug in ITK.   :-)


To put it short: The logic of the "IsInside" method
in the itkPolygonSpatialObject.txx file was flawed.

We have just committed a fix for it:
http://www.itk.org/cgi-bin/viewcvs.cgi/Code/SpatialObject/itkPolygonSpatialObject.txx?root=Insight&r1=1.24&r2=1.25&sortby=date

Also, for your convenience, we also added an example
on how to use the itkPolygonSpatialObject along with
the SpatialObjectToImageFilter in order to rasterize
a polygon.


You will find this new example in the file:

     Insight/Examples/Filtering/
               SpatialObjectToImage23.cxx

http://www.itk.org/cgi-bin/viewcvs.cgi/Examples/Filtering/SpatialObjectToImage3.cxx?root=Insight&sortby=date&view=log


Regarding the configuration of the polygon, please note
that you don't need to duplicate the first point.
(but it doesn't hurt do it either...)


If you look at the source code of the IsInside() method
in Insight/Code/SpatialObjects/itkPolygonSpatialObject.txx
lines :435-439 you will find:

   // if last point same as first, don't bother with it.
   if(start == last)
     {
     numpoints--;
     }

[this is before the changes committed for fixing the bug].


In other words, for the purpose of the IsInside() method,
this class is assuming that the polygon is always closed.

I guess the correctness of the assumption is arguable...
but that's what it does right now...


Please update your checkout of ITK, or apply the CVS diffs
in the files:

      itkPolygonSpatialObject.h
      itkPolygonSpatialObject.txx

and give it a try at running your test again.
(or run the new example in Examples/Filtering).

and let us know if you find any problem.


     Thanks


        Luis


----------------------
Aurelie Canale wrote:
> Hi Luis,
> 
> so as I said, I am using a combination of PolygonSpatialObject with 
> SpatialObjectToImageFilter. I am adding all the points of the polygon, 
> and at the end I add one again the first point, to close the polygon, I 
> check that the polygon is closed and it is the case. But at the end of 
> all the process, when I visualize the binary image, there is always a 
> kind of leak (see the attached image). Do you know where it could come 
> from? I saw that there was bug about the volume calculation of a polygon 
> in 2D, could it be linked to my issue?
> 
> Thanks for your help,
> Aurelie Canale
> 
> Luis Ibanez wrote:
> 
>>
>> 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!
> 
> 
> 
> ------------------------------------------------------------------------
>