Commit 4417ae3d authored by Stéphane Albert's avatar Stéphane Albert

ENH: Added HistogramModel::Percentile() to compute quantile percentage given...

ENH: Added HistogramModel::Percentile() to compute quantile percentage given intensity; Changed granularity of histogram (5x oversampling 256 RGBA values); Added mvdMath.h/.cxx to group math formulas (such as linear interpolation (lerp).
parent 62631e8a
......@@ -269,7 +269,7 @@ ColorDynamicsController
assert( imageModel->GetHistogramModel()!=NULL );
imageModel->GetSettings().DynamicsParam( 2 * channel ) =
imageModel->GetHistogramModel()->GetQuantile( channel, 0.01 * value );
imageModel->GetHistogramModel()->Quantile( channel, 0.01 * value );
emit ModelUpdated();
}
......@@ -286,7 +286,7 @@ ColorDynamicsController
assert( imageModel->GetHistogramModel()!=NULL );
imageModel->GetSettings().DynamicsParam( 2 * channel + 1 ) =
imageModel->GetHistogramModel()->GetQuantile( channel, 0.01 * value );
imageModel->GetHistogramModel()->Quantile( channel, 0.01 * value );
emit ModelUpdated();
}
......
......@@ -69,6 +69,74 @@ HistogramModel
{
}
/*******************************************************************************/
double
HistogramModel
::Percentile( CountType band, MeasurementType intensity, Bound bound ) const
{
// Get histogram of band.
Histogram::Pointer histogram( m_Histograms->GetNthElement( band ) );
// Contruct 1D measurement vector.
Histogram::MeasurementVectorType measurement( 1 );
measurement[ 0 ] = intensity;
// Get the index of measurement in 1D-histogram.
Histogram::IndexType index;
if( !histogram->GetIndex( measurement, index ) )
throw itk::RangeError( __FILE__, __LINE__ );
assert( Histogram::IndexType::GetIndexDimension()==1 );
// Min/max intensities of bin.
MeasurementType minI = histogram->GetBinMin( 0, index[ 0 ] );
MeasurementType maxI = histogram->GetBinMax( 0, index[ 0 ] );
// Min/max frequencies of
Histogram::FrequencyType frequency( histogram->GetFrequency( index ) );
// Initialize result.
double percent = frequency * ( intensity - minI ) / ( maxI - minI );
// Number of bins of histogram.
Histogram::SizeType::SizeValueType binCount = histogram->Size();
// Initialize bound indices.
assert( index[ 0 ]>=0 );
Histogram::SizeType::SizeValueType index0 = index[ 0 ];
Histogram::SizeType::SizeValueType i0 = 0;
Histogram::SizeType::SizeValueType iN = binCount;
switch( bound )
{
case BOUND_LOWER:
i0 = 0;
iN = index[ 0 ];
break;
case BOUND_UPPER:
i0 = index0 < binCount ? index0 + 1 : binCount;
iN = binCount;
break;
default:
assert( false && "Implement case statement of switch instruction." );
break;
};
// Traverse lower/upper bound.
Histogram::SizeType::SizeValueType i;
for( i=i0; i<iN; ++i )
percent += histogram->GetFrequency( i, 0 );
// Calculate frequency rate.
percent /= histogram->GetTotalFrequency();
// Return frequency rate.
return percent;
}
/*******************************************************************************/
void
HistogramModel
......
......@@ -53,6 +53,8 @@
/*****************************************************************************/
/* PRE-DECLARATION SECTION */
#define BINS_OVERSAMPLING_RATE 5
//
// External classes pre-declaration.
namespace
......@@ -84,6 +86,21 @@ class Monteverdi2_EXPORT HistogramModel :
//
// Public types.
public:
/** */
typedef
// itk::NumericTraits< T >::FloatType and
// itk::NumericTraits< T >::RealType do not depend on template
// parameter T. They are always typedef, respectively, as float
// and double.
//
// So, itk::NumericTraits< DefaultImageType::InternalPixelType
// >::RealType is equivalent to itk::NumericTraits< float
// >::RealType which is always an alias of double.
//
// This typedef is used for compatibility with
// itk::Histogram<>::MeasurementType.
itk::NumericTraits< DefaultImageType::InternalPixelType >::RealType
MeasurementType;
//
// Public methods.
......@@ -96,7 +113,13 @@ public:
virtual ~HistogramModel();
/** */
inline double GetQuantile( unsigned int band, double p ) const;
inline MeasurementType Quantile( CountType band, double p ) const;
/** */
double
Percentile( CountType band,
MeasurementType intensity,
Bound bound ) const;
/** */
inline DefaultImageType::PixelType GetMinPixel() const;
......@@ -131,21 +154,6 @@ protected:
//
// Private types.
private:
/** */
typedef
// itk::NumericTraits< T >::FloatType and
// itk::NumericTraits< T >::RealType do not depend on template
// parameter T. They are always typedef, respectively, as float
// and double.
//
// So, itk::NumericTraits< DefaultImageType::InternalPixelType
// >::RealType is equivalent to itk::NumericTraits< float
// >::RealType which is always an alias of double.
//
// This typedef is used for compatibility with
// itk::Histogram<>::MeasurementType.
itk::NumericTraits< DefaultImageType::InternalPixelType >::RealType
MeasurementType;
/** */
typedef itk::Statistics::Histogram< MeasurementType, 1 > Histogram;
......@@ -161,6 +169,7 @@ private:
template< typename TImageModel >
void template_BuildModel_M();
//
// Private attributes.
private:
......@@ -219,10 +228,10 @@ HistogramModel
/*******************************************************************************/
inline
double
HistogramModel::MeasurementType
HistogramModel
::GetQuantile( unsigned int band,
double p ) const
::Quantile( unsigned int band,
double p ) const
{
assert( band<m_Histograms->Size() );
......@@ -231,7 +240,6 @@ HistogramModel
/*******************************************************************************/
template< typename TImage >
inline
void
HistogramModel
::template_BuildModel_I()
......@@ -309,6 +317,7 @@ HistogramModel
// Setup histogram filter.
histogramFilter->GetFilter()->SetHistogramMin( m_MinPixel );
histogramFilter->GetFilter()->SetHistogramMax( m_MaxPixel );
histogramFilter->SetNumberOfBins( BINS_OVERSAMPLING_RATE * 256 );
histogramFilter->GetFilter()->SetSubSamplingRate( 1 );
// Go.
......@@ -329,7 +338,6 @@ HistogramModel
/*******************************************************************************/
template< typename TImageModel >
inline
void
HistogramModel
::template_BuildModel_M()
......@@ -403,6 +411,7 @@ HistogramModel
// Setup histogram filter.
histogramFilter->GetFilter()->SetHistogramMin( m_MinPixel );
histogramFilter->GetFilter()->SetHistogramMax( m_MaxPixel );
histogramFilter->GetFilter()->SetNumberOfBins( BINS_OVERSAMPLING_RATE * 256 );
histogramFilter->GetFilter()->SetSubSamplingRate( 1 );
// Go.
......
/*=========================================================================
Program: Monteverdi2
Language: C++
Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
See Copyright.txt for details.
Monteverdi2 is distributed under the CeCILL licence version 2. See
Licence_CeCILL_V2-en.txt or
http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt for more details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#include "mvdMath.h"
/*****************************************************************************/
/* INCLUDE SECTION */
//
// Qt includes (sorted by alphabetic order)
//// Must be included before system/custom includes.
//
// System includes (sorted by alphabetic order)
//
// ITK includes (sorted by alphabetic order)
//
// OTB includes (sorted by alphabetic order)
//
// Monteverdi includes (sorted by alphabetic order)
namespace mvd
{
/*
TRANSLATOR mvd::Math
Necessary for lupdate to be aware of C++ namespaces.
Context comment for translator.
*/
/*******************************************************************************/
} // end namespace 'mvd'
/*=========================================================================
Program: Monteverdi2
Language: C++
Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
See Copyright.txt for details.
Monteverdi2 is distributed under the CeCILL licence version 2. See
Licence_CeCILL_V2-en.txt or
http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt for more details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#ifndef __mvdMath_h
#define __mvdMath_h
//
// Configuration include.
//// Included at first position before any other ones.
#include "ConfigureMonteverdi2.h"
/*****************************************************************************/
/* INCLUDE SECTION */
#define USE_VNL 0
//
// Qt includes (sorted by alphabetic order)
//// Must be included before system/custom includes.
//
// System includes (sorted by alphabetic order)
//
// ITK includes (sorted by alphabetic order)
#if USE_VNL
#include "vnl_vector.h"
#endif
//
// OTB includes (sorted by alphabetic order)
//
// Monteverdi includes (sorted by alphabetic order)
#include "mvdTypes.h"
/*****************************************************************************/
/* PRE-DECLARATION SECTION */
//
// External classes pre-declaration.
namespace
{
}
namespace mvd
{
//
// Internal classes pre-declaration.
} // end of namespace 'mvd'.
/*****************************************************************************/
/* FUNCTIONS DECLARATION. */
namespace otb
{
} // end namespace 'otb'
namespace mvd
{
/**
*/
template< typename X, typename Y, typename K >
void
Lerp2( X& x, Y& y,
const K& k,
const X& x0, const Y& y0,
const X& x1, const Y& y1 );
/**
*/
template< typename X, typename Y >
const Y&
Lerp2( const X& x,
const X& x0, const Y& y0,
const X& x1, const Y& y1 );
/**
*/
#if USE_VNL
template< typename T, unsigned int N >
inline
const vnl_vector< T, N >&
Lerp( const T& k,
const vnl_vector< T, N >& v0,
const vnl_vector< T, N >& v1 );
#endif // USE_VNL
/**
*/
#if USE_VNL
template< typename T >
inline
const T&
Lerp( const T& k,
const T& x0, const T& x1,
const T& y0, const T& y1 );
#endif // USE_VNL
} // end namespace 'mvd'.
/*****************************************************************************/
/* INLINE SECTION */
namespace otb
{
} // end namespace 'otb'.
namespace mvd
{
/*******************************************************************************/
template< typename X, typename Y, typename K >
void
Lerp2( X& x, Y& y,
const K& k,
const X& x0, const Y& y0,
const X& x1, const Y& y1 )
{
const K& _1_minus_k( 1 - k );
x = k * x1 + _1_minus_k * x0;
y = k * y1 + _1_minus_k * y1;
}
/*******************************************************************************/
template< typename X, typename Y, typename K >
const Y&
Lerp2( const X& x,
const X& x0, const Y& y0,
const X& x1, const Y& y1 )
{
return y0 + (x - x0) * (y1 - y0) / (x1 - x0);
}
/*******************************************************************************/
#if USE_VNL
template< typename T >
inline
const vnl_vector< T, N >&
Lerp( const T& k,
const vnl_vector< T, N >& v0,
const vnl_vector< T, N >& v1 )
{
return k * v1 + ( T( 1 ) - k ) * v2;
}
#endif // USE_VNL
/*******************************************************************************/
#if USE_VNL
template< typename T >
inline
const T&
Lerp( const T& k,
const T& x0, const T& x1,
const T& y0, const T& y1 )
{
return Lerp( k, vnl_vector< T, 2 >( x0, x1 ), vnl_vector< T, 2 >( y0, y1 ) );
}
#endif // USE_VNL
} // end namespace 'mvd'
#endif // __mvdMath_h
......@@ -46,5 +46,10 @@ namespace mvd
*/
/*******************************************************************************/
const char*
BOUND_NAMES[ BOUND_COUNT ] = {
"BOUND_LOWER",
"BOUND_UPPER"
};
} // end namespace 'mvd'
......@@ -170,6 +170,25 @@ itkNumericTraitsGenericArrayDimensionsMacro( mvd::FixedArray, unsigned long long
namespace mvd
{
/*******************************************************************************/
/* Enumerations */
/**
*/
enum Bound
{
BOUND_LOWER = 0,
BOUND_UPPER,
BOUND_COUNT
};
/**
* Constant naming bound values.
*/
extern
const char*
BOUND_NAMES[ BOUND_COUNT ];
/*******************************************************************************/
/* Type definitions of scalar values. */
/**
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment