/*=========================================================================
  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 CLITKEXTRACTLUNGSFILTER_TXX
#define CLITKEXTRACTLUNGSFILTER_TXX

// clitk
#include "clitkImageCommon.h"
#include "clitkSetBackgroundImageFilter.h"
#include "clitkSegmentationFunctions.h"
#include "clitkAutoCropFilter.h"

// itk
#include "itkBinaryThresholdImageFilter.h"
#include "itkConnectedComponentImageFilter.h"
#include "itkRelabelComponentImageFilter.h"
#include "itkOtsuThresholdImageFilter.h"

//--------------------------------------------------------------------
template <class TInputImageType, class TMaskImageType>
clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
ExtractLungFilter():
  clitk::FilterBase(),
  itk::ImageToImageFilter<TInputImageType, TMaskImageType>()
{
  // Default global options
  this->SetNumberOfRequiredInputs(2);
  SetPatientMaskBackgroundValue(0);
  SetBackgroundValue(0); // Must be zero
  SetForegroundValue(1);

  // Step 1 default values
  SetUpperThreshold(-300);
  SetLowerThreshold(-1000);
  UseLowerThresholdOff();
  LabelParamType * p1 = new LabelParamType;
  p1->SetFirstKeep(1);
  p1->UseLastKeepOff();
  p1->AddLabelToRemove(2);
  SetLabelizeParameters1(p1);

  // Step 2 default values
  SetUpperThresholdForTrachea(-900);
  SetMultiplierForTrachea(5);
  SetThresholdStepSizeForTrachea(64);
  
  // Step 3 default values
  SetNumberOfHistogramBins(500);
  LabelParamType * p2 = new LabelParamType;
  p2->SetFirstKeep(1);
  p2->UseLastKeepOff();
  // p->AddLabelToRemove(2); // No label to remove by default
  SetLabelizeParameters2(p2);

  // Step 4 default values
  SetRadiusForTrachea(1);
  LabelParamType * p3 = new LabelParamType;
  p3->SetFirstKeep(1);
  p3->SetLastKeep(2);
  p3->UseLastKeepOff();
  SetLabelizeParameters3(p3);
}
//--------------------------------------------------------------------


//--------------------------------------------------------------------
template <class TInputImageType, class TMaskImageType>
void 
clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
SetInput(const TInputImageType * image) 
{
  this->SetNthInput(0, const_cast<TInputImageType *>(image));
}
//--------------------------------------------------------------------


//--------------------------------------------------------------------
template <class TInputImageType, class TMaskImageType>
void 
clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
SetInputPatientMask(TMaskImageType * image, MaskImagePixelType bg ) 
{
  this->SetNthInput(1, const_cast<TMaskImageType *>(image));
  SetPatientMaskBackgroundValue(bg);
}
//--------------------------------------------------------------------


//--------------------------------------------------------------------
template <class TInputImageType, class TMaskImageType>
void 
clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
AddSeed(InternalIndexType s) 
{ 
  m_Seeds.push_back(s);
}
//--------------------------------------------------------------------


//--------------------------------------------------------------------
template <class TInputImageType, class TMaskImageType>
template<class ArgsInfoType>
void 
clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
SetArgsInfo(ArgsInfoType mArgsInfo)
{
  SetVerboseOption_GGO(mArgsInfo);
  SetVerboseStep_GGO(mArgsInfo);
  SetWriteStep_GGO(mArgsInfo);
  SetVerboseWarningOff_GGO(mArgsInfo);

  SetUpperThreshold_GGO(mArgsInfo);
  SetLowerThreshold_GGO(mArgsInfo);
  SetLabelizeParameters1_GGO(mArgsInfo);
  if (!mArgsInfo.remove1_given) {
    GetLabelizeParameters1()->AddLabelToRemove(2); 
    if (GetVerboseOption()) VerboseOption("remove1", 2);
  }

  SetUpperThresholdForTrachea_GGO(mArgsInfo);
  SetMultiplierForTrachea_GGO(mArgsInfo);
  SetThresholdStepSizeForTrachea_GGO(mArgsInfo);
  AddSeed_GGO(mArgsInfo);

  SetNumberOfHistogramBins_GGO(mArgsInfo);
  SetLabelizeParameters2_GGO(mArgsInfo);

  SetRadiusForTrachea_GGO(mArgsInfo);
  SetLabelizeParameters3_GGO(mArgsInfo);
}
//--------------------------------------------------------------------


//--------------------------------------------------------------------
template <class TInputImageType, class TMaskImageType>
void 
clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
GenerateOutputInformation() 
{ 
  input   = dynamic_cast<const TInputImageType*>(itk::ProcessObject::GetInput(0));
  Superclass::GenerateOutputInformation();
//   MaskImagePointer output = this->GetOutput(0);

  // Get input pointers
  input   = dynamic_cast<const TInputImageType*>(itk::ProcessObject::GetInput(0));
  patient = dynamic_cast<const TMaskImageType*>(itk::ProcessObject::GetInput(1));

  // Check image
  if (!HasSameSizeAndSpacing<TInputImageType, TMaskImageType>(input, patient)) {
    this->SetLastError("* ERROR * the images (input and patient mask) must have the same size & spacing");
    return;
  }
  
  //--------------------------------------------------------------------
  //--------------------------------------------------------------------
  StartNewStep("Set background to initial image");
  working_input = SetBackground<TInputImageType, TMaskImageType>
    (input, patient, GetPatientMaskBackgroundValue(), -1000);
  StopCurrentStep<InputImageType>(working_input);

  //--------------------------------------------------------------------
  //--------------------------------------------------------------------
  StartNewStep("Remove Air");
  // Threshold to get air
  typedef itk::BinaryThresholdImageFilter<InputImageType, InternalImageType> InputBinarizeFilterType; 
  typename InputBinarizeFilterType::Pointer binarizeFilter=InputBinarizeFilterType::New();
  binarizeFilter->SetInput(working_input);
  if (m_UseLowerThreshold) binarizeFilter->SetLowerThreshold(m_LowerThreshold);
  binarizeFilter->SetUpperThreshold(m_UpperThreshold);
  binarizeFilter ->SetInsideValue(this->GetForegroundValue());
  binarizeFilter ->SetOutsideValue(this->GetBackgroundValue());
  binarizeFilter->Update();
  working_image = binarizeFilter->GetOutput();
  
  // Labelize and keep right labels
  working_image = Labelize<InternalImageType>(working_image, GetBackgroundValue(), true, GetMinimalComponentSize());
  working_image = RemoveLabels<InternalImageType>
    (working_image, GetBackgroundValue(), GetLabelizeParameters1()->GetLabelsToRemove());
  typename InternalImageType::Pointer air = KeepLabels<InternalImageType>
    (working_image, 
     GetBackgroundValue(), 
     GetForegroundValue(), 
     GetLabelizeParameters1()->GetFirstKeep(), 
     GetLabelizeParameters1()->GetLastKeep(), 
     GetLabelizeParameters1()->GetUseLastKeep());
 
  // Set Air to BG
  working_input = SetBackground<TInputImageType, InternalImageType>
    (working_input, air, this->GetForegroundValue(), this->GetBackgroundValue());

  // End
  StopCurrentStep<InputImageType>(working_input);

  //--------------------------------------------------------------------
  //--------------------------------------------------------------------
  StartNewStep("Find the trachea");
  if (m_Seeds.size() == 0) { // try to find seed
    // Search seed (parameters = UpperThresholdForTrachea)
    static const unsigned int Dim = InputImageType::ImageDimension;
    typename InternalImageType::RegionType sliceRegion = working_input->GetLargestPossibleRegion();
    typename InternalImageType::SizeType sliceRegionSize = sliceRegion.GetSize();
    typename InternalImageType::IndexType sliceRegionIndex = sliceRegion.GetIndex();
    sliceRegionIndex[Dim-1]=sliceRegionSize[Dim-1]-5;
    sliceRegionSize[Dim-1]=5;
    sliceRegion.SetSize(sliceRegionSize);
    sliceRegion.SetIndex(sliceRegionIndex);
    
    typedef  itk::ImageRegionConstIterator<InputImageType> IteratorType;
    IteratorType it(working_input, sliceRegion);
    it.GoToBegin();
    while (!it.IsAtEnd()) {
      if(it.Get() < GetUpperThresholdForTrachea() ) {
        AddSeed(it.GetIndex());
      }
      ++it;
    }
  }
  
  if (m_Seeds.size() != 0) {
    // Explosion controlled region growing
    typedef clitk::ExplosionControlledThresholdConnectedImageFilter<InputImageType, InternalImageType> ImageFilterType;
    typename ImageFilterType::Pointer f= ImageFilterType::New();
    f->SetInput(working_input);
    f->SetVerbose(false);
    f->SetLower(-2000);
    f->SetUpper(GetUpperThresholdForTrachea());
    f->SetMinimumLowerThreshold(-2000);
    f->SetMaximumUpperThreshold(0);
    f->SetAdaptLowerBorder(false);
    f->SetAdaptUpperBorder(true);
    f->SetMinimumSize(5000); 
    f->SetReplaceValue(1);
    f->SetMultiplier(GetMultiplierForTrachea());
    f->SetThresholdStepSize(GetThresholdStepSizeForTrachea());
    f->SetMinimumThresholdStepSize(1);
    for(unsigned int i=0; i<m_Seeds.size();i++) {
      //    std::cout<<"Adding seed " <<m_Seeds[i]<<"..."<<std::endl;
      f->AddSeed(m_Seeds[i]);
    }  
    f->Update();
    trachea_tmp = f->GetOutput();
    // Set output
    StopCurrentStep<InternalImageType>(trachea_tmp);

  }
  else { // Trachea not found
    this->SetWarning("* WARNING * No seed found for trachea.");
    // Set output
    StopCurrentStep();
  }

  //--------------------------------------------------------------------
  //--------------------------------------------------------------------
  StartNewStep("Extract the lung with Otsu filter");
  // Automated Otsu thresholding and relabeling
  typedef itk::OtsuThresholdImageFilter<InputImageType,InternalImageType> OtsuThresholdImageFilterType;
  typename OtsuThresholdImageFilterType::Pointer otsuFilter=OtsuThresholdImageFilterType::New();
  otsuFilter->SetInput(working_input);
  otsuFilter->SetNumberOfHistogramBins(GetNumberOfHistogramBins());
  otsuFilter->SetInsideValue(this->GetForegroundValue());
  otsuFilter->SetOutsideValue(this->GetBackgroundValue());
  otsuFilter->Update();
  working_image =  otsuFilter->GetOutput();

  // Set output
  StopCurrentStep<InternalImageType>(working_image);

  //--------------------------------------------------------------------
  //--------------------------------------------------------------------
  StartNewStep("Select labels");
  // Keep right labels
  working_image = LabelizeAndSelectLabels<InternalImageType>
    (working_image, 
     GetBackgroundValue(), 
     GetForegroundValue(), 
     false, 
     GetMinimalComponentSize(), 
     GetLabelizeParameters2());

  // Set output
  StopCurrentStep<InternalImageType>(working_image);
  
  //--------------------------------------------------------------------
  //--------------------------------------------------------------------
  if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
    StartNewStep("Remove the trachea");
    // Set the trachea
    working_image = SetBackground<InternalImageType, InternalImageType>
      (working_image, trachea_tmp, 1, -1);
  
   // Dilate the trachea 
    static const unsigned int Dim = InputImageType::ImageDimension;
    typedef itk::BinaryBallStructuringElement<InternalPixelType, Dim> KernelType;
    KernelType structuringElement;
    structuringElement.SetRadius(GetRadiusForTrachea());
    structuringElement.CreateStructuringElement();
    typedef clitk::ConditionalBinaryDilateImageFilter<InternalImageType, InternalImageType, KernelType> ConditionalBinaryDilateImageFilterType;
    typename ConditionalBinaryDilateImageFilterType::Pointer dilateFilter = ConditionalBinaryDilateImageFilterType::New();
    dilateFilter->SetBoundaryToForeground(false);
    dilateFilter->SetKernel(structuringElement);
    dilateFilter->SetBackgroundValue (1);
    dilateFilter->SetForegroundValue (-1);
    dilateFilter->SetInput (working_image);
    dilateFilter->Update();
    working_image = dilateFilter->GetOutput();  
    
    // Set trachea with dilatation
    trachea_tmp = SetBackground<InternalImageType, InternalImageType>
      (trachea_tmp, working_image, -1, this->GetForegroundValue()); 

    // Remove the trachea
    working_image = SetBackground<InternalImageType, InternalImageType>
      (working_image, working_image, -1, this->GetBackgroundValue());  
    
    // Label
    working_image = LabelizeAndSelectLabels<InternalImageType>
      (working_image, 
       GetBackgroundValue(), 
       GetForegroundValue(), 
       false, 
       GetMinimalComponentSize(), 
       GetLabelizeParameters3());
    
    // Set output
    StopCurrentStep<InternalImageType>(working_image);
  }
  

  //--------------------------------------------------------------------
  //--------------------------------------------------------------------
  typedef clitk::AutoCropFilter<InternalImageType> CropFilterType;
  typename CropFilterType::Pointer cropFilter = CropFilterType::New();
  if (m_Seeds.size() != 0) { // if ==0 ->no trachea found
    StartNewStep("Croping trachea");
    cropFilter->SetInput(trachea_tmp);
    cropFilter->Update(); // Needed
    typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
    typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
    caster->SetInput(cropFilter->GetOutput());
    caster->Update();   
    trachea = caster->GetOutput();
    StopCurrentStep<MaskImageType>(trachea);  
  }

  //--------------------------------------------------------------------
  //--------------------------------------------------------------------
  StartNewStep("Croping lung");
  cropFilter = CropFilterType::New(); // Needed to reset pipeline
  cropFilter->SetInput(working_image);
  cropFilter->Update();   
  working_image = cropFilter->GetOutput();
  StopCurrentStep<InternalImageType>(working_image);

  //--------------------------------------------------------------------
  //--------------------------------------------------------------------
  StartNewStep("Separate Left/Right lungs");
  // Initial label
  working_image = Labelize<InternalImageType>(working_image, 
                                              GetBackgroundValue(), 
                                              false, 
                                              GetMinimalComponentSize());

  // Count the labels
  typedef itk::StatisticsImageFilter<InternalImageType> StatisticsImageFilterType;
  typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
  statisticsImageFilter->SetInput(working_image);
  statisticsImageFilter->Update();
  unsigned int initialNumberOfLabels = statisticsImageFilter->GetMaximum();
  working_image = statisticsImageFilter->GetOutput();	
 
  // Decompose the first label
  static const unsigned int Dim = InputImageType::ImageDimension;
  if (initialNumberOfLabels<2) {
    // Structuring element radius
    typename InputImageType::SizeType radius;
    for (unsigned int i=0;i<Dim;i++) radius[i]=1;
    typedef clitk::DecomposeAndReconstructImageFilter<InternalImageType,InternalImageType> DecomposeAndReconstructFilterType;
    typename DecomposeAndReconstructFilterType::Pointer decomposeAndReconstructFilter=DecomposeAndReconstructFilterType::New();
    decomposeAndReconstructFilter->SetInput(working_image);
    decomposeAndReconstructFilter->SetVerbose(false);
    decomposeAndReconstructFilter->SetRadius(radius);
    decomposeAndReconstructFilter->SetMaximumNumberOfLabels(2);
    decomposeAndReconstructFilter->SetMinimumObjectSize(this->GetMinimalComponentSize());
    decomposeAndReconstructFilter->SetMinimumNumberOfIterations(1);
    decomposeAndReconstructFilter->SetBackgroundValue(this->GetBackgroundValue());
    decomposeAndReconstructFilter->SetForegroundValue(this->GetForegroundValue());
    decomposeAndReconstructFilter->SetFullyConnected(true);
    decomposeAndReconstructFilter->SetNumberOfNewLabels(1);
    decomposeAndReconstructFilter->Update();
    working_image = decomposeAndReconstructFilter->GetOutput();      
  }

  // Retain labels (lungs)
  typedef itk::ThresholdImageFilter<InternalImageType> ThresholdImageFilterType;
  typename ThresholdImageFilterType::Pointer thresholdFilter = ThresholdImageFilterType::New();
  thresholdFilter->SetInput(working_image);
  thresholdFilter->ThresholdAbove(2);
  thresholdFilter->SetOutsideValue(this->GetBackgroundValue());
  thresholdFilter->Update();
  working_image = thresholdFilter->GetOutput();
  StopCurrentStep<InternalImageType> (working_image);
}
//--------------------------------------------------------------------


//--------------------------------------------------------------------
template <class TInputImageType, class TMaskImageType>
void 
clitk::ExtractLungFilter<TInputImageType, TMaskImageType>::
GenerateData() {
  // Final Cast 
  typedef itk::CastImageFilter<InternalImageType, MaskImageType> CastImageFilterType;
  typename CastImageFilterType::Pointer caster= CastImageFilterType::New();
  caster->SetInput(working_image);
  caster->Update();
  // Set output
  //this->SetNthOutput(0, caster->GetOutput()); // -> no because redo filter otherwise
  this->GraftOutput(caster->GetOutput());
  return;
}
//--------------------------------------------------------------------

  
#endif //#define CLITKBOOLEANOPERATORLABELIMAGEFILTER_TXX