[Insight-users] fastest way to get image voxels(pixels) within an itkMesh ?
gregthom
gregthom99 at yahoo.com
Mon Jul 6 06:33:37 EDT 2009
Thanks Luis
I am making progress in my project. !
Best wishes
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
>>
>>
>
>
--
View this message in context: http://www.nabble.com/fastest-way-to-get-image-voxels%28pixels%29-within-an-itkMesh---tp24199751p24353129.html
Sent from the ITK - Users mailing list archive at Nabble.com.
More information about the Insight-users
mailing list