[Insight-users] fastest way to get image voxels(pixels) within an itkMesh ?
Luis Ibanez
luis.ibanez at kitware.com
Sun Jul 5 18:10:47 EDT 2009
Hi Gregthom,
If the image (B) and (C) are not using the same
grid then you need to use the interpolator.
The code that you need is the following (a mutation of
the GetValue() method from MeanSquaresImageToImageMetric)
Hereafter:
FixedImage == image (B)
MovingImage == image (C)
{
FixedImageConstPointer fixedImage = this->m_FixedImage;
FixedImagePointer outputImage = FixedImageType::New();
outputImage->CopyInformation( fixedImage );
outputImage->SetRegions( fixedImage->GetBufferedRegion() );
outputImage->Allocate();
outputImage->FillBuffer( 0 );
typedef itk::ImageRegionConstIteratorWithIndex<
FixedImageType> FixedIteratorType;
typedef itk::LinearInterpolateImageFunction< MovingImageType >
InterpolatorType;
InterpolatorType::Pointer interpolator = InterpolatorType::New();
interpolator->SetInputImage( movingImage );
FixedIteratorType ti( fixedImage, fixedImage->GetBufferedRegion() );
FixedIteratorType to( outputImage, outputImage->GetBufferedRegion());
ti.GoToBegin();
t0.GoToBegin();
InputPointType inputPoint;
while( !ti.IsAtEnd() )
{
fixedImage->TransformIndexToPhysicalPoint(
ti.GetIndex(), inputPoint );
if( interpolator->IsInsideBuffer( inputPoint ) )
{
const RealType movingValue =
interpolator->Evaluate( inputPoint );
to.Set( movingValue );
}
++ti;
++to;
}
}
Regards,
Luis
------------------
gregthom wrote:
> Hi Luis,
>
> Thank you for taking some time to look into this.
>
> Yes you describe the problem accurately.
>
>
>>>Can we assume that (B) and (C) share the same image grid ?
>
>
> Answer: not necessarily.
>
>
>>>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).
>
>
> I need the in the end an image in the grid of (B) with the values of C as
> you thought.
> Can you suggest how to obtain this efficiently ? I will start to work on the
> algorithm in the
> pseudocode, in fact I had already been doing something similar.
>
>
> Hope to hear from you
> best wishes
>
>
>
>
>
>
>
> Luis Ibanez wrote:
>
>>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
>>>>
>>>>
>>>
>>>
>>_____________________________________
>>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