[Insight-users] (no subject)

[Dženan Zukić]

More efficient would be to construct a closed (capped) polygonal cylinder,
transform it into binary mask (
http://www.itk.org/Doxygen/html/classitk_1_1TriangleMeshToBinaryImageFilter.html
or
its VTK equivalent) and then iterate that mask (some subregion, like you
already do) and check mask's intensity. 0=outside, 1=inside, 0<x<1=>edge

Here is an example conversion:

////voxelization using ITK implementation (bugged)
//VisualizingImageType::Pointer getMask()
//{
//    MeshType::Pointer mesh=makeCylinder();
//    typedef
itk::TriangleMeshToBinaryImageFilter<MeshType,VisualizingImageType>
MeshFilterType;
//    MeshFilterType::Pointer meshFilter = MeshFilterType::New();
//    meshFilter->SetInfoImage(visualizingImage);
//    meshFilter->SetInput(mesh);
//    meshFilter->Update();
//    mask=meshFilter->GetOutput();
//    return mask;
//}

////voxelization using VTK implementation (workaround)
VisualizingImageType::Pointer getMaskFromImageInformation(
                                    VisualizingImageType::RegionType region,
                                    VisualizingImageType::SpacingType
spacing,
                                    VisualizingImageType::PointType origin,
                                    VisualizingImageType::DirectionType
direction)
{
    VisualizingImageType::Pointer whiteITK=VisualizingImageType::New();
    whiteITK->SetRegions(region);
    whiteITK->Allocate();
    whiteITK->FillBuffer(1);

    typedef itk::ImageToVTKImageFilter<VisualizingImageType>
itkVtkConverter;
    itkVtkConverter::Pointer conv=itkVtkConverter::New();
    conv->SetInput(whiteITK);
    conv->Update();
    vtkSmartPointer<vtkImageData> whiteVTK=conv->GetOutput();

    vtkSmartPointer<vtkPolyData> poly=makeCylinder();

    vtkSmartPointer<vtkMatrix4x4> Mitk=vtkSmartPointer<vtkMatrix4x4>::New();
    for (int i=0; i<3; i++)
        for (int k=0; k<3; k++)
            Mitk->SetElement(i,k, direction(i,k));

    vtkSmartPointer<vtkMatrix4x4> Ms=vtkSmartPointer<vtkMatrix4x4>::New();
    Ms->SetElement(0,0,spacing[0]);
    Ms->SetElement(1,1,spacing[1]);
    Ms->SetElement(2,2,spacing[2]);

    vtkMatrix4x4::Multiply4x4(Mitk,Ms,Mitk);
    for (int i=0; i<3; i++)
        Mitk->SetElement(i,3, origin[i]);
    Mitk->Invert();

    vtkSmartPointer<vtkTransform> t=vtkSmartPointer<vtkTransform>::New();
    t->SetMatrix(Mitk);
    vtkSmartPointer<vtkTransformFilter>
tf=vtkSmartPointer<vtkTransformFilter>::New();
    tf->SetInputData(poly);
    tf->SetTransform(t);
    tf->Update();
    poly->SetPoints(tf->GetOutput()->GetPoints());

    vtkSmartPointer<vtkPolyDataToImageStencil> voxelizer =
        vtkSmartPointer<vtkPolyDataToImageStencil>::New();
    voxelizer->SetInputData(poly);
    voxelizer->SetOutputWholeExtent(whiteVTK->GetExtent());
    voxelizer->Update();

    vtkImageStencil *stencil=vtkImageStencil::New(); //crashes with smart
pointer
    //vtkSmartPointer<vtkImageStencil>
stencil=vtkSmartPointer<vtkImageStencil>::New();
    stencil->SetInputData(conv->GetOutput());
    stencil->SetStencilConnection(voxelizer->GetOutputPort());
    stencil->SetBackgroundValue(0);
    stencil->Update();
    whiteVTK=stencil->GetOutput();

    typedef itk::VTKImageToImageFilter<VisualizingImageType>
vtk2itkConverter;
    vtk2itkConverter::Pointer conv2=vtk2itkConverter::New();
    conv2->SetInput(whiteVTK);
    conv2->Update();
    mask=conv2->GetOutput();
    mask->DisconnectPipeline();
    mask->SetDirection(direction);
    mask->SetSpacing(spacing);
    mask->SetOrigin(origin);
    return mask;
}


On Thu, Aug 1, 2013 at 12:09 PM, Arindam Bhattacharya
<arindamb86 at gmail.com>wrote:

> hello all,
> Given a 3D  scalar volume.
> A particular grid point say P
> A direction say V1 ( a unit vector)
> A direction say V2 ( an orthogonal  unit vector to v1)
>
> I wish to find all grid points which are within a 'cylinder' defined by
> the point P the direction V1
> a length L and a breadth B along v2
>
> I have added a link
> <https://dl.dropboxusercontent.com/u/2989703/Note%20Aug%201%2C%202013.pdf> to
> show how it might look in  2D
>
> i am currently doing this by first finding all the neighbors of the point
> P in a region RxRxR
> and going through each point and finding the projection on the normal
> given by V1
> and the orthogonal plane given by v2 which is very slow,
> is there a faster approach ?
>
> thanks
> --
> Arindam Bhattacharya
> Graduate Student
>
>
> _____________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Kitware offers ITK Training Courses, for more information visit:
> http://www.kitware.com/products/protraining.php
>
> 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
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20130801/31187fc3/attachment.htm>