Authors belong to:
- University of LYON http://www.universite-lyon.fr/
- - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr
+ - Léon Bérard cancer center http://www.centreleonberard.fr
- CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
This software is distributed WITHOUT ANY WARRANTY; without even
- BSD See included LICENSE.txt file
- CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
-======================================================================-====*/
+===========================================================================**/
#ifndef clitkMotionMaskGenericFilter_txx
#define clitkMotionMaskGenericFilter_txx
*
===================================================*/
+#include "itkLabelImageToShapeLabelMapFilter.h"
+#include "itkShapeLabelObject.h"
namespace clitk
{
if (m_Verbose) std::cout<<"Reading the air feature image..."<<std::endl;
featureReader->Update();
air=featureReader->GetOutput();
+
+ //---------------------------------
+ // Pad
+ //---------------------------------
+ if(m_ArgsInfo.pad_flag) {
+ typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
+ IteratorType it(air, air->GetLargestPossibleRegion());
+ typename InternalImageType::IndexType index;
+ while(!it.IsAtEnd()) {
+ index=it.GetIndex();
+
+ //Pad borders of all dimensions but 2 (cranio-caudal)
+ for (unsigned int i=0; i<Dimension; i++){
+ if(i==2)
+ continue;
+ if(index[i]==air->GetLargestPossibleRegion().GetIndex()[i]
+ || index[i]==(unsigned int)air->GetLargestPossibleRegion().GetIndex()[i]+ (unsigned int) air->GetLargestPossibleRegion().GetSize()[i]-1)
+ it.Set(1);
+ }
+ ++it;
+ }
+ }
} else {
if (m_Verbose) std::cout<<"Extracting the air feature image..."<<std::endl;
//---------------------------------
if (m_Verbose) std::cout<<"Binarizing the image using thresholds "<<m_ArgsInfo.lowerThresholdAir_arg
<<", "<<m_ArgsInfo.upperThresholdAir_arg<<"..."<<std::endl;
binarizeFilter->Update();
+ air = binarizeFilter->GetOutput();
+
+ //---------------------------------
+ // Pad
+ //---------------------------------
+ if(m_ArgsInfo.pad_flag) {
+ typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
+ IteratorType it(air, air->GetLargestPossibleRegion());
+ typename InternalImageType::IndexType index;
+ while(!it.IsAtEnd()) {
+ index=it.GetIndex();
+
+ //Pad borders of all dimensions but 2 (cranio-caudal)
+ for (unsigned int i=0; i<Dimension; i++){
+ if(i==2)
+ continue;
+ if(index[i]==air->GetLargestPossibleRegion().GetIndex()[i]
+ || index[i]==(unsigned int)air->GetLargestPossibleRegion().GetIndex()[i]+ (unsigned int) air->GetLargestPossibleRegion().GetSize()[i]-1)
+ it.Set(binarizeFilter->GetInsideValue());
+ }
+ ++it;
+ }
+ }
//---------------------------------
// Remove lungs (in place)
//---------------------------------
typedef itk::ImageRegionIterator<InternalImageType> IteratorType;
- IteratorType itAir(binarizeFilter->GetOutput(), binarizeFilter->GetOutput()->GetLargestPossibleRegion());
+ IteratorType itAir(air, binarizeFilter->GetOutput()->GetLargestPossibleRegion());
IteratorType itLungs(lungs, binarizeFilter->GetOutput()->GetLargestPossibleRegion()); //lungs padded, used air region
itAir.GoToBegin();
itLungs.GoToBegin();
if (m_Verbose) std::cout<<"Mirroring the entire image along the CC axis..."<<std::endl;
mirrorPadImageFilter->Update();
air=mirrorPadImageFilter->GetOutput();
- writeImage<InternalImageType>(air,"/home/srit/tmp/air.mhd");
-
- //---------------------------------
- // Pad
- //---------------------------------
- if(m_ArgsInfo.pad_flag) {
- typedef itk::ImageRegionIteratorWithIndex<InternalImageType> IteratorType;
- IteratorType it(air, air->GetLargestPossibleRegion());
- typename InternalImageType::IndexType index;
- while(!it.IsAtEnd()) {
- index=it.GetIndex();
- for (unsigned int i=0; i<Dimension; i++)
- if(index[i]==air->GetLargestPossibleRegion().GetIndex()[i]
- || index[i]==(unsigned int)air->GetLargestPossibleRegion().GetIndex()[i]+ (unsigned int) air->GetLargestPossibleRegion().GetSize()[i]-1)
- it.Set(0);
- ++it;
- }
- }
+ //writeImage<InternalImageType>(air,"/home/srit/tmp/air.mhd");
return air;
}
connectFilter->SetBackgroundValue(0);
connectFilter->SetFullyConnected(true);
if (m_Verbose) std::cout<<"Labeling the connected components..."<<std::endl;
+ connectFilter->Update();
+ if (m_Verbose) std::cout<<"found "<< connectFilter->GetObjectCount() << std::endl;
//---------------------------------
// Sort the labels according to size
typename RelabelFilterType::Pointer relabelFilter=RelabelFilterType::New();
relabelFilter->SetInput(connectFilter->GetOutput());
if (m_Verbose) std::cout<<"Sorting the labels..."<<std::endl;
- // writeImage<InternalImageType> (relabelFilter->GetOutput(), "/home/jef/tmp/labels.mhd");
+ // writeImage<InternalImageType> (relabelFilter->GetOutput(), "/home/vdelmon/tmp/labels.mhd");
//---------------------------------
// Keep the label
//-------------------------------------------------------------------
template <unsigned int Dimension, class PixelType>
typename itk::Image<MotionMaskGenericFilter::InternalPixelType, Dimension>::Pointer
-MotionMaskGenericFilter::InitializeEllips( typename itk::Vector<double,Dimension> center, typename itk::Image<InternalPixelType,Dimension>::Pointer bones_low )
+MotionMaskGenericFilter::InitializeEllips( typename itk::Vector<double,Dimension> center, typename itk::Image<InternalPixelType,Dimension>::Pointer bones_low, typename itk::Image<InternalPixelType,Dimension>::Pointer lungs_low )
{
// ImageTypes
std::cout<<"=========================================="<<std::endl;
}
- //---------------------------------
- // Offset from center
- //---------------------------------
- typename itk::Vector<double,Dimension> offset;
- if(m_ArgsInfo.offset_given== Dimension) {
- for(unsigned int i=0; i<Dimension; i++)
- offset[i]=m_ArgsInfo.offset_arg[i];
- } else {
- offset.Fill(0.);
- offset[1]=-50;
+ itk::Vector<double,Dimension> centerEllips;
+ typename itk::Vector<double, Dimension> axes;
+ if (m_ArgsInfo.ellipseAutoDetect_flag) {
+ if(m_Verbose) {
+ std::cout<<"Auto-detecting intial ellipse..."<<std::endl;
+ }
+
+ // compute statistics of the (mirroed) lung region in order to guess the params of the ellipse
+ typedef itk::LabelImageToShapeLabelMapFilter<InternalImageType> LabelImageToShapeLabelMapFilter;
+ typename LabelImageToShapeLabelMapFilter::Pointer label_filter = LabelImageToShapeLabelMapFilter::New();
+ label_filter->SetInput(lungs_low);
+ label_filter->Update();
+
+ typename InternalImageType::SpacingType spacing = lungs_low->GetSpacing();
+ typedef typename LabelImageToShapeLabelMapFilter::OutputImageType::LabelObjectType LabelType;
+ const LabelType* label = dynamic_cast<LabelType*>(label_filter->GetOutput()->GetNthLabelObject(0));
+
+
+ // try to guess ideal ellipse axes from the lungs' bounding box
+#if ITK_VERSION_MAJOR >= 4
+ typename LabelType::RegionType lung_bbox = label->GetBoundingBox();
+#else
+ typename LabelType::RegionType lung_bbox = label->GetRegion();
+#endif
+ axes[0] = 0.95*lung_bbox.GetSize()[0]*spacing[0]/2;
+ axes[1] = 0.3*lung_bbox.GetSize()[1]*spacing[1]/2;
+ axes[2] = 0.95*lung_bbox.GetSize()[2]*spacing[2]/2;
+
+ // ellipse's center is the center of the lungs' bounding box
+ typename InternalImageType::PointType origin = lungs_low->GetOrigin();
+ centerEllips[0] = origin[0] + (lung_bbox.GetIndex()[0] + lung_bbox.GetSize()[0]/2)*spacing[0];
+ centerEllips[1] = origin[1] + (lung_bbox.GetIndex()[1] + lung_bbox.GetSize()[1]/2)*spacing[1];
+ centerEllips[2] = origin[2] + (lung_bbox.GetIndex()[2] + lung_bbox.GetSize()[2]/2)*spacing[2];
+
+ if(m_Verbose) {
+ std::cout << "Ellipse center at " << centerEllips << std::endl;
+ std::cout << "Ellipse axes as " << axes << std::endl;
+ std::cout << "Lung's bounding box (index,size) " << lung_bbox.GetIndex() << lung_bbox.GetSize() << std::endl;
+ }
}
- itk::Vector<double,Dimension> centerEllips=center+offset;
- if (m_Verbose) {
- std::cout<<"Placing the center of the initial ellipsoide at (mm) "<<std::endl;
- std::cout<<centerEllips[0];
- for (unsigned int i=1; i<Dimension; i++)
- std::cout<<" "<<centerEllips[i];
- std::cout<<std::endl;
+ else {
+ //---------------------------------
+ // Offset from center
+ //---------------------------------
+ typename itk::Vector<double,Dimension> offset;
+ if(m_ArgsInfo.offset_given== Dimension) {
+ for(unsigned int i=0; i<Dimension; i++)
+ offset[i]=m_ArgsInfo.offset_arg[i];
+ } else {
+ offset.Fill(0.);
+ offset[1]=-50;
+ }
+ centerEllips=center+offset;
+ if (m_Verbose) {
+ std::cout<<"Placing the center of the initial ellipsoide at (mm) "<<std::endl;
+ std::cout<<centerEllips[0];
+ for (unsigned int i=1; i<Dimension; i++)
+ std::cout<<" "<<centerEllips[i];
+ std::cout<<std::endl;
+ }
+
+ // Axes
+ if (m_ArgsInfo.axes_given==Dimension)
+ for(unsigned int i=0; i<Dimension; i++)
+ axes[i]=m_ArgsInfo.axes_arg[i];
+ else {
+ axes[0]=100;
+ axes[1]=30;
+ axes[2]=150;
+ }
}
//---------------------------------
ellips->Allocate();
ellips->FillBuffer(0);
- // Axes
- typename itk::Vector<double, Dimension> axes;
- if (m_ArgsInfo.axes_given==Dimension)
- for(unsigned int i=0; i<Dimension; i++)
- axes[i]=m_ArgsInfo.axes_arg[i];
- else {
- axes[0]=100;
- axes[1]=30;
- axes[2]=150;
- }
-
// Draw an ellips
IteratorType itEllips(ellips, ellips->GetLargestPossibleRegion());
itEllips.GoToBegin();
//----------------------------------------------------------------------------------------------------
// Ellipsoide initialization
//----------------------------------------------------------------------------------------------------
- typename InternalImageType::Pointer m_Ellips= InitializeEllips<Dimension, PixelType>(center,bones_low);
+ typename InternalImageType::Pointer m_Ellips= InitializeEllips<Dimension, PixelType>(center,bones_low,lungs_low);
//----------------------------------------------------------------------------------------------------
// Grow ellips
distanceMapImageFilter->SetInsideIsPositive(false);
if (m_Verbose) std::cout<<"Calculating the distance map..."<<std::endl;
distanceMapImageFilter->Update();
+ if (m_ArgsInfo.writeDistMap_given) {
+ writeImage<LevelSetImageType>(distanceMapImageFilter->GetOutput(), m_ArgsInfo.writeDistMap_arg, m_Verbose);
+
+ }
//---------------------------------
// Grow while monitoring detection point
distanceMapImageFilter->SetInsideIsPositive(false);
if (m_Verbose) std::cout<<"Calculating distance map..."<<std::endl;
distanceMapImageFilter->Update();
+ if (m_ArgsInfo.writeDistMap_given) {
+ writeImage<LevelSetImageType>(distanceMapImageFilter->GetOutput(), m_ArgsInfo.writeDistMap_arg, m_Verbose);
+
+ }
//---------------------------------
// Grow while monitoring lung volume coverage
git://git.creatis.insa-lyon.fr / clitk.git / blobdiff