#include "itkOrientImageFilter.h"
#include "itkSpatialOrientationAdapter.h"
#include "itkImageDuplicator.h"
+#include "itkRelabelComponentImageFilter.h"
+
#include <fcntl.h>
//--------------------------------------------------------------------
SetForegroundValue(1);
SetMinimalComponentSize(100);
VerboseRegionGrowingFlagOff();
+ RemoveSmallLabelBeforeSeparationFlagOn();
// Step 1 default values
SetUpperThreshold(-300);
TracheaVolumeMustBeCheckedFlagOn();
SetNumSlices(50);
SetMaxElongation(0.5);
- SetSeedPreProcessingThreshold(-400);
-
+ SetSeedPreProcessingThreshold(-400);
// Step 3 default values
SetNumberOfHistogramBins(500);
template <class ImageType>
void
clitk::ExtractLungFilter<ImageType>::
-AddSeed(InternalIndexType s)
+//AddSeed(InternalIndexType s)
+AddSeed(InputImagePointType s)
+{
+ m_SeedsInMM.push_back(s);
+}
+//--------------------------------------------------------------------
+
+
+//--------------------------------------------------------------------
+template <class ImageType>
+void
+clitk::ExtractLungFilter<ImageType>::
+AddSeedInPixels(InternalIndexType s)
{
m_Seeds.push_back(s);
}
StartNewStep("Search for the trachea");
SearchForTrachea();
PrintMemory(GetVerboseMemoryFlag(), "After SearchForTrachea");
+ if (m_Seeds.empty()) {
+ clitkExceptionMacro("No seeds for trachea... Aborting.");
+ }
//--------------------------------------------------------------------
//--------------------------------------------------------------------
//--------------------------------------------------------------------
//--------------------------------------------------------------------
StartNewStep("Separate Left/Right lungs");
- PrintMemory(GetVerboseMemoryFlag(), "Before Separate");
+ PrintMemory(GetVerboseMemoryFlag(), "Before Separate");
// Initial label
working_mask = Labelize<MaskImageType>(working_mask,
GetBackgroundValue(),
statisticsImageFilter->SetInput(working_mask);
statisticsImageFilter->Update();
unsigned int initialNumberOfLabels = statisticsImageFilter->GetMaximum();
- working_mask = statisticsImageFilter->GetOutput();
-
+ working_mask = statisticsImageFilter->GetOutput();
PrintMemory(GetVerboseMemoryFlag(), "After count label");
+
+ // If already 2 labels, but a too big differences, remove the
+ // smalest one (sometimes appends with the stomach
+ if (initialNumberOfLabels >1) {
+ if (GetRemoveSmallLabelBeforeSeparationFlag()) {
+ typedef itk::RelabelComponentImageFilter<MaskImageType, MaskImageType> RelabelFilterType;
+ typename RelabelFilterType::Pointer relabelFilter = RelabelFilterType::New();
+ relabelFilter->SetInput(working_mask);
+ relabelFilter->SetMinimumObjectSize(10);
+ relabelFilter->Update();
+ const std::vector<float> & a = relabelFilter->GetSizeOfObjectsInPhysicalUnits();
+ std::vector<MaskImagePixelType> remove_label;
+ for(unsigned int i=1; i<a.size(); i++) {
+ if (a[i] < 0.5*a[0]) { // more than 0.5 difference
+ remove_label.push_back(i+1); // label zero is BG
+ }
+ }
+ working_mask =
+ clitk::RemoveLabels<MaskImageType>(working_mask, GetBackgroundValue(), remove_label);
+ statisticsImageFilter->SetInput(working_mask);
+ statisticsImageFilter->Update();
+ initialNumberOfLabels = statisticsImageFilter->GetMaximum();
+ }
+ }
// Decompose the first label
if (initialNumberOfLabels<2) {
it.GoToBegin();
while (!it.IsAtEnd()) {
if(it.Get() < GetUpperThresholdForTrachea() ) {
- AddSeed(it.GetIndex());
+ AddSeedInPixels(it.GetIndex());
// DD(it.GetIndex());
}
++it;
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;
SearchForTracheaSeed2(int numberOfSlices)
{
if (m_Seeds.size() == 0) { // try to find seed (if not zero, it is given by user)
- if (GetVerboseFlag())
+ if (GetVerboseRegionGrowingFlag())
std::cout << "SearchForTracheaSeed2(" << numberOfSlices << ", " << GetMaxElongation() << ")" << std::endl;
typedef unsigned char MaskPixelType;
opening->Update();
typename SlicerFilterType::Pointer slicer = SlicerFilterType::New();
+ slicer->SetDirectionCollapseToIdentity();
slicer->SetInput(opening->GetOutput());
// label result
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();
+ unsigned int nlables = label_map->GetNumberOfLabelObjects();
+ for (unsigned int j = 0; j < nlables; j++) {
+ shape = label_map->GetNthLabelObject(j);
+ if (shape->Size() > 150 && shape->Size() <= max_size) {
+ double e = 1 - 1/shape->GetElongation();
//double area = 1 - r->Area() ;
if (e < max_elongation) {
nshapes++;
p2[2] = prev_e_centre[2];
double mag = (p2 - p1).GetNorm();
- if (GetVerboseFlag()) {
+ if (GetVerboseRegionGrowingFlag()) {
cout.precision(3);
cout << index[2] << ": ";
cout << "region(" << max_e_centre[0] << ", " << max_e_centre[1] << ", " << max_e_centre[2] << "); ";
prev_e_centre= max_e_centre;
}
+ else {
+ if (GetVerboseRegionGrowingFlag()) {
+ cout << "No shapes found at slice " << index[2] << std::endl;
+ }
+ }
}
size_t longest = 0;
}
}
- if (GetVerboseFlag())
- std::cout << "seed at: " << trachea_centre << std::endl;
- m_Seeds.push_back(trachea_centre);
+ if (longest > 0) {
+ if (GetVerboseRegionGrowingFlag())
+ std::cout << "seed at: " << trachea_centre << std::endl;
+ m_Seeds.push_back(trachea_centre);
+ }
}
return (m_Seeds.size() != 0);
f->SetUpper(GetUpperThresholdForTrachea());
f->SetMinimumLowerThreshold(-2000);
// f->SetMaximumUpperThreshold(0); // MAYBE TO CHANGE ???
- f->SetMaximumUpperThreshold(-700); // MAYBE TO CHANGE ???
+ f->SetMaximumUpperThreshold(-300); // MAYBE TO CHANGE ???
f->SetAdaptLowerBorder(false);
f->SetAdaptUpperBorder(true);
f->SetMinimumSize(5000);
f->SetVerbose(GetVerboseRegionGrowingFlag());
for(unsigned int i=0; i<m_Seeds.size();i++) {
f->AddSeed(m_Seeds[i]);
- // DD(m_Seeds[i]);
+ //DD(m_Seeds[i]);
}
f->Update();
PrintMemory(GetVerboseMemoryFlag(), "After RG update");
// compute trachea volume
// if volume not plausible -> skip more slices and restart
+ // If initial seed, convert from mm to pixels
+ if (m_SeedsInMM.size() > 0) {
+ for(unsigned int i=0; i<m_SeedsInMM.size(); i++) {
+ InputImageIndexType index;
+ working_input->TransformPhysicalPointToIndex(m_SeedsInMM[i], index);
+ m_Seeds.push_back(index);
+ }
+ }
+
bool has_seed;
bool stop = false;
double volume = 0.0;
std::cout << "\t Found trachea with volume " << volume << " cc." << std::endl;
}
}
- else {
- if (GetVerboseStepFlag()) {
- std::cout << "\t The volume of the trachea (" << volume
- << " cc) seems not correct. I skip some slices (" << skip << ") and restart to find seeds."
- << std::endl;
- }
- skip += 5;
- stop = false;
- // empty the list of seed
- m_Seeds.clear();
+ else
+ if (GetTracheaSeedAlgorithm() == 0) {
+ if (GetVerboseStepFlag()) {
+ std::cout << "\t The volume of the trachea (" << volume
+ << " cc) seems not correct. I skip some slices (" << skip << ") and restart to find seeds."
+ << std::endl;
+ }
+ skip += 5;
+ stop = false;
+ // empty the list of seed
+ m_Seeds.clear();
}
if (skip > 0.5 * working_input->GetLargestPossibleRegion().GetSize()[2]) {
// we want to skip more than a half of the image, it is probably a bug
git://git.creatis.insa-lyon.fr / clitk.git / blobdiff