[Insight-users] fastest way to get image voxels(pixels) within an itkMesh ?
Luis Ibanez
luis.ibanez at kitware.com
Sat Jul 4 16:00:32 EDT 2009
Hi Gregthom,
So you have:
A) a Mesh
B) a binary image that you get from rasterizing (A)
C) a grayscale image
You want to get the grayscales values of (C) that correspond
to positions of (B) that are inside of (A).
Can we assume that (B) and (C) share the same image grid ?
That is,
Do the images B and C have the same origin, spacing, orientation
and number of pixels ?
If yes,
then you could do this faster by using the MaskImageFilter,
and passing (B) and (C) as its inputs.
If no,
You will have to use code similar to what is used in the
ImageMetrics, when collecting the values of the moving image.
In pseudo code:
Run Mesh rasterization filter
Create interpolator and connect it to (C)
For all pixels in the mask (B)
if pixel is ON
convert index to physical point
Evaluate interpolator at point
It is not clear from your email,
if what you want at the end is an image in the grid of (B) with
the values of (C), or whether you want to simply have a list
(e.g. std::vector<> with the list of values from (C).
Please clarify what kind of final output you want/need.
Thanks
Luis
----------------
gregthom wrote:
> 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
>>
>>
>
>
More information about the Insight-users
mailing list