[Insight-users] fastest way to get image voxels(pixels) within an itkMesh ?

[gregthom]

Hello all,

Since my last post, I have been able to concoct a piece of code that does
what I originally wanted. i'e

take a mesh, binarize it, and then give me all voxels(physical) within my
binary volume. 
I want to use these physical coordinates to interpolate in another 3D volume
with scalar values, i.e for all voxels in my binary volume, what is the
scalar value interpolated from the 3D volume ? My main worry is about speed.
For the first part, is this really the fastest way ? And what about the
interpolation part ? how to do that ?

Thank you


Here is my working code
**********************************************************************************************************************************************
  typedef unsigned char  BinaryPixelType;
  typedef unsigned short PixelType;
  //typedef float PixelType;
  typedef itk::Image<BinaryPixelType, 3>      BinaryImageType;
  typedef itk::Image<PixelType, 3> ImageType;
  typedef itk::TriangleMeshToBinaryImageFilter<MeshType, BinaryImageType>
TriangleMeshToBinaryImageFilterType; 
  TriangleMeshToBinaryImageFilterType::Pointer triangleMeshToImage =
TriangleMeshToBinaryImageFilterType::New();
  typedef itk::MeshSpatialObject<MeshType> MeshSpatialObjectType;
  
  //convert matlab meshS (i.e meshS.vertices, meshS.triangles) to itkmesh
  MeshType::Pointer mymesh =  meshSstructToItkMesh(prhs[0]);   
  triangleMeshToImage->SetInput(mymesh);
  BinaryImageType::SizeType size; 
  size[0] = (unsigned int)((double *)(mxGetPr(prhs[3])))[0];  
  size[1] = (unsigned int)((double *)(mxGetPr(prhs[3])))[1];  
  size[2] = (unsigned int)((double *)(mxGetPr(prhs[3])))[2];  
  

  triangleMeshToImage->SetSpacing(spc);
  triangleMeshToImage->SetOrigin(org);
  triangleMeshToImage->SetSize (size);

  //typedef itk::ImageFileWriter<BinaryImageType> WriterType;
  //WriterType::Pointer writer = WriterType::New();
  //writer->SetFileName("testMeshBinaryImage.mha");
  //writer->SetInput(triangleMeshToImage->GetOutput());

  try
  {     
      triangleMeshToImage->Update();
      //writer->Update();
  }
  catch( itk::ExceptionObject & excep )
  {
    std::cerr << "Exception caught !" << std::endl;
    std::cerr << excep << std::endl;
  }

  //create meshtoimage filter
  BinaryImageType::ValueType  m_MaskValue =
itk::NumericTraits<BinaryImageType::ValueType>::One;  
  typedef  itk::ImageRegionConstIteratorWithIndex<BinaryImageType>
MaskIteratorType;  
  MaskIteratorType mit( triangleMeshToImage->GetOutput(),
triangleMeshToImage->GetOutput()->GetBufferedRegion() );   
  BinaryImageType::IndexType index;
  typedef itk::PointSet< double, 3 > PointSetType;
  typedef PointSetType::PointType InputPointType;

  // lets get number of voxels
  int nvoxels;
  for (mit.GoToBegin(); !mit.IsAtEnd(); ++mit)
  {
     index = mit.GetIndex();
     InputPointType iPoint;
     triangleMeshToImage->GetOutput()->TransformIndexToPhysicalPoint( index,
iPoint );  
     if (mit.Get() == m_MaskValue)
     {        
	 nvoxels++;         
     }
     else
     {
        
     }      
  }//

  //create output variable for result
  mwSize  dims[] = {nvoxels,3};  
  plhs[0]        = mxCreateNumericArray( 2, dims, mxDOUBLE_CLASS, mxREAL );
  double *pr     = mxGetPr( plhs[0] );

  // second pass thru mask to get only values within the mask
  unsigned int tcount;
  for (mit.GoToBegin(); !mit.IsAtEnd(); ++mit)
  {
     index = mit.GetIndex();
     InputPointType iPoint;
     triangleMeshToImage->GetOutput()->TransformIndexToPhysicalPoint( index,
iPoint );  
     if (mit.Get() == m_MaskValue)
     {        
	 //mexPrintf("Point in  mesh: [%d %f %f
%f]\n",tcount,iPoint[0],iPoint[1],iPoint[2]); 
         pr[0*nvoxels + tcount] = (double)iPoint[0];
         pr[1*nvoxels + tcount] = (double)iPoint[1];
         pr[2*nvoxels + tcount] = (double)iPoint[2];
  
         tcount++;    
     }
     else
     {
        
     }      
  }//
  





Dan Mueller-2 wrote:
> 
> Hi Greg,
> 
> You could try rasterizing your mesh, then iterating over the image
> considering only values under the raster mask. Further, you could
> constrict the iteration to the bounding box of the raster mask. If
> your mesh consists of triangles, you can use
> itkTriangleMeshToImageFilter to create a binary image mask.
> 
> If rasterizing the mesh is already too slow, I'm not sure what else
> you could do...
> 
> Some more information regarding your specific problem might also help
> discover solutions "outside the box".
> 
> Hope this helps.
> 
> Regards, Dan
> 
> 2009/6/25 gregthom <gregthom99 at yahoo.com>:
>>
>> Greetings
>>
>> This is more of a design help question. I have an itkMesh and and
>> itkImage
>> in the same frame of reference. I need to get the image values of all
>> image
>> voxels(pixels) bounded by the itkMesh. Anybody have any thoughts on how
>> to
>> achieve this ? What is the fastest way to achieve this ?
>>
>> Thank you
>>
>>
>> greg
>> --
>> View this message in context:
>> http://www.nabble.com/fastest-way-to-get-image-voxels%28pixels%29-within-an-itkMesh---tp24199751p24199751.html
>> Sent from the ITK - Users mailing list archive at Nabble.com.
>>
>> _____________________________________
>> Powered by www.kitware.com
>>
>> Visit other Kitware open-source projects at
>> http://www.kitware.com/opensource/opensource.html
>>
>> Please keep messages on-topic and check the ITK FAQ at:
>> http://www.itk.org/Wiki/ITK_FAQ
>>
>> Follow this link to subscribe/unsubscribe:
>> http://www.itk.org/mailman/listinfo/insight-users
>>
> _____________________________________
> Powered by www.kitware.com
> 
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
> 
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
> 
> Follow this link to subscribe/unsubscribe:
> http://www.itk.org/mailman/listinfo/insight-users
> 
> 

-- 
View this message in context: http://www.nabble.com/fastest-way-to-get-image-voxels%28pixels%29-within-an-itkMesh---tp24199751p24336657.html
Sent from the ITK - Users mailing list archive at Nabble.com.