<div dir="ltr"><div class="gmail_extra"><div><div><div><div><div><div><div>Hi guys,<br><br></div>i played around with your suggestions, but i'm still in deep water (and sinking).<br><br></div>By following the Lehmann-Beare paper The watershed transform in ITK - discussion<br>
and new developments i tried to process my dataset (this is a subset of the full dataset <a href="https://www.dropbox.com/s/h4g3dhx0s9e0k8k/FibreSeparator.tar.gz" target="_blank">https://www.dropbox.com/s/h4g3dhx0s9e0k8k/FibreSeparator.tar.gz</a>
, 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.<br><br></div>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.<br><br></div>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.<br><br></div>Thanks again,<br></div>Best,<br></div>Emiliano<br><br><br>#include <iostream><br>#include <stdio.h><br>#include <stdlib.h><br>#include <list><br>#include <ctime><br>#include <vector><br>
#include <string><br>#include <fstream><br>#include <math.h><br><br>#include "vtkImageReader.h"<br>#include "vtkImageData.h"<br>#include "vtkSmartPointer.h"<br>#include "vtkExtractVOI.h"<br>
<br>#include "itkImage.h"<br>#include "itkImageFileWriter.h"<br>#include "itkImageFileReader.h"<br>#include "itkMedianImageFilter.h"<br>#include "itkBinaryThresholdImageFilter.h"<br>
<br>#include "itkVTKImageToImageFilter.h"<br>#include "itkImageToVTKImageFilter.h"<br><br>#include "itkHessian3DToVesselnessMeasureImageFilter.h"<br>#include "itkHessianRecursiveGaussianImageFilter.h"<br>
#include "itkBinaryBallStructuringElement.h"<br><br>#include <itkOpeningByReconstructionImageFilter.h><br>#include "itkDanielssonDistanceMapImageFilter.h"<br>#include "itkInvertIntensityImageFilter.h"<br>
<br>using namespace std;<br>using namespace itk;<br><br>const unsigned int Dimension = 3;<br>typedef short InputPixelType;<br>typedef short OutputPixelType;<br>typedef itk::Image< InputPixelType, Dimension > InputImageType;<br>
typedef itk::Image< OutputPixelType, Dimension > OutputImageType;<br>typedef itk::Image< long unsigned int, Dimension > LabelImageType;<br><br>typedef VTKImageToImageFilter<OutputImageType> VTKImageToImageType;<br>
<br>void loadVolume(string filename, int maxX, int maxY, int maxZ, vtkImageReader *imgReader){<br> imgReader->SetFileDimensionality(3);<br> imgReader->SetFileName(filename.c_str());<br> imgReader->SetNumberOfScalarComponents(1);<br>
imgReader->SetDataExtent(0,maxX-1, 0,maxY-1, 0,maxZ-1);<br> imgReader->SetDataByteOrderToLittleEndian();<br> imgReader->Update();<br>}<br><br>void convertVTKtoITK(vtkImageData *image, VTKImageToImageType::Pointer convert){<br>
vtkImageData *myImage = (vtkImageData*)image;<br> convert->SetInput(myImage);<br> convert->Update();<br>}<br><br>int main(int argc, char* argv[]){<br><br> string fileName;<br> fileName = argv[1];<br>
<br> vtkImageReader *reader = vtkImageReader::New();<br> VTKImageToImageType::Pointer convert = VTKImageToImageType::New();<br><br> /*<br> * LOAD THE GIVEN RAW VOLUME IN VTK<br> */<br> loadVolume(fileName,800,800,576,reader);<br>
<br> /*<br> * EXTRACT SUBVOLUME<br> */<br> vtkSmartPointer<vtkExtractVOI> subvolumeA = vtkSmartPointer<vtkExtractVOI>::New();<br> subvolumeA->SetInputData(reader->GetOutput(0));<br> subvolumeA->SetVOI(100,450,100,450,100,450);<br>
subvolumeA->Update();<br> cout << "Subvolume extracted" << endl;<br><br> /*<br> * CONVERT FROM VTK TO ITK FILE FORMAT<br> */<br> convert->SetInput(subvolumeA->GetOutput());<br>
convert->Update();<br> cout << "Converted" << endl;<br><br> typedef itk::HessianRecursiveGaussianImageFilter< OutputImageType ><br> HessianFilterType;<br> HessianFilterType::Pointer hessianFilter = HessianFilterType::New();<br>
hessianFilter->SetInput( convert->GetOutput() );<br> hessianFilter->SetSigma(4.0);<br> hessianFilter->Update();<br><br> typedef itk::Hessian3DToVesselnessMeasureImageFilter< OutputPixelType ><br>
VesselnessMeasureFilterType;<br> VesselnessMeasureFilterType::Pointer vesselnessFilter = VesselnessMeasureFilterType::New();<br> vesselnessFilter->SetInput( hessianFilter->GetOutput() );<br> vesselnessFilter->SetAlpha1(0.5);<br>
vesselnessFilter->SetAlpha2(2.0);<br> vesselnessFilter->Update();<br> cout << "Frangi vesselness extracted" << endl;<br><br> typedef itk::MedianImageFilter< OutputImageType,OutputImageType > FilterType;<br>
FilterType::Pointer medianFilter = FilterType::New();<br> FilterType::InputSizeType radius;<br> radius.Fill(2.5);<br> medianFilter->SetRadius(radius);<br> medianFilter->SetInput( vesselnessFilter->GetOutput() );<br>
medianFilter->Update();<br> cout << "Median Filtered" << endl;<br><br> typedef itk::BinaryThresholdImageFilter <OutputImageType, OutputImageType> BinaryThresholdImageFilterType;<br>
BinaryThresholdImageFilterType::Pointer thresholdFilter = BinaryThresholdImageFilterType::New();<br> thresholdFilter->SetInput(medianFilter->GetOutput());<br> thresholdFilter->SetLowerThreshold(3.0f);<br>
thresholdFilter->SetUpperThreshold(255.0f);<br> thresholdFilter->SetInsideValue(255);<br> thresholdFilter->SetOutsideValue(0);<br> thresholdFilter->Update();<br> cout << "Binary Threshold" << endl;<br>
<br> typedef BinaryBallStructuringElement<short, 3 > StructuringElementType;<br> StructuringElementType structuringElement;<br> structuringElement.SetRadius(3.0);<br> structuringElement.CreateStructuringElement();<br>
<br> typedef OpeningByReconstructionImageFilter<OutputImageType,OutputImageType,StructuringElementType> OpeningFilter;<br> OpeningFilter::Pointer newOpeningFilter = OpeningFilter::New();<br> newOpeningFilter->SetInput(thresholdFilter->GetOutput());<br>
newOpeningFilter->SetKernel(structuringElement);<br> newOpeningFilter->Update();<br> cout << "Opening by reconstruction" << endl;<br><br> typedef InvertIntensityImageFilter <OutputImageType> InvertIntensityImageFilterType;<br>
InvertIntensityImageFilterType::Pointer invertIntensityFilter = InvertIntensityImageFilterType::New();<br> invertIntensityFilter->SetInput(newOpeningFilter->GetOutput());<br>// invertIntensityFilter->SetMaximum(255);<br>
invertIntensityFilter->Update();<br><br> typedef DanielssonDistanceMapImageFilter<OutputImageType, OutputImageType> DistanceFilterType;<br> DistanceFilterType::Pointer sddm = DistanceFilterType::New();<br>
// sddm->SetSquaredDistance(true);<br> sddm->SetInput(invertIntensityFilter->GetOutput());<br> sddm->Update();<br><br> invertIntensityFilter->SetInput(sddm->GetOutput());<br>// invertIntensityFilter->SetMaximum(255);<br>
invertIntensityFilter->Update();<br><br><br> /*<br> * SAVE TO FILE<br> */<br> typedef itk::ImageFileWriter< OutputImageType > WriterType;<br> WriterType::Pointer writer = WriterType::New();<br>
writer->SetInput( invertIntensityFilter->GetOutput() );<br> writer->SetFileName( "testAnalyzed.mhd" );<br><br> try{<br> writer->Update();<br> }<br> catch( itk::ExceptionObject & error ){<br>
std::cerr << "Error: " << error << std::endl;<br> return EXIT_FAILURE;<br> }<br><br> cout << "Volume Saved" << endl;<br><br>}<br><br><div class="gmail_quote">
2014-03-14 13:39 GMT+02:00 Dan Mueller <span dir="ltr"><<a href="mailto:dan.muel@gmail.com" target="_blank">dan.muel@gmail.com</a>></span>:<br><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<p dir="ltr">Hi,</p>
<p dir="ltr">You might try using markers computed from the regional maxima of a distance map (where inside is positive)...</p>
<p dir="ltr">Sharing your dataset may also help community members suggest other options.</p>
<p dir="ltr">HTH</p>
<p dir="ltr">Cheers, Dan</p><div class=""><div class="h5">
<div class="gmail_quote">On 14 Mar 2014 19:38, "Emiliano Pastorelli" <<a href="mailto:emiliano.pastorelli@gmail.com" target="_blank">emiliano.pastorelli@gmail.com</a>> wrote:<br type="attribution"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div dir="ltr"><div><div><div><div>Hi guys,<br><br></div>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.<br>
<br></div>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?<br>
<br></div>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?<br>
<br></div><div>Thanks and sorry for the possible naivety of some of my questions!<br><br></div><div>Best,<br></div><div>Emiliano<br></div></div><div class="gmail_extra"><br><br><div class="gmail_quote">2014-03-14 9:52 GMT+02:00 Emiliano Pastorelli <span dir="ltr"><<a href="mailto:emiliano.pastorelli@gmail.com" target="_blank">emiliano.pastorelli@gmail.com</a>></span>:<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">Hi guys,<br>
<br>
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.<br>
<br>
Thanks, i'll get back to you if I have further problems!<br>
<br>
Best,<br>
Emiliano<br>
<br>
Il 14/03/14 03:13, Bradley Lowekamp ha scritto:<div><div><br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
Hello,<br>
<br>
You may also find this SimpleITK IPython Notebook which used the this suggested watershed approach to separate some coins:<br>
<a href="http://simpleitk.github.io/SimpleITK-Notebooks/32_Watersheds_Segmentation.html" target="_blank">http://simpleitk.github.io/<u></u>SimpleITK-Notebooks/32_<u></u>Watersheds_Segmentation.html</a><br>
<br>
Also SimpleITK from Python may be a useful way to easily explore the many algorithms and infinite number of combinations available from ITK.<br>
<br>
Brad<br>
<br>
On Mar 13, 2014, at 5:37 PM, Dan Mueller <<a href="mailto:dan.muel@gmail.com" target="_blank">dan.muel@gmail.com</a>> wrote:<br>
<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
Hi Emiliano,<br>
<br>
One possibility might be to use the mathematical watershed algorithm<br>
to binary segment the fibres and separate touching labels.<br>
<br>
See here:<br>
<a href="http://www.insight-journal.org/browse/publication/92" target="_blank">http://www.insight-journal.<u></u>org/browse/publication/92</a><br>
<br>
and here:<br>
<a href="http://www.vincent-net.com/luc/papers/96semstats_morpho_topics.pdf" target="_blank">http://www.vincent-net.com/<u></u>luc/papers/96semstats_morpho_<u></u>topics.pdf</a> (pg 44)<br>
<br>
This would essentially reduce your problem to find a good "marker" for<br>
each fibre. In your case this should not be too difficult, as the<br>
skeleton of each fibre should suffice (i.e. the skeletons will not<br>
overlap).<br>
<br>
Good luck.<br>
<br>
Cheers, Dan<br>
<br>
On 13 March 2014 23:16, Emiliano Pastorelli<br>
<<a href="mailto:emiliano.pastorelli@gmail.com" target="_blank">emiliano.pastorelli@gmail.com</a><u></u>> wrote:<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
Hi All,<br>
<br>
i'm quite new to ITK, and i'm trying to manipulate a 3d volume image<br>
representing some fibres floating in a block of concrete. Due to the<br>
scanning inaccuracies, some of them actually seem to "touch", also when they<br>
wouldn't be touching in the real dataset.<br>
<br>
I tried to separate them in every way, through the Frangi Vesselness<br>
filters, and another bunch of attempts, but i don't seem to find a way to do<br>
it reasonably quickly and smartly.<br>
<br>
<a href="http://www.kyb3.org/images/touchingFibres.png" target="_blank">http://www.kyb3.org/images/<u></u>touchingFibres.png</a><br>
<br>
in the picture, exactly in the center you see an example of how the touching<br>
area usually looks like (the dataset is already divide in labeled regions as<br>
label map as well).<br>
I am pretty sure that some filter can do that, but for example i was taking<br>
a look at the connectedthresholdimagefilter and i couldn't really figure out<br>
if that's what i need, and in that case how i shall use the seeds (my<br>
software shall become flexible enough to deal with really similar in shapes<br>
but different in distribution, alignment and contact points datasets).<br>
<br>
Does anybody have any idea that might spare me another couple of endless<br>
days like the last ones?<br>
<br>
Thanks in advance,<br>
Emiliano<br>
______________________________<u></u>_________________<br>
Community mailing list<br>
<a href="mailto:Community@itk.org" target="_blank">Community@itk.org</a><br>
<a href="http://public.kitware.com/cgi-bin/mailman/listinfo/community" target="_blank">http://public.kitware.com/cgi-<u></u>bin/mailman/listinfo/community</a><br>
</blockquote>
______________________________<u></u>_________________<br>
Community mailing list<br>
<a href="mailto:Community@itk.org" target="_blank">Community@itk.org</a><br>
<a href="http://public.kitware.com/cgi-bin/mailman/listinfo/community" target="_blank">http://public.kitware.com/cgi-<u></u>bin/mailman/listinfo/community</a><br>
</blockquote></blockquote>
<br>
</div></div></blockquote></div><br></div>
</blockquote></div>
</div></div></blockquote></div><br></div></div>