clitkImageCommon.cxx
clitkTransformUtilities.cxx
clitkIO.cxx
+ clitkGateAsciiImageIO.cxx
+ clitkGateAsciiImageIOFactory.cxx
clitkVoxImageIO.cxx
clitkVoxImageIOFactory.cxx
clitkVfImageIO.cxx
vvImageWriter.cxx
)
-IF(UNIX)
- SET(clitkCommon_SRC
- ${clitkCommon_SRC}
- clitkGateAsciiImageIO.cxx
- clitkGateAsciiImageIOFactory.cxx
- )
-ENDIF(UNIX)
-
### Declare clitkCommon library
ADD_LIBRARY(clitkCommon STATIC ${clitkCommon_SRC})
===========================================================================**/
// std include
-#include <regex.h>
#include <cstdio>
#include <sstream>
#include <iostream>
std::ostream& operator<<(std::ostream& os, const clitk::GateAsciiImageIO::GateAsciiHeader& header)
{
- os << "matrix_size=[" << header.matrix_size[0] << "," << header.matrix_size[1] << "," << header.matrix_size[2] << "]" << endl;
- os << "resolution=[" << header.resolution[0] << "," << header.resolution[1] << "," << header.resolution[2] << "]" << endl;
- os << "voxel_size=[" << header.voxel_size[0] << "," << header.voxel_size[1] << "," << header.voxel_size[2] << "]" << endl;
- os << "nb_value=" << header.nb_value << endl;
- return os;
+ os << "matrix_size=[" << header.matrix_size[0] << "," << header.matrix_size[1] << "," << header.matrix_size[2] << "]" << endl;
+ os << "resolution=[" << header.resolution[0] << "," << header.resolution[1] << "," << header.resolution[2] << "]" << endl;
+ os << "voxel_size=[" << header.voxel_size[0] << "," << header.voxel_size[1] << "," << header.voxel_size[2] << "]" << endl;
+ os << "nb_value=" << header.nb_value << endl;
+ return os;
}
//--------------------------------------------------------------------
// Read Image Information
void clitk::GateAsciiImageIO::ReadImageInformation()
{
- FILE* handle = fopen(m_FileName.c_str(),"r");
- if (!handle) {
- itkGenericExceptionMacro(<< "Could not open file (for reading): " << m_FileName);
- return;
- }
-
- GateAsciiHeader header;
- if (!ReadHeader(handle,header)) {
- itkGenericExceptionMacro(<< "Could not read header: " << m_FileName);
- fclose(handle);
- return;
- }
+ FILE* handle = fopen(m_FileName.c_str(),"r");
+ if (!handle) {
+ itkGenericExceptionMacro(<< "Could not open file (for reading): " << m_FileName);
+ return;
+ }
+
+ GateAsciiHeader header;
+ if (!ReadHeader(handle,header)) {
+ itkGenericExceptionMacro(<< "Could not read header: " << m_FileName);
fclose(handle);
-
- int real_length = -1;
- double real_spacing = 0;
- for (int kk=0; kk<3; kk++) {
- if (header.resolution[kk]>1) {
- if (real_length>0) {
- itkGenericExceptionMacro(<< "Could not image dimension: " << m_FileName);
- return;
- }
- real_length = header.resolution[kk];
- real_spacing = header.voxel_size[kk];
- }
+ return;
+ }
+ fclose(handle);
+
+ int real_length = -1;
+ double real_spacing = 0;
+ for (int kk=0; kk<3; kk++) {
+ if (header.resolution[kk]>1) {
+ if (real_length>0) {
+ itkGenericExceptionMacro(<< "Could not image dimension: " << m_FileName);
+ return;
+ }
+ real_length = header.resolution[kk];
+ real_spacing = header.voxel_size[kk];
}
- assert(real_length == header.nb_value);
-
- // Set image information
- SetNumberOfDimensions(2);
- SetDimensions(0, real_length);
- SetDimensions(1, 1);
- SetSpacing(0, real_spacing);
- SetSpacing(1, 1);
- SetOrigin(0, 0);
- SetOrigin(1, 0);
- SetComponentType(itk::ImageIOBase::DOUBLE);
+ }
+ assert(real_length == header.nb_value);
+
+ // Set image information
+ SetNumberOfDimensions(2);
+ SetDimensions(0, real_length);
+ SetDimensions(1, 1);
+ SetSpacing(0, real_spacing);
+ SetSpacing(1, 1);
+ SetOrigin(0, 0);
+ SetOrigin(1, 0);
+ SetComponentType(itk::ImageIOBase::DOUBLE);
}
//--------------------------------------------------------------------
bool clitk::GateAsciiImageIO::CanReadFile(const char* FileNameToRead)
{
- std::string filename(FileNameToRead);
+ std::string filename(FileNameToRead);
- { // check extension
- std::string filenameext = GetExtension(filename);
- if (filenameext != "txt") return false;
- }
+ {
+ // check extension
+ std::string filenameext = GetExtension(filename);
+ if (filenameext != "txt") return false;
+ }
- { // check header
- FILE* handle = fopen(filename.c_str(),"r");
- if (!handle) return false;
+ {
+ // check header
+ FILE* handle = fopen(filename.c_str(),"r");
+ if (!handle) return false;
- GateAsciiHeader header;
- if (!ReadHeader(handle,header)) { fclose(handle); return false; }
- fclose(handle);
+ GateAsciiHeader header;
+ if (!ReadHeader(handle,header)) {
+ fclose(handle);
+ return false;
}
+ fclose(handle);
+ }
- return true;
+ return true;
}
//--------------------------------------------------------------------
// Read Line in file
bool clitk::GateAsciiImageIO::ReadLine(FILE* handle, std::string& line)
{
- std::stringstream stream;
- while (true)
- {
- char item;
- if (ferror(handle)) return false;
- if (fread(&item,1,1,handle) != 1) return false;
-
- if (item=='\n' or feof(handle)) {
- line = stream.str();
- return true;
- }
-
- stream << item;
+ std::stringstream stream;
+ while (true) {
+ char item;
+ if (ferror(handle)) return false;
+ if (fread(&item,1,1,handle) != 1) return false;
+
+ if (item=='\n' || feof(handle)) {
+ line = stream.str();
+ return true;
}
-}
-std::string ExtractMatch(const std::string& base, const regmatch_t& match)
-{
- return base.substr(match.rm_so,match.rm_eo-match.rm_so);
+ stream << item;
+ }
}
template <typename T>
T ConvertFromString(const std::string& value)
{
- std::stringstream stream;
- stream << value;
- T converted;
- stream >> converted;
- return converted;
+ std::stringstream stream;
+ stream << value;
+ T converted;
+ stream >> converted;
+ return converted;
+}
+
+bool
+clitk::GateAsciiImageIO::FindRegularExpressionNextLine(itksys::RegularExpression ®, std::string &s, FILE* handle)
+{
+ std::string line;
+ if(!ReadLine(handle,line)) return false;
+ if(!reg.compile(s.c_str())) return false;
+ return reg.find(line.c_str());
}
//--------------------------------------------------------------------
// Read Image Header
bool clitk::GateAsciiImageIO::ReadHeader(FILE* handle, GateAsciiHeader& header)
{
- assert(handle);
- std::string line;
-
- regex_t re_comment;
- regex_t re_matrix_size;
- regex_t re_resol;
- regex_t re_voxel_size;
- regex_t re_nb_value;
- regmatch_t matches[4];
-
- { // build regex
- assert(regcomp(&re_comment,"^#.+$",REG_EXTENDED|REG_NEWLINE) == 0);
- assert(regcomp(&re_matrix_size,"^# +Matrix *Size *= +\\(([0-9]+\\.?[0-9]*),([0-9]+\\.?[0-9]*),([0-9]+\\.?[0-9]*)\\)$",REG_EXTENDED|REG_NEWLINE) == 0);
- assert(regcomp(&re_resol,"^# +Resol *= +\\(([0-9]+),([0-9]+),([0-9]+)\\)$",REG_EXTENDED|REG_NEWLINE) == 0);
- assert(regcomp(&re_voxel_size,"^# +Voxel *Size *= +\\(([0-9]+\\.?[0-9]*),([0-9]+\\.?[0-9]*),([0-9]+\\.?[0-9]*)\\)$",REG_EXTENDED|REG_NEWLINE) == 0);
- assert(regcomp(&re_nb_value,"^# +nbVal *= +([0-9]+)$",REG_EXTENDED|REG_NEWLINE) == 0);
- }
+ assert(handle);
+
+ std::string regexstr[6] = {
+ "^#.+$",
+ "^# +Matrix *Size *= +\\(([0-9]+\\.?[0-9]*),([0-9]+\\.?[0-9]*),([0-9]+\\.?[0-9]*)\\)$",
+ "^# +Resol *= +\\(([0-9]+),([0-9]+),([0-9]+)\\)$",
+ "^# +Voxel *Size *= +\\(([0-9]+\\.?[0-9]*),([0-9]+\\.?[0-9]*),([0-9]+\\.?[0-9]*)\\)$",
+ "^# +nbVal *= +([0-9]+)$"
+ };
- if (!ReadLine(handle,line)) return false;
- if (regexec(&re_comment,line.c_str(),1,matches,0) == REG_NOMATCH) return false;
+ itksys::RegularExpression regex;
- if (!ReadLine(handle,line)) return false;
- if (regexec(&re_matrix_size,line.c_str(),4,matches,0) == REG_NOMATCH) return false;
- header.matrix_size[0] = ConvertFromString<double>(ExtractMatch(line,matches[1]));
- header.matrix_size[1] = ConvertFromString<double>(ExtractMatch(line,matches[2]));
- header.matrix_size[2] = ConvertFromString<double>(ExtractMatch(line,matches[3]));
+ if(!FindRegularExpressionNextLine(regex, regexstr[0], handle)) return false;
- if (!ReadLine(handle,line)) return false;
- if (regexec(&re_resol,line.c_str(),4,matches,0) == REG_NOMATCH) return false;
- header.resolution[0] = ConvertFromString<int>(ExtractMatch(line,matches[1]));
- header.resolution[1] = ConvertFromString<int>(ExtractMatch(line,matches[2]));
- header.resolution[2] = ConvertFromString<int>(ExtractMatch(line,matches[3]));
+ if(!FindRegularExpressionNextLine(regex, regexstr[1], handle)) return false;
+ header.matrix_size[0] = ConvertFromString<double>(regex.match(1));
+ header.matrix_size[1] = ConvertFromString<double>(regex.match(2));
+ header.matrix_size[2] = ConvertFromString<double>(regex.match(3));
- if (!ReadLine(handle,line)) return false;
- if (regexec(&re_voxel_size,line.c_str(),4,matches,0) == REG_NOMATCH) return false;
- header.voxel_size[0] = ConvertFromString<double>(ExtractMatch(line,matches[1]));
- header.voxel_size[1] = ConvertFromString<double>(ExtractMatch(line,matches[2]));
- header.voxel_size[2] = ConvertFromString<double>(ExtractMatch(line,matches[3]));
+ if(!FindRegularExpressionNextLine(regex, regexstr[2], handle)) return false;
+ header.resolution[0] = ConvertFromString<int>(regex.match(1));
+ header.resolution[1] = ConvertFromString<int>(regex.match(2));
+ header.resolution[2] = ConvertFromString<int>(regex.match(3));
- if (!ReadLine(handle,line)) return false;
- if (regexec(&re_nb_value,line.c_str(),2,matches,0) == REG_NOMATCH) return false;
- header.nb_value = ConvertFromString<int>(ExtractMatch(line,matches[1]));
+ if(!FindRegularExpressionNextLine(regex, regexstr[3], handle)) return false;
+ header.voxel_size[0] = ConvertFromString<double>(regex.match(1));
+ header.voxel_size[1] = ConvertFromString<double>(regex.match(2));
+ header.voxel_size[2] = ConvertFromString<double>(regex.match(3));
- if (!ReadLine(handle,line)) return false;
- if (regexec(&re_comment,line.c_str(),1,matches,0) == REG_NOMATCH) return false;
+ if(!FindRegularExpressionNextLine(regex, regexstr[4], handle)) return false;
+ header.nb_value = ConvertFromString<int>(regex.match(1));
- return true;
+ if(!FindRegularExpressionNextLine(regex, regexstr[0], handle)) return false;
+
+ return true;
}
//--------------------------------------------------------------------
// Read Image Content
void clitk::GateAsciiImageIO::Read(void* abstract_buffer)
{
- FILE* handle = fopen(m_FileName.c_str(),"r");
- if (!handle) {
- itkGenericExceptionMacro(<< "Could not open file (for reading): " << m_FileName);
- return;
- }
-
- GateAsciiHeader header;
- if (!ReadHeader(handle,header)) {
- itkGenericExceptionMacro(<< "Could not read header: " << m_FileName);
- fclose(handle);
- return;
- }
-
- {
- double* buffer = static_cast<double*>(abstract_buffer);
- int read_count = 0;
- while (true) {
- std::string line;
- if (!ReadLine(handle,line)) break;
- *buffer = ConvertFromString<double>(line);
- read_count++;
- buffer++;
- }
- assert(read_count == header.nb_value);
+ FILE* handle = fopen(m_FileName.c_str(),"r");
+ if (!handle) {
+ itkGenericExceptionMacro(<< "Could not open file (for reading): " << m_FileName);
+ return;
+ }
+
+ GateAsciiHeader header;
+ if (!ReadHeader(handle,header)) {
+ itkGenericExceptionMacro(<< "Could not read header: " << m_FileName);
+ fclose(handle);
+ return;
+ }
+
+ {
+ double* buffer = static_cast<double*>(abstract_buffer);
+ int read_count = 0;
+ while (true) {
+ std::string line;
+ if (!ReadLine(handle,line)) break;
+ *buffer = ConvertFromString<double>(line);
+ read_count++;
+ buffer++;
}
+ assert(read_count == header.nb_value);
+ }
- fclose(handle);
+ fclose(handle);
}
//--------------------------------------------------------------------
bool clitk::GateAsciiImageIO::CanWriteFile(const char* FileNameToWrite)
{
- if (GetExtension(std::string(FileNameToWrite)) != "txt") return false;
- return true;
+ if (GetExtension(std::string(FileNameToWrite)) != "txt") return false;
+ return true;
}
void clitk::GateAsciiImageIO::WriteImageInformation()
{
- cout << GetNumberOfDimensions() << endl;
+ cout << GetNumberOfDimensions() << endl;
}
bool clitk::GateAsciiImageIO::SupportsDimension(unsigned long dim)
{
- if (dim==2) return true;
- return false;
+ if (dim==2) return true;
+ return false;
}
//--------------------------------------------------------------------
// Write Image
void clitk::GateAsciiImageIO::Write(const void* abstract_buffer)
{
- const unsigned long nb_value = GetDimensions(0)*GetDimensions(1);
- std::stringstream stream;
- stream << "######################" << endl;
- stream << "# Matrix Size= (" << GetSpacing(0)*GetDimensions(0) << "," << GetSpacing(1)*GetDimensions(1) << ",1)" << endl;
- stream << "# Resol = (" << GetDimensions(0) << "," << GetDimensions(1) << ",1)" << endl;
- stream << "# VoxelSize = (" << GetSpacing(0) << "," << GetSpacing(1) << ",1)" << endl;
- stream << "# nbVal = " << nb_value << endl;
- stream << "######################" << endl;
-
- const double* buffer = static_cast<const double*>(abstract_buffer);
- for (unsigned long kk=0; kk<nb_value; kk++) { stream << buffer[kk] << endl; }
-
- FILE* handle = fopen(m_FileName.c_str(),"w");
- if (!handle) {
- itkGenericExceptionMacro(<< "Could not open file (for writing): " << m_FileName);
- return;
- }
-
- fwrite(stream.str().c_str(),1,stream.str().size(),handle);
-
- fclose(handle);
+ const unsigned long nb_value = GetDimensions(0)*GetDimensions(1);
+ std::stringstream stream;
+ stream << "######################" << endl;
+ stream << "# Matrix Size= (" << GetSpacing(0)*GetDimensions(0) << "," << GetSpacing(1)*GetDimensions(1) << ",1)" << endl;
+ stream << "# Resol = (" << GetDimensions(0) << "," << GetDimensions(1) << ",1)" << endl;
+ stream << "# VoxelSize = (" << GetSpacing(0) << "," << GetSpacing(1) << ",1)" << endl;
+ stream << "# nbVal = " << nb_value << endl;
+ stream << "######################" << endl;
+
+ const double* buffer = static_cast<const double*>(abstract_buffer);
+ for (unsigned long kk=0; kk<nb_value; kk++) {
+ stream << buffer[kk] << endl;
+ }
+
+ FILE* handle = fopen(m_FileName.c_str(),"w");
+ if (!handle) {
+ itkGenericExceptionMacro(<< "Could not open file (for writing): " << m_FileName);
+ return;
+ }
+
+ fwrite(stream.str().c_str(),1,stream.str().size(),handle);
+
+ fclose(handle);
}
/*=========================================================================
Program: vv http://www.creatis.insa-lyon.fr/rio/vv
- Authors belong to:
+ 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
#define CLITKGATEASCIIIMAGEIO_H
// itk include
-#include "itkImageIOBase.h"
+#include <itkImageIOBase.h>
+#include <itksys/RegularExpression.hxx>
#if defined (_MSC_VER) && (_MSC_VER < 1600)
//SR: taken from
#include <stdint.h>
#endif
-namespace clitk {
-
- //====================================================================
- // Class for reading gate ascii Image file format
- class GateAsciiImageIO: public itk::ImageIOBase
- {
- public:
- /** Standard class typedefs. */
- typedef GateAsciiImageIO Self;
- typedef itk::ImageIOBase Superclass;
- typedef itk::SmartPointer<Self> Pointer;
- typedef signed short int PixelType;
-
- struct GateAsciiHeader {
- double matrix_size[3];
- int resolution[3];
- double voxel_size[3];
- int nb_value;
- };
-
- GateAsciiImageIO():Superclass() {;}
-
- /** Method for creation through the object factory. */
- itkNewMacro(Self);
-
- /** Run-time type information (and related methods). */
- itkTypeMacro(GateAsciiImageIO, ImageIOBase);
-
- /*-------- This part of the interface deals with reading data. ------ */
- virtual void ReadImageInformation();
- virtual bool CanReadFile( const char* FileNameToRead );
- virtual void Read(void * buffer);
-
- /*-------- This part of the interfaces deals with writing data. ----- */
- virtual void WriteImageInformation();
- virtual bool CanWriteFile(const char* filename);
- virtual void Write(const void* buffer);
-
- virtual bool SupportsDimension(unsigned long dim);
-
- protected:
-
- static bool ReadHeader(FILE* handle, GateAsciiHeader& header);
- static bool ReadLine(FILE* handle, std::string& line);
-
- }; // end class GateAsciiImageIO
+namespace clitk
+{
+
+//====================================================================
+// Class for reading gate ascii Image file format
+class GateAsciiImageIO: public itk::ImageIOBase
+{
+public:
+ /** Standard class typedefs. */
+ typedef GateAsciiImageIO Self;
+ typedef itk::ImageIOBase Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef signed short int PixelType;
+
+ struct GateAsciiHeader {
+ double matrix_size[3];
+ int resolution[3];
+ double voxel_size[3];
+ int nb_value;
+ };
+
+ GateAsciiImageIO():Superclass() {}
+
+ /** Method for creation through the object factory. */
+ itkNewMacro(Self);
+
+ /** Run-time type information (and related methods). */
+ itkTypeMacro(GateAsciiImageIO, ImageIOBase);
+
+ /*-------- This part of the interface deals with reading data. ------ */
+ virtual void ReadImageInformation();
+ virtual bool CanReadFile( const char* FileNameToRead );
+ virtual void Read(void * buffer);
+
+ /*-------- This part of the interfaces deals with writing data. ----- */
+ virtual void WriteImageInformation();
+ virtual bool CanWriteFile(const char* filename);
+ virtual void Write(const void* buffer);
+
+ virtual bool SupportsDimension(unsigned long dim);
+
+protected:
+ static bool ReadHeader(FILE* handle, GateAsciiHeader& header);
+ static bool ReadLine(FILE* handle, std::string& line);
+ static bool FindRegularExpressionNextLine(itksys::RegularExpression ®, std::string &s, FILE* handle);
+
+}; // end class GateAsciiImageIO
} // end namespace
// explicit template instantiation
// Register factories
void clitk::RegisterClitkFactories()
{
-#ifdef unix
clitk::GateAsciiImageIOFactory::RegisterOneFactory();
-#endif
clitk::DicomRTDoseIOFactory::RegisterOneFactory();
#if ITK_VERSION_MAJOR <= 3
itk::ImageIOFactory::RegisterBuiltInFactories();
{
echo "$phase_file -> Resampling..."
clitkResampleImage -i $mask_dir_tmp/patient_$phase_nb.mhd -o $mask_dir_tmp/patient_$phase_nb.mhd --spacing $resample_spacing --interp $resample_algo
+ clitkResampleImage -i $mask_dir_tmp/patient_mask_$phase_nb.mhd -o $mask_dir_tmp/patient_mask_$phase_nb.mhd --spacing $resample_spacing --interp $resample_algo
clitkResampleImage -i $mask_dir_tmp/lungs_$phase_nb.mhd -o $mask_dir_tmp/lungs_$phase_nb.mhd --like $mask_dir_tmp/patient_$phase_nb.mhd
}
section "Step 2 : find trachea"
-option "skipslices" - "Number of slices to skip before searching seed" int no default="0"
-option "upperThresholdForTrachea" - "Initial upper threshold for trachea" double no default="-900"
-option "multiplierForTrachea" - "Multiplier for the region growing" double no default="5"
-option "thresholdStepSizeForTrachea" - "Threshold step size" int no default="64"
-option "seed" - "Index of the trachea seed point (in pixel, not in mm)" int no multiple
-option "doNotCheckTracheaVolume" - "If set, do not check the trachea volume" flag off
-option "verboseRG" - "Verbose RegionGrowing" flag off
+option "type" - "trachea finding algorithm (0: original; 1: LOLA11)" int no default="0"
+option "skipslices" - "0: Number of slices to skip before searching seed" int no default="0"
+option "upperThresholdForTrachea" - "0: Initial upper threshold for trachea" double no default="-900"
+option "multiplierForTrachea" - "0: Multiplier for the region growing" double no default="5"
+option "thresholdStepSizeForTrachea" - "0: Threshold step size" int no default="64"
+option "seed" - "0,1: Index of the trachea seed point (in pixel, not in mm)" int no multiple
+option "doNotCheckTracheaVolume" - "0,1: If set, do not check the trachea volume" flag off
+option "verboseRG" - "0,1: Verbose RegionGrowing" flag off
+option "maxElongation" - "1: Maximum allowed elongation of candidate regions for the trachea" float no default="0.5"
+option "numSlices" - "1: Number of slices to search for trachea" int no default="50"
+option "seedPreProcessingThreshold" - "1: Threshold for the pre-processing step of the algorithm" int no default="-400"
section "Step 3 : auto extract lung"
void SetLabelizeParameters1(LabelParamType * a) { m_LabelizeParameters1 = a; }
itkGetConstMacro(LabelizeParameters1, LabelParamType*);
+
+ itkSetMacro(TracheaSeedAlgorithm, int);
+ itkGetConstMacro(TracheaSeedAlgorithm, int);
// Step 2 options FindTrachea
itkSetMacro(UpperThresholdForTrachea, InputImagePixelType);
itkSetMacro(ThresholdStepSizeForTrachea, InputImagePixelType);
itkGetConstMacro(ThresholdStepSizeForTrachea, InputImagePixelType);
+ // options FindTrachea2
+ itkSetMacro(NumSlices, int);
+ itkGetConstMacro(NumSlices, int);
+ itkSetMacro(MaxElongation, double);
+ itkGetConstMacro(MaxElongation, double);
+ itkSetMacro(SeedPreProcessingThreshold, int);
+ itkGetConstMacro(SeedPreProcessingThreshold, int);
+
void AddSeed(InternalIndexType s);
std::vector<InternalIndexType> & GetSeeds() { return m_Seeds; }
itkSetMacro(AutoCrop, bool);
itkGetConstMacro(AutoCrop, bool);
itkBooleanMacro(AutoCrop);
-
+
protected:
ExtractLungFilter();
virtual ~ExtractLungFilter() {}
LabelParamType* m_LabelizeParameters1;
// Step 2
+ int m_TracheaSeedAlgorithm;
InputImagePixelType m_UpperThresholdForTrachea;
InputImagePixelType m_ThresholdStepSizeForTrachea;
double m_MultiplierForTrachea;
int m_NumberOfSlicesToSkipBeforeSearchingSeed;
bool m_TracheaVolumeMustBeCheckedFlag;
bool m_VerboseRegionGrowingFlag;
+ int m_NumSlices;
+ double m_MaxElongation;
+ int m_SeedPreProcessingThreshold;
// Step 3
int m_NumberOfHistogramBins;
// Functions for trachea extraction
bool SearchForTracheaSeed(int skip);
+ bool SearchForTracheaSeed2(int numberOfSlices);
void SearchForTrachea();
void TracheaRegionGrowing();
double ComputeTracheaVolume();
#include "itkBinaryMorphologicalOpeningImageFilter.h"
#include "itkBinaryMorphologicalClosingImageFilter.h"
#include "itkConstantPadImageFilter.h"
+#include <itkBinaryBallStructuringElement.h>
+#include "itkStatisticsLabelObject.h"
+#include "itkLabelMap.h"
+#include "itkLabelImageToShapeLabelMapFilter.h"
+#include "itkLabelMapToLabelImageFilter.h"
+#include "itkExtractImageFilter.h"
+#include "itkOrientImageFilter.h"
+#include "itkSpatialOrientationAdapter.h"
+#include "itkImageDuplicator.h"
+#include <fcntl.h>
//--------------------------------------------------------------------
template <class ImageType>
SetLabelizeParameters1(p1);
// Step 2 default values
- SetUpperThresholdForTrachea(-900);
+ SetTracheaSeedAlgorithm(0);
+ SetUpperThresholdForTrachea(-700);
SetMultiplierForTrachea(5);
SetThresholdStepSizeForTrachea(64);
SetNumberOfSlicesToSkipBeforeSearchingSeed(0);
TracheaVolumeMustBeCheckedFlagOn();
+ SetNumSlices(50);
+ SetMaxElongation(0.5);
+ SetSeedPreProcessingThreshold(-400);
+
// Step 3 default values
SetNumberOfHistogramBins(500);
}
//--------------------------------------------------------------------
+
+bool is_orientation_superior(itk::SpatialOrientation::ValidCoordinateOrientationFlags orientation)
+{
+ itk::SpatialOrientation::CoordinateTerms sup = itk::SpatialOrientation::ITK_COORDINATE_Superior;
+ std::cout << "orientation: " << std::hex << orientation << "; superior: " << std::hex << sup << std::endl;
+ std::cout << std::dec;
+
+ bool primary = (orientation & 0x0000ff) == sup;
+ bool secondary = ((orientation & 0x00ff00) >> 8) == sup;
+ bool tertiary = ((orientation & 0xff0000) >> 16) == sup;
+ return primary || secondary || tertiary;
+}
+
+//--------------------------------------------------------------------
+template <class ImageType>
+bool
+clitk::ExtractLungFilter<ImageType>::
+SearchForTracheaSeed2(int numberOfSlices)
+{
+ if (m_Seeds.size() == 0) { // try to find seed (if not zero, it is given by user)
+ if (GetVerboseFlag())
+ std::cout << "SearchForTracheaSeed2(" << numberOfSlices << ", " << GetMaxElongation() << ")" << std::endl;
+
+ typedef unsigned char MaskPixelType;
+ typedef itk::Image<MaskPixelType, ImageType::ImageDimension> MaskImageType;
+ typedef itk::Image<typename MaskImageType::PixelType, 2> MaskImageType2D;
+ typedef itk::BinaryThresholdImageFilter<ImageType, MaskImageType> ThresholdFilterType;
+ typedef itk::BinaryBallStructuringElement<MaskPixelType, MaskImageType::ImageDimension> KernelType;
+ typedef itk::BinaryMorphologicalClosingImageFilter<MaskImageType, MaskImageType, KernelType> ClosingFilterType;
+ typedef itk::BinaryMorphologicalOpeningImageFilter<MaskImageType, MaskImageType, KernelType> OpeningFilterType;
+ typedef itk::ExtractImageFilter<MaskImageType, MaskImageType2D> SlicerFilterType;
+ typedef itk::ConnectedComponentImageFilter<MaskImageType2D, MaskImageType2D> LabelFilterType;
+ typedef itk::ShapeLabelObject<MaskPixelType, MaskImageType2D::ImageDimension> ShapeLabelType;
+ typedef itk::LabelMap<ShapeLabelType> LabelImageType;
+ typedef itk::LabelImageToShapeLabelMapFilter<MaskImageType2D, LabelImageType> ImageToLabelMapFilterType;
+ typedef itk::LabelMapToLabelImageFilter<LabelImageType, MaskImageType2D> LabelMapToImageFilterType;
+ typedef itk::ImageFileWriter<MaskImageType2D> FileWriterType;
+
+ // threshold to isolate airawys and lungs
+ typename ThresholdFilterType::Pointer threshold = ThresholdFilterType::New();
+ threshold->SetLowerThreshold(-2000);
+ threshold->SetUpperThreshold(GetSeedPreProcessingThreshold());
+ threshold->SetInput(working_input);
+ threshold->Update();
+
+ KernelType kernel_closing, kernel_opening;
+
+ // remove small noise
+ typename OpeningFilterType::Pointer opening = OpeningFilterType::New();
+ kernel_opening.SetRadius(1);
+ opening->SetKernel(kernel_opening);
+ opening->SetInput(threshold->GetOutput());
+ opening->Update();
+
+ typename SlicerFilterType::Pointer slicer = SlicerFilterType::New();
+ slicer->SetInput(opening->GetOutput());
+
+ // label result
+ typename LabelFilterType::Pointer label_filter = LabelFilterType::New();
+ label_filter->SetInput(slicer->GetOutput());
+
+ // extract shape information from labels
+ typename ImageToLabelMapFilterType::Pointer label_to_map_filter = ImageToLabelMapFilterType::New();
+ label_to_map_filter->SetInput(label_filter->GetOutput());
+
+ typename LabelMapToImageFilterType::Pointer map_to_label_filter = LabelMapToImageFilterType::New();
+ typename FileWriterType::Pointer writer = FileWriterType::New();
+
+ typename ImageType::IndexType index;
+ typename ImageType::RegionType region = working_input->GetLargestPossibleRegion();
+ typename ImageType::SizeType size = region.GetSize();
+ typename ImageType::SpacingType spacing = working_input->GetSpacing();
+ typename ImageType::PointType origin = working_input->GetOrigin();
+
+ int nslices = min(numberOfSlices, size[2]);
+ int start = 0, increment = 1;
+ itk::SpatialOrientationAdapter orientation;
+ typename ImageType::DirectionType dir = working_input->GetDirection();
+ if (!is_orientation_superior(orientation.FromDirectionCosines(dir))) {
+ start = size[2]-1;
+ increment = -1;
+ }
+
+ typename MaskImageType::PointType image_centre;
+ image_centre[0] = size[0]/2;
+ image_centre[1] = size[1]/2;
+ image_centre[2] = 0;
+
+ typedef InternalIndexType SeedType;
+ SeedType trachea_centre, shape_centre, max_e_centre, prev_e_centre;
+ typedef std::list<SeedType> PointListType;
+ typedef std::list<PointListType> SequenceListType;
+ PointListType* current_sequence = NULL;
+ SequenceListType sequence_list;
+
+ prev_e_centre.Fill(0);
+ std::ostringstream file_name;
+ index[0] = index[1] = 0;
+ size[0] = size[1] = 512;
+ size[2] = 0;
+ while (nslices--) {
+ index[2] = start;
+ start += increment;
+
+ region.SetIndex(index);
+ region.SetSize(size);
+ slicer->SetExtractionRegion(region);
+ slicer->Update();
+ label_filter->SetInput(slicer->GetOutput());
+ label_filter->Update();
+
+ label_to_map_filter->SetInput(label_filter->GetOutput());
+ label_to_map_filter->Update();
+ typename LabelImageType::Pointer label_map = label_to_map_filter->GetOutput();
+
+ if (GetWriteStepFlag()) {
+ map_to_label_filter->SetInput(label_map);
+ writer->SetInput(map_to_label_filter->GetOutput());
+ file_name.str("");
+ file_name << "labels_";
+ file_name.width(3);
+ file_name.fill('0');
+ file_name << index[2] << ".mhd";
+ writer->SetFileName(file_name.str().c_str());
+ writer->Update();
+ }
+
+ typename LabelImageType::LabelObjectContainerType shapes_map = label_map->GetLabelObjectContainer();
+ typename LabelImageType::LabelObjectContainerType::const_iterator s;
+ typename ShapeLabelType::Pointer shape, max_e_shape;
+ double max_elongation = GetMaxElongation();
+ double max_size = size[0]*size[1]/128;
+ double max_e = 0;
+ int nshapes = 0;
+ for (s = shapes_map.begin(); s != shapes_map.end(); s++) {
+ shape = s->second;
+ if (shape->GetSize() > 150 && shape->GetSize() <= max_size) {
+ double e = 1 - 1/shape->GetBinaryElongation();
+ //double area = 1 - r->Area() ;
+ if (e < max_elongation) {
+ nshapes++;
+ shape_centre[0] = (shape->GetCentroid()[0] - origin[0])/spacing[0];
+ shape_centre[1] = (shape->GetCentroid()[1] - origin[1])/spacing[1];
+ shape_centre[2] = index[2];
+ //double d = 1 - (shape_centre - image_centre).Magnitude()/max_dist;
+ double dx = shape_centre[0] - image_centre[0];
+ double d = 1 - dx*2/size[0];
+ e = e + d;
+ if (e > max_e)
+ {
+ max_e = e;
+ max_e_shape = shape;
+ max_e_centre = shape_centre;
+ }
+ }
+ }
+ }
+
+ if (nshapes > 0)
+ {
+ itk::Point<typename SeedType::IndexValueType, ImageType::ImageDimension> p1, p2;
+ p1[0] = max_e_centre[0];
+ p1[1] = max_e_centre[1];
+ p1[2] = max_e_centre[2];
+
+ p2[0] = prev_e_centre[0];
+ p2[1] = prev_e_centre[1];
+ p2[2] = prev_e_centre[2];
+
+ double mag = (p2 - p1).GetNorm();
+ if (GetVerboseFlag()) {
+ cout.precision(3);
+ cout << index[2] << ": ";
+ cout << "region(" << max_e_centre[0] << ", " << max_e_centre[1] << ", " << max_e_centre[2] << "); ";
+ cout << "prev_region(" << prev_e_centre[0] << ", " << prev_e_centre[1] << ", " << prev_e_centre[2] << "); ";
+ cout << "mag(" << mag << "); " << endl;
+ }
+
+ if (mag > 5)
+ {
+ PointListType point_list;
+ point_list.push_back(max_e_centre);
+ sequence_list.push_back(point_list);
+ current_sequence = &sequence_list.back();
+ }
+ else if (current_sequence)
+ current_sequence->push_back(max_e_centre);
+
+ prev_e_centre= max_e_centre;
+ }
+ }
+
+ size_t longest = 0;
+ for (typename SequenceListType::iterator s = sequence_list.begin(); s != sequence_list.end(); s++)
+ {
+ if (s->size() > longest)
+ {
+ longest = s->size();
+ trachea_centre = s->front();
+ }
+ }
+
+ if (GetVerboseFlag())
+ std::cout << "seed at: " << trachea_centre << std::endl;
+ m_Seeds.push_back(trachea_centre);
+ }
+ return (m_Seeds.size() != 0);
+}
+//--------------------------------------------------------------------
+
//--------------------------------------------------------------------
template <class ImageType>
void
f->SetUpper(GetUpperThresholdForTrachea());
f->SetMinimumLowerThreshold(-2000);
// f->SetMaximumUpperThreshold(0); // MAYBE TO CHANGE ???
- f->SetMaximumUpperThreshold(-800); // MAYBE TO CHANGE ???
+ f->SetMaximumUpperThreshold(-700); // MAYBE TO CHANGE ???
f->SetAdaptLowerBorder(false);
f->SetAdaptUpperBorder(true);
f->SetMinimumSize(5000);
// compute trachea volume
// if volume not plausible -> skip more slices and restart
+ bool has_seed;
bool stop = false;
double volume = 0.0;
int skip = GetNumberOfSlicesToSkipBeforeSearchingSeed();
while (!stop) {
stop = true;
- if (SearchForTracheaSeed(skip)) {
+ if (GetTracheaSeedAlgorithm() == 0)
+ has_seed = SearchForTracheaSeed(skip);
+ else
+ has_seed = SearchForTracheaSeed2(GetNumSlices());
+
+ if (has_seed) {
TracheaRegionGrowing();
volume = ComputeTracheaVolume()/1000; // assume mm3, so divide by 1000 to get cc
if (GetWriteStepFlag()) {
if (GetTracheaVolumeMustBeCheckedFlag()) {
if ((volume > 10) && (volume < 65 )) { // depend on image size ...
// Typical volume 22.59 cm 3 (± 7.69 cm 3 ) [ Leader 2004 ]
- if (GetVerboseStepFlag()) {
+ if (GetVerboseStepFlag())
+ {
std::cout << "\t Found trachea with volume " << volume << " cc." << std::endl;
}
}
f->SetUpperThresholdForTrachea(mArgsInfo.upperThresholdForTrachea_arg);
f->SetMultiplierForTrachea(mArgsInfo.multiplierForTrachea_arg);
f->SetThresholdStepSizeForTrachea(mArgsInfo.thresholdStepSizeForTrachea_arg);
+ f->SetTracheaSeedAlgorithm(mArgsInfo.type_arg);
+ f->SetNumSlices(mArgsInfo.numSlices_arg);
+ f->SetMaxElongation(mArgsInfo.maxElongation_arg);
+ f->SetSeedPreProcessingThreshold(mArgsInfo.seedPreProcessingThreshold_arg);
typename FilterType::InputImageIndexType s;
if (mArgsInfo.seed_given) {
version "1.2"
purpose "Extract data from a Dicom wave file to a text file"
option "config" - "config file" string no
-option "InputFile" i "Dicom inputfile name" string no
-option "OutputFile" o "Text outputfile name" string no
+option "InputFile" i "Dicom inputfile name" string yes
+option "OutputFile" o "Text outputfile name" string yes
option "Verbose" v "verbose" int no
</property>
</widget>
</item>
- <item row="5" column="0" colspan="3">
+ <item row="5" column="0" colspan="4">
<widget class="QLabel" name="valueFusionnedLabel">
<property name="sizePolicy">
<sizepolicy hsizetype="Preferred" vsizetype="Fixed">
mPreset = 0;
mOverlayColor = 130;
- mFusionOpacity = 70;
- mFusionThresOpacity = 0;
+ mFusionOpacity = 30;
+ mFusionThresOpacity = 1;
mFusionColorMap = 3;
mFusionWindow = 1000;
mFusionLevel = 1000;
git://git.creatis.insa-lyon.fr / clitk.git / commitdiff