text: /*=========================================================================
   *
   *  Copyright Insight Software Consortium
   *
   *  Licensed under the Apache License, Version 2.0 (the "License");
   *  you may not use this file except in compliance with the License.
   *  You may obtain a copy of the License at
   *
   *         http://www.apache.org/licenses/LICENSE-2.0.txt
   *
   *  Unless required by applicable law or agreed to in writing, software
   *  distributed under the License is distributed on an "AS IS" BASIS,
   *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   *  See the License for the specific language governing permissions and
   *  limitations under the License.
   *
   *=========================================================================*/
  #ifndef __itkImageToHistogramFilter_hxx
  #define __itkImageToHistogramFilter_hxx
  
  #include "itkImageToHistogramFilter.h"
  #include "itkImageRegionConstIterator.h"
  
  namespace itk
  {
  namespace Statistics
  {
  template< typename TImage >
  ImageToHistogramFilter< TImage >
  ::ImageToHistogramFilter()
  {
    this->SetNumberOfRequiredInputs(1);
    this->SetNumberOfRequiredOutputs(1);
  
    this->ProcessObject::SetNthOutput( 0, this->MakeOutput(0) );
  
    // same default values as in the HistogramGenerator
  
    typename SimpleDataObjectDecorator<HistogramMeasurementType>::Pointer marginalScale =
      SimpleDataObjectDecorator<HistogramMeasurementType>::New();
    marginalScale->Set(100);
    this->ProcessObject::SetInput( "MarginalScale", marginalScale );
  
    SimpleDataObjectDecorator<bool>::Pointer autoMinMax =
      SimpleDataObjectDecorator<bool>::New();
    if( typeid(ValueType) == typeid(signed char) || typeid(ValueType) == typeid(unsigned char) )
      {
      autoMinMax->Set(false);
      }
    else
      {
      autoMinMax->Set(true);
      }
     this->ProcessObject::SetInput( "AutoMinimumMaximum", autoMinMax );
  }
  
  template< typename TImage >
  DataObject::Pointer
  ImageToHistogramFilter< TImage >
  ::MakeOutput( DataObjectPointerArraySizeType itkNotUsed(idx) )
  {
    return HistogramType::New().GetPointer();
  }
  
  template< typename TImage >
  const typename ImageToHistogramFilter< TImage >::HistogramType *
  ImageToHistogramFilter< TImage >
  ::GetOutput() const
  {
    const HistogramType *output =
      itkDynamicCastInDebugMode< const HistogramType * >( this->ProcessObject::GetOutput(0) );
  
    return output;
  }
  
  template< typename TImage >
  typename ImageToHistogramFilter< TImage >::HistogramType *
  ImageToHistogramFilter< TImage >
  ::GetOutput()
  {
    HistogramType *output =
      static_cast< HistogramType * >( this->ProcessObject::GetOutput(0) );
  
    return output;
  }
  
  template< typename TImage >
  void
  ImageToHistogramFilter< TImage >
  ::GraftOutput(DataObject *graft)
  {
    DataObject *output =
      const_cast< HistogramType * >( this->GetOutput() );
  
    // Call Histogram to copy meta-information, and the container
    output->Graft(graft);
  }
  
  
  template< typename TImage >
  void
  ImageToHistogramFilter< TImage >
  ::BeforeThreadedGenerateData()
  {
    // find the actual number of threads
    long nbOfThreads = this->GetNumberOfThreads();
    if ( itk::MultiThreader::GetGlobalMaximumNumberOfThreads() != 0 )
      {
      nbOfThreads = vnl_math_min( this->GetNumberOfThreads(), itk::MultiThreader::GetGlobalMaximumNumberOfThreads() );
      }
    // number of threads can be constrained by the region size, so call the
    // SplitRequestedRegion
    // to get the real number of threads which will be used
    RegionType splitRegion;  // dummy region - just to call the following method
    nbOfThreads = this->SplitRequestedRegion(0, nbOfThreads, splitRegion);
  
    // and allocate one histogram per thread
    m_Histograms.resize(nbOfThreads);
    m_Minimums.resize(nbOfThreads);
    m_Maximums.resize(nbOfThreads);
    m_Barrier = Barrier::New();
    m_Barrier->Initialize(nbOfThreads);
  }
  
  
  template< typename TImage >
  void
  ImageToHistogramFilter< TImage >
  ::ThreadedGenerateData(const RegionType & inputRegionForThread, ThreadIdType threadId)
  {
    long nbOfPixels = inputRegionForThread.GetNumberOfPixels();
    if( this->GetAutoMinimumMaximumInput() && this->GetAutoMinimumMaximum() )
      {
      // we'll have to iterate over all the pixels 2 times
      nbOfPixels *= 2;
      }
    ProgressReporter progress( this, threadId, nbOfPixels );
  
    if( threadId == 0 )
      {
      // just use the main one
      m_Histograms[threadId] = this->GetOutput();
      }
    else
      {
      m_Histograms[threadId] = HistogramType::New();
      }
    HistogramType * hist = m_Histograms[threadId];
    hist->SetClipBinsAtEnds(true);
  
    // the parameter needed to initialize the histogram
    unsigned int nbOfComponents = this->GetInput()->GetNumberOfComponentsPerPixel();
    HistogramSizeType size( nbOfComponents );
    HistogramMeasurementVectorType min( nbOfComponents );
    HistogramMeasurementVectorType max( nbOfComponents );
    if( this->GetHistogramSizeInput() )
      {
      // user provided value
      size = this->GetHistogramSize();
      }
    else
      {
      // use a default value, which must be computed at run time for the VectorImage
      size.Fill(256);
      }
  
    if( this->GetAutoMinimumMaximumInput() && this->GetAutoMinimumMaximum() )
      {
      // we have to compute the minimum and maximum values
      this->ThreadedComputeMinimumAndMaximum( inputRegionForThread, threadId, progress );
  
      // wait for the other threads to complete their part
      m_Barrier->Wait();
  
      // a non multithreaded part
      if( threadId == 0 )
        {
        min = m_Minimums[0];
        max = m_Maximums[0];
        for( unsigned int t=1; t<m_Minimums.size(); t++ )
          {
          for( unsigned int i=0; i<nbOfComponents; i++ )
            {
            min[i] = std::min( min[i], m_Minimums[t][i] );
            max[i] = std::max( max[i], m_Maximums[t][i] );
            }
          }
        this->ApplyMarginalScale( min, max, size );
        // store the values so they can be retreived by the other threads
        m_Minimums[0] = min;
        m_Maximums[0] = max;
        }
  
      // wait for all the threads to complete
      m_Barrier->Wait();
  
      // and get the values computed in the main thread
      min = m_Minimums[0];
      max = m_Maximums[0];
      }
    else
      {
      if( this->GetHistogramBinMinimumInput() )
        {
        min = this->GetHistogramBinMinimum();
        }
      else
        {
        min.Fill( NumericTraits<ValueType>::NonpositiveMin() - 0.5 );
        }
      if( this->GetHistogramBinMaximumInput() )
        {
        max = this->GetHistogramBinMaximum();
        }
      else
        {
        max.Fill( NumericTraits<ValueType>::max() + 0.5 );
        // this->ApplyMarginalScale( min, max, size );
        }
      }
  
    // finally, initialize the histogram
    hist->SetMeasurementVectorSize( nbOfComponents );
    hist->Initialize( size, min, max );
  
    // now fill the histograms
    this->ThreadedComputeHistogram( inputRegionForThread, threadId, progress );
  }
  
  
  template< typename TImage >
  void
  ImageToHistogramFilter< TImage >
  ::AfterThreadedGenerateData()
  {
    // group the results in the output histogram
    HistogramType * hist = m_Histograms[0];
    for( unsigned int i=1; i<m_Histograms.size(); i++ )
      {
      typedef typename HistogramType::ConstIterator         HistogramIterator;
      typedef typename HistogramType::AbsoluteFrequencyType AbsoluteFrequencyType;
  
      HistogramIterator hit = m_Histograms[i]->Begin();
      HistogramIterator end = m_Histograms[i]->End();
      while ( hit != end )
        {
        hist->IncreaseFrequencyOfMeasurement( hit.GetMeasurementVector(), hit.GetFrequency() );
        ++hit;
        }
      }
  
    // and drop the temporary histograms
    m_Histograms.clear();
    m_Minimums.clear();
    m_Maximums.clear();
    m_Barrier = NULL;
  }
  
  
  template< typename TImage >
  void
  ImageToHistogramFilter< TImage >
  ::ThreadedComputeMinimumAndMaximum(const RegionType & inputRegionForThread, ThreadIdType threadId, ProgressReporter & progress )
  {
    unsigned int nbOfComponents = this->GetInput()->GetNumberOfComponentsPerPixel();
    HistogramMeasurementVectorType min( nbOfComponents );
    HistogramMeasurementVectorType max( nbOfComponents );
  
    ImageRegionConstIterator< TImage > inputIt( this->GetInput(), inputRegionForThread );
    inputIt.GoToBegin();
    HistogramMeasurementVectorType m( nbOfComponents );
  
    min.Fill( NumericTraits<ValueType>::max() );
    max.Fill( NumericTraits<ValueType>::NonpositiveMin() );
    while ( !inputIt.IsAtEnd() )
      {
      const PixelType & p = inputIt.Get();
      NumericTraits<PixelType>::AssignToArray( p, m );
      for( unsigned int i=0; i<nbOfComponents; i++ )
        {
        min[i] = std::min( m[i], min[i] );
        max[i] = std::max( m[i], max[i] );
        }
      progress.CompletedPixel();  // potential exception thrown here
      ++inputIt;
      }
    m_Minimums[threadId] = min;
    m_Maximums[threadId] = max;
  }
  
  template< typename TImage >
  void
  ImageToHistogramFilter< TImage >
  ::ThreadedComputeHistogram(const RegionType & inputRegionForThread, ThreadIdType threadId, ProgressReporter & progress )
  {
    unsigned int nbOfComponents = this->GetInput()->GetNumberOfComponentsPerPixel();
    ImageRegionConstIterator< TImage > inputIt( this->GetInput(), inputRegionForThread );
    inputIt.GoToBegin();
    HistogramMeasurementVectorType m( nbOfComponents );
  
    while ( !inputIt.IsAtEnd() )
      {
      const PixelType & p = inputIt.Get();
      NumericTraits<PixelType>::AssignToArray( p, m );
      m_Histograms[threadId]->IncreaseFrequencyOfMeasurement( m, 1 );
      ++inputIt;
      progress.CompletedPixel();  // potential exception thrown here
      }
  }
  
  template< typename TImage >
  void
  ImageToHistogramFilter< TImage >
  ::ApplyMarginalScale( HistogramMeasurementVectorType & min, HistogramMeasurementVectorType & max, HistogramSizeType & size )
  {
    unsigned int nbOfComponents = this->GetInput()->GetNumberOfComponentsPerPixel();
    bool clipHistograms = true;
    for ( unsigned int i = 0; i < nbOfComponents; i++ )
      {
      if ( !NumericTraits< HistogramMeasurementType >::is_integer )
        {
        HistogramMeasurementType marginalScale = this->GetMarginalScale();
        const double margin =
          ( static_cast< HistogramMeasurementType >( max[i] - min[i] )
            / static_cast< HistogramMeasurementType >( size[i] ) )
          / static_cast< HistogramMeasurementType >( marginalScale );
  
        // Now we check if the max[i] value can be increased by
        // the margin value without saturating the capacity of the
        // HistogramMeasurementType
        if ( ( NumericTraits< HistogramMeasurementType >::max() - max[i] ) > margin )
          {
          max[i] = static_cast< HistogramMeasurementType >( max[i] + margin );
          }
        else
          {
          // an overflow would occur if we add 'margin' to the max
          // therefore we just compromise in setting max = max.
          // Histogram measurement type would force the clipping the max
          // value.
          // Therefore we must call the following to include the max value:
          clipHistograms = false;
          // The above function is okay since here we are within the
          // autoMinMax
          // computation and clearly the user intended to include min and max.
          }
        }
      else
        {
        // max[i] = SafeAssign(max[i] + NumericTraits<MeasurementType>::One);
        // if ( max[i] <= max[i] )
        if(max[i] <
            (static_cast<ValueType>
            (NumericTraits<HistogramMeasurementType>::max()) -
            NumericTraits<ValueType>::OneValue()))
          {
          max[i] = static_cast<HistogramMeasurementType>
            (max[i] + NumericTraits<ValueType>::OneValue());
          }
        else
          {
          // an overflow would have occurred, therefore set max to max
          // Histogram measurement type would force the clipping the max
          // value.
          // Therefore we must call the following to include the max value:
          clipHistograms = false;
          // The above function is okay since here we are within the
          // autoMinMax
          // computation and clearly the user intended to include min and max.
          }
        }
      }
    if( clipHistograms == false )
      {
      for( unsigned int i=0; i<m_Histograms.size(); i++ )
        {
        m_Histograms[i]->SetClipBinsAtEnds(false);
        }
      }
  }
  
  template< typename TImage >
  void
  ImageToHistogramFilter< TImage >
  ::PrintSelf(std::ostream & os, Indent indent) const
  {
    Superclass::PrintSelf(os, indent);
    // m_HistogramBinMinimum
    os << indent << "HistogramBinMinimum: " << this->GetHistogramBinMinimumInput() << std::endl;
    // m_HistogramBinMaximum
    os << indent << "HistogramBinMaximum: " << this->GetHistogramBinMaximumInput() << std::endl;
    // m_MarginalScale
    os << indent << "MarginalScale: " << this->GetMarginalScaleInput() << std::endl;
    // m_AutoMinimumMaximum
    os << indent << "AutoMinimumMaximum: " << this->GetAutoMinimumMaximumInput() << std::endl;
    // m_HistogramSize
    os << indent << "HistogramSize: " << this->GetHistogramSizeInput() << std::endl;
  }
  } // end of namespace Statistics
  } // end of namespace itk
  
  #endif
  