[Insight-users] Re: [Insight-developers] Filter Update Woes

[Karthik Krishnan]

To answer my own question, I found a bug in 
itkImageToParametricSpaceFilter.txx. The filter creates new point data 
for the mesh during every update call wiping out existing point data. 
The bug was resolved by replacing:

template <class TInputImage, class TOutputMesh>
void
ImageToParametricSpaceFilter<TInputImage,TOutputMesh>
::GenerateOutputInformation()
{
  OutputMeshPointer          mesh       = this->GetOutput();
  PointsContainerPointer     points     = mesh->GetPoints();
  PointDataContainerPointer pointData  = PointDataContainer::New();
}

with

template <class TInputImage, class TOutputMesh>
void
ImageToParametricSpaceFilter<TInputImage,TOutputMesh>
::GenerateOutputInformation()
{
  OutputMeshPointer          mesh       = this->GetOutput();
  PointsContainerPointer     points     = mesh->GetPoints();
 
  PointDataContainerPointer  pointData;
  if( mesh->GetPointData() )
    {
    pointData = mesh->GetPointData();
    }
  else
    { // Create
    pointData  = PointDataContainer::New();
    }
 }

Thanks
Regards
Karthik

Karthik Krishnan wrote:

> Hi,
>
> I am having a curious problem with updates. I've posted a minimal 
> example.
> The pipeline is as follows:
>
> Image  -> CastImageFilter -> ImageToParametricSpaceFilter -> 
> InteriorExteriorMeshFilter
>
>
> In the example, I create a 3x3x3 image. The cast image filter really 
> does nothing. It just spits the input as the output.  The 
> ImageToParametricSpaceFilter takes the three images (each of them 
> being the same input image).
> This InteriorExteriorMeshFilter uses a SphereSpatialFunction, whose 
> radius is so large that all points are inside.
>
> The output is fine if Update is called once.
> as in m_SpatialFunctionFilter->Update();
>
> If Update is called twice in succession as in
> as in m_SpatialFunctionFilter->Update(); 
> m_SpatialFunctionFilter->Update();
> the mesh point data gets garbled.       
> If I remove the filter, everything works fine.
>
> Any suggestions ?
>
> Thank you
> Regards
> Karthik
>
>------------------------------------------------------------------------
>
>#include "itkInteriorExteriorMeshFilter.h"
>#include "itkSphereSpatialFunction.h"
>#include "itkImageToParametricSpaceFilter.h"
>#include "itkImage.h"
>#include "itkMesh.h"
>#include "itkCastImageFilter.h"
>
>int main()
>{
>  // Create a simple image 3x3x3 image
>  typedef itk::Image< float, 3 > ImageType;
>  ImageType::Pointer image = ImageType::New();
>  ImageType::IndexType start;
>  ImageType::SizeType  size;
>  size.Fill(3);
>  start[0] =   0; 
>  start[1] =   0;  
>  start[2] =   0;
>  ImageType::RegionType region;
>  region.SetSize( size );
>  region.SetIndex( start );
>  image->SetRegions( region );
>  image->Allocate();
>
>  ImageType::IndexType pixelIndex;
>  ImageType::PixelType   pixelValue = 0;
>  for( unsigned int i = 0; i<size[0]; i++)
>    {
>    for( unsigned int j=0; j< size[1]; j++)
>      {
>      for( unsigned int k=0; k< size[1]; k++)
>        {
>        pixelIndex[0] = i;
>        pixelIndex[1] = j;
>        pixelIndex[2] = k;
>        image->SetPixel(   pixelIndex,   ++pixelValue  );
>        }
>      }
>    }
>  
>
>  // Use a cast image filter, or for that matter any other filter
>  // The filter here just passes the input to the output
>  typedef itk::CastImageFilter< ImageType, ImageType > CastImageFilterType;
>  CastImageFilterType::Pointer caster1 = CastImageFilterType::New();
>  CastImageFilterType::Pointer caster2 = CastImageFilterType::New();
>  CastImageFilterType::Pointer caster3 = CastImageFilterType::New();
>  caster1->SetInput(image);
>  caster2->SetInput(image);
>  caster3->SetInput(image);
>  
>
>  // Image to Parametric Space filter
>  typedef   itk::Point<float,3>                   MeshPointDataType;
>  typedef   itk::Mesh< MeshPointDataType, 3 >     MeshType;
>  typedef   itk::Mesh< MeshType::PointType, 3 >   ImageSpaceMeshType;
>  typedef   itk::ImageToParametricSpaceFilter< ImageType, MeshType >
>                                                  ParametricSpaceFilterType;
>  ParametricSpaceFilterType::Pointer m_ParametricSpace = ParametricSpaceFilterType::New();
>  m_ParametricSpace->SetInput( 0, caster1->GetOutput() );
>  m_ParametricSpace->SetInput( 1, caster2->GetOutput() );
>  m_ParametricSpace->SetInput( 2, caster3->GetOutput() );
>  
>  
>  // Map selected points in parametric space to image space.. 
>  typedef   itk::SphereSpatialFunction< 
>                                MeshType::PointDimension,
>                                MeshType::PointType >  SphereSpatialFunctionType;
>  typedef  SphereSpatialFunctionType          SpatialFunctionType;
>  typedef itk::InteriorExteriorMeshFilter<
>                                        MeshType,
>                                        MeshType,
>                                        SpatialFunctionType  >   
>                                                   SpatialFunctionFilterType;
>  SpatialFunctionFilterType::Pointer m_SpatialFunctionFilter 
>                                          = SpatialFunctionFilterType::New();
>  m_SpatialFunctionFilter->SetInput( m_ParametricSpace->GetOutput() );
>  SpatialFunctionType::Pointer m_SpatialFunction = SpatialFunctionType::New();
>  m_SpatialFunction->SetRadius( 100 );
>  m_SpatialFunctionFilter->SetSpatialFunction( m_SpatialFunction );
>
>
>  // Update and print out output....
>  {
>    m_SpatialFunctionFilter->Update();
>    MeshType::PointDataContainerPointer  pointData   = m_ParametricSpace->GetOutput()->GetPointData();
>    MeshType::PointsContainerPointer     points      = m_ParametricSpace->GetOutput()->GetPoints();
>    MeshType::PointsContainerIterator    pointIt     = points->Begin();
>    MeshType::PointDataContainerIterator pointDataIt = pointData->Begin();
>    while (pointIt != points->End() )
>    {
>      std::cout << "Image space point: " << pointDataIt.Value() 
>        << " mapped to Parametric space location: " << pointIt.Value() << std::endl;
>      ++pointIt;
>      ++pointDataIt;
>    }
>  }
>
>  // Update again and print out output....
>  {
>    m_SpatialFunctionFilter->Update();
>    MeshType::PointDataContainerPointer  pointData   = m_ParametricSpace->GetOutput()->GetPointData();
>    MeshType::PointsContainerPointer     points      = m_ParametricSpace->GetOutput()->GetPoints();
>    MeshType::PointsContainerIterator    pointIt     = points->Begin();
>    MeshType::PointDataContainerIterator pointDataIt = pointData->Begin();
>    while (pointIt != points->End() )
>    {
>      std::cout << "Image space point: " << pointDataIt.Value() 
>        << " mapped to Parametric space location: " << pointIt.Value() << std::endl;
>      ++pointIt;
>      ++pointDataIt;
>    }
>  }
>
>  return EXIT_SUCCESS;
>}
>  
>
>------------------------------------------------------------------------
>
>_______________________________________________
>Insight-developers mailing list
>Insight-developers at itk.org
>http://www.itk.org/mailman/listinfo/insight-developers
>  
>