From: dsarrut Date: Thu, 21 Apr 2011 05:05:39 +0000 (+0200) Subject: merge cvs -> git X-Git-Tag: v1.2.0~24^2~7 X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=34f45c6507a3605351e265f4d5eb5b4bb31a1fa6;p=clitk.git merge cvs -> git --- diff --git a/.gitignore b/.gitignore index 4b17a88..8453ef6 100644 --- a/.gitignore +++ b/.gitignore @@ -13,3 +13,5 @@ CMakeCache.txt build/* */.vimrc _* +#*# +notes.org diff --git a/itk/clitkMeshToBinaryImageFilter.h b/itk/clitkMeshToBinaryImageFilter.h new file mode 100644 index 0000000..b681a98 --- /dev/null +++ b/itk/clitkMeshToBinaryImageFilter.h @@ -0,0 +1,94 @@ +/*========================================================================= + Program: vv http://www.creatis.insa-lyon.fr/rio/vv + + Authors belong to: + - University of LYON http://www.universite-lyon.fr/ + - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr + - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the copyright notices for more information. + + It is distributed under dual licence + + - BSD See included LICENSE.txt file + - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html + ======================================================================-====*/ + +#ifndef CLITKMESHTOBINARYIMAGEFILTER_H +#define CLITKMESHTOBINARYIMAGEFILTER_H + +// clitk +#include "clitkCommon.h" + +namespace clitk { + + /* -------------------------------------------------------------------- + Convert a mesh composed of a list of 2D contours, into a 3D binray + image. + -------------------------------------------------------------------- */ + + template + class ITK_EXPORT MeshToBinaryImageFilter: public itk::ImageSource + { + + public: + /** Standard class typedefs. */ + typedef itk::ImageSource Superclass; + typedef MeshToBinaryImageFilter Self; + typedef itk::SmartPointer Pointer; + typedef itk::SmartPointer ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(MeshToBinaryImageFilter, Superclass); + + /** ImageDimension constants */ + itkStaticConstMacro(ImageDimension, unsigned int, ImageType::ImageDimension); + typedef itk::Image FloatImageType; + + /** Some convenient typedefs. */ + typedef typename ImageType::ConstPointer ImageConstPointer; + typedef typename ImageType::Pointer ImagePointer; + typedef typename ImageType::RegionType RegionType; + typedef typename ImageType::PixelType PixelType; + typedef typename ImageType::SpacingType SpacingType; + typedef typename ImageType::SizeType SizeType; + typedef typename ImageType::PointType PointType; + typedef itk::Image SliceType; + + /** Input : initial image and object */ + itkSetMacro(Mesh, vtkSmartPointer); + itkGetConstMacro(Mesh, vtkSmartPointer); + + itkSetMacro(LikeImage, ImagePointer); + itkGetConstMacro(LikeImage, ImagePointer); + + protected: + MeshToBinaryImageFilter(); + virtual ~MeshToBinaryImageFilter() {} + + virtual void GenerateOutputInformation(); + virtual void GenerateData(); + + ImagePointer m_LikeImage; + vtkSmartPointer m_Mesh; + + private: + MeshToBinaryImageFilter(const Self&); //purposely not implemented + void operator=(const Self&); //purposely not implemented + + }; // end class + //-------------------------------------------------------------------- + +} // end namespace clitk +//-------------------------------------------------------------------- + +#ifndef ITK_MANUAL_INSTANTIATION +#include "clitkMeshToBinaryImageFilter.txx" +#endif + +#endif diff --git a/itk/clitkMeshToBinaryImageFilter.txx b/itk/clitkMeshToBinaryImageFilter.txx new file mode 100644 index 0000000..db98970 --- /dev/null +++ b/itk/clitkMeshToBinaryImageFilter.txx @@ -0,0 +1,144 @@ +/*========================================================================= + Program: vv http://www.creatis.insa-lyon.fr/rio/vv + + Authors belong to: + - University of LYON http://www.universite-lyon.fr/ + - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr + - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the copyright notices for more information. + + It is distributed under dual licence + + - BSD See included LICENSE.txt file + - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html + ======================================================================-====*/ + +#include +#include +#include +#include + +#include "itkVTKImageImport.h" +#include "vtkImageExport.h" +#include "vtkImageData.h" + + + +//-------------------------------------------------------------------- +template +clitk::MeshToBinaryImageFilter:: +MeshToBinaryImageFilter():itk::ImageSource() +{ +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::MeshToBinaryImageFilter:: +GenerateOutputInformation() +{ + ImagePointer output = this->GetOutput(0); + + // Set the region to output + typename ImageType::RegionType m_Region = m_LikeImage->GetLargestPossibleRegion(); + typename ImageType::SizeType size = m_Region.GetSize(); + size[0]++; + size[1]++; + size[2]++; + m_Region.SetSize(size); + output->SetLargestPossibleRegion(m_Region); + output->SetRequestedRegion(m_Region); + output->SetBufferedRegion(m_Region); + output->SetRegions(m_Region); + output->Allocate(); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::MeshToBinaryImageFilter:: +GenerateData() +{ + // GO + vtkSmartPointer binary_image=vtkSmartPointer::New(); + binary_image->SetScalarTypeToUnsignedChar(); + + // Set spacing + PointType samp_origin = m_LikeImage->GetOrigin(); + SpacingType spacing=m_LikeImage->GetSpacing(); + double * spacing2 = new double[3]; + spacing2[0] = spacing[0]; + spacing2[1] = spacing[1]; + spacing2[2] = spacing[2]; + binary_image->SetSpacing(spacing2); + + // Set origin + /// Put the origin on a voxel to avoid small skips + binary_image->SetOrigin(samp_origin[0], samp_origin[1], samp_origin[2]); + + // Compute image bounds + binary_image->SetExtent(0,m_LikeImage->GetLargestPossibleRegion().GetSize()[0], + 0,m_LikeImage->GetLargestPossibleRegion().GetSize()[1], + 0,m_LikeImage->GetLargestPossibleRegion().GetSize()[2] + ); + + // Allocate data + binary_image->AllocateScalars(); + memset(binary_image->GetScalarPointer(),0, + binary_image->GetDimensions()[0]*binary_image->GetDimensions()[1]*binary_image->GetDimensions()[2]*sizeof(unsigned char)); + + // Image stencil + vtkSmartPointer sts=vtkSmartPointer::New(); + //The following line is extremely important + //http://www.nabble.com/Bug-in-vtkPolyDataToImageStencil--td23368312.html#a23370933 + sts->SetTolerance(0); + sts->SetInformationInput(binary_image); + + // Extrusion + vtkSmartPointer extrude=vtkSmartPointer::New(); + extrude->SetInput(m_Mesh); + // We extrude in the -slice_spacing direction to respect the FOCAL convention + extrude->SetVector(0, 0, -m_LikeImage->GetSpacing()[2]); + sts->SetInput(extrude->GetOutput()); + + // Stencil + vtkSmartPointer stencil=vtkSmartPointer::New(); + stencil->SetStencil(sts->GetOutput()); + stencil->SetInput(binary_image); + + // Convert VTK to ITK + vtkImageExport * m_Exporter = vtkImageExport::New(); + typedef itk::VTKImageImport< ImageType > ImporterFilterType; + typename ImporterFilterType::Pointer m_Importer = ImporterFilterType::New(); + + m_Importer->SetUpdateInformationCallback( m_Exporter->GetUpdateInformationCallback()); + m_Importer->SetPipelineModifiedCallback( m_Exporter->GetPipelineModifiedCallback()); + m_Importer->SetWholeExtentCallback( m_Exporter->GetWholeExtentCallback()); + m_Importer->SetSpacingCallback( m_Exporter->GetSpacingCallback()); + m_Importer->SetOriginCallback( m_Exporter->GetOriginCallback()); + m_Importer->SetScalarTypeCallback( m_Exporter->GetScalarTypeCallback()); + m_Importer->SetNumberOfComponentsCallback( m_Exporter->GetNumberOfComponentsCallback()); + m_Importer->SetPropagateUpdateExtentCallback( m_Exporter->GetPropagateUpdateExtentCallback()); + m_Importer->SetUpdateDataCallback( m_Exporter->GetUpdateDataCallback()); + m_Importer->SetDataExtentCallback( m_Exporter->GetDataExtentCallback()); + m_Importer->SetBufferPointerCallback( m_Exporter->GetBufferPointerCallback()); + m_Importer->SetCallbackUserData( m_Exporter->GetCallbackUserData()); + + m_Exporter->SetInput( stencil->GetOutput() ); + m_Importer->Update(); + + writeImage(m_Importer->GetOutput(), "f.mhd"); + + this->SetNthOutput(0, m_Importer->GetOutput()); + + return; +} +//-------------------------------------------------------------------- + diff --git a/itk/clitkSegmentationUtils.h b/itk/clitkSegmentationUtils.h index fcc7f28..5aa9b9a 100644 --- a/itk/clitkSegmentationUtils.h +++ b/itk/clitkSegmentationUtils.h @@ -3,7 +3,7 @@ Authors belong to: - University of LYON http://www.universite-lyon.fr/ - - Léon Bérard cancer center http://www.centreleonberard.fr + - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr This software is distributed WITHOUT ANY WARRANTY; without even @@ -14,7 +14,7 @@ - BSD See included LICENSE.txt file - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html - ===========================================================================**/ + ======================================================================-====*/ #ifndef CLITKSEGMENTATIONUTILS_H #define CLITKSEGMENTATIONUTILS_H @@ -209,6 +209,11 @@ namespace clitk { ComputeCentroids(const ImageType * image, typename ImageType::PixelType BG, std::vector & centroids); + template + void + ComputeCentroids2(const ImageType * image, + typename ImageType::PixelType BG, + std::vector & centroids); //-------------------------------------------------------------------- @@ -258,14 +263,20 @@ namespace clitk { typedef std::map MapPoint2DType; typedef std::vector VectorPoint3DType; + typedef std::vector VectorPoint2DType; + public: static void Convert2DTo3D(const PointType2D & p2D, const ImageType * image, const int slice, PointType3D & p3D); - static void Convert2DTo3DList(const MapPoint2DType & map, + static void Convert2DMapTo3DList(const MapPoint2DType & map, const ImageType * image, VectorPoint3DType & list); + static void Convert2DListTo3DList(const VectorPoint2DType & p, + int slice, + const ImageType * image, + VectorPoint3DType & list); }; //-------------------------------------------------------------------- diff --git a/itk/clitkSegmentationUtils.txx b/itk/clitkSegmentationUtils.txx index d2e1758..7acec20 100644 --- a/itk/clitkSegmentationUtils.txx +++ b/itk/clitkSegmentationUtils.txx @@ -14,7 +14,7 @@ - BSD See included LICENSE.txt file - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html - ===========================================================================**/ + ======================================================================-====*/ // clitk #include "clitkSetBackgroundImageFilter.h" @@ -516,6 +516,47 @@ namespace clitk { //-------------------------------------------------------------------- + //-------------------------------------------------------------------- + template + void + ComputeCentroids2(const ImageType * image, + typename ImageType::PixelType BG, + std::vector & centroids) + { + typedef long LabelType; + static const unsigned int Dim = ImageType::ImageDimension; + typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType; + typedef itk::LabelMap< LabelObjectType > LabelMapType; + typedef itk::LabelImageToLabelMapFilter ImageToMapFilterType; + typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New(); + typedef itk::ShapeLabelMapFilter ShapeFilterType; + typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New(); + imageToLabelFilter->SetBackgroundValue(BG); + imageToLabelFilter->SetInput(image); + statFilter->SetInput(imageToLabelFilter->GetOutput()); + statFilter->Update(); + typename LabelMapType::Pointer labelMap = statFilter->GetOutput(); + + centroids.clear(); + typename ImageType::PointType dummy; + centroids.push_back(dummy); // label 0 -> no centroid, use dummy point + for(uint i=1; iGetNumberOfLabelObjects()+1; i++) { + centroids.push_back(labelMap->GetLabelObject(i)->GetCentroid()); + } + + for(uint i=1; iGetNumberOfLabelObjects()+1; i++) { + DD(labelMap->GetLabelObject(i)->GetBinaryPrincipalAxes()); + DD(labelMap->GetLabelObject(i)->GetBinaryFlatness()); + DD(labelMap->GetLabelObject(i)->GetRoundness ()); + + // search for the point on the boundary alog PA + + } + + } + //-------------------------------------------------------------------- + + //-------------------------------------------------------------------- template void @@ -558,7 +599,7 @@ namespace clitk { //-------------------------------------------------------------------- template void - PointsUtils::Convert2DTo3DList(const MapPoint2DType & map, + PointsUtils::Convert2DMapTo3DList(const MapPoint2DType & map, const ImageType * image, VectorPoint3DType & list) { @@ -572,6 +613,24 @@ namespace clitk { } //-------------------------------------------------------------------- + + //-------------------------------------------------------------------- + template + void + PointsUtils::Convert2DListTo3DList(const VectorPoint2DType & p2D, + int slice, + const ImageType * image, + VectorPoint3DType & list) + { + for(uint i=0; i void @@ -940,7 +999,7 @@ namespace clitk { } // Convert 2D points in slice into 3D points - clitk::PointsUtils::Convert2DTo3DList(position2D, input, A); + clitk::PointsUtils::Convert2DMapTo3DList(position2D, input, A); // Create additional point just right to the previous ones, on the // given lineDirection, in order to create a horizontal/vertical line. diff --git a/itk/clitkSliceBySliceRelativePositionFilter.h b/itk/clitkSliceBySliceRelativePositionFilter.h index 35b3b8d..53c5a87 100644 --- a/itk/clitkSliceBySliceRelativePositionFilter.h +++ b/itk/clitkSliceBySliceRelativePositionFilter.h @@ -3,7 +3,7 @@ Authors belong to: - University of LYON http://www.universite-lyon.fr/ - - Léon Bérard cancer center http://www.centreleonberard.fr + - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr This software is distributed WITHOUT ANY WARRANTY; without even @@ -14,7 +14,7 @@ - BSD See included LICENSE.txt file - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html - ===========================================================================**/ + ======================================================================-====*/ #ifndef CLITKSLICEBYSLICERELATIVEPOSITIONFILTER_H #define CLITKSLICEBYSLICERELATIVEPOSITIONFILTER_H @@ -86,6 +86,17 @@ namespace clitk { itkSetMacro(UseASingleObjectConnectedComponentBySliceFlag, bool); itkBooleanMacro(UseASingleObjectConnectedComponentBySliceFlag); + itkGetConstMacro(CCLSelectionFlag, bool); + itkSetMacro(CCLSelectionFlag, bool); + itkBooleanMacro(CCLSelectionFlag); + itkGetConstMacro(CCLSelectionDimension, int); + itkSetMacro(CCLSelectionDimension, int); + itkGetConstMacro(CCLSelectionDirection, int); + itkSetMacro(CCLSelectionDirection, int); + itkGetConstMacro(CCLSelectionIgnoreSingleCCLFlag, bool); + itkSetMacro(CCLSelectionIgnoreSingleCCLFlag, bool); + itkBooleanMacro(CCLSelectionIgnoreSingleCCLFlag); + protected: SliceBySliceRelativePositionFilter(); virtual ~SliceBySliceRelativePositionFilter() {} @@ -102,6 +113,10 @@ namespace clitk { int m_Direction; bool m_IgnoreEmptySliceObjectFlag; bool m_UseASingleObjectConnectedComponentBySliceFlag; + bool m_CCLSelectionFlag; + int m_CCLSelectionDimension; + int m_CCLSelectionDirection; + bool m_CCLSelectionIgnoreSingleCCLFlag; private: SliceBySliceRelativePositionFilter(const Self&); //purposely not implemented diff --git a/itk/clitkSliceBySliceRelativePositionFilter.txx b/itk/clitkSliceBySliceRelativePositionFilter.txx index 9d355bf..76b0fec 100644 --- a/itk/clitkSliceBySliceRelativePositionFilter.txx +++ b/itk/clitkSliceBySliceRelativePositionFilter.txx @@ -3,7 +3,7 @@ Authors belong to: - University of LYON http://www.universite-lyon.fr/ - - Léon Bérard cancer center http://www.centreleonberard.fr + - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr This software is distributed WITHOUT ANY WARRANTY; without even @@ -14,7 +14,7 @@ - BSD See included LICENSE.txt file - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html - ===========================================================================**/ + ======================================================================-====*/ // clitk #include "clitkSegmentationUtils.h" @@ -37,6 +37,10 @@ SliceBySliceRelativePositionFilter(): this->VerboseStepFlagOff(); this->WriteStepFlagOff(); this->SetCombineWithOrFlag(false); + CCLSelectionFlagOn(); + SetCCLSelectionDimension(0); + SetCCLSelectionDirection(1); + CCLSelectionIgnoreSingleCCLFlagOff(); } //-------------------------------------------------------------------- @@ -187,6 +191,39 @@ GenerateOutputInformation() mObjectSlices[i] = KeepLabels(mObjectSlices[i], 0, 1, 1, 1, true); } + // Select a single according to a position if more than one CCL + if (GetCCLSelectionFlag()) { + // if several CCL, choose the most extrema according a direction, + // if not -> should we consider this slice ? + if (nb<2) { + if (GetCCLSelectionIgnoreSingleCCLFlag()) { + mObjectSlices[i] = SetBackground(mObjectSlices[i], mObjectSlices[i], + 1, this->GetBackgroundValue(), + true); + } + } + int dim = GetCCLSelectionDimension(); + int direction = GetCCLSelectionDirection(); + std::vector centroids; + ComputeCentroids(mObjectSlices[i], this->GetBackgroundValue(), centroids); + uint index=1; + for(uint j=1; j centroids[index][dim]) index = j; + } + else { + if (centroids[j][dim] < centroids[index][dim]) index = j; + } + } + for(uint v=1; v(mObjectSlices[i], mObjectSlices[i], + (char)v, this->GetBackgroundValue(), + true); + } + } + } + // Relative position typedef clitk::AddRelativePositionConstraintToLabelImageFilter RelPosFilterType; typename RelPosFilterType::Pointer relPosFilter = RelPosFilterType::New(); @@ -216,6 +253,15 @@ GenerateOutputInformation() mInputSlices[i] = KeepLabels(mInputSlices[i], 0, 1, 1, 1, true); } + /* + // Select unique CC according to the most in a given direction + if (GetUniqueConnectedComponentBySliceAccordingToADirection()) { + int nb; + mInputSlices[i] = LabelizeAndCountNumberOfObjects(mInputSlices[i], 0, true, 1, nb); + std::vector & centroids; + ComputeCentroids + } + */ } } diff --git a/itk/itkVTKImageToImageFilter.h b/itk/itkVTKImageToImageFilter.h index 0d1204f..7d0eb97 100644 --- a/itk/itkVTKImageToImageFilter.h +++ b/itk/itkVTKImageToImageFilter.h @@ -3,7 +3,7 @@ Authors belong to: - University of LYON http://www.universite-lyon.fr/ - - Léon Bérard cancer center http://www.centreleonberard.fr + - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr This software is distributed WITHOUT ANY WARRANTY; without even @@ -14,94 +14,95 @@ - BSD See included LICENSE.txt file - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html -===========================================================================**/ -#ifndef __itkVTKImageToImageFilter_h -#define __itkVTKImageToImageFilter_h -#include "itkVTKImageImport.h" -#include "vtkImageExport.h" -#include "vtkImageData.h" - -#ifndef vtkFloatingPointType -#define vtkFloatingPointType float -#endif - -namespace itk -{ - -/** \class VTKImageToImageFilter - * \brief Converts a VTK image into an ITK image and plugs a - * vtk data pipeline to an ITK datapipeline. - * - * This class puts together an itkVTKImageImporter and a vtkImageExporter. - * It takes care of the details related to the connection of ITK and VTK - * pipelines. The User will perceive this filter as an adaptor to which - * a vtkImage can be plugged as input and an itk::Image is produced as - * output. - * - * \ingroup ImageFilters - */ -template -class ITK_EXPORT VTKImageToImageFilter : public ProcessObject -{ -public: - /** Standard class typedefs. */ - typedef VTKImageToImageFilter Self; - typedef ProcessObject Superclass; - typedef SmartPointer Pointer; - typedef SmartPointer ConstPointer; - - /** Method for creation through the object factory. */ - itkNewMacro(Self); - - /** Run-time type information (and related methods). */ - itkTypeMacro(VTKImageToImageFilter, ProcessObject); - - /** Some typedefs. */ - typedef TOutputImage OutputImageType; - typedef typename OutputImageType::ConstPointer OutputImagePointer; - typedef VTKImageImport< OutputImageType > ImporterFilterType; - typedef typename ImporterFilterType::Pointer ImporterFilterPointer; - - /** Get the output in the form of a vtkImage. - This call is delegated to the internal vtkImageImporter filter */ - const OutputImageType * GetOutput() const; - - /** Set the input in the form of a vtkImageData */ - void SetInput( vtkImageData * ); - - /** Return the internal VTK image exporter filter. - This is intended to facilitate users the access - to methods in the exporter */ - vtkImageExport * GetExporter() const; - - /** Return the internal ITK image importer filter. - This is intended to facilitate users the access - to methods in the importer */ - ImporterFilterType * GetImporter() const; - - /** This call delegate the update to the importer */ - void Update(); - -protected: - VTKImageToImageFilter(); - virtual ~VTKImageToImageFilter(); - -private: - VTKImageToImageFilter(const Self&); //purposely not implemented - void operator=(const Self&); //purposely not implemented - - ImporterFilterPointer m_Importer; - vtkImageExport * m_Exporter; - -}; - -} // end namespace itk - -#ifndef ITK_MANUAL_INSTANTIATION -#include "itkVTKImageToImageFilter.txx" -#endif - -#endif - - - +======================================================================-====*/ +#ifndef __itkVTKImageToImageFilter_h +#define __itkVTKImageToImageFilter_h +#include "itkVTKImageImport.h" +#include "vtkImageExport.h" +#include "vtkImageData.h" + +#ifndef vtkFloatingPointType +#define vtkFloatingPointType float +#endif + +namespace itk +{ + +/** \class VTKImageToImageFilter + * \brief Converts a VTK image into an ITK image and plugs a + * vtk data pipeline to an ITK datapipeline. + * + * This class puts together an itkVTKImageImporter and a vtkImageExporter. + * It takes care of the details related to the connection of ITK and VTK + * pipelines. The User will perceive this filter as an adaptor to which + * a vtkImage can be plugged as input and an itk::Image is produced as + * output. + * + * \ingroup ImageFilters + */ +template +class ITK_EXPORT VTKImageToImageFilter : public ProcessObject +{ +public: + /** Standard class typedefs. */ + typedef VTKImageToImageFilter Self; + typedef ProcessObject Superclass; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(VTKImageToImageFilter, ProcessObject); + + /** Some typedefs. */ + typedef TOutputImage OutputImageType; + typedef typename OutputImageType::ConstPointer OutputImagePointer; + typedef VTKImageImport< OutputImageType > ImporterFilterType; + typedef typename ImporterFilterType::Pointer ImporterFilterPointer; + + /** Get the output in the form of a vtkImage. + This call is delegated to the internal vtkImageImporter filter */ + // const + OutputImageType * GetOutput() const; + + /** Set the input in the form of a vtkImageData */ + void SetInput( vtkImageData * ); + + /** Return the internal VTK image exporter filter. + This is intended to facilitate users the access + to methods in the exporter */ + vtkImageExport * GetExporter() const; + + /** Return the internal ITK image importer filter. + This is intended to facilitate users the access + to methods in the importer */ + ImporterFilterType * GetImporter() const; + + /** This call delegate the update to the importer */ + void Update(); + +protected: + VTKImageToImageFilter(); + virtual ~VTKImageToImageFilter(); + +private: + VTKImageToImageFilter(const Self&); //purposely not implemented + void operator=(const Self&); //purposely not implemented + + ImporterFilterPointer m_Importer; + vtkImageExport * m_Exporter; + +}; + +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkVTKImageToImageFilter.txx" +#endif + +#endif + + + diff --git a/itk/itkVTKImageToImageFilter.txx b/itk/itkVTKImageToImageFilter.txx index 27474df..27a6595 100644 --- a/itk/itkVTKImageToImageFilter.txx +++ b/itk/itkVTKImageToImageFilter.txx @@ -3,7 +3,7 @@ Authors belong to: - University of LYON http://www.universite-lyon.fr/ - - Léon Bérard cancer center http://www.centreleonberard.fr + - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr This software is distributed WITHOUT ANY WARRANTY; without even @@ -14,7 +14,7 @@ - BSD See included LICENSE.txt file - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html -===========================================================================**/ +======================================================================-====*/ #ifndef _itkVTKImageToImageFilter_txx #define _itkVTKImageToImageFilter_txx #include "itkVTKImageToImageFilter.h" @@ -86,7 +86,8 @@ VTKImageToImageFilter * Get an itk::Image as output */ template -const typename VTKImageToImageFilter::OutputImageType * +//const +typename VTKImageToImageFilter::OutputImageType * VTKImageToImageFilter ::GetOutput() const { diff --git a/notes.org~ b/notes.org~ new file mode 100644 index 0000000..ec58a0b --- /dev/null +++ b/notes.org~ @@ -0,0 +1,12 @@ + +* GIT emacs command +C-x v d +.gitignore -> contains ignored files +commit (local) +push (repository) +pull + +* Test in temp +git clone http://www.creatis.insa-lyon.fr/~dsarrut/clitk3.pub.git clitk3 +need itk4 (git) ? + diff --git a/segmentation/CMakeLists.txt b/segmentation/CMakeLists.txt index b3e22d2..75250cd 100644 --- a/segmentation/CMakeLists.txt +++ b/segmentation/CMakeLists.txt @@ -43,9 +43,9 @@ IF(CLITK_BUILD_SEGMENTATION) ADD_EXECUTABLE(clitkExtractMediastinum clitkExtractMediastinum.cxx ${clitkExtractMediastinum_GGO_C}) TARGET_LINK_LIBRARIES(clitkExtractMediastinum clitkCommon clitkSegmentationGgoLib ${ITK_LIBRARIES}) - # WRAP_GGO(clitkExtractLymphStations_GGO_C clitkExtractLymphStations.ggo) - # ADD_EXECUTABLE(clitkExtractLymphStations clitkExtractLymphStations.cxx clitkFilterWithAnatomicalFeatureDatabaseManagement.cxx ${clitkExtractLymphStations_GGO_C}) - # TARGET_LINK_LIBRARIES(clitkExtractLymphStations clitkSegmentationGgoLib clitkCommon ${ITK_LIBRARIES} ITKStatistics) + WRAP_GGO(clitkExtractLymphStations_GGO_C clitkExtractLymphStations.ggo) + ADD_EXECUTABLE(clitkExtractLymphStations clitkExtractLymphStations.cxx clitkFilterWithAnatomicalFeatureDatabaseManagement.cxx ${clitkExtractLymphStations_GGO_C}) + TARGET_LINK_LIBRARIES(clitkExtractLymphStations clitkSegmentationGgoLib clitkCommon ITKIO ITKStatistics vtkHybrid) WRAP_GGO(clitkMorphoMath_GGO_C clitkMorphoMath.ggo) ADD_EXECUTABLE(clitkMorphoMath clitkMorphoMath.cxx ${clitkMorphoMath_GGO_C}) @@ -83,6 +83,10 @@ IF(CLITK_BUILD_SEGMENTATION) ADD_EXECUTABLE(clitkFillImageRegion clitkFillImageRegion.cxx clitkFillImageRegionGenericFilter.cxx ${clitkFillImageRegion_GGO_C}) TARGET_LINK_LIBRARIES(clitkFillImageRegion clitkCommon ${ITK_LIBRARIES}) + WRAP_GGO(clitkTestFilter_GGO_C clitkTestFilter.ggo) + ADD_EXECUTABLE(clitkTestFilter clitkTestFilter.cxx ${clitkTestFilter_GGO_C}) + TARGET_LINK_LIBRARIES(clitkTestFilter clitkSegmentationGgoLib clitkCommon vtkHybrid ITKIO) + ENDIF(CLITK_BUILD_SEGMENTATION) # ADD_EXECUTABLE(ScalarImageMarkovRandomField1 ScalarImageMarkovRandomField1.cxx) diff --git a/segmentation/clitkExtractLymphStation_2RL.txx b/segmentation/clitkExtractLymphStation_2RL.txx new file mode 100644 index 0000000..baa5730 --- /dev/null +++ b/segmentation/clitkExtractLymphStation_2RL.txx @@ -0,0 +1,853 @@ + + +// vtk +#include +#include +#include + +// clitk +#include "clitkMeshToBinaryImageFilter.h" + +// itk +#include + +//-------------------------------------------------------------------- +template +class comparePointsX { +public: + bool operator() (PointType i, PointType j) { return (i[0] +class comparePointsWithAngle { +public: + bool operator() (PairType i, PairType j) { return (i.second < j.second); } +}; +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void HypercubeCorners(std::vector > & out) { + std::vector > previous; + HypercubeCorners(previous); + out.resize(previous.size()*2); + for(uint i=0; i p; + if (i +void HypercubeCorners<1>(std::vector > & out) { + out.resize(2); + out[0][0] = 0; + out[1][0] = 1; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void ComputeImageBoundariesCoordinates(typename ImageType::Pointer image, + std::vector & bounds) +{ + // Get image max/min coordinates + const uint dim=ImageType::ImageDimension; + typedef typename ImageType::PointType PointType; + typedef typename ImageType::IndexType IndexType; + PointType min_c, max_c; + IndexType min_i, max_i; + min_i = image->GetLargestPossibleRegion().GetIndex(); + for(uint i=0; iGetLargestPossibleRegion().GetSize()[i] + min_i[i]; + image->TransformIndexToPhysicalPoint(min_i, min_c); + image->TransformIndexToPhysicalPoint(max_i, max_c); + + // Get corners coordinates + HypercubeCorners(bounds); + for(uint i=0; i +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_2RL_SetDefaultValues() +{ + SetFuzzyThreshold("2RL", "CommonCarotidArtery", 0.7); + SetFuzzyThreshold("2RL", "BrachioCephalicTrunk", 0.7); + SetFuzzyThreshold("2RL", "BrachioCephalicVein", 0.3); + SetFuzzyThreshold("2RL", "Aorta", 0.7); + SetFuzzyThreshold("2RL", "SubclavianArteryRight", 0.5); + SetFuzzyThreshold("2RL", "SubclavianArteryLeft", 0.8); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_2RL_SI_Limits() +{ + // Apex of the chest or Sternum & Carina. + StartNewStep("[Station 2RL] Inf/Sup limits with Sternum and TopOfAorticArch/CaudalMarginOfLeftBrachiocephalicVein"); + + /* Rod says: "For the inferior border, unlike in Atlas – UM, there + is now a difference between 2R and 2L. 2R stops at the + intersection of the caudal margin of the innominate vein with the + trachea. 2L extends less inferiorly to the superior border of the + aortic arch." */ + + /* Get TopOfAorticArch and CaudalMarginOfLeftBrachiocephalicVein + - TopOfAorticArch -> can be obtain from Aorta -> most sup part. + - CaudalMarginOfLeftBrachiocephalicVein -> must inf part of BrachicephalicVein + */ + MaskImagePointer Aorta = GetAFDB()->template GetImage("Aorta"); + MaskImagePointType p = Aorta->GetOrigin(); // initialise to avoid warning + clitk::FindExtremaPointInAGivenDirection(Aorta, GetBackgroundValue(), 2, false, p); + double TopOfAorticArchZ=p[2]; + GetAFDB()->SetDouble("TopOfAorticArchZ", TopOfAorticArchZ); + + MaskImagePointer BrachioCephalicVein = GetAFDB()->template GetImage("BrachioCephalicVein"); + clitk::FindExtremaPointInAGivenDirection(BrachioCephalicVein, GetBackgroundValue(), 2, true, p); + double CaudalMarginOfLeftBrachiocephalicVeinZ=p[2]; + GetAFDB()->SetDouble("CaudalMarginOfLeftBrachiocephalicVeinZ", CaudalMarginOfLeftBrachiocephalicVeinZ); + + // First, cut on the most inferior part. Add one slice because this + // inf slice should not be included. + double inf = std::min(CaudalMarginOfLeftBrachiocephalicVeinZ, TopOfAorticArchZ) + m_Working_Support->GetSpacing()[2]; + + // Get Sternum and search for the upper position + MaskImagePointer Sternum = GetAFDB()->template GetImage("Sternum"); + + // Search most sup point + clitk::FindExtremaPointInAGivenDirection(Sternum, GetBackgroundValue(), 2, false, p); + double m_SternumZ = p[2]; + // +Sternum->GetSpacing()[2]; // One more slice, because it is below this point // NOT HERE ?? WHY DIFFERENT FROM 3A ?? + + //* Crop support : + m_Working_Support = + clitk::CropImageAlongOneAxis(m_Working_Support, 2, + inf, m_SternumZ, true, + GetBackgroundValue()); + + StopCurrentStep(m_Working_Support); + m_ListOfStations["2R"] = m_Working_Support; + m_ListOfStations["2L"] = m_Working_Support; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_2RL_Ant_Limits() +{ + // ----------------------------------------------------- + /* Rod says: "The anterior border, as with the Atlas – UM, is + posterior to the vessels (right subclavian vein, left + brachiocephalic vein, right brachiocephalic vein, left subclavian + artery, left common carotid artery and brachiocephalic trunk). + These vessels are not included in the nodal station. The anterior + border is drawn to the midpoint of the vessel and an imaginary + line joins the middle of these vessels. Between the vessels, + station 2 is in contact with station 3a." */ + + // ----------------------------------------------------- + StartNewStep("[Station 2RL] Ant limits with CommonCarotidArtery"); + + // Get CommonCarotidArtery + MaskImagePointer CommonCarotidArtery = GetAFDB()->template GetImage("CommonCarotidArtery"); + + // Remove Ant to CommonCarotidArtery + typedef SliceBySliceRelativePositionFilter SliceRelPosFilterType; + typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New(); + sliceRelPosFilter->VerboseStepFlagOff(); + sliceRelPosFilter->WriteStepFlagOff(); + sliceRelPosFilter->SetInput(m_Working_Support); + sliceRelPosFilter->SetInputObject(CommonCarotidArtery); + sliceRelPosFilter->SetDirection(2); + sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("2RL", "CommonCarotidArtery")); + sliceRelPosFilter->AddOrientationTypeString("NotAntTo"); + sliceRelPosFilter->IntermediateSpacingFlagOn(); + sliceRelPosFilter->SetIntermediateSpacing(2); + sliceRelPosFilter->UniqueConnectedComponentBySliceOff(); + sliceRelPosFilter->UseASingleObjectConnectedComponentBySliceFlagOff(); + sliceRelPosFilter->AutoCropFlagOn(); + sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn(); + sliceRelPosFilter->RemoveObjectFlagOff(); + sliceRelPosFilter->Update(); + m_Working_Support = sliceRelPosFilter->GetOutput(); + + // End + StopCurrentStep(m_Working_Support); + m_ListOfStations["2R"] = m_Working_Support; + m_ListOfStations["2L"] = m_Working_Support; + + // ----------------------------------------------------- + // Remove Ant to H line from the Ant most part of the + // CommonCarotidArtery until we reach the first slice of + // BrachioCephalicTrunk + StartNewStep("[Station 2RL] Ant limits with CommonCarotidArtery, H line"); + + // First, find the first slice of BrachioCephalicTrunk + MaskImagePointer BrachioCephalicTrunk = GetAFDB()->template GetImage("BrachioCephalicTrunk"); + MaskImagePointType p = BrachioCephalicTrunk->GetOrigin(); // initialise to avoid warning + clitk::FindExtremaPointInAGivenDirection(BrachioCephalicTrunk, GetBackgroundValue(), 2, false, p); + double TopOfBrachioCephalicTrunkZ=p[2] + BrachioCephalicTrunk->GetSpacing()[2]; // Add one slice + + // Remove CommonCarotidArtery below this point + CommonCarotidArtery = clitk::CropImageBelow(CommonCarotidArtery, 2, TopOfBrachioCephalicTrunkZ, true, GetBackgroundValue()); + + // Find most Ant points + std::vector ccaAntPositionA; + std::vector ccaAntPositionB; + clitk::SliceBySliceBuildLineSegmentAccordingToExtremaPosition(CommonCarotidArtery, + GetBackgroundValue(), 2, + 1, true, // Ant + 0, // Horizontal line + -3, // margin + ccaAntPositionA, + ccaAntPositionB); + // Cut ant to this line + clitk::SliceBySliceSetBackgroundFromLineSeparation(m_Working_Support, + ccaAntPositionA, + ccaAntPositionB, + GetBackgroundValue(), 1, 10); + + // End + StopCurrentStep(m_Working_Support); + m_ListOfStations["2R"] = m_Working_Support; + m_ListOfStations["2L"] = m_Working_Support; + + // ----------------------------------------------------- + // Ant limit with the BrachioCephalicTrunk + StartNewStep("[Station 2RL] Ant limits with BrachioCephalicTrunk line"); + + // Remove Ant to BrachioCephalicTrunk + m_Working_Support = + clitk::SliceBySliceRelativePosition(m_Working_Support, BrachioCephalicTrunk, 2, + GetFuzzyThreshold("2RL", "BrachioCephalicTrunk"), "NotAntTo", false, 2, true, false); + // End + StopCurrentStep(m_Working_Support); + m_ListOfStations["2R"] = m_Working_Support; + m_ListOfStations["2L"] = m_Working_Support; + + // ----------------------------------------------------- + // Ant limit with the BrachioCephalicTrunk H line + StartNewStep("[Station 2RL] Ant limits with BrachioCephalicTrunk, Horizontal line"); + + // Find most Ant points + std::vector bctAntPositionA; + std::vector bctAntPositionB; + clitk::SliceBySliceBuildLineSegmentAccordingToExtremaPosition(BrachioCephalicTrunk, + GetBackgroundValue(), 2, + 1, true, // Ant + 0, // Horizontal line + -1, // margin + bctAntPositionA, + bctAntPositionB); + // Cut ant to this line + clitk::SliceBySliceSetBackgroundFromLineSeparation(m_Working_Support, + bctAntPositionA, + bctAntPositionB, + GetBackgroundValue(), 1, 10); + // End + StopCurrentStep(m_Working_Support); + m_ListOfStations["2R"] = m_Working_Support; + m_ListOfStations["2L"] = m_Working_Support; + + // ----------------------------------------------------- + // Ant limit with the BrachioCephalicVein + StartNewStep("[Station 2RL] Ant limits with BrachioCephalicVein"); + + // Remove Ant to BrachioCephalicVein + MaskImagePointer BrachioCephalicVein = GetAFDB()->template GetImage("BrachioCephalicVein"); + m_Working_Support = + clitk::SliceBySliceRelativePosition(m_Working_Support, BrachioCephalicVein, 2, + GetFuzzyThreshold("2RL", "BrachioCephalicVein"), "NotAntTo", false, 2, true, false); + // End + StopCurrentStep(m_Working_Support); + m_ListOfStations["2R"] = m_Working_Support; + m_ListOfStations["2L"] = m_Working_Support; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_2RL_Ant_Limits2() +{ + // ----------------------------------------------------- + /* Rod says: "The anterior border, as with the Atlas – UM, is + posterior to the vessels (right subclavian vein, left + brachiocephalic vein, right brachiocephalic vein, left subclavian + artery, left common carotid artery and brachiocephalic trunk). + These vessels are not included in the nodal station. The anterior + border is drawn to the midpoint of the vessel and an imaginary + line joins the middle of these vessels. Between the vessels, + station 2 is in contact with station 3a." */ + + // ----------------------------------------------------- + StartNewStep("[Station 2RL] Ant limits with vessels centroids"); + + /* Here, we consider the vessels form a kind of anterior barrier. We + link all vessels centroids and remove what is post to it. + - select the list of structure + vessel1 = BrachioCephalicTrunk + vessel2 = BrachioCephalicVein (warning several CCL, keep most at Right) + vessel3 = CommonCarotidArtery + vessel4 = SubclavianArtery + other = Thyroid + other = Aorta + - crop images as needed + - slice by slice, choose the CCL and extract centroids + - slice by slice, sort according to polar angle wrt Trachea centroid. + - slice by slice, link centroids and close contour + - remove outside this contour + - merge with support + */ + + // Read structures + MaskImagePointer BrachioCephalicTrunk = GetAFDB()->template GetImage("BrachioCephalicTrunk"); + MaskImagePointer BrachioCephalicVein = GetAFDB()->template GetImage("BrachioCephalicVein"); + MaskImagePointer CommonCarotidArtery = GetAFDB()->template GetImage("CommonCarotidArtery"); + MaskImagePointer SubclavianArtery = GetAFDB()->template GetImage("SubclavianArtery"); + MaskImagePointer Thyroid = GetAFDB()->template GetImage("Thyroid"); + MaskImagePointer Aorta = GetAFDB()->template GetImage("Aorta"); + MaskImagePointer Trachea = GetAFDB()->template GetImage("Trachea"); + + // Resize all structures like support + BrachioCephalicTrunk = + clitk::ResizeImageLike(BrachioCephalicTrunk, m_Working_Support, GetBackgroundValue()); + CommonCarotidArtery = + clitk::ResizeImageLike(CommonCarotidArtery, m_Working_Support, GetBackgroundValue()); + SubclavianArtery = + clitk::ResizeImageLike(SubclavianArtery, m_Working_Support, GetBackgroundValue()); + Thyroid = + clitk::ResizeImageLike(Thyroid, m_Working_Support, GetBackgroundValue()); + Aorta = + clitk::ResizeImageLike(Aorta, m_Working_Support, GetBackgroundValue()); + BrachioCephalicVein = + clitk::ResizeImageLike(BrachioCephalicVein, m_Working_Support, GetBackgroundValue()); + Trachea = + clitk::ResizeImageLike(Trachea, m_Working_Support, GetBackgroundValue()); + + // Extract slices + std::vector slices_BrachioCephalicTrunk; + clitk::ExtractSlices(BrachioCephalicTrunk, 2, slices_BrachioCephalicTrunk); + std::vector slices_BrachioCephalicVein; + clitk::ExtractSlices(BrachioCephalicVein, 2, slices_BrachioCephalicVein); + std::vector slices_CommonCarotidArtery; + clitk::ExtractSlices(CommonCarotidArtery, 2, slices_CommonCarotidArtery); + std::vector slices_SubclavianArtery; + clitk::ExtractSlices(SubclavianArtery, 2, slices_SubclavianArtery); + std::vector slices_Thyroid; + clitk::ExtractSlices(Thyroid, 2, slices_Thyroid); + std::vector slices_Aorta; + clitk::ExtractSlices(Aorta, 2, slices_Aorta); + std::vector slices_Trachea; + clitk::ExtractSlices(Trachea, 2, slices_Trachea); + uint n= slices_BrachioCephalicTrunk.size(); + + // Get the boundaries of one slice + std::vector bounds; + ComputeImageBoundariesCoordinates(slices_BrachioCephalicTrunk[0], bounds); + + // For all slices, for all structures, find the centroid and build the contour + // List of 3D points (for debug) + std::vector p3D; + + vtkSmartPointer append = vtkSmartPointer::New(); + for(uint i=0; i(slices_CommonCarotidArtery[i], + GetBackgroundValue(), true, 1); + slices_SubclavianArtery[i] = Labelize(slices_SubclavianArtery[i], + GetBackgroundValue(), true, 1); + slices_BrachioCephalicTrunk[i] = Labelize(slices_BrachioCephalicTrunk[i], + GetBackgroundValue(), true, 1); + slices_BrachioCephalicVein[i] = Labelize(slices_BrachioCephalicVein[i], + GetBackgroundValue(), true, 1); + slices_Thyroid[i] = Labelize(slices_Thyroid[i], + GetBackgroundValue(), true, 1); + slices_Aorta[i] = Labelize(slices_Aorta[i], + GetBackgroundValue(), true, 1); + + // Search centroids + std::vector points2D; + std::vector centroids1; + std::vector centroids2; + std::vector centroids3; + std::vector centroids4; + std::vector centroids5; + std::vector centroids6; + ComputeCentroids(slices_CommonCarotidArtery[i], GetBackgroundValue(), centroids1); + ComputeCentroids(slices_SubclavianArtery[i], GetBackgroundValue(), centroids2); + ComputeCentroids(slices_BrachioCephalicTrunk[i], GetBackgroundValue(), centroids3); + ComputeCentroids(slices_Thyroid[i], GetBackgroundValue(), centroids4); + ComputeCentroids(slices_Aorta[i], GetBackgroundValue(), centroids5); + ComputeCentroids(slices_BrachioCephalicVein[i], GetBackgroundValue(), centroids6); + + // BrachioCephalicVein -> when it is separated into two CCL, we + // only consider the most at Right one + if (centroids6.size() > 2) { + if (centroids6[1][0] < centroids6[2][0]) centroids6.erase(centroids6.begin()+2); + else centroids6.erase(centroids6.begin()+1); + } + + // BrachioCephalicVein -> when SubclavianArtery has 2 CCL + // (BrachioCephalicTrunk is divided) -> forget BrachioCephalicVein + if ((centroids3.size() ==1) && (centroids2.size() > 2)) { + centroids6.clear(); + } + + for(uint j=1; j centroids_trachea; + ComputeCentroids(slices_Trachea[i], GetBackgroundValue(), centroids_trachea); + typedef std::pair PointAngleType; + std::vector angles; + for(uint j=0; j0) angle = atan(y/x); + if ((x<0) && (y>=0)) angle = atan(y/x)+M_PI; + if ((x<0) && (y<0)) angle = atan(y/x)-M_PI; + if (x==0) { + if (y>0) angle = M_PI/2.0; + if (y<0) angle = -M_PI/2.0; + if (y==0) angle = 0; + } + angle = clitk::rad2deg(angle); + // Angle is [-180;180] wrt the X axis. We change the X axis to + // be the vertical line, because we want to sort from Right to + // Left from Post to Ant. + if (angle>0) + angle = (270-angle); + if (angle<0) { + angle = -angle-90; + if (angle<0) angle = 360-angle; + } + PointAngleType a; + a.first = points2D[j]; + a.second = angle; + angles.push_back(a); + } + + // Do nothing if less than 2 points --> n + if (points2D.size() < 3) { //continue; + continue; + } + + // Sort points2D according to polar angles + std::sort(angles.begin(), angles.end(), comparePointsWithAngle()); + for(uint j=0; j toadd; + uint index = 0; + double dmax = 5; + while (indexdmax) { + + MaskSlicePointType b; + b[0] = a[0]+(c[0]-a[0])/2.0; + b[1] = a[1]+(c[1]-a[1])/2.0; + + // Compute distance to trachea centroid + MaskSlicePointType m = centroids_trachea[1]; + double da = m.EuclideanDistanceTo(a); + double db = m.EuclideanDistanceTo(b); + //double dc = m.EuclideanDistanceTo(c); + + // Mean distance, find point on the line from trachea centroid + // to b + double alpha = (da+db)/2.0; + MaskSlicePointType v; + double n = sqrt( pow(b[0]-m[0], 2) + pow(b[1]-m[1], 2)); + v[0] = (b[0]-m[0])/n; + v[1] = (b[1]-m[1])/n; + MaskSlicePointType r; + r[0] = m[0]+alpha*v[0]; + r[1] = m[1]+alpha*v[1]; + points2D.insert(points2D.begin()+index+1, r); + } + else { + index++; + } + } + // DDV(points2D, points2D.size()); + + // Add some points to close the contour + // - H line towards Right + MaskSlicePointType p = points2D[0]; + p[0] = bounds[0][0]; + points2D.insert(points2D.begin(), p); + // - corner Right/Post + p = bounds[0]; + points2D.insert(points2D.begin(), p); + // - H line towards Left + p = points2D.back(); + p[0] = bounds[2][0]; + points2D.push_back(p); + // - corner Right/Post + p = bounds[2]; + points2D.push_back(p); + // Close contour with the first point + points2D.push_back(points2D[0]); + // DDV(points2D, points2D.size()); + + // Build 3D points from the 2D points + std::vector points3D; + clitk::PointsUtils::Convert2DListTo3DList(points2D, i, m_Working_Support, points3D); + for(uint x=0; x mesh = Build3DMeshFrom2DContour(points3D); + append->AddInput(mesh); + } + + // DEBUG: write points3D + clitk::WriteListOfLandmarks(p3D, "vessels-centroids.txt"); + + // Build the final 3D mesh form the list 2D mesh + append->Update(); + vtkSmartPointer mesh = append->GetOutput(); + + // Debug, write the mesh + /* + vtkSmartPointer w = vtkSmartPointer::New(); + w->SetInput(mesh); + w->SetFileName("bidon.vtk"); + w->Write(); + */ + + // Compute a single binary 3D image from the list of contours + clitk::MeshToBinaryImageFilter::Pointer filter = + clitk::MeshToBinaryImageFilter::New(); + filter->SetMesh(mesh); + filter->SetLikeImage(m_Working_Support); + filter->Update(); + MaskImagePointer binarizedContour = filter->GetOutput(); + + // Inverse binary mask if needed. We test a point that we know must be in FG. If it is not, inverse + ImagePointType p = p3D[2]; // This is the first centroid of the first slice + p[1] += 50; // 50 mm Post from this point must be kept + ImageIndexType index; + binarizedContour->TransformPhysicalPointToIndex(p, index); + bool isInside = (binarizedContour->GetPixel(index) != GetBackgroundValue()); + + // remove from support + typedef clitk::BooleanOperatorLabelImageFilter BoolFilterType; + typename BoolFilterType::Pointer boolFilter = BoolFilterType::New(); + boolFilter->InPlaceOn(); + boolFilter->SetInput1(m_Working_Support); + boolFilter->SetInput2(binarizedContour); + boolFilter->SetBackgroundValue1(GetBackgroundValue()); + boolFilter->SetBackgroundValue2(GetBackgroundValue()); + if (isInside) + boolFilter->SetOperationType(BoolFilterType::And); + else + boolFilter->SetOperationType(BoolFilterType::AndNot); + boolFilter->Update(); + m_Working_Support = boolFilter->GetOutput(); + + // End + StopCurrentStep(m_Working_Support); + m_ListOfStations["2R"] = m_Working_Support; + m_ListOfStations["2L"] = m_Working_Support; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_2RL_Post_Limits() +{ + StartNewStep("[Station 2RL] Post limits with post wall of Trachea"); + + // Get Trachea + MaskImagePointer Trachea = GetAFDB()->template GetImage("Trachea"); + + // Resize like the current support (to have the same number of slices) + Trachea = clitk::ResizeImageLike(Trachea, m_Working_Support, GetBackgroundValue()); + + // Find extrema post positions + std::vector tracheaPostPositionsA; + std::vector tracheaPostPositionsB; + clitk::SliceBySliceBuildLineSegmentAccordingToExtremaPosition(Trachea, + GetBackgroundValue(), 2, + 1, false, // Post + 0, // Horizontal line + 1, + tracheaPostPositionsA, + tracheaPostPositionsB); + // Cut post to this line + clitk::SliceBySliceSetBackgroundFromLineSeparation(m_Working_Support, + tracheaPostPositionsA, + tracheaPostPositionsB, + GetBackgroundValue(), 1, -10); + // END + StopCurrentStep(m_Working_Support); + m_ListOfStations["2R"] = m_Working_Support; + m_ListOfStations["2L"] = m_Working_Support; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +// Build a vtk mesh from a list of slice number/closed-contours +template +vtkSmartPointer +clitk::ExtractLymphStationsFilter:: +Build3DMeshFrom2DContour(const std::vector & points) +{ + // create a contour, polydata. + vtkSmartPointer mesh = vtkSmartPointer::New(); + mesh->Allocate(); //for cell structures + mesh->SetPoints(vtkPoints::New()); + vtkIdType ids[2]; + int point_number = points.size(); + for (unsigned int i=0; iGetPoints()->InsertNextPoint(points[i][0],points[i][1],points[i][2]); + ids[0]=i; + ids[1]=(ids[0]+1)%point_number; //0-1,1-2,...,n-1-0 + mesh->GetLines()->InsertNextCell(2,ids); + } + // Return + return mesh; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_2RL_LR_Limits() +{ + // --------------------------------------------------------------------------- + StartNewStep("[Station 2RL] Left/Right limits with Aorta"); + MaskImagePointer Aorta = GetAFDB()->template GetImage("Aorta"); + // DD(GetFuzzyThreshold("2RL", "BrachioCephalicVein")); + m_Working_Support = + clitk::SliceBySliceRelativePosition(m_Working_Support, Aorta, 2, + GetFuzzyThreshold("2RL", "Aorta"), + "RightTo", false, 2, true, false); + // END + StopCurrentStep(m_Working_Support); + m_ListOfStations["2R"] = m_Working_Support; + m_ListOfStations["2L"] = m_Working_Support; + + // --------------------------------------------------------------------------- + StartNewStep("[Station 2RL] Left/Right limits with SubclavianArtery (Right)"); + + // SliceBySliceRelativePosition + select CCL most at Right + MaskImagePointer SubclavianArtery = GetAFDB()->template GetImage("SubclavianArtery"); + typedef SliceBySliceRelativePositionFilter SliceRelPosFilterType; + typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New(); + sliceRelPosFilter->VerboseStepFlagOff(); + sliceRelPosFilter->WriteStepFlagOff(); + sliceRelPosFilter->SetInput(m_Working_Support); + sliceRelPosFilter->SetInputObject(SubclavianArtery); + sliceRelPosFilter->SetDirection(2); + sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("2RL", "SubclavianArteryRight")); + sliceRelPosFilter->AddOrientationTypeString("NotRightTo"); + sliceRelPosFilter->IntermediateSpacingFlagOn(); + sliceRelPosFilter->SetIntermediateSpacing(2); + sliceRelPosFilter->UniqueConnectedComponentBySliceOff(); + sliceRelPosFilter->UseASingleObjectConnectedComponentBySliceFlagOff(); + + sliceRelPosFilter->CCLSelectionFlagOn(); // select one CCL by slice + sliceRelPosFilter->SetCCLSelectionDimension(0); // select according to X (0) axis + sliceRelPosFilter->SetCCLSelectionDirection(-1); // select most at Right + sliceRelPosFilter->CCLSelectionIgnoreSingleCCLFlagOn(); // ignore if only one CCL + + sliceRelPosFilter->AutoCropFlagOn(); + sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn(); + sliceRelPosFilter->RemoveObjectFlagOff(); + sliceRelPosFilter->Update(); + m_Working_Support = sliceRelPosFilter->GetOutput(); + + // END + StopCurrentStep(m_Working_Support); + m_ListOfStations["2R"] = m_Working_Support; + m_ListOfStations["2L"] = m_Working_Support; + + + // --------------------------------------------------------------------------- + StartNewStep("[Station 2RL] Left/Right limits with SubclavianArtery (Left)"); + + // SliceBySliceRelativePosition + select CCL most at Right + sliceRelPosFilter = SliceRelPosFilterType::New(); + sliceRelPosFilter->VerboseStepFlagOff(); + sliceRelPosFilter->WriteStepFlagOff(); + sliceRelPosFilter->SetInput(m_Working_Support); + sliceRelPosFilter->SetInputObject(SubclavianArtery); + sliceRelPosFilter->SetDirection(2); + sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("2RL", "SubclavianArteryLeft")); + sliceRelPosFilter->AddOrientationTypeString("NotLeftTo"); + sliceRelPosFilter->IntermediateSpacingFlagOn(); + sliceRelPosFilter->SetIntermediateSpacing(2); + sliceRelPosFilter->UniqueConnectedComponentBySliceOff(); + sliceRelPosFilter->UseASingleObjectConnectedComponentBySliceFlagOff(); + + sliceRelPosFilter->CCLSelectionFlagOn(); // select one CCL by slice + sliceRelPosFilter->SetCCLSelectionDimension(0); // select according to X (0) axis + sliceRelPosFilter->SetCCLSelectionDirection(+1); // select most at Left + sliceRelPosFilter->CCLSelectionIgnoreSingleCCLFlagOff(); // do not ignore if only one CCL + + sliceRelPosFilter->AutoCropFlagOn(); + sliceRelPosFilter->IgnoreEmptySliceObjectFlagOn(); + sliceRelPosFilter->RemoveObjectFlagOff(); + sliceRelPosFilter->Update(); + m_Working_Support = sliceRelPosFilter->GetOutput(); + + // END + StopCurrentStep(m_Working_Support); + m_ListOfStations["2R"] = m_Working_Support; + m_ListOfStations["2L"] = m_Working_Support; +} +//-------------------------------------------------------------------- + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_2RL_Remove_Structures() +{ + Remove_Structures("2RL", "BrachioCephalicVein"); + Remove_Structures("2RL", "CommonCarotidArtery"); + Remove_Structures("2RL", "SubclavianArtery"); + Remove_Structures("2RL", "Thyroid"); + Remove_Structures("2RL", "Aorta"); + + // END + StopCurrentStep(m_Working_Support); + m_ListOfStations["2R"] = m_Working_Support; + m_ListOfStations["2L"] = m_Working_Support; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_2RL_SeparateRL() +{ + // --------------------------------------------------------------------------- + StartNewStep("[Station 2RL] Separate 2R/2L according to Trachea"); + + /*Rod says: + + "For station 2 there is a shift, dividing 2R from 2L, from midline + to the left paratracheal border." + + Algo: + - Consider Trachea SliceBySlice + - find extrema at Left + - add margins towards Right + - remove what is at Left of this line + */ + + // Get Trachea + MaskImagePointer Trachea = GetAFDB()->template GetImage("Trachea"); + + // Resize like the current support (to have the same number of slices) + Trachea = clitk::ResizeImageLike(Trachea, m_Working_Support, GetBackgroundValue()); + + // Find extrema post positions + std::vector tracheaLeftPositionsA; + std::vector tracheaLeftPositionsB; + clitk::SliceBySliceBuildLineSegmentAccordingToExtremaPosition(Trachea, + GetBackgroundValue(), 2, + 0, false, // Left + 1, // Vertical line + 1, // margins + tracheaLeftPositionsA, + tracheaLeftPositionsB); + // Copy support for R and L + typedef itk::ImageDuplicator DuplicatorType; + DuplicatorType::Pointer duplicator = DuplicatorType::New(); + duplicator->SetInputImage(m_Working_Support); + duplicator->Update(); + MaskImageType::Pointer m_Working_Support2 = duplicator->GetOutput(); + + // Cut post to this line for Right part + clitk::SliceBySliceSetBackgroundFromLineSeparation(m_Working_Support, + tracheaLeftPositionsA, + tracheaLeftPositionsB, + GetBackgroundValue(), 0, -10); + writeImage(m_Working_Support, "R.mhd"); + + // Cut post to this line for Left part + clitk::SliceBySliceSetBackgroundFromLineSeparation(m_Working_Support2, + tracheaLeftPositionsA, + tracheaLeftPositionsB, + GetBackgroundValue(), 0, +10); + writeImage(m_Working_Support2, "L.mhd"); + + // END + StopCurrentStep(m_Working_Support); + m_ListOfStations["2R"] = m_Working_Support; + m_ListOfStations["2L"] = m_Working_Support2; +} +//-------------------------------------------------------------------- diff --git a/segmentation/clitkExtractLymphStation_3A.txx b/segmentation/clitkExtractLymphStation_3A.txx index 0446913..a55071c 100644 --- a/segmentation/clitkExtractLymphStation_3A.txx +++ b/segmentation/clitkExtractLymphStation_3A.txx @@ -1,23 +1,3 @@ -/*========================================================================= - Program: vv http://www.creatis.insa-lyon.fr/rio/vv - - Authors belong to: - - University of LYON http://www.universite-lyon.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 - the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR - PURPOSE. See the copyright notices for more information. - - It is distributed under dual licence - - - BSD See included LICENSE.txt file - - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html -===========================================================================*/ - -#include -#include //-------------------------------------------------------------------- template @@ -25,6 +5,8 @@ void clitk::ExtractLymphStationsFilter:: ExtractStation_3A_SetDefaultValues() { + SetFuzzyThreshold("3A", "Sternum", 0.5); + SetFuzzyThreshold("3A", "SubclavianArtery", 0.5); } //-------------------------------------------------------------------- @@ -73,7 +55,27 @@ ExtractStation_3A_Ant_Limits() MaskImagePointer Sternum = GetAFDB()->template GetImage("Sternum"); m_Working_Support = clitk::SliceBySliceRelativePosition(m_Working_Support, Sternum, 2, - 0.5, "PostTo", + GetFuzzyThreshold("3A", "Sternum"), "PostTo", + false, 3, true, false); + StopCurrentStep(m_Working_Support); + m_ListOfStations["3A"] = m_Working_Support; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_3A_Post_Limits() +{ + StartNewStep("[Station 3A] Post limits with SubclavianArtery"); + + // Get Sternum, keep posterior part. + MaskImagePointer SubclavianArtery = GetAFDB()->template GetImage("SubclavianArtery"); + m_Working_Support = + clitk::SliceBySliceRelativePosition(m_Working_Support, SubclavianArtery, 2, + GetFuzzyThreshold("3A", "SubclavianArtery"), "AntTo", false, 3, true, false); StopCurrentStep(m_Working_Support); m_ListOfStations["3A"] = m_Working_Support; diff --git a/segmentation/clitkExtractLymphStation_3P.txx b/segmentation/clitkExtractLymphStation_3P.txx index 7f1d5bf..1c9fdc0 100644 --- a/segmentation/clitkExtractLymphStation_3P.txx +++ b/segmentation/clitkExtractLymphStation_3P.txx @@ -1,23 +1,3 @@ -/*========================================================================= - Program: vv http://www.creatis.insa-lyon.fr/rio/vv - - Authors belong to: - - University of LYON http://www.universite-lyon.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 - the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR - PURPOSE. See the copyright notices for more information. - - It is distributed under dual licence - - - BSD See included LICENSE.txt file - - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html -===========================================================================*/ - -#include -#include //-------------------------------------------------------------------- template @@ -82,13 +62,13 @@ ExtractStation_3P_Remove_Structures() StartNewStep("[Station 3P] Remove some structures."); - Remove_Structures("Aorta"); - Remove_Structures("VertebralBody"); - Remove_Structures("SubclavianArtery"); - Remove_Structures("Esophagus"); - Remove_Structures("Azygousvein"); - Remove_Structures("Thyroid"); - Remove_Structures("VertebralArtery"); + Remove_Structures("3P", "Aorta"); + Remove_Structures("3P", "VertebralBody"); + Remove_Structures("3P", "SubclavianArtery"); + Remove_Structures("3P", "Esophagus"); + Remove_Structures("3P", "Azygousvein"); + Remove_Structures("3P", "Thyroid"); + Remove_Structures("3P", "VertebralArtery"); StopCurrentStep(m_Working_Support); m_ListOfStations["3P"] = m_Working_Support; diff --git a/segmentation/clitkExtractLymphStation_4RL.txx b/segmentation/clitkExtractLymphStation_4RL.txx index 5aad362..462ae9d 100644 --- a/segmentation/clitkExtractLymphStation_4RL.txx +++ b/segmentation/clitkExtractLymphStation_4RL.txx @@ -1,20 +1,3 @@ -/*========================================================================= - Program: vv http://www.creatis.insa-lyon.fr/rio/vv - - Authors belong to: - - University of LYON http://www.universite-lyon.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 - the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR - PURPOSE. See the copyright notices for more information. - - It is distributed under dual licence - - - BSD See included LICENSE.txt file - - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html -===========================================================================*/ //-------------------------------------------------------------------- template void diff --git a/segmentation/clitkExtractLymphStation_7.txx b/segmentation/clitkExtractLymphStation_7.txx index c581166..7bb4643 100644 --- a/segmentation/clitkExtractLymphStation_7.txx +++ b/segmentation/clitkExtractLymphStation_7.txx @@ -1,20 +1,3 @@ -/*========================================================================= - Program: vv http://www.creatis.insa-lyon.fr/rio/vv - - Authors belong to: - - University of LYON http://www.universite-lyon.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 - the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR - PURPOSE. See the copyright notices for more information. - - It is distributed under dual licence - - - BSD See included LICENSE.txt file - - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html -===========================================================================*/ //-------------------------------------------------------------------- template @@ -22,38 +5,12 @@ void clitk::ExtractLymphStationsFilter:: ExtractStation_7_SetDefaultValues() { - SetFuzzyThresholdForS7("Bronchi", 0.1); - SetFuzzyThresholdForS7("LeftSuperiorPulmonaryVein", 0.3); - SetFuzzyThresholdForS7("RightSuperiorPulmonaryVein", 0.2); - SetFuzzyThresholdForS7("RightPulmonaryArtery", 0.3); - SetFuzzyThresholdForS7("LeftPulmonaryArtery", 0.5); - SetFuzzyThresholdForS7("SVC", 0.2); -} -//-------------------------------------------------------------------- - -//-------------------------------------------------------------------- -template -void -clitk::ExtractLymphStationsFilter:: -SetFuzzyThresholdForS7(std::string tag, double value) -{ - m_FuzzyThresholdForS7[tag] = value; -} -//-------------------------------------------------------------------- - - -//-------------------------------------------------------------------- -template -double -clitk::ExtractLymphStationsFilter:: -GetFuzzyThresholdForS7(std::string tag) -{ - if (m_FuzzyThresholdForS7.find(tag) != m_FuzzyThresholdForS7.end()) { - return m_FuzzyThresholdForS7[tag]; - } - else { - clitkExceptionMacro("Could not find options "+tag+" in the m_FuzzyThresholdForS7 list"); - } + SetFuzzyThreshold("7", "Bronchi", 0.1); + SetFuzzyThreshold("7", "LeftSuperiorPulmonaryVein", 0.3); + SetFuzzyThreshold("7", "RightSuperiorPulmonaryVein", 0.2); + SetFuzzyThreshold("7", "RightPulmonaryArtery", 0.3); + SetFuzzyThreshold("7", "LeftPulmonaryArtery", 0.5); + SetFuzzyThreshold("7", "SVC", 0.2); } //-------------------------------------------------------------------- @@ -175,7 +132,7 @@ ExtractStation_7_RL_Limits() m_Working_Support = clitk::SliceBySliceRelativePosition(m_Working_Support, m_LeftBronchus, 2, - GetFuzzyThresholdForS7("Bronchi"), "RightTo", + GetFuzzyThreshold("7", "Bronchi"), "RightTo", false, 3, false); StopCurrentStep(m_Working_Support); @@ -184,7 +141,7 @@ ExtractStation_7_RL_Limits() StartNewStep("[Station7] Limits with bronchus : LeftTo the right bronchus"); m_Working_Support = clitk::SliceBySliceRelativePosition(m_Working_Support, m_RightBronchus, 2, - GetFuzzyThresholdForS7("Bronchi"), "LeftTo", + GetFuzzyThreshold("7", "Bronchi"), "LeftTo", false, 3, false); StopCurrentStep(m_Working_Support); @@ -198,7 +155,7 @@ ExtractStation_7_RL_Limits() sliceRelPosFilter->SetInput(m_Working_Support); sliceRelPosFilter->SetInputObject(LeftSuperiorPulmonaryVein); sliceRelPosFilter->SetDirection(2); - sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThresholdForS7("LeftSuperiorPulmonaryVein")); + sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("7", "LeftSuperiorPulmonaryVein")); sliceRelPosFilter->AddOrientationTypeString("NotLeftTo"); sliceRelPosFilter->AddOrientationTypeString("NotAntTo"); sliceRelPosFilter->SetIntermediateSpacingFlag(true); @@ -223,7 +180,7 @@ ExtractStation_7_RL_Limits() sliceRelPosFilter->SetInput(m_Working_Support); sliceRelPosFilter->SetInputObject(RightSuperiorPulmonaryVein); sliceRelPosFilter->SetDirection(2); - sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThresholdForS7("RightSuperiorPulmonaryVein")); + sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("7", "RightSuperiorPulmonaryVein")); sliceRelPosFilter->AddOrientationTypeString("NotRightTo"); sliceRelPosFilter->AddOrientationTypeString("NotAntTo"); sliceRelPosFilter->AddOrientationTypeString("NotPostTo"); @@ -248,7 +205,7 @@ ExtractStation_7_RL_Limits() sliceRelPosFilter->SetInput(m_Working_Support); sliceRelPosFilter->SetInputObject(RightPulmonaryArtery); sliceRelPosFilter->SetDirection(2); - sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThresholdForS7("RightPulmonaryArtery")); + sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("7", "RightPulmonaryArtery")); sliceRelPosFilter->AddOrientationTypeString("NotAntTo"); sliceRelPosFilter->SetIntermediateSpacingFlag(true); sliceRelPosFilter->SetIntermediateSpacing(3); @@ -266,7 +223,7 @@ ExtractStation_7_RL_Limits() sliceRelPosFilter->SetInput(m_Working_Support); sliceRelPosFilter->SetInputObject(LeftPulmonaryArtery); sliceRelPosFilter->SetDirection(2); - sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThresholdForS7("LeftPulmonaryArtery")); + sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("7", "LeftPulmonaryArtery")); sliceRelPosFilter->AddOrientationTypeString("NotAntTo"); sliceRelPosFilter->SetIntermediateSpacingFlag(true); sliceRelPosFilter->SetIntermediateSpacing(3); @@ -283,7 +240,7 @@ ExtractStation_7_RL_Limits() sliceRelPosFilter->SetInput(m_Working_Support); sliceRelPosFilter->SetInputObject(SVC); sliceRelPosFilter->SetDirection(2); - sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThresholdForS7("SVC")); + sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("7", "SVC")); sliceRelPosFilter->AddOrientationTypeString("NotRightTo"); sliceRelPosFilter->AddOrientationTypeString("NotAntTo"); sliceRelPosFilter->SetIntermediateSpacingFlag(true); @@ -386,14 +343,14 @@ ExtractStation_7_Remove_Structures() //-------------------------------------------------------------------- StartNewStep("[Station7] remove some structures"); - Remove_Structures("AzygousVein"); - Remove_Structures("Aorta"); - Remove_Structures("Esophagus"); - Remove_Structures("RightPulmonaryArtery"); - Remove_Structures("LeftPulmonaryArtery"); - Remove_Structures("LeftSuperiorPulmonaryVein"); - Remove_Structures("PulmonaryTrunk"); - Remove_Structures("VertebralBody"); + Remove_Structures("7", "AzygousVein"); + Remove_Structures("7", "Aorta"); + Remove_Structures("7", "Esophagus"); + Remove_Structures("7", "RightPulmonaryArtery"); + Remove_Structures("7", "LeftPulmonaryArtery"); + Remove_Structures("7", "LeftSuperiorPulmonaryVein"); + Remove_Structures("7", "PulmonaryTrunk"); + Remove_Structures("7", "VertebralBody"); // END StopCurrentStep(m_Working_Support); diff --git a/segmentation/clitkExtractLymphStation_8.txx b/segmentation/clitkExtractLymphStation_8.txx index 75c0129..bb1ce79 100644 --- a/segmentation/clitkExtractLymphStation_8.txx +++ b/segmentation/clitkExtractLymphStation_8.txx @@ -1,20 +1,3 @@ -/*========================================================================= - Program: vv http://www.creatis.insa-lyon.fr/rio/vv - - Authors belong to: - - University of LYON http://www.universite-lyon.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 - the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR - PURPOSE. See the copyright notices for more information. - - It is distributed under dual licence - - - BSD See included LICENSE.txt file - - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html -===========================================================================*/ #include #include @@ -31,7 +14,7 @@ ExtractStation_8_SetDefaultValues() SetEsophagusDiltationForAnt(p); p[0] = 5; p[1] = 10; p[2] = 1; SetEsophagusDiltationForRight(p); - SetFuzzyThresholdForS8(0.5); + SetFuzzyThreshold("8", "Esophagus", 0.5); SetInjectedThresholdForS8(150); } //-------------------------------------------------------------------- @@ -180,9 +163,9 @@ ExtractStation_8_Post_Limits() // Convert 2D points in slice into 3D points std::vector vertebralAntPositions; - clitk::PointsUtils::Convert2DTo3DList(vertebralAntPositionBySlice, - VertebralBody, - vertebralAntPositions); + clitk::PointsUtils::Convert2DMapTo3DList(vertebralAntPositionBySlice, + VertebralBody, + vertebralAntPositions); // DEBUG : write list of points clitk::WriteListOfLandmarks(vertebralAntPositions, @@ -485,7 +468,7 @@ ExtractStation_8_Ant_Inf_Limits() relPosFilter->UniqueConnectedComponentBySliceOff(); relPosFilter->SetIntermediateSpacing(3); relPosFilter->IntermediateSpacingFlagOn(); - relPosFilter->SetFuzzyThreshold(GetFuzzyThresholdForS8()); + relPosFilter->SetFuzzyThreshold(GetFuzzyThreshold("8", "Esophagus")); relPosFilter->RemoveObjectFlagOff(); // Do not exclude here because it is dilated relPosFilter->CombineWithOrFlagOff(); // NO ! relPosFilter->IgnoreEmptySliceObjectFlagOn(); @@ -1148,8 +1131,8 @@ ExtractStation_8_Remove_Structures() //-------------------------------------------------------------------- StartNewStep("[Station8] remove some structures"); - Remove_Structures("Aorta"); - Remove_Structures("Esophagus"); + Remove_Structures("8", "Aorta"); + Remove_Structures("8", "Esophagus"); // END StopCurrentStep(m_Working_Support); diff --git a/segmentation/clitkExtractLymphStations.cxx b/segmentation/clitkExtractLymphStations.cxx index 6ddd818..ec41a4e 100644 --- a/segmentation/clitkExtractLymphStations.cxx +++ b/segmentation/clitkExtractLymphStations.cxx @@ -3,7 +3,7 @@ Authors belong to: - University of LYON http://www.universite-lyon.fr/ - - Léon Bérard cancer center http://www.centreleonberard.fr + - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr This software is distributed WITHOUT ANY WARRANTY; without even @@ -14,7 +14,7 @@ - BSD See included LICENSE.txt file - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html -===========================================================================**/ +======================================================================-====*/ // clitk #include "clitkExtractLymphStations_ggo.h" diff --git a/segmentation/clitkExtractLymphStations.ggo b/segmentation/clitkExtractLymphStations.ggo index 57cf24a..d789824 100644 --- a/segmentation/clitkExtractLymphStations.ggo +++ b/segmentation/clitkExtractLymphStations.ggo @@ -20,16 +20,28 @@ option "output" o "Output lungs mask filename" string no option "station" - "Force to compute station even if already exist in the DB" string no multiple section "Options for Station 8" -option "maxAntSpine" - "Distance max to anterior part of the spine in mm" double no default="10" -option "esophagusDilatationForAnt" - "Dilatation of esophagus, in mm, for 'anterior' limits (default=15,2,1)" double no multiple +option "maxAntSpine" - "Distance max to anterior part of the spine in mm" double no default="10" +option "esophagusDilatationForAnt" - "Dilatation of esophagus, in mm, for 'anterior' limits (default=15,2,1)" double no multiple option "esophagusDilatationForRight" - "Dilatation of esophagus, in mm, for 'right' limits (default=5,10,1)" double no multiple -option "fuzzyThresholdForS8" - "Threshold for 'Post' to dilated Esophagus" double default="0.5" no -option "injectedThresholdForS8" - "Threshold injected areas in the ct" double default="150" no +option "tS8_Esophagus" - "Threshold for 'Post' to dilated Esophagus" double default="0.5" no +option "injectedThresholdForS8" - "Threshold injected areas in the ct" double default="150" no section "Options for Station 7" -option "tS7_Bronchi" - "Threshold for Left/Right bronchi" double default="0.1" no -option "tS7_LeftSuperiorPulmonaryVein" - "Threshold for LeftSuperiorPulmonaryVein" double default="0.3" no +option "tS7_Bronchi" - "Threshold for Left/Right bronchi" double default="0.1" no +option "tS7_LeftSuperiorPulmonaryVein" - "Threshold for LeftSuperiorPulmonaryVein" double default="0.3" no option "tS7_RightSuperiorPulmonaryVein" - "Threshold for RightSuperiorPulmonaryVein" double default="0.2" no -option "tS7_RightPulmonaryArtery" - "Threshold for RightPulmonaryArtery" double default="0.3" no -option "tS7_LeftPulmonaryArtery" - "Threshold for LeftPulmonaryArtery (NOT USEFUL YET)" double default="0.5" no -option "tS7_SVC" - "Threshold for SVC" double default="0.2" no \ No newline at end of file +option "tS7_RightPulmonaryArtery" - "Threshold for RightPulmonaryArtery" double default="0.3" no +option "tS7_LeftPulmonaryArtery" - "Threshold for LeftPulmonaryArtery (NOT USEFUL YET)" double default="0.5" no +option "tS7_SVC" - "Threshold for SVC" double default="0.2" no + +section "Options for Station 3A" +option "tS3A_Sternum" - "Threshold for Post to Sternum" double default="0.5" no +option "tS3A_SubclavianArtery" - "Threshold for Ant to SubclavianArtery" double default="0.2" no + +section "Options for Station 2RL" +option "tS2RL_CommonCarotidArtery" - "Threshold for NotAntTo to CommonCarotidArtery" double default="0.7" no +option "tS2RL_BrachioCephalicTrunk" - "Threshold for NotAntTo to BrachioCephalicTrunk" double default="0.7" no +option "tS2RL_BrachioCephalicVein" - "Threshold for NotAntTo to BrachioCephalicVein" double default="0.3" no +option "tS2RL_Aorta" - "Threshold for RightTo to Aorta" double default="0.7" no +option "tS2RL_SubclavianArteryRight" - "Threshold for NotLeft to SubclavianArteryRight" double default="0.5" no +option "tS2RL_SubclavianArteryLeft" - "Threshold for NotRight to SubclavianArteryLeft" double default="0.8" no diff --git a/segmentation/clitkExtractLymphStationsFilter.h b/segmentation/clitkExtractLymphStationsFilter.h index 6f81423..8d160b1 100644 --- a/segmentation/clitkExtractLymphStationsFilter.h +++ b/segmentation/clitkExtractLymphStationsFilter.h @@ -3,7 +3,7 @@ Authors belong to: - University of LYON http://www.universite-lyon.fr/ - - Léon Bérard cancer center http://www.centreleonberard.fr + - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr This software is distributed WITHOUT ANY WARRANTY; without even @@ -14,7 +14,7 @@ - BSD See included LICENSE.txt file - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html - ===========================================================================**/ + ======================================================================-====*/ #ifndef CLITKEXTRACTLYMPHSTATIONSFILTER_H #define CLITKEXTRACTLYMPHSTATIONSFILTER_H @@ -23,6 +23,9 @@ #include "clitkFilterBase.h" #include "clitkFilterWithAnatomicalFeatureDatabaseManagement.h" +// vtk +#include + namespace clitk { //-------------------------------------------------------------------- @@ -72,6 +75,7 @@ namespace clitk { typedef itk::Image MaskSliceType; typedef typename MaskSliceType::Pointer MaskSlicePointer; + typedef typename MaskSliceType::PointType MaskSlicePointType; /** ImageDimension constants */ itkStaticConstMacro(ImageDimension, unsigned int, ImageType::ImageDimension); @@ -89,19 +93,16 @@ namespace clitk { itkGetConstMacro(EsophagusDiltationForAnt, MaskImagePointType); itkSetMacro(EsophagusDiltationForRight, MaskImagePointType); itkGetConstMacro(EsophagusDiltationForRight, MaskImagePointType); - itkSetMacro(FuzzyThresholdForS8, double); - itkGetConstMacro(FuzzyThresholdForS8, double); - itkSetMacro(InjectedThresholdForS8, double); itkGetConstMacro(InjectedThresholdForS8, double); // Station 7 - void SetFuzzyThresholdForS7(std::string tag, double value); - double GetFuzzyThresholdForS7(std::string tag); // All stations bool GetComputeStation(std::string s); void AddComputeStation(std::string station) ; + void SetFuzzyThreshold(std::string station, std::string tag, double value); + double GetFuzzyThreshold(std::string station, std::string tag); protected: ExtractLymphStationsFilter(); @@ -120,14 +121,17 @@ namespace clitk { std::map m_ComputeStationMap; bool CheckForStation(std::string station); - void Remove_Structures(std::string s); + void Remove_Structures(std::string station, std::string s); + + // Global parameters + typedef std::map FuzzyThresholdByStructureType; + std::map m_FuzzyThreshold; // Station 8 double m_DistanceMaxToAnteriorPartOfTheSpine; double m_DiaphragmInferiorLimit; double m_CarinaZ; double m_OriginOfRightMiddleLobeBronchusZ; - double m_FuzzyThresholdForS8; double m_InjectedThresholdForS8; MaskImagePointer m_Esophagus; MaskImagePointType m_EsophagusDiltationForAnt; @@ -159,11 +163,24 @@ namespace clitk { void ExtractStation_3P_LR_sup_Limits_2(); void ExtractStation_3P_LR_inf_Limits(); + // Station 2RL + void ExtractStation_2RL(); + void ExtractStation_2RL_SetDefaultValues(); + void ExtractStation_2RL_SI_Limits(); + void ExtractStation_2RL_Ant_Limits(); + void ExtractStation_2RL_Ant_Limits2(); + void ExtractStation_2RL_Post_Limits(); + void ExtractStation_2RL_LR_Limits(); + void ExtractStation_2RL_Remove_Structures(); + void ExtractStation_2RL_SeparateRL(); + vtkSmartPointer Build3DMeshFrom2DContour(const std::vector & points); + // Station 3A void ExtractStation_3A(); void ExtractStation_3A_SetDefaultValues(); void ExtractStation_3A_SI_Limits(); void ExtractStation_3A_Ant_Limits(); + void ExtractStation_3A_Post_Limits(); // Station 7 void ExtractStation_7(); @@ -173,7 +190,6 @@ namespace clitk { void ExtractStation_7_Posterior_Limits(); void ExtractStation_7_Remove_Structures(); MaskImagePointer m_Working_Trachea; - std::map m_FuzzyThresholdForS7; MaskImagePointer m_LeftBronchus; MaskImagePointer m_RightBronchus; typedef std::vector ListOfPointsType; @@ -212,6 +228,7 @@ namespace clitk { #include "clitkExtractLymphStationsFilter.txx" #include "clitkExtractLymphStation_8.txx" #include "clitkExtractLymphStation_3P.txx" +#include "clitkExtractLymphStation_2RL.txx" #include "clitkExtractLymphStation_3A.txx" #include "clitkExtractLymphStation_7.txx" #include "clitkExtractLymphStation_4RL.txx" diff --git a/segmentation/clitkExtractLymphStationsFilter.old.h b/segmentation/clitkExtractLymphStationsFilter.old.h new file mode 100644 index 0000000..6be4764 --- /dev/null +++ b/segmentation/clitkExtractLymphStationsFilter.old.h @@ -0,0 +1,124 @@ +/*========================================================================= + Program: vv http://www.creatis.insa-lyon.fr/rio/vv + + Authors belong to: + - University of LYON http://www.universite-lyon.fr/ + - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr + - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the copyright notices for more information. + + It is distributed under dual licence + + - BSD See included LICENSE.txt file + - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html + ======================================================================-====*/ + +#ifndef CLITKEXTRACTLYMPHSTATIONSFILTER_H +#define CLITKEXTRACTLYMPHSTATIONSFILTER_H + +// clitk +#include "clitkFilterBase.h" +#include "clitkFilterWithAnatomicalFeatureDatabaseManagement.h" + +namespace clitk { + + //-------------------------------------------------------------------- + /* + Try to extract the LymphStations part of a thorax CT. + Inputs : + - Patient label image + - Lungs label image + - Bones label image + */ + //-------------------------------------------------------------------- + + template + class ITK_EXPORT ExtractLymphStationsFilter: + public virtual clitk::FilterBase, + public clitk::FilterWithAnatomicalFeatureDatabaseManagement, + public itk::ImageToImageFilter > + { + + public: + /** Standard class typedefs. */ + typedef itk::ImageToImageFilter > Superclass; + typedef ExtractLymphStationsFilter Self; + typedef itk::SmartPointer Pointer; + typedef itk::SmartPointer ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(ExtractLymphStationsFilter, InPlaceImageFilter); + + /** Some convenient typedefs. */ + typedef TImageType ImageType; + typedef typename ImageType::ConstPointer ImageConstPointer; + typedef typename ImageType::Pointer ImagePointer; + typedef typename ImageType::RegionType ImageRegionType; + typedef typename ImageType::PixelType ImagePixelType; + typedef typename ImageType::SizeType ImageSizeType; + typedef typename ImageType::IndexType ImageIndexType; + typedef typename ImageType::PointType ImagePointType; + + typedef uchar MaskImagePixelType; + typedef itk::Image MaskImageType; + typedef typename MaskImageType::Pointer MaskImagePointer; + + /** Connect inputs */ + // void SetInput(const TImageType * image); + + /** ImageDimension constants */ + itkStaticConstMacro(ImageDimension, unsigned int, TImageType::ImageDimension); + FILTERBASE_INIT; + + itkGetConstMacro(BackgroundValue, ImagePixelType); + itkGetConstMacro(ForegroundValue, ImagePixelType); + itkSetMacro(BackgroundValue, ImagePixelType); + itkSetMacro(ForegroundValue, ImagePixelType); + + itkSetMacro(FuzzyThreshold, double); + itkGetConstMacro(FuzzyThreshold, double); + + protected: + ExtractLymphStationsFilter(); + virtual ~ExtractLymphStationsFilter() {} + + virtual void GenerateOutputInformation(); + virtual void GenerateInputRequestedRegion(); + virtual void GenerateData(); + + MaskImagePointer m_mediastinum; + MaskImagePointer m_trachea; + MaskImagePointer m_working_mediastinum; + MaskImagePointer m_working_trachea; + MaskImagePointer m_output; + + ImagePixelType m_BackgroundValue; + ImagePixelType m_ForegroundValue; + double m_CarinaZPositionInMM; + bool m_CarinaZPositionInMMIsSet; + double m_MiddleLobeBronchusZPositionInMM; + + double m_IntermediateSpacing; + double m_FuzzyThreshold; + + private: + ExtractLymphStationsFilter(const Self&); //purposely not implemented + void operator=(const Self&); //purposely not implemented + + }; // end class + //-------------------------------------------------------------------- + +} // end namespace clitk +//-------------------------------------------------------------------- + +#ifndef ITK_MANUAL_INSTANTIATION +#include "clitkExtractLymphStationsFilter.txx" +#endif + +#endif diff --git a/segmentation/clitkExtractLymphStationsFilter.old.txx b/segmentation/clitkExtractLymphStationsFilter.old.txx new file mode 100644 index 0000000..c95be39 --- /dev/null +++ b/segmentation/clitkExtractLymphStationsFilter.old.txx @@ -0,0 +1,395 @@ +/*========================================================================= + Program: vv http://www.creatis.insa-lyon.fr/rio/vv + + Authors belong to: + - University of LYON http://www.universite-lyon.fr/ + - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr + - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the copyright notices for more information. + + It is distributed under dual licence + + - BSD See included LICENSE.txt file + - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html + ======================================================================-====*/ + +#ifndef CLITKEXTRACTLYMPHSTATIONSFILTER_TXX +#define CLITKEXTRACTLYMPHSTATIONSFILTER_TXX + +// clitk +#include "clitkCommon.h" +#include "clitkExtractLymphStationsFilter.h" +#include "clitkAddRelativePositionConstraintToLabelImageFilter.h" +#include "clitkSegmentationUtils.h" +#include "clitkAutoCropFilter.h" +#include "clitkSegmentationUtils.h" +#include "clitkSliceBySliceRelativePositionFilter.h" + +// itk +#include +#include +#include +#include +#include +#include +#include + +// itk ENST +#include "RelativePositionPropImageFilter.h" + +//-------------------------------------------------------------------- +template +clitk::ExtractLymphStationsFilter:: +ExtractLymphStationsFilter(): + clitk::FilterBase(), + clitk::FilterWithAnatomicalFeatureDatabaseManagement(), + itk::ImageToImageFilter() +{ + this->SetNumberOfRequiredInputs(1); + SetBackgroundValue(0); + SetForegroundValue(1); + SetFuzzyThreshold(0.5); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +GenerateOutputInformation() { + // Superclass::GenerateOutputInformation(); + + // Get input + LoadAFDB(); + m_mediastinum = GetAFDB()->template GetImage ("mediastinum"); + m_trachea = GetAFDB()->template GetImage ("trachea"); + ImagePointType carina; + GetAFDB()->GetPoint3D("carina", carina); + m_CarinaZPositionInMM = carina[2]; + DD(m_CarinaZPositionInMM); + ImagePointType mlb; + GetAFDB()->GetPoint3D("middleLobeBronchus", mlb); + m_MiddleLobeBronchusZPositionInMM = mlb[2]; + DD(m_MiddleLobeBronchusZPositionInMM); + + // ---------------------------------------------------------------- + // ---------------------------------------------------------------- + // Superior limit = carina + // Inferior limit = origin right middle lobe bronchus + StartNewStep("Inf/Sup mediastinum limits with carina/bronchus"); + m_working_mediastinum = + clitk::CropImageAlongOneAxis(m_mediastinum, 2, + m_MiddleLobeBronchusZPositionInMM, + m_CarinaZPositionInMM, true, + GetBackgroundValue()); + StopCurrentStep(m_working_mediastinum); + m_output = m_working_mediastinum; + + // ---------------------------------------------------------------- + // ---------------------------------------------------------------- + // Superior limit = carina + // Inferior limit = origin right middle lobe bronchus + StartNewStep("Inf/Sup trachea limits with carina/bronchus"); + m_working_trachea = + clitk::CropImageAlongOneAxis(m_trachea, 2, + m_MiddleLobeBronchusZPositionInMM, + m_CarinaZPositionInMM, true, + GetBackgroundValue()); + StopCurrentStep(m_working_trachea); + + // ---------------------------------------------------------------- + // ---------------------------------------------------------------- + // Separate trachea in two CCL + StartNewStep("Separate trachea under carina"); + + // Labelize and consider two main labels + m_working_trachea = Labelize(m_working_trachea, 0, true, 1); + + // Carina position must at the first slice that separate the two main bronchus (not superiorly) + // Check that upper slice is composed of at least two labels + typedef itk::ImageSliceIteratorWithIndex SliceIteratorType; + SliceIteratorType iter(m_working_trachea, m_working_trachea->GetLargestPossibleRegion()); + iter.SetFirstDirection(0); + iter.SetSecondDirection(1); + iter.GoToReverseBegin(); // Start from the end (because image is IS not SI) + int maxLabel=0; + while (!iter.IsAtReverseEndOfSlice()) { + while (!iter.IsAtReverseEndOfLine()) { + // DD(iter.GetIndex()); + if (iter.Get() > maxLabel) maxLabel = iter.Get(); + --iter; + } + iter.PreviousLine(); + } + if (maxLabel < 2) { + clitkExceptionMacro("First slice form Carina does not seems to seperate the two main bronchus. Abort"); + } + + // Compute centroid of both parts to identify the left from the right bronchus + typedef long LabelType; + static const unsigned int Dim = ImageType::ImageDimension; + typedef itk::ShapeLabelObject< LabelType, Dim > LabelObjectType; + typedef itk::LabelMap< LabelObjectType > LabelMapType; + typedef itk::LabelImageToLabelMapFilter ImageToMapFilterType; + typename ImageToMapFilterType::Pointer imageToLabelFilter = ImageToMapFilterType::New(); + typedef itk::ShapeLabelMapFilter ShapeFilterType; + typename ShapeFilterType::Pointer statFilter = ShapeFilterType::New(); + imageToLabelFilter->SetBackgroundValue(GetBackgroundValue()); + imageToLabelFilter->SetInput(m_working_trachea); + statFilter->SetInput(imageToLabelFilter->GetOutput()); + statFilter->Update(); + typename LabelMapType::Pointer labelMap = statFilter->GetOutput(); + + ImagePointType C1 = labelMap->GetLabelObject(1)->GetCentroid(); + ImagePointType C2 = labelMap->GetLabelObject(2)->GetCentroid(); + + ImagePixelType leftLabel; + ImagePixelType rightLabel; + if (C1[0] < C2[0]) { leftLabel = 1; rightLabel = 2; } + else { leftLabel = 2; rightLabel = 1; } + DD(leftLabel); + DD(rightLabel); + + StopCurrentStep(m_working_trachea); + + //----------------------------------------------------- + StartNewStep("Left limits with bronchus (slice by slice)"); + // Select LeftLabel (set one label to Backgroundvalue) + MaskImagePointer leftBronchus = + SetBackground(m_working_trachea, m_working_trachea, + rightLabel, GetBackgroundValue()); + MaskImagePointer rightBronchus = + SetBackground(m_working_trachea, m_working_trachea, + leftLabel, GetBackgroundValue()); + writeImage(leftBronchus, "left.mhd"); + writeImage(rightBronchus, "right.mhd"); + + m_working_mediastinum = + clitk::SliceBySliceRelativePosition(m_working_mediastinum, + leftBronchus, 2, + GetFuzzyThreshold(), "RightTo", + true, 4); + m_working_mediastinum = + clitk::SliceBySliceRelativePosition(m_working_mediastinum, + rightBronchus, + 2, GetFuzzyThreshold(), "LeftTo", + true, 4); + m_working_mediastinum = + clitk::SliceBySliceRelativePosition(m_working_mediastinum, + leftBronchus, + 2, GetFuzzyThreshold(), "AntTo", + true, 4, true); // NOT + m_working_mediastinum = + clitk::SliceBySliceRelativePosition(m_working_mediastinum, + rightBronchus, + 2, GetFuzzyThreshold(), "AntTo", + true, 4, true); + m_working_mediastinum = + clitk::SliceBySliceRelativePosition(m_working_mediastinum, + leftBronchus, + 2, GetFuzzyThreshold(), "PostTo", + true, 4, true); + m_working_mediastinum = + clitk::SliceBySliceRelativePosition(m_working_mediastinum, + rightBronchus, + 2, GetFuzzyThreshold(), "PostTo", + true, 4, true); + m_output = m_working_mediastinum; + StopCurrentStep(m_output); + + //----------------------------------------------------- + //----------------------------------------------------- + StartNewStep("Posterior limits"); + + // Left Bronchus slices + typedef clitk::ExtractSliceFilter ExtractSliceFilterType; + typedef typename ExtractSliceFilterType::SliceType SliceType; + typename ExtractSliceFilterType::Pointer + extractSliceFilter = ExtractSliceFilterType::New(); + extractSliceFilter->SetInput(leftBronchus); + extractSliceFilter->SetDirection(2); + extractSliceFilter->Update(); + std::vector leftBronchusSlices; + extractSliceFilter->GetOutputSlices(leftBronchusSlices); + + // Right Bronchus slices + extractSliceFilter = ExtractSliceFilterType::New(); + extractSliceFilter->SetInput(rightBronchus); + extractSliceFilter->SetDirection(2); + extractSliceFilter->Update(); + std::vector rightBronchusSlices ; + extractSliceFilter->GetOutputSlices(rightBronchusSlices); + + assert(leftBronchusSlices.size() == rightBronchusSlices.size()); + + std::vector leftPoints; + std::vector rightPoints; + for(uint i=0; i(leftBronchusSlices[i], 0, true, 10); + leftBronchusSlices[i] = KeepLabels(leftBronchusSlices[i], + GetBackgroundValue(), + GetForegroundValue(), 1, 1, true); + rightBronchusSlices[i] = Labelize(rightBronchusSlices[i], 0, true, 10); + rightBronchusSlices[i] = KeepLabels(rightBronchusSlices[i], + GetBackgroundValue(), + GetForegroundValue(), 1, 1, true); + double distance_max_from_center_point = 15; + + // ------- Find point in left Bronchus ------- + // find right most point in left = rightMost + SliceType::PointType a; + SliceType::PointType rightMost = + clitk::FindExtremaPointInAGivenDirection(leftBronchusSlices[i], + GetBackgroundValue(), + 0, false, a, 0); + // find post most point in left, not far away from rightMost + SliceType::PointType p = + clitk::FindExtremaPointInAGivenDirection(leftBronchusSlices[i], + GetBackgroundValue(), + 1, false, rightMost, + distance_max_from_center_point); + MaskImageType::PointType pp; + pp[0] = p[0]; pp[1] = p[1]; + pp[2] = i*leftBronchus->GetSpacing()[2] + leftBronchus->GetOrigin()[2]; + leftPoints.push_back(pp); + + // ------- Find point in right Bronchus ------- + // find left most point in right = leftMost + SliceType::PointType leftMost = + clitk::FindExtremaPointInAGivenDirection(rightBronchusSlices[i], + GetBackgroundValue(), + 0, true, a, 0); + // find post most point in left, not far away from leftMost + p = clitk::FindExtremaPointInAGivenDirection(rightBronchusSlices[i], + GetBackgroundValue(), + 1, false, leftMost, + distance_max_from_center_point); + pp[0] = p[0]; pp[1] = p[1]; + pp[2] = i*rightBronchus->GetSpacing()[2] + rightBronchus->GetOrigin()[2]; + rightPoints.push_back(pp); + } + + // DEBUG + std::ofstream osl; + openFileForWriting(osl, "left.txt"); + osl << "LANDMARKS1" << std::endl; + std::ofstream osr; + openFileForWriting(osr, "right.txt"); + osr << "LANDMARKS1" << std::endl; + for(uint i=0; i to accelerate, start with formula, when change sign -> stop and fill + */ + // typedef itk::ImageSliceIteratorWithIndex SliceIteratorType; + iter = SliceIteratorType(m_working_mediastinum, m_working_mediastinum->GetLargestPossibleRegion()); + iter.SetFirstDirection(0); + iter.SetSecondDirection(1); + iter.GoToBegin(); + int i=0; + MaskImageType::PointType A; + MaskImageType::PointType B; + MaskImageType::PointType C; + while (!iter.IsAtEnd()) { + A = leftPoints[i]; + B = rightPoints[i]; + C = A; + C[1] -= 10; // I know I must keep this point + double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]); + bool isPositive = s<0; + while (!iter.IsAtEndOfSlice()) { + while (!iter.IsAtEndOfLine()) { + // Very slow, I know ... but image should be very small + m_working_mediastinum->TransformIndexToPhysicalPoint(iter.GetIndex(), C); + double s = (B[0] - A[0]) * (C[1] - A[1]) - (B[1] - A[1]) * (C[0] - A[0]); + if (s == 0) iter.Set(2); + if (isPositive) { + if (s > 0) iter.Set(GetBackgroundValue()); + } + else { + if (s < 0) iter.Set(GetBackgroundValue()); + } + ++iter; + } + iter.NextLine(); + } + iter.NextSlice(); + ++i; + } + + //----------------------------------------------------- + //----------------------------------------------------- + // StartNewStep("Anterior limits"); + + + // MISSING FROM NOW + + // Station 4R, Station 4L, the right pulmonary artery, and/or the + // left superior pulmonary vein + + + //----------------------------------------------------- + //----------------------------------------------------- + // ALSO SUBSTRACT ARTERY/VEIN (initially in the support) + + + // Set output image information (required) + MaskImagePointer outputImage = this->GetOutput(0); + outputImage->SetRegions(m_working_mediastinum->GetLargestPossibleRegion()); + outputImage->SetOrigin(m_working_mediastinum->GetOrigin()); + outputImage->SetRequestedRegion(m_working_mediastinum->GetLargestPossibleRegion()); + +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +GenerateInputRequestedRegion() { + DD("GenerateInputRequestedRegion"); + // Call default + // Superclass::GenerateInputRequestedRegion(); + // Following needed because output region can be greater than input (trachea) + //ImagePointer mediastinum = dynamic_cast(itk::ProcessObject::GetInput(0)); + //ImagePointer trachea = dynamic_cast(itk::ProcessObject::GetInput(1)); + DD("set reg"); + m_mediastinum->SetRequestedRegion(m_mediastinum->GetLargestPossibleRegion()); + m_trachea->SetRequestedRegion(m_trachea->GetLargestPossibleRegion()); + DD("end"); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +GenerateData() { + DD("GenerateData"); + // Final Step -> graft output (if SetNthOutput => redo) + this->GraftOutput(m_output); +} +//-------------------------------------------------------------------- + + +#endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX diff --git a/segmentation/clitkExtractLymphStationsFilter.txx b/segmentation/clitkExtractLymphStationsFilter.txx index 801aeeb..1823e99 100644 --- a/segmentation/clitkExtractLymphStationsFilter.txx +++ b/segmentation/clitkExtractLymphStationsFilter.txx @@ -3,7 +3,7 @@ Authors belong to: - University of LYON http://www.universite-lyon.fr/ - - Léon Bérard cancer center http://www.centreleonberard.fr + - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr This software is distributed WITHOUT ANY WARRANTY; without even @@ -14,7 +14,7 @@ - BSD See included LICENSE.txt file - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html - ===========================================================================**/ + ======================================================================-====*/ #ifndef CLITKEXTRACTLYMPHSTATIONSFILTER_TXX #define CLITKEXTRACTLYMPHSTATIONSFILTER_TXX @@ -55,6 +55,7 @@ ExtractLymphStationsFilter(): // Default values ExtractStation_8_SetDefaultValues(); ExtractStation_3P_SetDefaultValues(); + ExtractStation_2RL_SetDefaultValues(); ExtractStation_3A_SetDefaultValues(); ExtractStation_7_SetDefaultValues(); } @@ -83,6 +84,12 @@ GenerateOutputInformation() { ExtractStation_3P(); StopSubStep(); + // Extract Station2RL + StartNewStep("Station 2RL"); + StartSubStep(); + ExtractStation_2RL(); + StopSubStep(); + // Extract Station3A StartNewStep("Station 3A"); StartSubStep(); @@ -102,20 +109,6 @@ GenerateOutputInformation() { //ExtractStation_4RL(); StopSubStep(); } - - - // - // typedef clitk::BooleanOperatorLabelImageFilter BFilter; - //BFilter::Pointer merge = BFilter::New(); - // writeImage(m_Output, "ouput.mhd"); - //writeImage(m_Working_Support, "ws.mhd"); - /*merge->SetInput1(m_Station7); - merge->SetInput2(m_Station4RL); // support - merge->SetOperationType(BFilter::AndNot); CHANGE OPERATOR - merge->SetForegroundValue(4); - merge->Update(); - m_Output = merge->GetOutput(); - */ } //-------------------------------------------------------------------- @@ -275,6 +268,7 @@ ExtractStation_3A() if (CheckForStation("3A")) { ExtractStation_3A_SI_Limits(); ExtractStation_3A_Ant_Limits(); + ExtractStation_3A_Post_Limits(); // Store image filenames into AFDB writeImage(m_ListOfStations["3A"], "seg/Station3A.mhd"); GetAFDB()->SetImageFilename("Station3A", "seg/Station3A.mhd"); @@ -284,6 +278,32 @@ ExtractStation_3A() //-------------------------------------------------------------------- +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +ExtractStation_2RL() +{ + if (CheckForStation("2RL")) { + ExtractStation_2RL_SI_Limits(); + ExtractStation_2RL_Post_Limits(); + ExtractStation_2RL_Ant_Limits2(); + ExtractStation_2RL_Ant_Limits(); + ExtractStation_2RL_LR_Limits(); + ExtractStation_2RL_Remove_Structures(); + ExtractStation_2RL_SeparateRL(); + + // Store image filenames into AFDB + writeImage(m_ListOfStations["2R"], "seg/Station2R.mhd"); + writeImage(m_ListOfStations["2L"], "seg/Station2L.mhd"); + GetAFDB()->SetImageFilename("Station2R", "seg/Station2R.mhd"); + GetAFDB()->SetImageFilename("Station2L", "seg/Station2L.mhd"); + WriteAFDB(); + } +} +//-------------------------------------------------------------------- + + //-------------------------------------------------------------------- template void @@ -305,10 +325,10 @@ ExtractStation_4RL() { template void clitk::ExtractLymphStationsFilter:: -Remove_Structures(std::string s) +Remove_Structures(std::string station, std::string s) { try { - StartNewStep("[Station7] Remove "+s); + StartNewStep("[Station"+station+"] Remove "+s); MaskImagePointer Structure = GetAFDB()->template GetImage(s); clitk::AndNot(m_Working_Support, Structure, GetBackgroundValue()); } @@ -319,5 +339,38 @@ Remove_Structures(std::string s) //-------------------------------------------------------------------- +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +SetFuzzyThreshold(std::string station, std::string tag, double value) +{ + m_FuzzyThreshold[station][tag] = value; +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +double +clitk::ExtractLymphStationsFilter:: +GetFuzzyThreshold(std::string station, std::string tag) +{ + if (m_FuzzyThreshold.find(station) == m_FuzzyThreshold.end()) { + clitkExceptionMacro("Could not find options for station "+station+" in the list (while searching for tag "+tag+")."); + return 0.0; + } + + if (m_FuzzyThreshold[station].find(tag) == m_FuzzyThreshold[station].end()) { + clitkExceptionMacro("Could not find options "+tag+" in the list of FuzzyThreshold for station "+station+"."); + return 0.0; + } + + return m_FuzzyThreshold[station][tag]; +} +//-------------------------------------------------------------------- + + + #endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX diff --git a/segmentation/clitkExtractLymphStationsFilter.txx.save b/segmentation/clitkExtractLymphStationsFilter.txx.save new file mode 100644 index 0000000..a7b3244 --- /dev/null +++ b/segmentation/clitkExtractLymphStationsFilter.txx.save @@ -0,0 +1,289 @@ +/*========================================================================= + Program: vv http://www.creatis.insa-lyon.fr/rio/vv + + Authors belong to: + - University of LYON http://www.universite-lyon.fr/ + - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr + - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the copyright notices for more information. + + It is distributed under dual licence + + - BSD See included LICENSE.txt file + - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html + ======================================================================-====*/ + +#ifndef CLITKEXTRACTLYMPHSTATIONSFILTER_TXX +#define CLITKEXTRACTLYMPHSTATIONSFILTER_TXX + +// clitk +#include "clitkCommon.h" +#include "clitkExtractLymphStationsFilter.h" +#include "clitkAddRelativePositionConstraintToLabelImageFilter.h" +#include "clitkSegmentationUtils.h" +#include "clitkAutoCropFilter.h" +#include "clitkSegmentationUtils.h" +#include "clitkSliceBySliceRelativePositionFilter.h" + +// itk +#include +#include +#include +#include +#include +#include + +// itk ENST +#include "RelativePositionPropImageFilter.h" + +//-------------------------------------------------------------------- +template +clitk::ExtractLymphStationsFilter:: +ExtractLymphStationsFilter(): + clitk::FilterBase(), + clitk::FilterWithAnatomicalFeatureDatabaseManagement(), + itk::ImageToImageFilter() +{ + this->SetNumberOfRequiredInputs(1); + SetBackgroundValue(0); + SetForegroundValue(1); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +GenerateOutputInformation() { + // Superclass::GenerateOutputInformation(); + + // Get input + LoadAFDB(); + m_mediastinum = GetAFDB()->template GetImage ("mediastinum"); + m_trachea = GetAFDB()->template GetImage ("trachea"); + ImagePointType carina; + GetAFDB()->GetPoint3D("carina", carina); + DD(carina); + m_CarinaZPositionInMM = carina[2]; + DD(m_CarinaZPositionInMM); + ImagePointType mlb; + GetAFDB()->GetPoint3D("middleLobeBronchus", mlb); + m_MiddleLobeBronchusZPositionInMM = mlb[2]; + DD(m_MiddleLobeBronchusZPositionInMM); + + // ---------------------------------------------------------------- + // ---------------------------------------------------------------- + // Superior limit = carina + // Inferior limit = origin right middle lobe bronchus + StartNewStep("Inf/Sup mediastinum limits with carina/bronchus"); + ImageRegionType region = m_mediastinum->GetLargestPossibleRegion(); DD(region); + ImageSizeType size = region.GetSize(); + ImageIndexType index = region.GetIndex(); + DD(m_CarinaZPositionInMM); + DD(m_MiddleLobeBronchusZPositionInMM); + index[2] = floor((m_MiddleLobeBronchusZPositionInMM - m_mediastinum->GetOrigin()[2]) / m_mediastinum->GetSpacing()[2]); + size[2] = 1+ceil((m_CarinaZPositionInMM-m_MiddleLobeBronchusZPositionInMM) / m_mediastinum->GetSpacing()[2]); // +1 to + region.SetSize(size); + region.SetIndex(index); + DD(region); + typedef itk::RegionOfInterestImageFilter CropFilterType; + typename CropFilterType::Pointer cropFilter = CropFilterType::New(); + cropFilter->SetInput(m_mediastinum); + cropFilter->SetRegionOfInterest(region); + cropFilter->Update(); + m_working_image = cropFilter->GetOutput(); + // Auto Crop (because following rel pos is faster) + m_working_image = clitk::AutoCrop(m_working_image, 0); + StopCurrentStep(m_working_image); + m_output = m_working_image; + + // ---------------------------------------------------------------- + // ---------------------------------------------------------------- + // Separate trachea in two CCL + StartNewStep("Separate trachea under carina"); + // DD(region); + ImageRegionType trachea_region = m_trachea->GetLargestPossibleRegion(); + for(int i=0; i<3; i++) { + index[i] = floor(((index[i]*m_mediastinum->GetSpacing()[i])+m_mediastinum->GetOrigin()[i] + -m_trachea->GetOrigin()[i])/m_trachea->GetSpacing()[i]); + // DD(index[i]); + size[i] = ceil((size[i]*m_mediastinum->GetSpacing()[i])/m_trachea->GetSpacing()[i]); + // DD(size[i]); + if (index[i] < 0) { + size[i] += index[i]; + index[i] = 0; + } + if (size[i]+index[i] > (trachea_region.GetSize()[i] + trachea_region.GetIndex()[i])) { + size[i] = trachea_region.GetSize()[i] + trachea_region.GetIndex()[i] - index[i]; + } + } + // DD(index); + // DD(size); + region.SetIndex(index); + region.SetSize(size); + // typedef itk::RegionOfInterestImageFilter CropFilterType; + // typename CropFilterType::Pointer + cropFilter = CropFilterType::New(); + // m_trachea.Print(std::cout); + cropFilter->SetInput(m_trachea); + cropFilter->SetRegionOfInterest(region); + cropFilter->Update(); + m_working_trachea = cropFilter->GetOutput(); + + // Labelize and consider two main labels + m_working_trachea = Labelize(m_working_trachea, 0, true, 1); + + // Detect wich label is at Left + typedef itk::ImageSliceConstIteratorWithIndex SliceIteratorType; + SliceIteratorType iter(m_working_trachea, m_working_trachea->GetLargestPossibleRegion()); + iter.SetFirstDirection(0); + iter.SetSecondDirection(1); + iter.GoToBegin(); + bool stop = false; + ImagePixelType leftLabel; + ImagePixelType rightLabel; + while (!stop) { + if (iter.Get() != GetBackgroundValue()) { + // DD(iter.GetIndex()); + // DD((int)iter.Get()); + leftLabel = iter.Get(); + stop = true; + } + ++iter; + } + if (leftLabel == 1) rightLabel = 2; + else rightLabel = 1; + DD((int)leftLabel); + DD((int)rightLabel); + + // End step + StopCurrentStep(m_working_trachea); + + //----------------------------------------------------- + /* DD("TEST Skeleton"); + typedef itk::BinaryThinningImageFilter SkeletonFilterType; + typename SkeletonFilterType::Pointer skeletonFilter = SkeletonFilterType::New(); + skeletonFilter->SetInput(m_working_trachea); + skeletonFilter->Update(); + writeImage(skeletonFilter->GetOutput(), "skel.mhd"); + writeImage(skeletonFilter->GetThinning(), "skel2.mhd"); + */ + + //----------------------------------------------------- + StartNewStep("Left limits with bronchus (slice by slice)"); + // Select LeftLabel (set right label to 0) + MaskImagePointer temp = SetBackground(m_working_trachea, m_working_trachea, rightLabel, 0); + writeImage(temp, "temp1.mhd"); + + m_working_image = + clitk::SliceBySliceRelativePosition(m_working_image, + temp, + 2, GetFuzzyThreshold(), "RightTo"); + /* + typedef clitk::SliceBySliceRelativePositionFilter SliceRelPosFilterType; + typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New(); + sliceRelPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId()); + sliceRelPosFilter->VerboseStepOn(); + sliceRelPosFilter->WriteStepOn(); + sliceRelPosFilter->SetInput(m_working_image); + sliceRelPosFilter->SetInputObject(temp); + sliceRelPosFilter->SetDirection(2); + sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThreshold()); + sliceRelPosFilter->SetOrientationTypeString("RightTo"); + sliceRelPosFilter->Update(); + m_working_image = sliceRelPosFilter->GetOutput();*/ + writeImage(m_working_image, "afterslicebyslice.mhd"); + + + typedef clitk::SliceBySliceRelativePositionFilter SliceRelPosFilterType; + typename SliceRelPosFilterType::Pointer sliceRelPosFilter = SliceRelPosFilterType::New(); + + + //----------------------------------------------------- + StartNewStep("Right limits with bronchus (slice by slice)"); + // Select LeftLabel (set right label to 0) + temp = SetBackground(m_working_trachea, m_working_trachea, leftLabel, 0); + writeImage(temp, "temp2.mhd"); + + sliceRelPosFilter = SliceRelPosFilterType::New(); + sliceRelPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId()); + sliceRelPosFilter->VerboseStepOff(); + sliceRelPosFilter->WriteStepOff(); + sliceRelPosFilter->SetInput(m_working_image); + sliceRelPosFilter->SetInputObject(temp); + sliceRelPosFilter->SetDirection(2); + sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThreshold()); + sliceRelPosFilter->SetOrientationTypeString("LeftTo"); + sliceRelPosFilter->Update(); + m_working_image = sliceRelPosFilter->GetOutput(); + writeImage(m_working_image, "afterslicebyslice.mhd"); + DD("end"); + m_output = m_working_image; + StopCurrentStep(m_output); + + //----------------------------------------------------- + StartNewStep("Trial Post position according to trachea"); + // Post: do not extend past the post wall of the 2 main bronchi + sliceRelPosFilter = SliceRelPosFilterType::New(); + sliceRelPosFilter->SetCurrentStepBaseId(this->GetCurrentStepId()); + sliceRelPosFilter->VerboseStepOn(); + sliceRelPosFilter->WriteStepOn(); + sliceRelPosFilter->SetInput(m_output); + sliceRelPosFilter->SetInputObject(m_trachea); + sliceRelPosFilter->SetDirection(2); + sliceRelPosFilter->SetFuzzyThreshold(GetFuzzyThreshold()); + // It means "not PostTo" (different from AntTo) + sliceRelPosFilter->NotFlagOn(); + //sliceRelPosFilter->SetOrientationType(SliceRelPosFilterType::RelPosFilterType::PostTo); + sliceRelPosFilter->SetOrientationTypeString("PostTo"); + sliceRelPosFilter->Update(); + m_output = sliceRelPosFilter->GetOutput(); + writeImage(m_output, "postrelpos.mhd"); + + + // Set output image information (required) + MaskImagePointer outputImage = this->GetOutput(0); + outputImage->SetRegions(m_working_image->GetLargestPossibleRegion()); + outputImage->SetOrigin(m_working_image->GetOrigin()); + outputImage->SetRequestedRegion(m_working_image->GetLargestPossibleRegion()); + DD("end2"); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +GenerateInputRequestedRegion() { + DD("GenerateInputRequestedRegion"); + // Call default + // Superclass::GenerateInputRequestedRegion(); + // Following needed because output region can be greater than input (trachea) + //ImagePointer mediastinum = dynamic_cast(itk::ProcessObject::GetInput(0)); + //ImagePointer trachea = dynamic_cast(itk::ProcessObject::GetInput(1)); + DD("set reg"); + m_mediastinum->SetRequestedRegion(m_mediastinum->GetLargestPossibleRegion()); + m_trachea->SetRequestedRegion(m_trachea->GetLargestPossibleRegion()); + DD("end"); +} +//-------------------------------------------------------------------- + + +//-------------------------------------------------------------------- +template +void +clitk::ExtractLymphStationsFilter:: +GenerateData() { + DD("GenerateData"); + // Final Step -> graft output (if SetNthOutput => redo) + this->GraftOutput(m_output); +} +//-------------------------------------------------------------------- + + +#endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX diff --git a/segmentation/clitkExtractLymphStationsGenericFilter.h b/segmentation/clitkExtractLymphStationsGenericFilter.h index 389e68d..afeb66b 100644 --- a/segmentation/clitkExtractLymphStationsGenericFilter.h +++ b/segmentation/clitkExtractLymphStationsGenericFilter.h @@ -3,7 +3,7 @@ Authors belong to: - University of LYON http://www.universite-lyon.fr/ - - Léon Bérard cancer center http://www.centreleonberard.fr + - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr This software is distributed WITHOUT ANY WARRANTY; without even @@ -14,7 +14,7 @@ - BSD See included LICENSE.txt file - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html -===========================================================================**/ +======================================================================-====*/ #ifndef CLITKEXTRACTLYMPHSTATIONSSGENERICFILTER_H #define CLITKEXTRACTLYMPHSTATIONSSGENERICFILTER_H diff --git a/segmentation/clitkExtractLymphStationsGenericFilter.txx b/segmentation/clitkExtractLymphStationsGenericFilter.txx index 977cf93..ffa6c3c 100644 --- a/segmentation/clitkExtractLymphStationsGenericFilter.txx +++ b/segmentation/clitkExtractLymphStationsGenericFilter.txx @@ -3,7 +3,7 @@ Authors belong to: - University of LYON http://www.universite-lyon.fr/ - - Léon Bérard cancer center http://www.centreleonberard.fr + - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr This software is distributed WITHOUT ANY WARRANTY; without even @@ -14,7 +14,7 @@ - BSD See included LICENSE.txt file - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html - ===========================================================================**/ + ======================================================================-====*/ #ifndef CLITKEXTRACTLYMPHSTATIONSSGENERICFILTER_TXX #define CLITKEXTRACTLYMPHSTATIONSSGENERICFILTER_TXX @@ -68,8 +68,10 @@ SetOptionsFromArgsInfoToFilter(FilterType * f) f->SetWriteStepFlag(mArgsInfo.writeStep_flag); f->SetVerboseMemoryFlag(mArgsInfo.verboseMemory_flag); f->SetAFDBFilename(mArgsInfo.afdb_arg); + + // Station 8 f->SetDistanceMaxToAnteriorPartOfTheSpine(mArgsInfo.maxAntSpine_arg); - f->SetFuzzyThresholdForS8(mArgsInfo.fuzzyThresholdForS8_arg); + f->SetFuzzyThreshold("8", "Esophagus", mArgsInfo.tS8_Esophagus_arg); f->SetInjectedThresholdForS8(mArgsInfo.injectedThresholdForS8_arg); // Check multiple options for radius dilatation @@ -110,12 +112,24 @@ SetOptionsFromArgsInfoToFilter(FilterType * f) f->AddComputeStation(mArgsInfo.station_arg[i]); // Station 7 - f->SetFuzzyThresholdForS7("Bronchi", mArgsInfo.tS7_Bronchi_arg); - f->SetFuzzyThresholdForS7("LeftSuperiorPulmonaryVein", mArgsInfo.tS7_LeftSuperiorPulmonaryVein_arg); - f->SetFuzzyThresholdForS7("RightSuperiorPulmonaryVein", mArgsInfo.tS7_RightSuperiorPulmonaryVein_arg); - f->SetFuzzyThresholdForS7("RightPulmonaryArtery", mArgsInfo.tS7_RightPulmonaryArtery_arg); - f->SetFuzzyThresholdForS7("LeftPulmonaryArtery", mArgsInfo.tS7_LeftPulmonaryArtery_arg); - f->SetFuzzyThresholdForS7("SVC", mArgsInfo.tS7_SVC_arg); + f->SetFuzzyThreshold("7", "Bronchi", mArgsInfo.tS7_Bronchi_arg); + f->SetFuzzyThreshold("7", "LeftSuperiorPulmonaryVein", mArgsInfo.tS7_LeftSuperiorPulmonaryVein_arg); + f->SetFuzzyThreshold("7", "RightSuperiorPulmonaryVein", mArgsInfo.tS7_RightSuperiorPulmonaryVein_arg); + f->SetFuzzyThreshold("7", "RightPulmonaryArtery", mArgsInfo.tS7_RightPulmonaryArtery_arg); + f->SetFuzzyThreshold("7", "LeftPulmonaryArtery", mArgsInfo.tS7_LeftPulmonaryArtery_arg); + f->SetFuzzyThreshold("7", "SVC", mArgsInfo.tS7_SVC_arg); + + // Station 3A + f->SetFuzzyThreshold("3A", "Sternum", mArgsInfo.tS3A_Sternum_arg); + f->SetFuzzyThreshold("3A", "SubclavianArtery", mArgsInfo.tS3A_SubclavianArtery_arg); + + // Station 2RL + f->SetFuzzyThreshold("2RL", "CommonCarotidArtery", mArgsInfo.tS2RL_CommonCarotidArtery_arg); + f->SetFuzzyThreshold("2RL", "BrachioCephalicTrunk", mArgsInfo.tS2RL_BrachioCephalicTrunk_arg); + f->SetFuzzyThreshold("2RL", "BrachioCephalicVein", mArgsInfo.tS2RL_BrachioCephalicVein_arg); + f->SetFuzzyThreshold("2RL", "Aorta", mArgsInfo.tS2RL_Aorta_arg); + f->SetFuzzyThreshold("2RL", "SubclavianArteryLeft", mArgsInfo.tS2RL_SubclavianArteryLeft_arg); + f->SetFuzzyThreshold("2RL", "SubclavianArteryRight", mArgsInfo.tS2RL_SubclavianArteryRight_arg); } //-------------------------------------------------------------------- diff --git a/segmentation/clitkExtractPatientFilter.h b/segmentation/clitkExtractPatientFilter.h index 8ff1943..b612fd3 100644 --- a/segmentation/clitkExtractPatientFilter.h +++ b/segmentation/clitkExtractPatientFilter.h @@ -3,7 +3,7 @@ Authors belong to: - University of LYON http://www.universite-lyon.fr/ - - Léon Bérard cancer center http://www.centreleonberard.fr + - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr This software is distributed WITHOUT ANY WARRANTY; without even @@ -14,13 +14,14 @@ - BSD See included LICENSE.txt file - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html - ===========================================================================**/ + ======================================================================-====*/ #ifndef CLITKEXTRACTPATIENTFILTER_H #define CLITKEXTRACTPATIENTFILTER_H #include "clitkFilterBase.h" #include "clitkFilterWithAnatomicalFeatureDatabaseManagement.h" +#include "clitkSegmentationUtils.h" namespace clitk { diff --git a/segmentation/clitkReconstructThroughDilationImageFilter.txx b/segmentation/clitkReconstructThroughDilationImageFilter.txx old mode 100755 new mode 100644 index 5522766..613f10b --- a/segmentation/clitkReconstructThroughDilationImageFilter.txx +++ b/segmentation/clitkReconstructThroughDilationImageFilter.txx @@ -3,7 +3,7 @@ Authors belong to: - University of LYON http://www.universite-lyon.fr/ - - Léon Bérard cancer center http://www.centreleonberard.fr + - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr This software is distributed WITHOUT ANY WARRANTY; without even @@ -14,20 +14,10 @@ - BSD See included LICENSE.txt file - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html -===========================================================================**/ +======================================================================-====*/ #ifndef clitkReconstructThroughDilationImageFilter_txx #define clitkReconstructThroughDilationImageFilter_txx -/* ================================================= - * @file clitkReconstructThroughDilationImageFilter.txx - * @author - * @date - * - * @brief - * - ===================================================*/ - - namespace clitk { diff --git a/segmentation/clitkTestArtery.cxx b/segmentation/clitkTestArtery.cxx new file mode 100644 index 0000000..8e36523 --- /dev/null +++ b/segmentation/clitkTestArtery.cxx @@ -0,0 +1,39 @@ +/*========================================================================= + Program: vv http://www.creatis.insa-lyon.fr/rio/vv + + Authors belong to: + - University of LYON http://www.universite-lyon.fr/ + - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr + - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the copyright notices for more information. + + It is distributed under dual licence + + - BSD See included LICENSE.txt file + - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html +======================================================================-====*/ + +// clitk +#include "clitkCommon.h" + +//-------------------------------------------------------------------- +int main(int argc, char * argv[]) +{ + // Typedef + typedef uchar MaskImagePixelType; + typedef short ImagePixelType; + typedef itk::Image ImageType; + typedef itk::Image MaskImageType; + + // Input + + + + + + return EXIT_SUCCESS; +} // This is the end, my friend +//-------------------------------------------------------------------- diff --git a/segmentation/clitkTestFilter.cxx b/segmentation/clitkTestFilter.cxx new file mode 100644 index 0000000..64716ed --- /dev/null +++ b/segmentation/clitkTestFilter.cxx @@ -0,0 +1,523 @@ + +/*========================================================================= + Program: vv http://www.creatis.insa-lyon.fr/rio/vv + + Authors belong to: + - University of LYON http://www.universite-lyon.fr/ + - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr + - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the copyright notices for more information. + + It is distributed under dual licence + + - BSD See included LICENSE.txt file + - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html + ======================================================================-====*/ + +// clitk +#include "clitkTestFilter_ggo.h" +#include "clitkImageCommon.h" +#include "clitkBooleanOperatorLabelImageFilter.h" +#include "clitkAutoCropFilter.h" +//#include "clitkRelativePositionConstraintLabelImageFilter.h" +#include "clitkResampleImageWithOptionsFilter.h" +#include "clitkAddRelativePositionConstraintToLabelImageFilter.h" +#include "clitkExtractLungFilter.h" +#include "clitkExtractPatientFilter.h" +#include "clitkExtractMediastinumFilter.h" +//#include "clitkTestStation7.h" +#include "clitkSegmentationUtils.h" +#include "clitkMorphoMathFilter.h" + +// ITK ENST +#include "RelativePositionPropImageFilter.h" + +// VTK +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + +template +class comparePointsX { +public: + bool operator() (PointType i, PointType j) { return (i[0] +void HypercubeCorners(std::vector > & out) { + std::vector > previous; + HypercubeCorners(previous); + out.resize(previous.size()*2); + for(uint i=0; i p; + if (i +void HypercubeCorners<1>(std::vector > & out) { + out.resize(2); + out[0][0] = 0; + out[1][0] = 1; +} + + +//-------------------------------------------------------------------- +int main(int argc, char * argv[]) { + + // Init command line + GGO(clitkTestFilter, args_info); + + // Image types + //typedef unsigned char PixelType; + typedef short PixelType; + static const int Dim=3; + typedef itk::Image InputImageType; + + // Read inputs + InputImageType::Pointer input1; + InputImageType::Pointer input2; + InputImageType::Pointer input3; + if (args_info.input1_given) input1 = clitk::readImage(args_info.input1_arg, true); + if (args_info.input2_given) input2 = clitk::readImage(args_info.input2_arg, true); + if (args_info.input3_given) input3 = clitk::readImage(args_info.input3_arg, true); + + // Declare output + InputImageType::Pointer output; + + //-------------------------------------------------------------------- + // Filter test BooleanOperatorLabelImageFilter + if (0) { + typedef clitk::BooleanOperatorLabelImageFilter FilterType; + FilterType::Pointer filter = FilterType::New(); + filter->SetInput1(input1); + filter->SetInput2(input2); + output = clitk::NewImageLike(input1); + filter->GraftOutput(output); /// TO VERIFY !!! + filter->Update(); + filter->SetInput2(input3); + filter->Update(); + output = filter->GetOutput(); + clitk::writeImage(output, args_info.output_arg); + } + + //-------------------------------------------------------------------- + // Filter test AutoCropLabelImageFilter + if (0) { + typedef clitk::AutoCropFilter FilterType; + FilterType::Pointer filter = FilterType::New(); + filter->SetInput(input1); + filter->Update(); + output = filter->GetOutput(); + clitk::writeImage(output, args_info.output_arg); + } + + //-------------------------------------------------------------------- + // Filter test RelativePositionPropImageFilter + if (0) { + typedef itk::Image OutputImageType; + OutputImageType::Pointer outputf; + typedef itk::RelativePositionPropImageFilter FilterType; + FilterType::Pointer filter = FilterType::New(); + filter->SetInput(input1); + + filter->SetAlpha1(clitk::deg2rad(args_info.angle1_arg)); // xy plane + filter->SetAlpha2(clitk::deg2rad(args_info.angle2_arg)); + filter->SetK1(M_PI/2.0); // Opening parameter, default = pi/2 + filter->SetFast(true); + filter->SetRadius(2); + filter->SetVerboseProgress(true); + + /* A1 A2 + Left 0 0 + Right 180 0 + Ant 90 0 + Post -90 0 + Inf 0 90 + Sup 0 -90 + */ + + filter->Update(); + clitk::writeImage(filter->GetOutput(), args_info.output_arg); + } + + //-------------------------------------------------------------------- + // Filter test + if (0) { + typedef itk::Image OutputImageType; + typedef clitk::ResampleImageWithOptionsFilter FilterType; + FilterType::Pointer filter = FilterType::New(); + filter->SetInput(input1); + filter->SetOutputIsoSpacing(1); + //filter->SetNumberOfThreads(4); // auto + filter->SetGaussianFilteringEnabled(false); + filter->Update(); + clitk::writeImage(filter->GetOutput(), args_info.output_arg); + } + + //-------------------------------------------------------------------- + // Filter AddRelativePositionConstraintToLabelImageFilter test + if (0) { + /* + typedef clitk::AddRelativePositionConstraintToLabelImageFilter FilterType; + FilterType::Pointer filter = FilterType::New(); + filter->SetInput(input1); + filter->SetInputObject(input2); + filter->SetOrientationType(FilterType::LeftTo); + filter->SetIntermediateSpacing(5); + filter->SetFuzzyThreshold(0.5); + filter->VerboseStepOn(); + filter->WriteStepOff(); + filter->Update(); + + filter->SetInput(filter->GetOutput()); + filter->SetInputObject(input2); + filter->SetOrientationType(FilterType::RightTo); + filter->SetIntermediateSpacing(5); + filter->SetFuzzyThreshold(0.5); + filter->Update(); + + clitk::writeImage(filter->GetOutput(), args_info.output_arg); + */ + } + + //-------------------------------------------------------------------- + // Filter test ExtractPatientFilter + if (0) { + /* + typedef itk::Image OutputImageType; + typedef clitk::ExtractPatientFilter FilterType; + FilterType::Pointer filter = FilterType::New(); + filter->SetInput(input1); + filter->VerboseStepOn(); + filter->WriteStepOn(); + // options (default) + filter->SetUpperThreshold(-300); + filter->FinalOpenCloseOff(); // default=on (not rezally needed ?) + filter->Update(); + OutputImageType::Pointer output = filter->GetOutput(); + clitk::writeImage(output, args_info.output_arg); + */ + } + + //-------------------------------------------------------------------- + // Filter test ExtractLungsFilter + if (0) { + /* + typedef itk::Image OutputImageType; // to change into char + typedef clitk::ExtractLungFilter FilterType; + FilterType::Pointer filter = FilterType::New(); + // DD(filter->GetNumberOfSteps()); + filter->SetInput(input1); // CT + filter->SetInputPatientMask(input2, 0); // Patient mask and BG value + filter->VerboseStepOn(); + filter->WriteStepOn(); + // options (default) + //filter->SetBackgroundValue(0); + filter->SetUpperThreshold(-300); + // filter->SetMinimumComponentSize(100); + + filter->Update(); + OutputImageType::Pointer output = filter->GetOutput(); + clitk::writeImage(output, args_info.output_arg); + */ + } + + //-------------------------------------------------------------------- + // Filter test ExtractMediastinumFilter + if (0) { + /* + typedef clitk::ExtractMediastinumFilter FilterType; + FilterType::Pointer filter = FilterType::New(); + filter->SetInputPatientLabelImage(input1); + filter->SetInputLungLabelImage(input2, 0, 1, 2); // BG, FG Left Lung, FG Right Lung + filter->SetInputBonesLabelImage(input3); + filter->VerboseStepOn(); + filter->WriteStepOn(); + filter->Update(); + output = filter->GetOutput(); + clitk::writeImage(output, args_info.output_arg); + */ + } + + //-------------------------------------------------------------------- + // Test for auto register sub-task in a segmentation process + if (0) { + // ExtractLymphStation_7 * s7 = new ExtractLymphStation_7; + // s7->SetArgsInfo(args_info); + // GetParent->SetArgsInfo<> + //s7->StartSegmentation(); + } + + //-------------------------------------------------------------------- + // Test for biinary image from a contour set + if (0) { + DD("here"); + + // Type of a slice + typedef itk::Image SliceType; + + // Build the list of slices + std::vector slices; + clitk::ExtractSlices(input1, 2, slices); + DD(slices.size()); + + // HOW TO DO SEVERAL BY SLICE !!! not a map ? + + // Compute centroids 3D centroids by slices + int BG = 0; + std::map > centroids3D; + for(uint i=0; i(slices[i], BG, true, 1); + // ComputeCentroids + std::vector temp; + clitk::ComputeCentroids(slices[i], BG, temp); + for(uint j=1; j::Convert2DTo3D(temp[j], input1, i, a); + centroids3D[i].push_back(a); + } + } + + // REPRENDRE POUR TOUT FAIRE BY SLICE (pas de i); + + // Take a given slice i=29 + int index=29; + SliceType::Pointer slice = slices[index]; + std::vector points = centroids3D[index]; + + // Sort points from R to L + std::sort(points.begin(), points.end(), + comparePointsX()); + + // Slice corners (quel sens ?) + + // Compute min and max coordinates + const uint dim=3; + typedef InputImageType::PointType PointType; + typedef InputImageType::IndexType IndexType; + PointType min_c, max_c; + IndexType min_i, max_i; + min_i = input1->GetLargestPossibleRegion().GetIndex(); + for(uint i=0; iGetLargestPossibleRegion().GetSize()[i] + min_i[i]; + input1->TransformIndexToPhysicalPoint(min_i, min_c); + input1->TransformIndexToPhysicalPoint(max_i, max_c); + + // Compute the corners coordinates + std::vector > l; + HypercubeCorners(l); + for(uint i=0; i mesh = vtkSmartPointer::New(); + mesh->Allocate(); //for cell structures + mesh->SetPoints(vtkPoints::New()); + vtkIdType ids[2]; + int point_number = points.size(); + for (unsigned int i=0; iGetPoints()->InsertNextPoint(points[i][0],points[i][1],points[i][2]); + ids[0]=i; + ids[1]=(ids[0]+1)%point_number; //0-1,1-2,...,n-1-0 + mesh->GetLines()->InsertNextCell(2,ids); + } + + vtkSmartPointer w = vtkSmartPointer::New(); + w->SetInput(mesh); + w->SetFileName("bidon.vtk"); + w->Write(); + + // binarize + DD("binarize"); + double *bounds=mesh->GetBounds(); + DDV(bounds, 6); + + DD("create image"); + vtkSmartPointer binary_image=vtkSmartPointer::New(); + binary_image->SetScalarTypeToUnsignedChar(); + ///Use the smallest mask in which the mesh fits + // Add two voxels on each side to make sure the mesh fits + // double * samp_origin= + InputImageType::PointType samp_origin = input1->GetOrigin(); + // double * + InputImageType::SpacingType spacing=input1->GetSpacing(); + double * spacing2 = new double[3]; + spacing2[0] = spacing[0]; + spacing2[1] = spacing[1]; + spacing2[2] = spacing[2]; + DD(spacing2[0]); + binary_image->SetSpacing(spacing2); + /// Put the origin on a voxel to avoid small skips + DD(floor((bounds[0]-samp_origin[0])/spacing[0]-2)*spacing[0]); + DD(bounds[0]); + DD(samp_origin[0]); + DD(spacing[0]); + binary_image->SetOrigin(//samp_origin[0], samp_origin[1], samp_origin[2]); + floor((bounds[0]-samp_origin[0])/spacing[0]-2)*spacing[0]+samp_origin[0], + floor((bounds[2]-samp_origin[1])/spacing[1]-2)*spacing[1]+samp_origin[1], + floor((bounds[4]-samp_origin[2])/spacing[2]-2)*spacing[2]+samp_origin[2]); + double * origin=binary_image->GetOrigin(); + binary_image->SetExtent(0,ceil((bounds[1]-origin[0])/spacing[0]+4), // Joel used +4 here (?) + 0,ceil((bounds[3]-origin[1])/spacing[1]+4), + 0,ceil((bounds[5]-origin[2])/spacing[2]+4)); + binary_image->AllocateScalars(); + memset(binary_image->GetScalarPointer(),0, + binary_image->GetDimensions()[0]*binary_image->GetDimensions()[1]*binary_image->GetDimensions()[2]*sizeof(unsigned char)); + + vtkSmartPointer sts=vtkSmartPointer::New(); + //The following line is extremely important + //http://www.nabble.com/Bug-in-vtkPolyDataToImageStencil--td23368312.html#a23370933 + sts->SetTolerance(0); + sts->SetInformationInput(binary_image); + + bool extrude=true; + if (extrude) { + vtkSmartPointer extrude=vtkSmartPointer::New(); + extrude->SetInput(mesh); + + /// NO ????We extrude in the -slice_spacing direction to respect the FOCAL convention ??? + + extrude->SetVector(0, 0, input1->GetSpacing()[2]);//slice_spacing); // 2* ? yes use a single one + sts->SetInput(extrude->GetOutput()); + } else + sts->SetInput(mesh); + + DD("stencil"); + vtkSmartPointer stencil=vtkSmartPointer::New(); + stencil->SetStencil(sts->GetOutput()); + stencil->SetInput(binary_image); + stencil->Update(); + + DD("write"); + vtkSmartPointer ww = vtkSmartPointer::New(); + ww->SetInput(stencil->GetOutput()); + ww->SetFileName("binary.mhd"); + ww->Write(); + } + //-------------------------------------------------------------------- + + //-------------------------------------------------------------------- + // Test for vessels ReconstructionByDilatation + if (1) { + // Read input CT (already croped) + // treshold 3D + // erode n times (or in 2D ?) + // slices extract + // SBS + // - CCL (keep mask) + // - for all CCL, ReconstructionByDilatation in initial mask + // joint for output + + // input1 + DD("binarize") + int BG = 0; + int FG = 1; + typedef itk::Image MaskImageType; + typedef itk::BinaryThresholdImageFilter BinarizeFilterType; + BinarizeFilterType::Pointer binarizeFilter = BinarizeFilterType::New(); + binarizeFilter->SetInput(input1); + binarizeFilter->SetLowerThreshold(150); + binarizeFilter->SetInsideValue(FG); + binarizeFilter->SetOutsideValue(BG); + binarizeFilter->Update(); + MaskImageType::Pointer mask = binarizeFilter->GetOutput(); + clitk::writeImage(mask, "m.mhd"); + + // Extract slices + typedef itk::Image SliceType; + std::vector slices_mask; + clitk::ExtractSlices(mask, 2, slices_mask); + DD(slices_mask.size()); + + std::vector debug_eroded; + std::vector debug_labeled; + + // Loop + for(uint i=0; i::Pointer f = clitk::MorphoMathFilter::New(); + f->SetInput(slices_mask[i]); + f->SetBackgroundValue(BG); + f->SetForegroundValue(FG); + f->SetRadius(1); + f->SetOperationType(0); // Erode + f->Update(); + SliceType::Pointer eroded = f->GetOutput(); + debug_eroded.push_back(eroded); + + // CCL + DD("CCL"); + SliceType::Pointer labeled = Labelize(slices_mask[i], 0, false, 10); + debug_labeled.push_back(labeled); + + // ReconstructionByDilatation + + + } // end loop + + MaskImageType::Pointer eroded = clitk::JoinSlices(debug_eroded, mask, 2); + clitk::writeImage(eroded, "eroded.mhd"); + + MaskImageType::Pointer labeled = clitk::JoinSlices(debug_labeled, mask, 2); + clitk::writeImage(labeled, "labeled.mhd"); + + } + //-------------------------------------------------------------------- + + // This is the end my friend + return EXIT_SUCCESS; +}// end main +//-------------------------------------------------------------------- diff --git a/segmentation/clitkTestFilter.ggo b/segmentation/clitkTestFilter.ggo new file mode 100644 index 0000000..6ab0ebb --- /dev/null +++ b/segmentation/clitkTestFilter.ggo @@ -0,0 +1,16 @@ +#File clitkTestFilter.ggo +package "clitkTestFilter" +version "1.0" +purpose "Test a filter" + +option "config" - "Config file" string no +option "verbose" v "Verbose" flag off + +option "input1" i "Input 1 image filename" string no +option "input2" j "Input 2 image filename" string no +option "input3" k "Input 3 image filename" string no +option "output" o "Output image filename" string no + +option "angle1" a "First angle (degree)" float default = "0" no +option "angle2" b "Second angle (degree)" float default = "0" no + diff --git a/tools/clitkCropImage.ggo b/tools/clitkCropImage.ggo old mode 100755 new mode 100644 index 4931b5b..e8650de --- a/tools/clitkCropImage.ggo +++ b/tools/clitkCropImage.ggo @@ -15,9 +15,9 @@ option "input" i "Input image filename" string yes option "output" o "Output image filename" string yes section "Used determined crop" -option "boundingBox" b "Bounding box of the crop region" int no multiple -option "lower" l "Size of the lower crop region" int no multiple -option "upper" u "Size of the upper crop region" int no multiple +option "boundingBox" b "Bounding box of the crop region (in 2D: =x1,y2, x2,y2)" int no multiple +option "lower" l "Size of the lower crop region (multiple values)" int no multiple +option "upper" u "Size of the upper crop region (multiple values)" int no multiple option "origin" - "Set new origin to zero" flag off section "AutoCrop with BG value"