[ITK Community] Isolating 3d touching regions

[Emiliano Pastorelli]

Hi guys,

i played around with your suggestions, but i'm still in deep water (and
sinking).

By following the Lehmann-Beare paper The watershed transform in ITK -
discussion
and new developments i tried to process my dataset (this is a subset of the
full dataset https://www.dropbox.com/s/h4g3dhx0s9e0k8k/FibreSeparator.tar.gz,
the one i load in my software is raw while this is a mhd+raw, but
otherwise it's just part of it), but i get stuck even before reaching the
watershed part.

Until after the OpeningByReconstructionImageFilter filter, everything goes
reasonably well, i still have to give the binary lower threshold manually
(and i'd like to avoid that in a final version, but at the moment i can
stand it). Then it comes to the invert-distance-invert and the data goes
wild. When i try to visualize in paraview the obtained image, i have no
signs of fibres anymore, but only a blue uniform cube. I tried the
signedDanielsson distance and in that case my cube was colours with 3d
rings of colors.

I get the suggestion that Dan gave to get the maxima of distances of a
region to obtain two markers, but due to my still huge ignorance of ITK,
how would you suggest to proceed to obtain those? at the moment i feel
pretty lost even on where to start to get that.

Thanks again,
Best,
Emiliano


#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <list>
#include <ctime>
#include <vector>
#include <string>
#include <fstream>
#include <math.h>

#include "vtkImageReader.h"
#include "vtkImageData.h"
#include "vtkSmartPointer.h"
#include "vtkExtractVOI.h"

#include "itkImage.h"
#include "itkImageFileWriter.h"
#include "itkImageFileReader.h"
#include "itkMedianImageFilter.h"
#include "itkBinaryThresholdImageFilter.h"

#include "itkVTKImageToImageFilter.h"
#include "itkImageToVTKImageFilter.h"

#include "itkHessian3DToVesselnessMeasureImageFilter.h"
#include "itkHessianRecursiveGaussianImageFilter.h"
#include "itkBinaryBallStructuringElement.h"

#include <itkOpeningByReconstructionImageFilter.h>
#include "itkDanielssonDistanceMapImageFilter.h"
#include "itkInvertIntensityImageFilter.h"

using namespace std;
using namespace itk;

const unsigned int Dimension = 3;
typedef short InputPixelType;
typedef short OutputPixelType;
typedef itk::Image< InputPixelType, Dimension > InputImageType;
typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
typedef itk::Image< long unsigned int, Dimension > LabelImageType;

typedef VTKImageToImageFilter<OutputImageType> VTKImageToImageType;

void loadVolume(string filename, int maxX, int maxY, int maxZ,
vtkImageReader *imgReader){
    imgReader->SetFileDimensionality(3);
    imgReader->SetFileName(filename.c_str());
    imgReader->SetNumberOfScalarComponents(1);
    imgReader->SetDataExtent(0,maxX-1, 0,maxY-1, 0,maxZ-1);
    imgReader->SetDataByteOrderToLittleEndian();
    imgReader->Update();
}

void convertVTKtoITK(vtkImageData *image, VTKImageToImageType::Pointer
convert){
    vtkImageData *myImage = (vtkImageData*)image;
    convert->SetInput(myImage);
    convert->Update();
}

int main(int argc, char* argv[]){

    string fileName;
    fileName = argv[1];

    vtkImageReader *reader = vtkImageReader::New();
    VTKImageToImageType::Pointer convert = VTKImageToImageType::New();

    /*
     * LOAD THE GIVEN RAW VOLUME IN VTK
     */
    loadVolume(fileName,800,800,576,reader);

    /*
     * EXTRACT SUBVOLUME
     */
    vtkSmartPointer<vtkExtractVOI> subvolumeA =
vtkSmartPointer<vtkExtractVOI>::New();
    subvolumeA->SetInputData(reader->GetOutput(0));
    subvolumeA->SetVOI(100,450,100,450,100,450);
    subvolumeA->Update();
    cout << "Subvolume extracted" << endl;

    /*
     * CONVERT FROM VTK TO ITK FILE FORMAT
     */
    convert->SetInput(subvolumeA->GetOutput());
    convert->Update();
    cout << "Converted" << endl;

    typedef itk::HessianRecursiveGaussianImageFilter< OutputImageType >
    HessianFilterType;
    HessianFilterType::Pointer hessianFilter = HessianFilterType::New();
    hessianFilter->SetInput( convert->GetOutput() );
    hessianFilter->SetSigma(4.0);
    hessianFilter->Update();

    typedef itk::Hessian3DToVesselnessMeasureImageFilter< OutputPixelType >
    VesselnessMeasureFilterType;
    VesselnessMeasureFilterType::Pointer vesselnessFilter =
VesselnessMeasureFilterType::New();
    vesselnessFilter->SetInput( hessianFilter->GetOutput() );
    vesselnessFilter->SetAlpha1(0.5);
    vesselnessFilter->SetAlpha2(2.0);
    vesselnessFilter->Update();
    cout << "Frangi vesselness extracted" << endl;

    typedef itk::MedianImageFilter< OutputImageType,OutputImageType >
FilterType;
    FilterType::Pointer medianFilter = FilterType::New();
    FilterType::InputSizeType radius;
    radius.Fill(2.5);
    medianFilter->SetRadius(radius);
    medianFilter->SetInput( vesselnessFilter->GetOutput() );
    medianFilter->Update();
    cout << "Median Filtered" << endl;

    typedef itk::BinaryThresholdImageFilter <OutputImageType,
OutputImageType> BinaryThresholdImageFilterType;
    BinaryThresholdImageFilterType::Pointer thresholdFilter =
BinaryThresholdImageFilterType::New();
    thresholdFilter->SetInput(medianFilter->GetOutput());
    thresholdFilter->SetLowerThreshold(3.0f);
    thresholdFilter->SetUpperThreshold(255.0f);
    thresholdFilter->SetInsideValue(255);
    thresholdFilter->SetOutsideValue(0);
    thresholdFilter->Update();
    cout << "Binary Threshold" << endl;

    typedef BinaryBallStructuringElement<short, 3 > StructuringElementType;
    StructuringElementType structuringElement;
    structuringElement.SetRadius(3.0);
    structuringElement.CreateStructuringElement();

    typedef
OpeningByReconstructionImageFilter<OutputImageType,OutputImageType,StructuringElementType>
OpeningFilter;
    OpeningFilter::Pointer newOpeningFilter = OpeningFilter::New();
    newOpeningFilter->SetInput(thresholdFilter->GetOutput());
    newOpeningFilter->SetKernel(structuringElement);
    newOpeningFilter->Update();
    cout << "Opening by reconstruction" << endl;

    typedef InvertIntensityImageFilter <OutputImageType>
InvertIntensityImageFilterType;
    InvertIntensityImageFilterType::Pointer invertIntensityFilter =
InvertIntensityImageFilterType::New();
    invertIntensityFilter->SetInput(newOpeningFilter->GetOutput());
//    invertIntensityFilter->SetMaximum(255);
    invertIntensityFilter->Update();

    typedef DanielssonDistanceMapImageFilter<OutputImageType,
OutputImageType>  DistanceFilterType;
    DistanceFilterType::Pointer sddm = DistanceFilterType::New();
//    sddm->SetSquaredDistance(true);
    sddm->SetInput(invertIntensityFilter->GetOutput());
    sddm->Update();

    invertIntensityFilter->SetInput(sddm->GetOutput());
//    invertIntensityFilter->SetMaximum(255);
    invertIntensityFilter->Update();


    /*
     * SAVE TO FILE
     */
    typedef itk::ImageFileWriter< OutputImageType > WriterType;
    WriterType::Pointer writer = WriterType::New();
    writer->SetInput( invertIntensityFilter->GetOutput() );
    writer->SetFileName( "testAnalyzed.mhd" );

    try{
        writer->Update();
    }
    catch( itk::ExceptionObject & error ){
        std::cerr << "Error: " << error << std::endl;
        return EXIT_FAILURE;
    }

    cout << "Volume Saved" << endl;

}

2014-03-14 13:39 GMT+02:00 Dan Mueller <dan.muel at gmail.com>:

> Hi,
>
> You might try using markers computed from the regional maxima of a
> distance map (where inside is positive)...
>
> Sharing your dataset may also help community members suggest other options.
>
> HTH
>
> Cheers, Dan
> On 14 Mar 2014 19:38, "Emiliano Pastorelli" <emiliano.pastorelli at gmail.com>
> wrote:
>
>> Hi guys,
>>
>> I started taking a look at the watershed and might indeed do the work.
>> Due to my inexperience though, finding the marker seems a less trivial
>> problem than what Dan mentions as not so difficult.
>>
>> The problem is that somehow, i was trying to avoid skeletonization (i
>> will need later to compare my results to a similar dataset analyzed in a
>> different way that is centered on skeletonization). Do you happen to have
>> any idea on how to get those markers otherwise?
>>
>> And should I go for the skeletonization (should it be the only solution,
>> i might have to) wouldn't the skeleton of the touching fibres be detected a
>> single skeleton? Also, what ITK algorithm would you best suggest should it
>> be needed to skeletonize a dataset like this one?
>>
>> Thanks and sorry for the possible naivety of some of my questions!
>>
>> Best,
>> Emiliano
>>
>>
>> 2014-03-14 9:52 GMT+02:00 Emiliano Pastorelli <
>> emiliano.pastorelli at gmail.com>:
>>
>>> Hi guys,
>>>
>>> thanks, there's quite some material to chew on there, i'll take a look
>>> at it today in the office. The whole process for me needs to be somehow
>>> automatized from the beginning to the end to achieve what i need, so the
>>> user will only have to specify the fibres size and length and no manual
>>> thresholding.
>>>
>>> Thanks, i'll get back to you if I have further problems!
>>>
>>> Best,
>>> Emiliano
>>>
>>> Il 14/03/14 03:13, Bradley Lowekamp ha scritto:
>>>
>>>  Hello,
>>>>
>>>> You may also find this SimpleITK IPython Notebook which used the this
>>>> suggested watershed approach to separate some coins:
>>>> http://simpleitk.github.io/SimpleITK-Notebooks/32_
>>>> Watersheds_Segmentation.html
>>>>
>>>> Also SimpleITK from Python may be a useful way to easily explore the
>>>> many algorithms and infinite number of combinations available from ITK.
>>>>
>>>> Brad
>>>>
>>>> On Mar 13, 2014, at 5:37 PM, Dan Mueller <dan.muel at gmail.com> wrote:
>>>>
>>>>  Hi Emiliano,
>>>>>
>>>>> One possibility might be to use the mathematical watershed algorithm
>>>>> to binary segment the fibres and separate touching labels.
>>>>>
>>>>> See here:
>>>>> http://www.insight-journal.org/browse/publication/92
>>>>>
>>>>> and here:
>>>>> http://www.vincent-net.com/luc/papers/96semstats_morpho_topics.pdf(pg 44)
>>>>>
>>>>> This would essentially reduce your problem to find a good "marker" for
>>>>> each fibre. In your case this should not be too difficult, as the
>>>>> skeleton of each fibre should suffice (i.e. the skeletons will not
>>>>> overlap).
>>>>>
>>>>> Good luck.
>>>>>
>>>>> Cheers, Dan
>>>>>
>>>>> On 13 March 2014 23:16, Emiliano Pastorelli
>>>>> <emiliano.pastorelli at gmail.com> wrote:
>>>>>
>>>>>> Hi All,
>>>>>>
>>>>>> i'm quite new to ITK, and i'm trying to manipulate a 3d volume image
>>>>>> representing some fibres floating in a block of concrete. Due to the
>>>>>> scanning inaccuracies, some of them actually seem to "touch", also
>>>>>> when they
>>>>>> wouldn't be touching in the real dataset.
>>>>>>
>>>>>> I tried to separate them in every way, through the Frangi Vesselness
>>>>>> filters, and another bunch of attempts, but i don't seem to find a
>>>>>> way to do
>>>>>> it reasonably quickly and smartly.
>>>>>>
>>>>>> http://www.kyb3.org/images/touchingFibres.png
>>>>>>
>>>>>> in the picture, exactly in the center you see an example of how the
>>>>>> touching
>>>>>> area usually looks like (the dataset is already divide in labeled
>>>>>> regions as
>>>>>> label map as well).
>>>>>> I am pretty sure that some filter can do that, but for example i was
>>>>>> taking
>>>>>> a look at the connectedthresholdimagefilter and i couldn't really
>>>>>> figure out
>>>>>> if that's what i need, and in that case how i shall use the seeds (my
>>>>>> software shall become flexible enough to deal with really similar in
>>>>>> shapes
>>>>>> but different in distribution, alignment and contact points datasets).
>>>>>>
>>>>>> Does anybody have any idea that might spare me another couple of
>>>>>> endless
>>>>>> days like the last ones?
>>>>>>
>>>>>> Thanks in advance,
>>>>>> Emiliano
>>>>>> _______________________________________________
>>>>>> Community mailing list
>>>>>> Community at itk.org
>>>>>> http://public.kitware.com/cgi-bin/mailman/listinfo/community
>>>>>>
>>>>> _______________________________________________
>>>>> Community mailing list
>>>>> Community at itk.org
>>>>> http://public.kitware.com/cgi-bin/mailman/listinfo/community
>>>>>
>>>>
>>>
>>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/community/attachments/20140317/bf8d81cf/attachment-0002.html>