diff --git a/Modules/Applications/AppFiltering/app/CMakeLists.txt b/Modules/Applications/AppFiltering/app/CMakeLists.txt index e1743080394eb81864737e76e958d0a237857f3c..fdd85cc2c869229032e0f6d299bf7703487b223b 100644 --- a/Modules/Applications/AppFiltering/app/CMakeLists.txt +++ b/Modules/Applications/AppFiltering/app/CMakeLists.txt @@ -22,3 +22,8 @@ otb_create_application( NAME Smoothing SOURCES otbSmoothing.cxx LINK_LIBRARIES ${${otb-module}_LIBRARIES}) + +otb_create_application( + NAME ContrastEnhancement + SOURCES otbContrastEnhancement.cxx + LINK_LIBRARIES ${${otb-module}_LIBRARIES}) diff --git a/Modules/Applications/AppFiltering/app/otbContrastEnhancement.cxx b/Modules/Applications/AppFiltering/app/otbContrastEnhancement.cxx new file mode 100644 index 0000000000000000000000000000000000000000..977d4e75a469f86c3d611ab3f7343487df7f6054 --- /dev/null +++ b/Modules/Applications/AppFiltering/app/otbContrastEnhancement.cxx @@ -0,0 +1,988 @@ +/* + * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) + * + * This file is part of Orfeo Toolbox + * + * https://www.orfeo-toolbox.org/ + * + * 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 + * + * Unless required by applicable law or agreened 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. + */ + +#include "otbWrapperApplication.h" +#include "otbWrapperApplicationFactory.h" + +#include "otbVectorImageToImageListFilter.h" +#include "otbImageListToVectorImageFilter.h" +#include "otbStreamingStatisticsVectorImageFilter.h" +#include "otbStreamingStatisticsImageFilter.h" +#include "otbUnaryFunctorImageFilter.h" +#include "itkStreamingImageFilter.h" +#include "otbInPlacePassFilter.h" + +#include "otbComputeHistoFilter.h" +#include "otbComputeGainLutFilter.h" +#include "otbApplyGainFilter.h" +#include "otbImageFileWriter.h" +#include "itkImageRegionIterator.h" +#include <string> +#include "otbStreamingHistogramVectorImageFilter.h" + +namespace otb +{ + +namespace Wrapper +{ + +namespace Functor +{ + +class LuminanceOperator +{ +typedef FloatVectorImageType::PixelType OutPixel; +typedef FloatVectorImageType::PixelType InPixel; +public: + LuminanceOperator() {} + unsigned int GetOutputSize() + { + return 1; + } + virtual ~LuminanceOperator() { } + + OutPixel operator() ( InPixel input ) + { + OutPixel out(1); + out[0] = m_LumCoef[0] * input[m_Rgb[0]] + + m_LumCoef[1] * input[m_Rgb[1]] + + m_LumCoef[2] * input[m_Rgb[2]] ; + return out; + } // end operator () + + + void SetRgb( std::vector<int> rgb) + { + m_Rgb = rgb; + } + + void SetLumCoef(std::vector<float> lumCoef) + { + m_LumCoef = lumCoef; + } + std::vector<float> GetLumCoef() + { + return m_LumCoef; + } +private: + std::vector<int> m_Rgb; + std::vector<float> m_LumCoef; +}; // end of functor class MultiplyOperator + +} // end of functor + +class ContrastEnhancement : public Application +{ +public: + /** Standard class typedefs. */ + typedef ContrastEnhancement Self; + typedef Application Superclass; + typedef itk::SmartPointer < Self > Pointer; + typedef itk::SmartPointer < const Self > ConstPointer; + + typedef otb::VectorImage < unsigned int , 2 > HistogramType; + typedef otb::VectorImage < double , 2 > LutType; + + typedef FloatImageType::PixelType ImagePixelType; + + typedef otb::ComputeHistoFilter < FloatImageType , + HistogramType > + HistoFilterType; + typedef otb::ComputeGainLutFilter < HistogramType , + LutType > + GainLutFilterType; + typedef otb::ApplyGainFilter < FloatImageType , + LutType , FloatImageType > + ApplyFilterType; + typedef otb::ImageList < FloatImageType > ImageListType; + + typedef otb::VectorImageToImageListFilter < FloatVectorImageType, + ImageListType > + VectorToImageListFilterType; + + typedef otb::ImageListToVectorImageFilter < ImageListType, + FloatVectorImageType > + ImageListToVectorFilterType; + + typedef otb::StreamingStatisticsVectorImageFilter < FloatVectorImageType > + VectorStatsFilterType; + + typedef otb::StreamingStatisticsImageFilter < FloatImageType > + StatsFilterType; + + typedef otb::UnaryFunctorImageFilter < FloatVectorImageType , + FloatVectorImageType , Functor::LuminanceOperator > + LuminanceFunctorType; + + typedef itk::StreamingImageFilter < LutType , LutType > + StreamingImageFilterType; + + typedef otb::InPlacePassFilter < FloatImageType > BufferFilterType; + + typedef otb::StreamingHistogramVectorImageFilter < FloatVectorImageType > + HistoPersistentFilterType; + + /** Standard macro */ + itkNewMacro( Self ); + + itkTypeMacro( ContrastEnhancement , otb::Application ); + +private: + + void DoInit() override + { + SetName("ContrastEnhancement"); + SetDescription("This application is the implementation of the histogram " + "equalization algorithm. It can be used to enhance contrast in an image " + "or to reduce the dynamic of th image without loosing too much contrast. " + "It offers several options as a no data value, " + "a contrast limitation factor, a local version of the algorithm and " + "also a mode to equalized the luminance of the image."); + + // Documentation + SetDocName("Contrast Enhancement"); + SetDocLongDescription("This application is the implementation of the " + "histogram equalization algorithm. The idea of the algorithm is to use " + "the whole available dynamic. In order to do so it computes a histogram " + "over the image and then use the whole dynamic : meaning flattening the " + "histogram. That gives us gain for each bin that transform the original " + "histogram into the flat one. This gain is then apply on the original " + "image. Upon this coarse algorithm we added several option to allow " + "a finer result. First there is the limitation of the contrast. Many " + "ways can be used to do it, we choose to limit the contrast by modifying " + "the original histogram. To do so we clip the histogram at a given " + "height and redistribute equally among the bins the clipped population. " + "Then we add a local version of the algorithm. It is possible to apply " + "the algorithm on tiles of the image. That gives us gain depending on " + "the value of the pixel and its position in the image. In order to " + "smoothen the result we interpolate the gain between tiles."); + SetDocLimitations("None"); + SetDocAuthors("OTB-Team"); + SetDocSeeAlso(" "); + + AddDocTag(Tags::Filter); + + AddParameter(ParameterType_InputImage, "in", "Input Image"); + SetParameterDescription("in", "Input image."); + + AddParameter(ParameterType_OutputImage, "out", "Output Image"); + SetParameterDescription("out", "Output image."); + + AddParameter(ParameterType_Int , "bins" , "Number of bin"); + SetDefaultParameterInt("bins", 256); + SetParameterDescription("bins", + "Number of bin used to create the histogram"); + + AddParameter(ParameterType_Float , "hfact" , "Contrast Limitation"); + SetParameterDescription("hfact","This parameter will set the maximum " + "height accepted in a bin on the input image histogram. " + "The maximum height will be computed as hfact*eqHeight where eqHeight " + "is the height of the theoretical flat histogram. The higher hfact, the " + "higher the contrast."); + MandatoryOff("hfact"); + + AddParameter(ParameterType_Float , "nodata" , "Nodata Value"); + SetParameterDescription("nodata","If there is a value in the " + "image that has no visualization meaning, it can be ignored by the " + "algorithm."); + MandatoryOff("nodata"); + + AddParameter(ParameterType_Choice , "spatial" , "Spatial parameters " + "for the histogram computation"); + AddChoice( "spatial.local" , "Local" ); + SetParameterDescription("spatial.local" , "The histograms will be " + "computed on the each thumbnail. Each of the histogram will be " + "equalized and the corresponding gain will be interpolated."); + AddChoice( "spatial.global" , "Global" ); + SetParameterDescription("spatial.global" , "The histogram will be " + "computed on the whole image. The equalization will be done on " + "this single histogram."); + + + AddParameter(ParameterType_Int,"spatial.local.h" , + "Thumbnail height in pixel"); + AddParameter(ParameterType_Int,"spatial.local.w" , + "Thumbnail width in pixel"); + + AddParameter(ParameterType_Choice , "minmax" , "Minimum and maximum " + "definition"); + SetParameterDescription("minmax","Minimum and maximum value that will " + "bound the histogram."); + AddChoice( "minmax.auto" , "Automatic" ); + SetParameterDescription("minmax.auto" , "Minimum and maximum value will " + "be computed on the image (nodata value won't be taken " + "into account) . Each band will have a minimum and a maximum."); + AddParameter(ParameterType_Empty, "minmax.auto.global", "Global"); + SetParameterDescription("minmax.auto.global" , "Automatic" + "Min/max computation will result in the same minimum and maximum for " + "all the bands."); + AddChoice( "minmax.manuel" , "Manuel" ); + SetParameterDescription("minmax.auto","Minimum and maximum value will be " + "set by the user"); + AddParameter(ParameterType_Float , "minmax.manuel.min" , "Minimum"); + AddParameter(ParameterType_Float , "minmax.manuel.max" , "Maximum"); + MandatoryOff("minmax.manuel.min"); + MandatoryOff("minmax.manuel.max"); + + AddParameter(ParameterType_Choice , "mode" , "What to equalized"); + AddChoice( "mode.each" , "Channels" ); + SetParameterDescription( "mode.each" , + "Each channel is equalized independently" ); + AddChoice( "mode.lum" , "Luminance" ); + SetParameterDescription( "mode.lum" , + "The luminance is equalized and then a gain is applied " + "on each channels. This gain for each channels will depend on" + "the weight (coef) of the channel in the luminance." ); + AddParameter(ParameterType_Group , "mode.lum.red" , "Red Channel" ); + AddParameter(ParameterType_Int , "mode.lum.red.ch" , "Red Channel" ); + SetDefaultParameterInt("mode.lum.red.ch", 0 ); + AddParameter(ParameterType_Float , "mode.lum.red.coef" , + "Value for luminance computation" ); + SetDefaultParameterFloat("mode.lum.red.coef", 0.21 ); + + AddParameter(ParameterType_Group , "mode.lum.green" , "Green Channel" ); + AddParameter(ParameterType_Int , "mode.lum.green.ch" , "Greenen Channel" ); + SetDefaultParameterInt("mode.lum.green.ch", 1 ); + AddParameter(ParameterType_Float , "mode.lum.green.coef" , + "Value for luminance computation" ); + SetDefaultParameterFloat("mode.lum.green.coef", 0.71 ); + + AddParameter(ParameterType_Group , "mode.lum.blue" , "Blue Channel" ); + AddParameter(ParameterType_Int , "mode.lum.blue.ch" , "Blue Channel" ); + SetDefaultParameterInt("mode.lum.blue.ch", 2 ); + AddParameter(ParameterType_Float , "mode.lum.blue.coef" , + "Value for luminance computation" ); + SetDefaultParameterFloat("mode.lum.blue.coef", 0.08 ); + + SetDefaultParameterInt( "spatial.local.w" , 256 ); + SetDefaultParameterInt( "spatial.local.h" , 256 ); + + SetMinimumParameterIntValue("mode.lum.red.ch", 0); + SetMinimumParameterIntValue("mode.lum.green.ch", 0); + SetMinimumParameterIntValue("mode.lum.blue.ch", 0); + SetMinimumParameterIntValue("bins", 2); + SetMinimumParameterIntValue("spatial.local.h", 1); + SetMinimumParameterIntValue("spatial.local.w", 1); + + SetExampleComment( "Local contrast enhancement by luminance" , 0 ); + SetDocExampleParameterValue( "in" , "couleurs.tif" ); + SetDocExampleParameterValue( "out" , "equalizedcouleurs.tif float" ); + SetDocExampleParameterValue( "bins" , "256" ); + SetDocExampleParameterValue( "spatial.local.w" , "500" ); + SetDocExampleParameterValue( "spatial.local.h" , "500"); + SetDocExampleParameterValue( "mode" , "lum" ); + + AddRAMParameter(); + } + + void DoUpdateParameters() override + { + if ( HasValue("in") ) + { + FloatVectorImageType * inImage = GetParameterImage("in"); + FloatVectorImageType::RegionType::SizeType size; + size = inImage->GetLargestPossibleRegion().GetSize() ; + + // if ( !HasUserValue("spatial.local.w") ) + // SetParameterInt( "spatial.local.w" , size[0] ); + + // if ( !HasUserValue("spatial.local.h") ) + // SetParameterInt( "spatial.local.h" , size[1] ); + + if ( GetParameterString("spatial") == "local" && + HasValue("spatial.local.h") && HasValue("spatial.local.w") && + HasValue("bins") ) + CheckValidity(); + + + if ( !HasUserValue("nodata") && IsParameterEnabled("nodata")) + SetDefaultValue( inImage , "NODATA" ); + + if ( GetParameterString( "mode" ) == "lum" && + !HasUserValue("mode.lum.red.ch") && + !HasUserValue("mode.lum.green.ch") && + !HasUserValue("mode.lum.blue.ch") ) + SetDefaultValue( inImage , "RGB" ); + + // if ( HasUserValue("minmax.manuel.min") && + // HasUserValue("minmax.manuel.max") ) + // { + // if ( GetParameterFloat( "minmax.manuel.min" ) > + // GetParameterFloat( "minmax.manuel.max" ) ) + // { + // float temp = GetParameterFloat( "minmax.manuel.min" ); + // SetParameterFloat( "minmax.manuel.min" , + // GetParameterFloat( "minmax.manuel.max" )); + // SetParameterFloat( "minmax.manuel.max" , temp ); + // } + // else if ( GetParameterFloat( "minmax.manuel.min" ) == + // GetParameterFloat( "minmax.manuel.max" ) ) + // { + // std::ostringstream oss; + // oss<<"Warning minimum and maximum are equal."<<std::endl; + // otbAppLogINFO( << oss.str() ); + // } + // } + } + + if ( GetParameterString("minmax") == "manuel" ) + { + MandatoryOn("minmax.manuel.min"); + MandatoryOn("minmax.manuel.max"); + } + else if ( GetParameterString("minmax") == "auto" ) + { + MandatoryOff("minmax.manuel.min"); + MandatoryOff("minmax.manuel.max"); + } + } + + void DoExecute() override + { + m_MinMaxMode = GetParameterString("minmax"); + m_EqMode = GetParameterString("mode"); + m_SpatialMode = GetParameterString("spatial"); + FloatVectorImageType * inImage = GetParameterImage("in"); + WarningGlobalOrNot( inImage ); + WarningMinMax(); + LogInfo(); + ImageListType::Pointer outputImageList ( ImageListType::New() ); + m_VectorToImageListFilter = VectorToImageListFilterType::New() ; + m_VectorToImageListFilter->SetInput( inImage ); + m_VectorToImageListFilter->UpdateOutputInformation(); + ImageListType::Pointer inputImageList = + m_VectorToImageListFilter->GetOutput(); + unsigned int nbChannel = inImage->GetVectorLength (); + + if ( m_EqMode == "each") + { + // Each channel will be equalized + m_GainLutFilter.resize(nbChannel); + m_ApplyFilter.resize(nbChannel); + m_BufferFilter.resize(nbChannel); + m_StreamingFilter.resize(nbChannel); + PerBandEqualization( inImage , inputImageList , + nbChannel , outputImageList ); + } + else if ( m_EqMode == "lum") + { + std::vector<int> rgb( 3 , 0 ); + rgb[0] = GetParameterInt("mode.lum.red.ch"); + rgb[1] = GetParameterInt("mode.lum.green.ch"); + rgb[2] = GetParameterInt("mode.lum.blue.ch"); + ComputeLuminance( inImage , rgb ); + LuminanceEqualization( inputImageList , rgb , outputImageList ); + } + + m_ImageListToVectorFilterOut = ImageListToVectorFilterType::New() ; + m_ImageListToVectorFilterOut->SetInput(outputImageList); + SetParameterOutputImage( "out" , + m_ImageListToVectorFilterOut->GetOutput() ); + } + + // Look for default values in the image metadata + void SetDefaultValue( const FloatVectorImageType * inImage , + std::string what) + { + typedef ImageMetadataInterfaceBase ImageMetadataInterfaceType; + ImageMetadataInterfaceType::Pointer metadataInterface = + ImageMetadataInterfaceFactory + ::CreateIMI(inImage->GetMetaDataDictionary()); + if ( what == "NODATA" ) + { + std::vector<double> values; + std::vector<bool> flags; + + bool ret = metadataInterface->GetNoDataFlags(flags,values); + + if(ret && !values.empty() && !flags.empty() && flags[0]) + { + SetParameterFloat( "nodata" , static_cast<float>( values[0] ) ); + } + else + { + SetParameterFloat( "nodata" , 0 ); + } + } + else if ( what == "RGB" ) + { + std::vector<unsigned int> rgb = + metadataInterface->GetDefaultDisplay() ; + unsigned int m = inImage->GetVectorLength (); + SetParameterInt( "mode.lum.red.ch" , rgb[0] ); + SetParameterInt( "mode.lum.green.ch" , rgb[1] ); + SetParameterInt( "mode.lum.blue.ch" , rgb[2] ); + if( m < rgb[ 0 ] ) + { + SetParameterFloat ("mode.lum.red.coef" , 0.0 ); + SetParameterInt( "mode.lum.red.ch" , 0 ); + } + if( m < rgb[ 1 ] ) + { + SetParameterFloat ("mode.lum.green.coef" , 0.0 ); + SetParameterInt( "mode.lum.gre.ch" , 0 ); + } + if( m < rgb[ 2 ] ) + { + SetParameterFloat ("mode.lum.blue.coef" , 0.0 ); + SetParameterInt( "mode.lum.blue.ch" , 0 ); + } + } + } + + // Log info for the user + void LogInfo() + { + std::ostringstream oss; + oss << "The application has been launched with the following parameters :" + <<std::endl; + oss << "- number of bins : " << GetParameterInt("bins") << std::endl; + if ( HasValue("hfact") ) + { + oss << "- contrast limtaition factor : " + << GetParameterFloat("hfact") << std::endl; + } + else + { + oss << "- no contrast limitation factor" << std::endl; + } + oss << "- spatial parameters : " << m_SpatialMode ; + if ( m_SpatialMode == "local" ) + { + oss<< " with a thumbnail of " << m_ThumbSize[0] <<" X "<< m_ThumbSize[1] ; + } + oss << std::endl << "- equalisation of "; + if ( m_EqMode == "each" ) + { + oss << "each channel"; + } + else + { + oss << "the luminance"; + } + oss << std::endl << "- Min/Max parameters : "; + if ( m_MinMaxMode == "auto" ) + { + oss << "automatic"; + if ( IsParameterEnabled( "minmax.auto.global" ) ) + { + oss << " and global"; + } + } + else + { + oss << GetParameterFloat("minmax.manuel.min") << "/" << + GetParameterFloat("minmax.manuel.max"); + } + + otbAppLogINFO( << oss.str() ); + } + + // Force global computation if the thumbsize is equal to the largest region + void WarningGlobalOrNot( FloatVectorImageType * input) + { + auto size = input->GetLargestPossibleRegion().GetSize(); + if ( m_SpatialMode == "global" ) + { + m_ThumbSize[0] = size[0]; + m_ThumbSize[1] = size[1]; + } + else + { + m_ThumbSize[0] = GetParameterInt("spatial.local.w"); + m_ThumbSize[1] = GetParameterInt("spatial.local.h"); + if ( size[0] == m_ThumbSize[0] && size[1] == m_ThumbSize[1] ) + { + std::ostringstream oss; + oss<<"Warning you choose to compute the histogram with a local " + "method whereas you have selected the whole image for the thumbnail " + "size. In order to use less memory consider using the global option for " + "histogram computation."; + otbAppLogWARNING( << oss.str() ); + } + } + } + + // Check for min max validity + void WarningMinMax() + { + if ( m_MinMaxMode == "manuel" && + GetParameterFloat( "minmax.manuel.min" ) > + GetParameterFloat( "minmax.manuel.max" ) ) + { + std::ostringstream oss; + oss<<"The minimum (" << GetParameterFloat( "minmax.manuel.min" ) << + ") is superior to the maximum (" + << GetParameterFloat( "minmax.manuel.max" ) + << ") please correct this error or allow the application to compute " + "those parameters"; + otbAppLogFATAL( << oss.str() ) + } + } + + // Check if the image size is a multiple of the thumbnail size + void CheckValidity() + { + std::ostringstream oss; + long nbPixel = GetParameterInt("spatial.local.w") * + GetParameterInt("spatial.local.h"); + if ( nbPixel < 10 * GetParameterInt("bins")) + { + oss<<"Warning in parameters selection the thumbnail size is small " + "in comparison with the number of bin. Histogram may not have much sens. " + "For better result enlarge thumbnail size or reduce number of bin."; + otbAppLogINFO( << oss.str() ); + } + } + + // Compute min max from a vector image + void ComputeVectorMinMax( const FloatVectorImageType::Pointer inImage , + FloatVectorImageType::PixelType & max , + FloatVectorImageType::PixelType & min ) + { + if ( m_MinMaxMode == "manuel" ) + { + min.Fill( GetParameterFloat("minmax.manuel.min") ); + max.Fill( GetParameterFloat("minmax.manuel.max") ); + } + else + { + VectorStatsFilterType::Pointer + statFilter ( VectorStatsFilterType::New() ); + statFilter->SetIgnoreInfiniteValues(true); + if( IsParameterEnabled("nodata") ) + { + statFilter->SetIgnoreUserDefinedValue(true); + statFilter->SetUserIgnoredValue( GetParameterFloat("nodata") ); + } + statFilter->SetInput( inImage ); + AddProcess(statFilter->GetStreamer(), "Computing statistics"); + statFilter->Update(); + min = statFilter->GetMinimum(); + max = statFilter->GetMaximum(); + if ( IsParameterEnabled("minmax.auto.global") ) + { + float temp(min[0]); + for ( unsigned int i = 1 ; i < min.GetSize() ; i++ ) + { + temp = std::min(temp , min[i]); + } + min.Fill(temp); + temp = max[0]; + for ( unsigned int i = 1 ; i < max.GetSize() ; i++ ) + { + temp = std::max(temp , max[i]); + } + max.Fill(temp); + } + } + std::ostringstream oss; + oss<<"Minimum and maximum are for each channel : "; + if ( IsParameterEnabled("minmax.auto.global") || + m_MinMaxMode == "manuel" ) + { + oss<<std::endl<<min[0]<<" and "<<max[0]; + } + else + { + for ( unsigned int i = 0 ; i < min.GetSize() ; i++ ) + { + oss<<std::endl<<min[i]<<" and "<<max[i]; + } + } + otbAppLogINFO( << oss.str() ); + } + + // Prepare the first half of the pipe that is common to every methode of + // equalization + void SetUpPipeline( unsigned int channel , + const FloatImageType::Pointer input ) + { + m_GainLutFilter[channel] = GainLutFilterType::New(); + m_ApplyFilter[channel] = ApplyFilterType::New(); + m_StreamingFilter[channel] = StreamingImageFilterType::New(); + m_BufferFilter[channel] = BufferFilterType::New(); + m_BufferFilter[channel]->SetInput( input ); + m_GainLutFilter[channel]->SetInput ( m_Histogram[channel] ); + m_StreamingFilter[channel]->SetInput( + m_GainLutFilter[channel]->GetOutput() ); + m_ApplyFilter[channel]->SetInputImage ( + m_BufferFilter[channel]->GetOutput() ); + m_ApplyFilter[channel]->SetInputLut( + m_StreamingFilter[channel]->GetOutput() ); + } + + // Function corresponding to the "each" mode + void PerBandEqualization( const FloatVectorImageType::Pointer inImage , + const ImageListType::Pointer inputImageList , + const unsigned int nbChannel, + ImageListType::Pointer outputImageList ) + { + FloatVectorImageType::PixelType min(nbChannel) , max(nbChannel); + min.Fill(0); + max.Fill(0); + ComputeVectorMinMax( inImage , max , min ); + + if ( m_SpatialMode == "global" ) + PersistentComputation( inImage , nbChannel , max , min ); + else + { + float thresh (-1); + if ( HasValue("hfact") ) + { + thresh = GetParameterInt("hfact"); + } + + m_HistoFilter.resize(nbChannel); + m_Histogram.resize(nbChannel); + for (unsigned int channel = 0 ; channel < nbChannel ; channel++) + { + m_HistoFilter[channel] = HistoFilterType::New(); + m_HistoFilter[channel]->SetInput( + inputImageList->GetNthElement(channel) ); + SetHistoFilterParameter( m_HistoFilter[channel] , + min[channel] , + max[channel] , + GetParameterInt("bins") , + thresh ); + m_Histogram[channel] = m_HistoFilter[channel]->GetHistoOutput(); + } + } + + for ( unsigned int channel = 0 ; channel < nbChannel ; channel++ ) + { + SetUpPipeline( channel , inputImageList->GetNthElement(channel) ); + + if ( min[channel] == max[channel] ) + { + std::ostringstream oss; + oss<< "Channel " <<channel << " is constant : " << "min = " << + min[channel] << " and max = " << max[channel] ; + otbAppLogINFO( << oss.str() ); + m_BufferFilter[channel]->SetInput( + inputImageList->GetNthElement(channel) ); + outputImageList->PushBack( m_BufferFilter[channel]->GetOutput() ); + continue; + } + + SetGainLutFilterParameter( m_GainLutFilter[channel] , + min[channel] , + max[channel]); + SetApplyFilterParameter( m_ApplyFilter[channel] , + min[channel] , + max[channel]); + + outputImageList->PushBack( m_ApplyFilter[channel]->GetOutput() ); + } + } + + // Compute the luminance with user parameters + void ComputeLuminance( const FloatVectorImageType::Pointer inImage , + std::vector < int > rgb ) + { + // Retreive coeffs for each channel + std::vector < float > lumCoef( 3 , 0.0 ); + lumCoef[0] = GetParameterFloat("mode.lum.red.coef"); + lumCoef[1] = GetParameterFloat("mode.lum.green.coef"); + lumCoef[2] = GetParameterFloat("mode.lum.blue.coef"); + // Normalize those coeffs + float sum = std::accumulate( lumCoef.begin() , lumCoef.end() , 0.0 ); + assert(sum>0); + for (int i = 0 ; i<3 ; i++ ) + { + lumCoef[i] /= sum; + } + m_LuminanceFunctor = LuminanceFunctorType::New() ; + m_LuminanceFunctor->GetFunctor().SetRgb( rgb ); + m_LuminanceFunctor->GetFunctor().SetLumCoef( lumCoef ); + m_LuminanceFunctor->SetInput( inImage ); + m_LuminanceFunctor->UpdateOutputInformation(); + } + + // Equalize the luminance and apply the corresponding gain on each channel + // used to compute this luminance + void LuminanceEqualization( const ImageListType::Pointer inputImageList , + const std::vector < int > rgb , + ImageListType::Pointer outputImageList ) + { + m_GainLutFilter.resize( 1 , GainLutFilterType::New() ); + m_HistoFilter.resize( 1 , HistoFilterType::New() ); + m_StreamingFilter.resize( 1 , StreamingImageFilterType::New() ); + m_ApplyFilter.resize(3); + m_BufferFilter.resize(3); + FloatVectorImageType::PixelType min(1) , max(1); + ComputeVectorMinMax( m_LuminanceFunctor->GetOutput() , max , min ); + + if ( m_SpatialMode == "global" ) + PersistentComputation( m_LuminanceFunctor->GetOutput() , 1 , max , min ); + else + { + float thresh (-1); + if ( HasValue("hfact") ) + { + thresh = GetParameterInt("hfact"); + } + m_Histogram.resize(1); + m_HistoFilter[0] = HistoFilterType::New(); + m_LuminanceToImageListFilter = VectorToImageListFilterType::New(); + m_LuminanceToImageListFilter->SetInput( + m_LuminanceFunctor->GetOutput() ); + SetHistoFilterParameter( m_HistoFilter[0] , + min[0] , + max[0] , + GetParameterInt("bins") , + thresh ); + m_LuminanceToImageListFilter->UpdateOutputInformation(); + m_HistoFilter[0]->SetInput( + m_LuminanceToImageListFilter->GetOutput()->GetNthElement(0) ) ; + m_Histogram[0] = m_HistoFilter[0]->GetHistoOutput(); + } + + m_GainLutFilter[0]->SetInput ( m_Histogram[0] ); + m_StreamingFilter[0]->SetInput( m_GainLutFilter[0]->GetOutput() ); + + SetGainLutFilterParameter( m_GainLutFilter[0] , + min[0] , + max[0] ); + + for ( int channel = 0 ; channel < 3 ; channel++ ) + { + + m_BufferFilter[channel] = BufferFilterType::New(); + m_ApplyFilter[channel] = ApplyFilterType::New(); + SetApplyFilterParameter( m_ApplyFilter[channel] , + min[0] , + max[0] ); + + m_BufferFilter[channel]->SetInput( + inputImageList->GetNthElement( rgb[channel] ) ); + m_ApplyFilter[channel]->SetInputImage ( + m_BufferFilter[channel]->GetOutput() ); + m_ApplyFilter[channel]->SetInputLut( + m_StreamingFilter[0]->GetOutput() ); + + + outputImageList->PushBack( m_ApplyFilter[channel]->GetOutput() ); + } + } + + // Function that compute histograms with HistoPersistentFilterType + void PersistentComputation( const FloatVectorImageType::Pointer inImage , + const unsigned int nbChannel , + const FloatVectorImageType::PixelType & max , + const FloatVectorImageType::PixelType & min) + { + + HistoPersistentFilterType::Pointer histoPersistent ( + HistoPersistentFilterType::New()); + unsigned int nbBin ( GetParameterInt( "bins" ) ); + histoPersistent->SetInput( inImage ); + FloatVectorImageType::PixelType pixel( nbChannel ) , step( nbChannel ); + pixel.Fill( nbBin ); + histoPersistent->GetFilter()->SetNumberOfBins( pixel ); + if ( IsParameterEnabled("nodata") ) + { + histoPersistent->GetFilter()->NoDataFlagOn(); + histoPersistent->GetFilter()->SetNoDataValue( + GetParameterFloat( "nodata" ) ); + } + step = ( max - min ) / nbBin; + pixel = min - 0.5 * step; + histoPersistent->GetFilter()->SetHistogramMin(pixel); + pixel = max + 0.5 * step; + histoPersistent->GetFilter()->SetHistogramMax(pixel); + AddProcess(histoPersistent->GetStreamer(), "Computing histogram"); + histoPersistent->Update(); + HistoPersistentFilterType::HistogramListType * histoList = + histoPersistent->GetHistogramList(); + + Transfer( histoList ); + + if ( HasValue("hfact") ) + { + Threshold( histoList , nbBin ); + } + } + + // Threshold function that is normaly done in ComputeHistoFilter and here is + // used on the output of HistoPersistentFilterType. + void Threshold( HistoPersistentFilterType::HistogramListType * histoList , + unsigned int nbBin ) + { + for ( unsigned int j = 0 ; j < histoList->Size() ; j++ ) + { + unsigned int rest(0) , height ( static_cast<unsigned int>( + GetParameterFloat( "hfact" ) * + ( histoList->GetNthElement(j)->GetTotalFrequency() / nbBin ) ) ); + + HistogramType::IndexType zero; + HistogramType::Pointer & histoToThresh = m_Histogram[j]; + zero.Fill(0); + for ( unsigned int i = 0 ; i < nbBin ; i++ ) + { + if ( histoToThresh->GetPixel(zero)[i] > height ) + { + rest += histoToThresh->GetPixel(zero)[i] - height ; + histoToThresh->GetPixel(zero)[i] = height ; + } + } + height = rest / nbBin; + rest = rest % nbBin; + for ( unsigned int i = 0 ; i < nbBin ; i++ ) + { + histoToThresh->GetPixel(zero)[i] += height ; + if ( i > (nbBin - rest)/2 && i <= (nbBin - rest)/2 + rest ) + { + ++histoToThresh->GetPixel(zero)[i]; + } + } + } + } + + // Transfer the output of the HistoPersistentFilterType to the input type + // of ComputeGainLutFilter + void Transfer( HistoPersistentFilterType::HistogramListType * histoList ) + { + unsigned int nbBin( GetParameterInt( "bins" ) ); + + HistoPersistentFilterType::HistogramType::Pointer histo; + FloatImageType::SpacingType inputSpacing ( + GetParameterImage("in")->GetSignedSpacing() ); + FloatImageType::PointType inputOrigin ( + GetParameterImage("in")->GetOrigin() ); + + HistogramType::SpacingType histoSpacing ; + histoSpacing[0] = inputSpacing[0] * m_ThumbSize[0] ; + histoSpacing[1] = inputSpacing[1] * m_ThumbSize[1] ; + HistogramType::PointType histoOrigin ; + histoOrigin[0] = histoSpacing[0] / 2 + + inputOrigin[0] - inputSpacing[0] / 2 ; + histoOrigin[1] = histoSpacing[1] / 2 + + inputOrigin[1] - inputSpacing[1] / 2 ; + + for ( unsigned int i = 0 ; i < histoList->Size() ; i++ ) + { + HistogramType::Pointer histoVectorImage = + CreateHistogramFrom( histoList->GetNthElement( i ) , + nbBin , + histoSpacing , + histoOrigin ); + + m_Histogram.push_back( histoVectorImage ); + } + } + + // Creating a histogram (VectorImage) from an itkHistogram + HistogramType::Pointer CreateHistogramFrom( + HistoPersistentFilterType::HistogramType::Pointer histo , + unsigned int nbBin , + HistogramType::SpacingType histoSpacing , + HistogramType::PointType histoOrigin ) + { + HistogramType::SizeType sizeOne; + sizeOne.Fill( 1 ); + HistogramType::IndexType index; + index.Fill(0); + + HistogramType::PixelType histoPixel( nbBin ); + + HistogramType::Pointer histoVectorImage ( + HistogramType::New() ); + + histoVectorImage->SetVectorLength( nbBin ); + histoVectorImage->SetBufferedRegion( sizeOne ); + histoVectorImage->SetRequestedRegion( sizeOne ); + histoVectorImage->SetLargestPossibleRegion( sizeOne ); + histoVectorImage->SetSignedSpacing( histoSpacing ) ; + histoVectorImage->SetOrigin( histoOrigin ); + histoVectorImage->Allocate(); + for (unsigned int j = 0 ; j < nbBin ; j++ ) + { + histoPixel[j] = histo->GetFrequency( j ); + } + histoVectorImage->SetPixel( index , histoPixel ); + return histoVectorImage; + } + + // Set correct parameters for the ComputeHistoFilter + void SetHistoFilterParameter( HistoFilterType::Pointer histoFilter , + float min , + float max , + unsigned int nbBin , + float thresh = -1 ) + { + histoFilter->SetMin( min ); + histoFilter->SetMax( max ); + histoFilter->SetNbBin( nbBin ); + histoFilter->SetThumbSize( m_ThumbSize ); + histoFilter->SetThreshold( thresh ); + if ( IsParameterEnabled("nodata") ) + { + histoFilter->SetNoData( GetParameterFloat( "nodata" ) ); + histoFilter->SetNoDataFlag( true ); + } + } + + // Set correct parameters for the ComputeGainLutFilter + void SetGainLutFilterParameter( GainLutFilterType::Pointer gainLutFilter , + ImagePixelType min , + ImagePixelType max ) + { + gainLutFilter->SetMin( min ); + gainLutFilter->SetMax( max ); + gainLutFilter->SetNbPixel( m_ThumbSize[0]*m_ThumbSize[1] ); + } + + // Set correct parameters for the ApplyGainFilter + void SetApplyFilterParameter( ApplyFilterType::Pointer applyFilter , + ImagePixelType min , + ImagePixelType max ) + { + applyFilter->SetMin( min ); + applyFilter->SetMax( max ); + applyFilter->SetThumbSize( m_ThumbSize ); + if ( IsParameterEnabled("nodata") ) + { + applyFilter->SetNoData( GetParameterFloat( "nodata" ) ); + applyFilter->SetNoDataFlag( true ); + } + } + + std::string m_SpatialMode , m_MinMaxMode , m_EqMode ; + FloatImageType::SizeType m_ThumbSize; + ImageListToVectorFilterType::Pointer m_ImageListToVectorFilterOut; + LuminanceFunctorType::Pointer m_LuminanceFunctor; + VectorToImageListFilterType::Pointer m_LuminanceToImageListFilter; + VectorToImageListFilterType::Pointer m_VectorToImageListFilter; + std::vector < GainLutFilterType::Pointer > m_GainLutFilter; + std::vector < HistoFilterType::Pointer > m_HistoFilter; + std::vector < HistogramType::Pointer > m_Histogram; + std::vector < ApplyFilterType::Pointer > m_ApplyFilter; + std::vector < StreamingImageFilterType::Pointer > m_StreamingFilter; + std::vector < BufferFilterType::Pointer > m_BufferFilter; + +}; + + + +} //End namespace Wrapper +} //End namespace otb + +OTB_APPLICATION_EXPORT(otb::Wrapper::ContrastEnhancement) diff --git a/Modules/Applications/AppFiltering/otb-module.cmake b/Modules/Applications/AppFiltering/otb-module.cmake index 5e8fc45196bb23ef0d7e19e6c062b7d7a5d3a5a2..e5c818748dcf89c3083847334d93d96819405c75 100644 --- a/Modules/Applications/AppFiltering/otb-module.cmake +++ b/Modules/Applications/AppFiltering/otb-module.cmake @@ -26,6 +26,9 @@ otb_module(OTBAppFiltering OTBITK OTBApplicationEngine OTBImageBase + OTBContrast + OTBStatistics + OTBStreaming TEST_DEPENDS OTBTestKernel diff --git a/Modules/Applications/AppFiltering/test/CMakeLists.txt b/Modules/Applications/AppFiltering/test/CMakeLists.txt index 9c4a8026882033f886e31cc4d222c21945eb3f1c..f8195d22ce294a43281177f623e47c534de4927a 100644 --- a/Modules/Applications/AppFiltering/test/CMakeLists.txt +++ b/Modules/Applications/AppFiltering/test/CMakeLists.txt @@ -48,6 +48,61 @@ otb_test_application(NAME apTvUtSmoothingTest_OutXML VALID --compare-image ${NOTOL} ${BASELINE}/apTvUtSmoothingTest.tif ${TEMP}/apTvUtSmoothingTest_OutXML.tif) + +#----------- Contrast TESTS ---------------- + +otb_test_application(NAME apTvUtContrastTest_base + APP ContrastEnhancement + OPTIONS -in ${INPUTDATA}/QB_Suburb.png + -out ${TEMP}/apTvUtContrastTest_base.tif int16 + -bins 256 + -spatial.local.h 51 + -spatial.local.w 67 + VALID --compare-image ${NOTOL} + ${BASELINE}/apTvUtContrastTest_base.tif + ${TEMP}/apTvUtContrastTest_base.tif) + +otb_test_application(NAME apTvUtContrastTest_base_glob + APP ContrastEnhancement + OPTIONS -in ${INPUTDATA}/QB_Suburb.png + -out ${TEMP}/apTvUtContrastTest_base_glob.tif int16 + -bins 256 + -spatial global + -minmax manuel + -minmax.manuel.min 0 + -minmax.manuel.max 255 + VALID --compare-image ${NOTOL} + ${BASELINE}/apTvUtContrastTest_base_glob.tif + ${TEMP}/apTvUtContrastTest_base_glob.tif) + +otb_test_application(NAME apTvUtContrastTest_lum_glob + APP ContrastEnhancement + OPTIONS -in ${INPUTDATA}/anaglyphInput1.tif + -out ${TEMP}/apTvUtContrastTest_lum_glob.tif int16 + -bins 256 + -spatial global + -hfact 2.7 + -nodata 0 + -mode lum + VALID --compare-image ${NOTOL} + ${BASELINE}/apTvUtContrastTest_lum_glob.tif + ${TEMP}/apTvUtContrastTest_lum_glob.tif) + +otb_test_application(NAME apTvUtContrastTest_lum + APP ContrastEnhancement + OPTIONS -in ${INPUTDATA}/anaglyphInput1.tif + -out ${TEMP}/apTvUtContrastTest_lum.tif int16 + -bins 256 + -spatial.local.h 33 + -spatial.local.w 47 + -hfact 2.1 + -nodata 0 + -mode lum + VALID --compare-image ${NOTOL} + ${BASELINE}/apTvUtContrastTest_lum.tif + ${TEMP}/apTvUtContrastTest_lum.tif) + + diff --git a/Modules/Filtering/Contrast/CMakeLists.txt b/Modules/Filtering/Contrast/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..45a89f243f34d7d41633b84f7eb0132e6eb8d69c --- /dev/null +++ b/Modules/Filtering/Contrast/CMakeLists.txt @@ -0,0 +1,22 @@ +# +# Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) +# +# This file is part of Orfeo Toolbox +# +# https://www.orfeo-toolbox.org/ +# +# 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 +# +# 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. +# + +project(OTBContrast) +otb_module_impl() \ No newline at end of file diff --git a/Modules/Filtering/Contrast/include/otbApplyGainFilter.h b/Modules/Filtering/Contrast/include/otbApplyGainFilter.h new file mode 100644 index 0000000000000000000000000000000000000000..c735725382971d33b31273edcffcda0316184db0 --- /dev/null +++ b/Modules/Filtering/Contrast/include/otbApplyGainFilter.h @@ -0,0 +1,145 @@ +/* + * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) + * + * This file is part of Orfeo Toolbox + * + * https://www.orfeo-toolbox.org/ + * + * 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 + * + * 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 otbApplyGainFilter_h +#define otbApplyGainFilter_h + +#include "itkImageToImageFilter.h" +#include "otbImage.h" + +namespace otb +{ + +/** \class ApplyGainFilter + * \brief Apply gain on the input image with a bilineare interpolation + * + * This class implements the third part of the CLAHE algorithm. It's aim + * is to apply the computed gain with a bilineare interpolation. The gain + * is in a look up table, and the minimum and maximum asked by the filter + * should be the same as the one used to compute those look up table. + * + * \ingroup OTBContrast + */ + +template < class TInputImage , class TLut , class TOutputImage > +class ITK_EXPORT ApplyGainFilter : + public itk::ImageToImageFilter< TInputImage , TOutputImage > +{ +public : + /** typedef for standard classes. */ + + typedef TInputImage InputImageType; + typedef TOutputImage OutputImageType; + + typedef ApplyGainFilter Self; + typedef itk::ImageToImageFilter< InputImageType, OutputImageType > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef TLut LutType; + typedef typename InputImageType::InternalPixelType InputPixelType; + typedef typename OutputImageType::InternalPixelType OutputPixelType; + typedef typename OutputImageType::RegionType OutputImageRegionType; + + /** Method for creation through the object factory. */ + itkNewMacro(Self) + /** Run-time type information (and related methods). */ + itkTypeMacro(ComputeHistoFilter, ImageToImageFilter) + + /** Get/Set macro to get/set the nodata value */ + itkSetMacro(NoData, InputPixelType) + itkGetMacro(NoData, InputPixelType) + + /** Get/Set macro to get/set the nodata flag value */ + itkBooleanMacro(NoDataFlag) + itkGetMacro(NoDataFlag, bool) + itkSetMacro(NoDataFlag, bool) + + /** Get/Set macro to get/set the ThumbSizeFromSpacing flag value */ + itkBooleanMacro(ThumbSizeFromSpacing) + itkGetMacro(ThumbSizeFromSpacing, bool) + itkSetMacro(ThumbSizeFromSpacing, bool) + + /** Get/Set macro to get/set the thumbnail's size */ + itkSetMacro(ThumbSize, typename InputImageType::SizeType) + itkGetMacro(ThumbSize, typename InputImageType::SizeType) + + /** Get/Set macro to get/set the minimum value */ + itkSetMacro(Min, InputPixelType) + itkGetMacro(Min, InputPixelType) + + /** Get/Set macro to get/set the maximum value */ + itkSetMacro(Max, InputPixelType) + itkGetMacro(Max, InputPixelType) + + /** Set the input look up table*/ + void SetInputLut( const LutType * lut) ; + + /** Set the input image*/ + void SetInputImage( const InputImageType * input) ; + +protected : + ApplyGainFilter(); + ~ApplyGainFilter() override {} + void PrintSelf(std::ostream& os, itk::Indent indent) const override; + + /** Get the input image*/ + const InputImageType * GetInputImage() const; + + /** Get the input look up table*/ + const LutType * GetInputLut() const; + + void GenerateInputRequestedRegion() override; + + void BeforeThreadedGenerateData() override; + + void ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, + itk::ThreadIdType threadId) override; + void VerifyInputInformation() override {} ; + +private : + ApplyGainFilter(const Self &) = delete ; + void operator =(const Self&) = delete ; + + /** Bilinear interpolation of the gain between the different window.*/ + double InterpolateGain( typename LutType::ConstPointer gridLut , + unsigned int pixelValue , + typename InputImageType::IndexType index); + + InputPixelType m_NoData; + InputPixelType m_Min; + InputPixelType m_Max; + bool m_NoDataFlag; + bool m_ThumbSizeFromSpacing; + double m_Step; + typename LutType::SizeType m_LutSize; + typename InputImageType::SizeType m_ThumbSize; + + +}; + +} // End namespace otb + +#ifndef OTB_MANUAL_INSTANTIATION +#include "otbApplyGainFilter.txx" +#endif + + +#endif diff --git a/Modules/Filtering/Contrast/include/otbApplyGainFilter.txx b/Modules/Filtering/Contrast/include/otbApplyGainFilter.txx new file mode 100644 index 0000000000000000000000000000000000000000..28bd851c0dbc83959e497147b413fb4cc8ad224b --- /dev/null +++ b/Modules/Filtering/Contrast/include/otbApplyGainFilter.txx @@ -0,0 +1,224 @@ +/* + * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) + * + * This file is part of Orfeo Toolbox + * + * https://www.orfeo-toolbox.org/ + * + * 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 + * + * 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 otbApplyGainFilter_txx +#define otbApplyGainFilter_txx + +#include "otbApplyGainFilter.h" +#include "itkImageRegionIterator.h" +#include "itkImageRegionConstIteratorWithIndex.h" +#include "itkContinuousIndex.h" + +#include <limits> + +namespace otb +{ +template <class TInputImage , class TLut , class TOutputImage > +ApplyGainFilter < TInputImage , TLut , TOutputImage > +::ApplyGainFilter() +{ + this->SetNumberOfRequiredInputs(2); + m_Min = std::numeric_limits< InputPixelType >::quiet_NaN(); + m_Max = std::numeric_limits< InputPixelType >::quiet_NaN(); + m_NoData = std::numeric_limits< InputPixelType >::quiet_NaN(); + m_NoDataFlag = false; + m_ThumbSizeFromSpacing = true; + m_Step = -1; +} + +template <class TInputImage , class TLut , class TOutputImage > +void ApplyGainFilter < TInputImage , TLut , TOutputImage > +::SetInputImage(const InputImageType * input) +{ + // Process object is not const-correct so the const casting is required. + this->SetNthInput( 0 , const_cast<InputImageType *>( input ) ); +} + +template <class TInputImage , class TLut , class TOutputImage > +const TInputImage * ApplyGainFilter < TInputImage , TLut , TOutputImage > +::GetInputImage() const +{ + return static_cast< const InputImageType * > + ( this->itk::ProcessObject::GetInput(0) ); +} + +template <class TInputImage , class TLut , class TOutputImage > +void ApplyGainFilter < TInputImage , TLut , TOutputImage > +::SetInputLut(const LutType * lut) +{ + // Process object is not const-correct so the const casting is required. + this->SetNthInput( 1 , const_cast<LutType *>( lut ) ); +} + +template <class TInputImage , class TLut , class TOutputImage > +const TLut * ApplyGainFilter < TInputImage , TLut , TOutputImage > +::GetInputLut() const +{ + return static_cast< const LutType * > + ( this->itk::ProcessObject::GetInput(1) ); +} + +template <class TInputImage , class TLut , class TOutputImage > +void ApplyGainFilter < TInputImage , TLut , TOutputImage > +::GenerateInputRequestedRegion() +{ + typename InputImageType::Pointer input ( + const_cast<InputImageType *>( GetInputImage() ) ); + typename LutType::Pointer lut ( const_cast<LutType *>( GetInputLut() ) ); + typename OutputImageType::Pointer output ( this->GetOutput() ); + + lut->SetRequestedRegion( lut->GetLargestPossibleRegion() ); + input->SetRequestedRegion( output->GetRequestedRegion() ); + if ( input->GetRequestedRegion().GetNumberOfPixels() == 0 ) + { + input->SetRequestedRegionToLargestPossibleRegion(); + } +} + +template <class TInputImage , class TLut , class TOutputImage > +void ApplyGainFilter < TInputImage , TLut , TOutputImage > +::BeforeThreadedGenerateData() +{ + typename LutType::ConstPointer lut ( GetInputLut() ); + typename InputImageType::ConstPointer input ( GetInputImage() ); + if ( m_ThumbSizeFromSpacing ) + { + m_ThumbSize[0] = std::round( lut->GetSignedSpacing()[0] + / input->GetSignedSpacing()[0] ); + m_ThumbSize[1] = std::round( lut->GetSignedSpacing()[1] + / input->GetSignedSpacing()[1] ); + } + m_Step = static_cast<double>( m_Max - m_Min ) \ + / static_cast<double>( lut->GetVectorLength() - 1 ); +} + +template <class TInputImage , class TLut , class TOutputImage > +void ApplyGainFilter < TInputImage , TLut , TOutputImage > +::ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread , + itk::ThreadIdType itkNotUsed(threadId) ) +{ + assert( m_Step > 0 ); + // TODO error + // support progress methods/callbacks + // itk::ProgressReporter progress(this , threadId , + // outputRegionForThread.GetNumberOfPixels() ); + + typename InputImageType::ConstPointer input ( GetInputImage() ); + typename LutType::ConstPointer lut ( GetInputLut() ); + typename OutputImageType::Pointer output ( this->GetOutput() ); + typename InputImageType::RegionType inputRegionForThread( + outputRegionForThread ); + + + itk::ImageRegionConstIteratorWithIndex < InputImageType > it ( input , + inputRegionForThread ); + itk::ImageRegionIterator <OutputImageType > oit ( output , + outputRegionForThread ); + + unsigned int pixelLutValue(0); + double gain(0.0) , newValue(0); + InputPixelType currentPixel(0); + + for( it.GoToBegin() , oit.GoToBegin() ; !oit.IsAtEnd() || !it.IsAtEnd() ; + ++oit , ++it ) + { + currentPixel = it.Get(); + newValue = static_cast< double > ( currentPixel ); + if( !( ( currentPixel == m_NoData && m_NoDataFlag ) || + currentPixel > m_Max || currentPixel < m_Min ) ) + { + pixelLutValue = static_cast< unsigned int > ( + std::round( ( currentPixel - m_Min ) / m_Step ) ); + gain = InterpolateGain( lut , pixelLutValue , it.GetIndex() ); + newValue *= gain; + } + oit.Set( static_cast<OutputPixelType>( newValue ) ); + } + assert ( oit.IsAtEnd() && it.IsAtEnd() ); +} + +template <class TInputImage , class TLut , class TOutputImage > +double ApplyGainFilter < TInputImage , TLut , TOutputImage > +::InterpolateGain( typename LutType::ConstPointer gridLut, + unsigned int pixelLutValue , + typename InputImageType::IndexType index) +{ + typename InputImageType::PointType pixelPoint; + typename itk::ContinuousIndex< double , 2 > pixelIndex; + typename InputImageType::ConstPointer input ( GetInputImage() ); + typename LutType::ConstPointer lut ( GetInputLut() ); + input->TransformIndexToPhysicalPoint( index , pixelPoint ); + lut->TransformPhysicalPointToContinuousIndex( pixelPoint , pixelIndex ); + std::vector< typename LutType::IndexType > neighbors(4); + neighbors[0][0] = std::floor(pixelIndex[0]) ; + neighbors[0][1] = std::floor(pixelIndex[1]) ; + neighbors[1][0] = neighbors[0][0] + 1 ; + neighbors[1][1] = neighbors[0][1] ; + neighbors[2][0] = neighbors[0][0] ; + neighbors[2][1] = neighbors[0][1] + 1 ; + neighbors[3][0] = neighbors[0][0] + 1 ; + neighbors[3][1] = neighbors[0][1] + 1 ; + float gain(0.f) , w(0.f) , wtm(0.f); + typename LutType::IndexType maxIndex; + maxIndex[0] = lut->GetLargestPossibleRegion().GetSize()[0]; + maxIndex[1] = lut->GetLargestPossibleRegion().GetSize()[1]; + for ( auto i : neighbors ) + { + if ( i[0] < 0 || i[1] < 0 || i[0] >= maxIndex[0] || i[1] >= maxIndex[1] ) + continue; + if ( gridLut->GetPixel(i)[pixelLutValue] == -1 ) + continue; + wtm = ( 1 - std::abs( pixelIndex[0] - i[0] ) ) + * ( 1 - std::abs( pixelIndex[1] - i[1] ) ); + gain += gridLut->GetPixel(i)[pixelLutValue] * wtm; + w += wtm; + } + if ( w == 0 ) + { + w = 1; + gain = 1; + } + + return gain/w; +} + +/** + * Standard "PrintSelf" method + */ +template <class TInputImage , class TLut , class TOutputImage > +void ApplyGainFilter < TInputImage , TLut , TOutputImage > +::PrintSelf(std::ostream& os, itk::Indent indent) const +{ + Superclass::PrintSelf(os, indent); + os << indent << "Is no data activated : " << m_NoDataFlag << std::endl; + os << indent << "No Data : " << m_NoData << std::endl; + os << indent << "Minimum : " << m_Min << std::endl; + os << indent << "Maximum : " << m_Max << std::endl; + os << indent << "Step : " << m_Step << std::endl; + os << indent << "Look up table size : " << m_LutSize << std::endl; + os << indent << "Is ThumbSize from sapcing is activated : " + << m_NoDataFlag << std::endl; + os << indent << "Thumbnail size : " << m_ThumbSize << std::endl; +} + + +} // End namespace otb + +#endif diff --git a/Modules/Filtering/Contrast/include/otbCLHistogramEqualizationFilter.h b/Modules/Filtering/Contrast/include/otbCLHistogramEqualizationFilter.h new file mode 100644 index 0000000000000000000000000000000000000000..edd5167a9b901f069f0f63658f12255d8037a697 --- /dev/null +++ b/Modules/Filtering/Contrast/include/otbCLHistogramEqualizationFilter.h @@ -0,0 +1,186 @@ +/* + * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) + * + * This file is part of Orfeo Toolbox + * + * https://www.orfeo-toolbox.org/ + * + * 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 + * + * 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 otbCLHistogramEqualizationFilter_h +#define otbCLHistogramEqualizationFilter_h + +#include "itkImageToImageFilter.h" +#include "otbImage.h" +#include "otbVectorImage.h" +#include "otbComputeHistoFilter.h" +#include "otbComputeGainLutFilter.h" +#include "otbApplyGainFilter.h" +#include "itkStreamingImageFilter.h" +#include "otbInPlacePassFilter.h" + +namespace otb +{ + +/** \class CLHistogramEqualizationFilter + * \brief Implement CLAHE algorithm + * + * This class implement CLAHE algorithm. It is a composite filter that gathers + * the 3 filters (ComputeHisto, ComputeGainLut, ApplyGain) and pipes them with + * additional filters (InPlacePass and StreamingImage) in order to make + * streaming available. + * \ingroup OTBContrast + */ + +template < class TInputImage , class TOutputImage > +class ITK_EXPORT CLHistogramEqualizationFilter : + public itk::ImageToImageFilter< TInputImage , TOutputImage > +{ +public : + /** typedef for standard classes. */ + + typedef TInputImage InputImageType; + typedef TOutputImage OutputImageType; + + typedef CLHistogramEqualizationFilter Self; + typedef itk::ImageToImageFilter< InputImageType, OutputImageType > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef otb::VectorImage< unsigned int , 2 > HistogramType; + typedef otb::VectorImage< double , 2 > LutType; + + typedef itk::StreamingImageFilter< LutType , LutType > + StreamingImageFilter; + + typedef otb::InPlacePassFilter < InputImageType > BufferFilter; + + typedef otb::ComputeHistoFilter< InputImageType , HistogramType > + HistoFilter; + + typedef otb::ComputeGainLutFilter< HistogramType , LutType > + GainLutFilter; + + typedef otb::ApplyGainFilter< InputImageType , LutType , OutputImageType > + ApplyGainFilter; + + typedef typename InputImageType::PixelType InputPixelType; + typedef typename OutputImageType::RegionType OutputImageRegionType; + /** Method for creation through the object factory. */ + itkNewMacro(Self) + /** Run-time type information (and related methods). */ + itkTypeMacro(CLHistogramEqualizationFilter, ImageToImageFilter) + + itkGetMacro(Min , InputPixelType) + // itkSetMacro(Min, InputPixelType) + void SetMin( InputPixelType min ) + { + m_HistoFilter->SetMin(min); + m_GainLutFilter->SetMin(min); + m_ApplyGainFilter->SetMin(min); + m_Min = min; + }; + + itkGetMacro(Max , InputPixelType) + // itkSetMacro(Max, InputPixelType) + void SetMax( InputPixelType max ) + { + m_HistoFilter->SetMax(max); + m_GainLutFilter->SetMax(max); + m_ApplyGainFilter->SetMax(max); + m_Max = max; + }; + + itkGetMacro(NbBin , unsigned long) + // itkSetMacro(NbBin, unsigned long) + void SetNbBin( unsigned long bin ) + { + m_HistoFilter->SetNbBin(bin); + m_NbBin = bin; + }; + + itkGetMacro(ThumbSize , typename InputImageType::SizeType) + // itkSetMacro(ThumbSize, typename InputImageType::SizeType) + void SetThumbSize( typename InputImageType::SizeType size ) + { + m_HistoFilter->SetThumbSize(size); + m_ApplyGainFilter->SetThumbSize(size); + m_GainLutFilter->SetNbPixel(size[0]*size[1]); + m_ThumbSize = size; + }; + + itkGetMacro(Threshold , double) + // itkSetMacro(Threshold, double) + void SetThreshold( double t ) + { + m_HistoFilter->SetThreshold(t); + m_Threshold = t; + }; + + itkGetMacro(NoData , InputPixelType) + // itkSetMacro(NoData, InputPixelType) + void SetNoData( InputPixelType n ) + { + m_HistoFilter->SetNoData(n); + m_ApplyGainFilter->SetNoData(n); + m_NoData = n; + } + + + itkGetMacro(NoDataFlag , bool) + // itkSetMacro(NoDataFlag, bool) + void SetNoDataFlag( bool flag ) + { + m_HistoFilter->SetNoDataFlag(flag); + m_ApplyGainFilter->SetNoDataFlag(flag); + m_NoDataFlag = flag; + } + +protected : + CLHistogramEqualizationFilter(); + ~CLHistogramEqualizationFilter() override {} + + void PrintSelf(std::ostream& os, itk::Indent indent) const override ; + + void UpdateOutputInformation() override; + + void PropagateRequestedRegion( itk::DataObject * output ) override; + + void GenerateData() override; + +private : + CLHistogramEqualizationFilter(const Self &) = delete ; + void operator =(const Self&) = delete ; + + typename HistoFilter::Pointer m_HistoFilter; + typename GainLutFilter::Pointer m_GainLutFilter; + typename ApplyGainFilter::Pointer m_ApplyGainFilter; + typename StreamingImageFilter::Pointer m_StreamingImageFilter; + typename BufferFilter::Pointer m_BufferFilter; + InputPixelType m_Min , m_Max , m_NoData; + unsigned long m_NbBin; + typename InputImageType::SizeType m_ThumbSize; + double m_Threshold; + bool m_NoDataFlag; + +}; + +} // End namespace otb + +#ifndef OTB_MANUAL_INSTANTIATION +#include "otbCLHistogramEqualizationFilter.txx" +#endif + + +#endif diff --git a/Modules/Filtering/Contrast/include/otbCLHistogramEqualizationFilter.txx b/Modules/Filtering/Contrast/include/otbCLHistogramEqualizationFilter.txx new file mode 100644 index 0000000000000000000000000000000000000000..f5cbebd96ba57b459509f3434e8d26ff61497f5b --- /dev/null +++ b/Modules/Filtering/Contrast/include/otbCLHistogramEqualizationFilter.txx @@ -0,0 +1,102 @@ +/* + * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) + * + * This file is part of Orfeo Toolbox + * + * https://www.orfeo-toolbox.org/ + * + * 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 + * + * 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 otbCLHistogramEqualizationFilter_txx +#define otbCLHistogramEqualizationFilter_txx + +#include "otbCLHistogramEqualizationFilter.h" +#include "itkImageRegionIterator.h" + +#include <limits> + +namespace otb +{ +template < class TInputImage , class TOutputImage > +CLHistogramEqualizationFilter < TInputImage , TOutputImage > +::CLHistogramEqualizationFilter(): +m_HistoFilter( HistoFilter::New() ) , +m_GainLutFilter ( GainLutFilter::New() ) , +m_ApplyGainFilter ( ApplyGainFilter::New() ) , +m_StreamingImageFilter ( StreamingImageFilter::New() ) , +m_BufferFilter ( BufferFilter::New() ) +{ + m_Min = std::numeric_limits< InputPixelType >::quiet_NaN(); + m_Max = std::numeric_limits< InputPixelType >::quiet_NaN(); + m_NbBin = 256; + m_Threshold = -1; + m_NoDataFlag = false; + m_NoData = std::numeric_limits< InputPixelType >::quiet_NaN(); + m_ThumbSize.Fill(0); + m_GainLutFilter->SetInput( m_HistoFilter->GetHistoOutput() ); + m_StreamingImageFilter->SetInput( m_GainLutFilter->GetOutput() ); + m_ApplyGainFilter->SetInputLut( m_StreamingImageFilter->GetOutput() ); + m_ApplyGainFilter->SetInputImage( m_BufferFilter->GetOutput() ); +} + +template < class TInputImage , class TOutputImage > +void CLHistogramEqualizationFilter < TInputImage , TOutputImage > +::UpdateOutputInformation() +{ + const InputImageType * input = this->GetInput() ; + m_HistoFilter->SetInput( input ); + m_BufferFilter->SetInput( input ); + m_ApplyGainFilter->GetOutput()->UpdateOutputInformation(); + this->GetOutput()->CopyInformation( m_ApplyGainFilter->GetOutput() ); +} + +template < class TInputImage , class TOutputImage > +void CLHistogramEqualizationFilter < TInputImage , TOutputImage > +::PropagateRequestedRegion( itk::DataObject * output ) +{ + m_ApplyGainFilter->GetOutput()->SetRequestedRegion( static_cast<OutputImageType *>(output)->GetRequestedRegion() ); + m_ApplyGainFilter->GetOutput()->PropagateRequestedRegion(); +} + +template < class TInputImage , class TOutputImage > +void CLHistogramEqualizationFilter < TInputImage , TOutputImage > +::GenerateData() +{ + m_ApplyGainFilter->GraftOutput( this->GetOutput() ); + m_ApplyGainFilter->Update(); + this->GraftOutput( m_ApplyGainFilter->GetOutput() ); + +} + +/** + * Standard "PrintSelf" method + */ +template <class TInputImage , class TOutputImage > +void CLHistogramEqualizationFilter < TInputImage , TOutputImage > +::PrintSelf(std::ostream& os, itk::Indent indent) const +{ + Superclass::PrintSelf(os, indent); + os << indent << "Minimum : " << m_Min << std::endl; + os << indent << "Maximum : " << m_Max << std::endl; + os << indent << "Bin Number : " << m_NbBin << std::endl; + os << indent << "Thumbnail size : " << m_ThumbSize << std::endl; + os << indent << "Threshold value : " << m_Threshold << std::endl; + os << indent << "Is no data activated : " << m_NoDataFlag << std::endl; + os << indent << "No Data : " << m_NoData << std::endl; +} + + +} // End namespace otb + +#endif diff --git a/Modules/Filtering/Contrast/include/otbComputeGainLutFilter.h b/Modules/Filtering/Contrast/include/otbComputeGainLutFilter.h new file mode 100644 index 0000000000000000000000000000000000000000..431b6f1f0f5f55d0bc5c403a1df86e8f5f26e255 --- /dev/null +++ b/Modules/Filtering/Contrast/include/otbComputeGainLutFilter.h @@ -0,0 +1,128 @@ +/* + * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) + * + * This file is part of Orfeo Toolbox + * + * https://www.orfeo-toolbox.org/ + * + * 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 + * + * 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 otbComputeGainLutFilter_h +#define otbComputeGainLutFilter_h + +#include "itkImageToImageFilter.h" +#include "otbImage.h" + +namespace otb +{ + +/** \class ComputeGainLutFilter + * \brief Compute the gain for each pixel value from a histogram + * + * This class implements the second part of the CLAHE algorithm. It's aim + * is to compute a look up table filled with gain that need to be applied + * on the input to match the target histogram. To keep consistency with + * the other parts of the algorithm it needs the minimum and maximum value + * of the input image and also the theorical number of pixel per histogram. + * + * \ingroup OTBContrast + */ + +template <class TInputImage , class TOutputImage > +class ITK_EXPORT ComputeGainLutFilter : + public itk::ImageToImageFilter< TInputImage , TOutputImage > +{ +public: + + /** typedef for standard classes. */ + + typedef TInputImage InputImageType; + typedef TOutputImage OutputImageType; + + typedef ComputeGainLutFilter Self; + typedef itk::ImageToImageFilter< InputImageType, OutputImageType > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef typename InputImageType::PixelType HistoType; + typedef typename OutputImageType::PixelType LutType; + typedef typename OutputImageType::InternalPixelType OutputPixelType; + typedef typename OutputImageType::RegionType OutputImageRegionType; + + /** Method for creation through the object factory. */ + itkNewMacro(Self) + /** Run-time type information (and related methods). */ + itkTypeMacro(ComputeGainLutFilter, ImageToImageFilter) + + /** Get/Set macro to get/set the number of pixel by histogram */ + itkSetMacro(NbPixel, unsigned long) + itkGetMacro(NbPixel, unsigned long) + + /** Get/Set macro to get/set the minimum value */ + itkSetMacro(Min, double) + itkGetMacro(Min, double) + + /** Get/Set macro to get/set the maximum value */ + itkSetMacro(Max, double) + itkGetMacro(Max, double) + +protected: + ComputeGainLutFilter() ; + ~ComputeGainLutFilter() override {} + void PrintSelf(std::ostream& os, itk::Indent indent) const override; + + void BeforeThreadedGenerateData() override ; + + void ThreadedGenerateData( + const OutputImageRegionType & outputRegionForThread , + itk::ThreadIdType threadId) override ; + +private: + ComputeGainLutFilter(const Self &) = delete ; + void operator =(const Self&) = delete ; + + /** Post-process the look up tabe to get a gain instead of a simple value */ + OutputPixelType PostProcess( unsigned int countMapValue , + unsigned int countValue ) ; + + /** Equalized input histogram regarding the target and filling the + * corresponding look up table */ + void Equalized( const HistoType & inputHisto , + HistoType & targetHisto , + LutType & lut) ; + + /** Create target depending on the number of pixel in the input histogram */ + void CreateTarget( const HistoType & inputHisto , + HistoType & targetHisto ) ; + //TODO Give the opportunity to choose the histogram target + + /** Check whether the input histogram has enought pixel to be meaningful */ + bool IsValid(const HistoType & inputHisto ) ; + + double m_Min; + double m_Max; + double m_Step; + unsigned int m_NbBin; + unsigned long m_NbPixel; + +}; + + // End namespace otb +} + +#ifndef OTB_MANUAL_INSTANTIATION +#include "otbComputeGainLutFilter.txx" +#endif + +#endif diff --git a/Modules/Filtering/Contrast/include/otbComputeGainLutFilter.txx b/Modules/Filtering/Contrast/include/otbComputeGainLutFilter.txx new file mode 100644 index 0000000000000000000000000000000000000000..ccf1b857528fb6cd5d85512545a9eacad24fc895 --- /dev/null +++ b/Modules/Filtering/Contrast/include/otbComputeGainLutFilter.txx @@ -0,0 +1,191 @@ +/* + * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) + * + * This file is part of Orfeo Toolbox + * + * https://www.orfeo-toolbox.org/ + * + * 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 + * + * 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 otbComputeGainLutFilter_txx +#define otbComputeGainLutFilter_txx + +#include "otbComputeGainLutFilter.h" +#include "itkImageRegionIterator.h" + +#include <limits> +#include <numeric> + +namespace otb +{ + +template < class TInputImage, class TOutputImage > +ComputeGainLutFilter < TInputImage , TOutputImage > +::ComputeGainLutFilter() +{ + m_NbBin = 256; + m_NbPixel = 0; + m_Min = std::numeric_limits< double >::quiet_NaN(); + m_Max = std::numeric_limits< double >::quiet_NaN(); + m_Step = -1; +} + +template <class TInputImage , class TOutputImage > +void ComputeGainLutFilter <TInputImage , TOutputImage > +::BeforeThreadedGenerateData() +{ + m_NbBin = this->GetInput()->GetNumberOfComponentsPerPixel(); + m_Step = static_cast<double>( m_Max - m_Min ) \ + / static_cast<double>( m_NbBin -1 ); +} + +template <class TInputImage , class TOutputImage > +void ComputeGainLutFilter <TInputImage , TOutputImage > +::ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread, + itk::ThreadIdType itkNotUsed(threadId) ) +{ + assert( m_Step > 0 ); + assert( m_NbBin > 0 ); + // TODO error + // itk::ProgressReporter progress(this , threadId , + // outputRegionForThread.GetNumberOfPixels() ); + + typename InputImageType::ConstPointer input ( this->GetInput() ); + typename OutputImageType::Pointer output ( this->GetOutput() ); + + itk::ImageRegionConstIterator < InputImageType > it ( input , + outputRegionForThread ); + + itk::ImageRegionIterator <OutputImageType > oit ( output , + outputRegionForThread ); + + HistoType target; + target.SetSize( m_NbBin ); + + HistoType currentHisto; + + LutType lut; + lut.SetSize( m_NbBin ); + + for (it.GoToBegin() , oit.GoToBegin() ; !oit.IsAtEnd() || !it.IsAtEnd() ; + ++oit , ++it ) + { + currentHisto = it.Get(); + target.Fill(0); + lut.Fill(-1); + if ( IsValid( currentHisto ) ) + { + CreateTarget( currentHisto , target ); + Equalized( currentHisto , target , lut ); + } + oit.Set( lut ); + } + assert ( oit.IsAtEnd() && it.IsAtEnd() ); +} + +template <class TInputImage, class TOutputImage > +typename TOutputImage::InternalPixelType + ComputeGainLutFilter < TInputImage , TOutputImage > +::PostProcess( unsigned int countValue , + unsigned int countMapValue ) +{ + double denum ( countValue * m_Step + m_Min ); + if ( denum == 0 ) + return 0; + return static_cast< OutputPixelType > ( (countMapValue * m_Step + m_Min) + / denum ); +} + +template <class TInputImage, class TOutputImage > +void ComputeGainLutFilter < TInputImage , TOutputImage > +::Equalized( const HistoType & inputHisto , + HistoType & targetHisto , + LutType & lut) +{ + unsigned int countValue(0) , countMapValue(0) ; + lut[countValue] = 1; // Black stays black + ++countValue; + unsigned int countInput ( inputHisto[ 0 ] + inputHisto[ countValue ] ); + lut[m_NbBin - 1 ] = 1 ; // White stays white + unsigned int countTarget ( targetHisto[ countMapValue ] ); + + while ( ( countMapValue < m_NbBin ) && countValue < ( m_NbBin - 1 ) ) + { + if ( countInput > countTarget ) + { + ++countMapValue; + countTarget += targetHisto[ countMapValue ]; + } + else + { + lut[countValue] = PostProcess( countValue , countMapValue ); + ++countValue; + countInput += inputHisto[ countValue ]; + } + } + for (unsigned int i = 0 ; i < m_NbBin ; i++) + { + if (lut[i] == -1) + { + lut[i] = 1; + } + } +} + +template <class TInputImage, class TOutputImage > +void ComputeGainLutFilter < TInputImage , TOutputImage > +::CreateTarget( const HistoType & inputHisto , + HistoType & targetHisto ) +{ + unsigned int nbPixel(0); + for ( unsigned int i = 0 ; i < m_NbBin ; i++ ) + { + nbPixel += inputHisto[i]; + } + unsigned int rest ( nbPixel % m_NbBin ) , height ( nbPixel / m_NbBin ); + targetHisto.Fill(height); + for ( unsigned int i = 0 ; i < rest ; i++ ) + { + ++targetHisto[ ( m_NbBin - rest ) / 2 + i ]; + } +} + +template <class TInputImage, class TOutputImage > +bool ComputeGainLutFilter < TInputImage , TOutputImage > +::IsValid( const HistoType & inputHisto ) +{ + long acc = std::accumulate( &inputHisto[0] , + &inputHisto[ m_NbBin - 1 ] , + 0); + return acc >= (0.5 * m_NbPixel); +} + +/** + * Standard "PrintSelf" method + */ +template <class TInputImage , class TOutputImage > +void ComputeGainLutFilter < TInputImage , TOutputImage > +::PrintSelf(std::ostream& os, itk::Indent indent) const +{ + Superclass::PrintSelf(os, indent); + os << indent << "Minimum: " << m_Min << std::endl; + os << indent << "Maximum: " << m_Max << std::endl; + os << indent << "Step: " << m_Step << std::endl; + os << indent << "Number of bin: " << m_NbBin << std::endl; + os << indent << "Number of pixel by histogram: " << m_NbPixel << std::endl; +} + +} // End namespace otb + +#endif diff --git a/Modules/Filtering/Contrast/include/otbComputeHistoFilter.h b/Modules/Filtering/Contrast/include/otbComputeHistoFilter.h new file mode 100644 index 0000000000000000000000000000000000000000..2f54203b2458a43ebb39c8dcad41e94ae4569a03 --- /dev/null +++ b/Modules/Filtering/Contrast/include/otbComputeHistoFilter.h @@ -0,0 +1,160 @@ +/* + * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) + * + * This file is part of Orfeo Toolbox + * + * https://www.orfeo-toolbox.org/ + * + * 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 + * + * 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 otbComputeHistoFilter_h +#define otbComputeHistoFilter_h + +#include "itkImageToImageFilter.h" +#include "otbImage.h" +#include "itkImageRegionIterator.h" + +namespace otb +{ + +/** \class ComputeHistoFilter + * \brief Compute local histogram with several parameters + * + * This class implements the first part of the CLAHE algorithm. It aims + * to compute local histogram with several input parameters such as + * nodata value, threshold, thumbnail size and number of bin. Mandatory parameters are min + * and max value as it will be used in the histogram computation. + * + * \ingroup OTBContrast + */ + +template < class TInputImage , class TOutputImage > +class ITK_EXPORT ComputeHistoFilter : + public itk::ImageToImageFilter< TInputImage , TOutputImage > +{ +public: + /** typedef for standard classes. */ + + typedef TInputImage InputImageType; + typedef TOutputImage OutputImageType; + + typedef ComputeHistoFilter Self; + typedef itk::ImageToImageFilter< InputImageType, OutputImageType > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef typename InputImageType::InternalPixelType InputPixelType; + typedef typename InputImageType::IndexType IndexType; + typedef typename InputImageType::SizeType SizeType; + typedef typename OutputImageType::RegionType OutputImageRegionType; + + /** Method for creation through the object factory. */ + itkNewMacro(Self) + + /** Run-time type information (and related methods). */ + itkTypeMacro(ComputeHistoFilter, ImageToImageFilter) + + /** Get/Set macro to get/set the number of bin. Default value is 256 */ + itkSetMacro(NbBin, unsigned int) + itkGetMacro(NbBin, unsigned int) + + /** Get/Set macro to get/set the minimum value */ + itkSetMacro(Min, InputPixelType) + itkGetMacro(Min, InputPixelType) + + /** Get/Set macro to get/set the maximum value */ + itkSetMacro(Max, InputPixelType) + itkGetMacro(Max, InputPixelType) + + /** Get/Set macro to get/set the nodata value */ + itkSetMacro(NoData, InputPixelType) + itkGetMacro(NoData, InputPixelType) + + /** Get/Set macro to get/set the nodata flag value */ + itkBooleanMacro(NoDataFlag) + itkGetMacro(NoDataFlag, bool) + itkSetMacro(NoDataFlag, bool) + + /** Get/Set macro to get/set the thumbnail's size */ + itkSetMacro(ThumbSize, SizeType) + itkGetMacro(ThumbSize, SizeType) + + /** Get/Set macro to get/set the threshold parameter */ + itkSetMacro(Threshold , float) + itkGetMacro(Threshold , float) + + typename OutputImageType::Pointer GetHistoOutput(); + + virtual itk::ProcessObject::DataObjectPointer + MakeOutput(itk::ProcessObject::DataObjectPointerArraySizeType idx) + override; + + virtual itk::ProcessObject::DataObjectPointer + MakeOutput(const itk::ProcessObject::DataObjectIdentifierType &) + override; + + +protected: + ComputeHistoFilter(); + ~ComputeHistoFilter() override {} + void PrintSelf(std::ostream& os, itk::Indent indent) const override; + + void GenerateInputRequestedRegion() override; + + void GenerateOutputInformation() override; + + // Call BeforeThreadedGenerateData after getting the number of thread + void GenerateData() override; + + + void BeforeThreadedGenerateData() override; + + virtual void ThreadedGenerateData( + const OutputImageRegionType & outputRegionForThread , + itk::ThreadIdType threadId) override; + + void AfterThreadedGenerateData() override; + + void GenerateOutputRequestedRegion( itk::DataObject *output ) override; + +private: + ComputeHistoFilter(const Self &) = delete ; + void operator =(const Self&) = delete ; + + void SetRequestedRegion( itk::ImageBase<2> * image ) ; + + void ApplyThreshold( + typename itk::ImageRegionIterator < OutputImageType > oit , + unsigned int total ); + + std::vector< typename OutputImageType::PixelType > m_HistoThread; + InputPixelType m_Min; + InputPixelType m_Max; + InputPixelType m_NoData; + SizeType m_ThumbSize; + bool m_NoDataFlag; + double m_Step; + float m_Threshold; + unsigned int m_NbBin; + unsigned int m_ValidThreads; + +}; + +} // End namespace otb + +#ifndef OTB_MANUAL_INSTANTIATION +#include "otbComputeHistoFilter.txx" +#endif + +#endif diff --git a/Modules/Filtering/Contrast/include/otbComputeHistoFilter.txx b/Modules/Filtering/Contrast/include/otbComputeHistoFilter.txx new file mode 100644 index 0000000000000000000000000000000000000000..c29b223c648ba16b16caac2ec38431474898d233 --- /dev/null +++ b/Modules/Filtering/Contrast/include/otbComputeHistoFilter.txx @@ -0,0 +1,369 @@ +/* + * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) + * + * This file is part of Orfeo Toolbox + * + * https://www.orfeo-toolbox.org/ + * + * 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 + * + * 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 otbComputeHistoFilter_txx +#define otbComputeHistoFilter_txx + +#include "otbComputeHistoFilter.h" + +#include <limits> + +namespace otb +{ + +template <class TInputImage, class TOutputImage > +ComputeHistoFilter < TInputImage , TOutputImage > +::ComputeHistoFilter() +{ + this->SetNumberOfRequiredOutputs(2); + this->SetNthOutput( 0, this->MakeOutput(0) ); + this->SetNthOutput( 1, this->MakeOutput(1) ); + m_Min = std::numeric_limits< InputPixelType >::quiet_NaN(); + m_Max = std::numeric_limits< InputPixelType >::quiet_NaN(); + m_NoDataFlag = false; + m_NoData = std::numeric_limits< InputPixelType >::quiet_NaN(); + m_NbBin = 256; + m_Threshold = -1; + m_ThumbSize.Fill(0); + m_ValidThreads = 1; + m_Step = -1; +} + +template <class TInputImage, class TOutputImage > +itk::ProcessObject::DataObjectPointer +ComputeHistoFilter < TInputImage , TOutputImage > +::MakeOutput(itk::ProcessObject::DataObjectPointerArraySizeType idx) +{ + itk::DataObject::Pointer output; + + switch ( idx ) + { + case 0: + output = ( OutputImageType::New() ).GetPointer(); + break; + case 1: + output = ( OutputImageType::New() ).GetPointer(); + break; + default: + std::cerr << "No output " << idx << std::endl; + output = NULL; + break; + } + return output; +} + +template <class TInputImage, class TOutputImage > +itk::ProcessObject::DataObjectPointer +ComputeHistoFilter < TInputImage , TOutputImage > +::MakeOutput(const itk::ProcessObject::DataObjectIdentifierType & name) +{ + return Superclass::MakeOutput( name ); +} + +template <class TInputImage, class TOutputImage > +void ComputeHistoFilter < TInputImage , TOutputImage > +::GenerateInputRequestedRegion() +{ + typename Superclass::InputImagePointer inputPtr ( + const_cast<InputImageType *>( this->GetInput() ) ); + inputPtr->SetRequestedRegion( this->GetOutput()->GetRequestedRegion() ); + if ( inputPtr->GetRequestedRegion().GetNumberOfPixels() == 0 ) + { + inputPtr->SetRequestedRegionToLargestPossibleRegion(); + } +} + +template <class TInputImage, class TOutputImage > +void ComputeHistoFilter < TInputImage , TOutputImage > +::GenerateOutputInformation() +{ + Superclass::GenerateOutputInformation(); + typename InputImageType::ConstPointer input ( this->GetInput() ); + typename OutputImageType::Pointer output ( this->GetHistoOutput() ); + typename OutputImageType::Pointer outImage ( this->GetOutput() ); + + typename InputImageType::RegionType inputLargestRegion ( + input->GetLargestPossibleRegion()); + outImage->SetLargestPossibleRegion( inputLargestRegion ); + + typename OutputImageType::IndexType start; + typename OutputImageType::SizeType size; + + start.Fill(0); + + assert( m_ThumbSize[0] != 0 ); + assert( m_ThumbSize[1] != 0 ); + + //TODO throw error if 0 + + size[0] = std::ceil( inputLargestRegion.GetSize()[0] / + static_cast< double > ( m_ThumbSize[0] ) ); + size[1] = std::ceil( inputLargestRegion.GetSize()[1] / + static_cast< double > ( m_ThumbSize[1] ) ); + + typename OutputImageType::RegionType region; + region.SetSize(size); + region.SetIndex(start); + output->SetNumberOfComponentsPerPixel(m_NbBin); + output->SetLargestPossibleRegion(region); + typename InputImageType::SpacingType inputSpacing ( input->GetSignedSpacing() ); + typename InputImageType::PointType inputOrigin ( input->GetOrigin() ); + + typename OutputImageType::SpacingType histoSpacing ; + histoSpacing[0] = inputSpacing[0] * m_ThumbSize[0] ; + histoSpacing[1] = inputSpacing[1] * m_ThumbSize[1] ; + output->SetSignedSpacing( histoSpacing ) ; + + typename OutputImageType::PointType histoOrigin ; + histoOrigin[0] = histoSpacing[0] / 2 + inputOrigin[0] - inputSpacing[0] / 2 ; + histoOrigin[1] = histoSpacing[1] / 2 + inputOrigin[1] - inputSpacing[1] / 2 ; + output->SetOrigin( histoOrigin ); +} + +template <class TInputImage, class TOutputImage > +void ComputeHistoFilter < TInputImage , TOutputImage > +::GenerateOutputRequestedRegion( itk::DataObject * itkNotUsed(output) ) +{ + if ( GetHistoOutput()->GetRequestedRegion().GetNumberOfPixels() == 0 ) + { + GetHistoOutput()->SetRequestedRegionToLargestPossibleRegion(); + } + typename OutputImageType::Pointer outImage ( this->GetOutput() ); + SetRequestedRegion( outImage ); +} + +template <class TInputImage, class TOutputImage > +void ComputeHistoFilter < TInputImage , TOutputImage > +::GenerateData() +{ + this->AllocateOutputs(); + + // Set up the multithreaded processing + typename itk::ImageSource<OutputImageType>::ThreadStruct str; + str.Filter = this; + + // Get the output pointer + const OutputImageType * outputPtr ( this->GetOutput() ); + const itk::ImageRegionSplitterBase * splitter ( + this->GetImageRegionSplitter() ); + m_ValidThreads = + splitter->GetNumberOfSplits( outputPtr->GetRequestedRegion() , + this->GetNumberOfThreads() ); + + this->BeforeThreadedGenerateData(); + + this->GetMultiThreader()->SetNumberOfThreads( m_ValidThreads ); + this->GetMultiThreader()->SetSingleMethod(this->ThreaderCallback, &str); + // multithread the execution + this->GetMultiThreader()->SingleMethodExecute(); + + this->AfterThreadedGenerateData(); +} + +template <class TInputImage, class TOutputImage > +void ComputeHistoFilter < TInputImage , TOutputImage > +::BeforeThreadedGenerateData() +{ + // Initializing output + typename OutputImageType::Pointer output ( this->GetHistoOutput() ); + typename OutputImageType::PixelType zeroPixel(m_NbBin) ; + zeroPixel.Fill(0); + output->FillBuffer( zeroPixel ); + + // Initializing shared variable with thread number parameter + SizeType outSize ( output->GetRequestedRegion().GetSize() ); + m_HistoThread.resize( m_ValidThreads * outSize[0] * outSize[1] ); + m_HistoThread.shrink_to_fit(); + std::fill( m_HistoThread.begin() , m_HistoThread.end() , zeroPixel ); + + m_Step = static_cast<double>( m_Max - m_Min ) \ + / static_cast<double>( m_NbBin -1 ); +} + +template <class TInputImage, class TOutputImage > +void ComputeHistoFilter < TInputImage , TOutputImage > +::ThreadedGenerateData(const OutputImageRegionType & outputRegionForThread , + itk::ThreadIdType threadId ) +{ + assert(m_Step>0); + // TODO throw error + + // itk::ProgressReporter progress(this , threadId , + // outputRegionForThread.GetNumberOfPixels() ); + typename InputImageType::ConstPointer input ( this->GetInput() ); + typename OutputImageType::Pointer output ( this->GetHistoOutput() ); + + OutputImageRegionType histoRegion ( + GetHistoOutput()->GetRequestedRegion() ); + SizeType outSize ( histoRegion.GetSize() ); + IndexType outIndex ( histoRegion.GetIndex() ); + + typename InputImageType::RegionType region; + + unsigned int threadIndex ( threadId * outSize[0] * outSize[1] ) , pixel(0) ; + + for ( unsigned int nthHisto = 0 ; + nthHisto < outSize[0] * outSize[1] ; nthHisto++ ) + { + IndexType start; + start[0] = ( outIndex[0] + nthHisto % outSize[0] ) * m_ThumbSize[0]; + start[1] = ( outIndex[1] + nthHisto / outSize[0] ) * m_ThumbSize[1]; + region.SetSize( m_ThumbSize ); + region.SetIndex(start); + + if ( !region.Crop( outputRegionForThread ) ) + continue; + + typename itk::ImageRegionConstIterator < InputImageType > + it( input , region ); + InputPixelType currentPixel(0); + for ( it.GoToBegin() ; !it.IsAtEnd() ; ++it ) + { + currentPixel = it.Get(); + if( ( currentPixel == m_NoData && m_NoDataFlag ) || + currentPixel > m_Max || currentPixel < m_Min ) + continue; + + pixel = static_cast< unsigned int >( + std::round( ( currentPixel - m_Min ) / m_Step ) ); + ++m_HistoThread[threadIndex + nthHisto][pixel]; + } + + } +} + +template <class TInputImage, class TOutputImage > +void ComputeHistoFilter < TInputImage , TOutputImage > +::AfterThreadedGenerateData() +{ + typename OutputImageType::Pointer output ( this->GetHistoOutput() ); + typename itk::ImageRegionIterator < OutputImageType > + oit( output , output->GetRequestedRegion() ); + OutputImageRegionType histoRegion ( + GetHistoOutput()->GetRequestedRegion() ); + SizeType outSize ( histoRegion.GetSize() ); + IndexType outIndex ( histoRegion.GetIndex() ); + + unsigned int agreg(0) , total(0) ; + + for ( oit.GoToBegin() ; !oit.IsAtEnd() ; ++oit ) + { + total = 0; + + for ( unsigned int i = 0 ; i < m_NbBin ; i++ ) + { + agreg = 0; + + for ( unsigned int threadId = 0 ; threadId < m_ValidThreads ; threadId++ ) + { + agreg += m_HistoThread[ threadId * outSize[0] * outSize[1] + + ( oit.GetIndex()[1] - outIndex[1] ) * outSize[0] + + ( ( oit.GetIndex()[0] - outIndex[0] ) ) ][i] ; + } + oit.Get()[i] = agreg; + total += agreg; + } + if ( m_Threshold > 0 ) + ApplyThreshold( oit , total ); + } +} + +template <class TInputImage, class TOutputImage > +void ComputeHistoFilter < TInputImage , TOutputImage > +::ApplyThreshold( typename itk::ImageRegionIterator < OutputImageType > oit , + unsigned int total ) +{ + unsigned int rest(0); + unsigned int height ( static_cast<unsigned int>( + m_Threshold * ( total / m_NbBin ) ) ); + // Do i need to static_cast in a an assignation? + // Warning!!!! Need to handle out of bound int!!! + + for( unsigned int i = 0 ; i < m_NbBin ; i++ ) + { + if ( static_cast<unsigned int>( oit.Get()[i] ) > height ) + { + rest += oit.Get()[i] - height ; + oit.Get()[i] = height ; + } + } + height = rest / m_NbBin; + rest = rest % m_NbBin; + for( unsigned int i = 0 ; i < m_NbBin ; i++ ) + { + oit.Get()[i] += height ; + if ( i > (m_NbBin - rest)/2 && i <= (m_NbBin - rest)/2 + rest ) + { + ++oit.Get()[i]; + } + } +} + +template <class TInputImage, class TOutputImage > +typename TOutputImage::Pointer ComputeHistoFilter < TInputImage , TOutputImage > +::GetHistoOutput() +{ + assert( this->itk::ProcessObject::GetOutput( 1 ) ); + + return dynamic_cast< TOutputImage * >( + this->itk::ProcessObject::GetOutput(1) ); +} + +template <class TInputImage, class TOutputImage > +void ComputeHistoFilter < TInputImage , TOutputImage > +::SetRequestedRegion( itk::ImageBase<2> * image ) +{ + OutputImageRegionType histoRegion ( + GetHistoOutput()->GetRequestedRegion() ); + + IndexType start ; + start[0] = histoRegion.GetIndex()[0] * m_ThumbSize[0]; + start[1] = histoRegion.GetIndex()[1] * m_ThumbSize[1]; + + SizeType size ; + size[0] = histoRegion.GetSize()[0] * m_ThumbSize[0]; + size[1] = histoRegion.GetSize()[1] * m_ThumbSize[1]; + + typename OutputImageType::RegionType outputRequestedRegion; + outputRequestedRegion.SetIndex( start ); + outputRequestedRegion.SetSize( size ); + + outputRequestedRegion.Crop( image->GetLargestPossibleRegion() ); + image->SetRequestedRegion( outputRequestedRegion ); +} + +template <class TInputImage , class TOutputImage > +void ComputeHistoFilter < TInputImage , TOutputImage > +::PrintSelf(std::ostream& os, itk::Indent indent) const +{ + Superclass::PrintSelf(os, indent); + os << indent << "Is no data activated: " << m_NoDataFlag << std::endl; + os << indent << "No Data: " << m_NoData << std::endl; + os << indent << "Minimum: " << m_Min << std::endl; + os << indent << "Maximum: " << m_Max << std::endl; + os << indent << "Step: " << m_Step << std::endl; + os << indent << "Number of bin: " << m_NbBin << std::endl; + os << indent << "Thumbnail size: " << m_ThumbSize << std::endl; + os << indent << "Threshold value: " << m_Threshold << std::endl; +} + +} // End namespace otb + +#endif diff --git a/Modules/Filtering/Contrast/otb-module.cmake b/Modules/Filtering/Contrast/otb-module.cmake new file mode 100644 index 0000000000000000000000000000000000000000..58c014ca07a0434a46052b4eaba0806f6d24a986 --- /dev/null +++ b/Modules/Filtering/Contrast/otb-module.cmake @@ -0,0 +1,37 @@ +# +# Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) +# +# This file is part of Orfeo Toolbox +# +# https://www.orfeo-toolbox.org/ +# +# 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 +# +# 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. +# + +set(DOCUMENTATION "work in progress") + +otb_module(OTBContrast + DEPENDS + OTBImageManipulation + OTBITK + OTBCommon + OTBImageBase + + TEST_DEPENDS + OTBTestKernel + OTBImageBase + OTBImageIO + + DESCRIPTION + "${DOCUMENTATION}" +) diff --git a/Modules/Filtering/Contrast/test/CMakeLists.txt b/Modules/Filtering/Contrast/test/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..ee6950adaac2463eb8867f66454961453baa3d9b --- /dev/null +++ b/Modules/Filtering/Contrast/test/CMakeLists.txt @@ -0,0 +1,88 @@ +# +# Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) +# +# This file is part of Orfeo Toolbox +# +# https://www.orfeo-toolbox.org/ +# +# 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 +# +# 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. +# + +otb_module_test() + +set(OTBContrastTests +otbContrastTestDriver.cxx +otbComputeHistoFilterNew.cxx +otbComputeHistoFilter.cxx +otbApplyGainFilterNew.cxx +otbApplyGainFilter.cxx +otbComputeGainLutFilterNew.cxx +otbComputeGainLutFilter.cxx +otbCLHistogramEqualizationFilterNew.cxx +otbCLHistogramEqualizationFilter.cxx +otbHelperCLAHE.cxx +) + +add_executable(otbContrastTestDriver ${OTBContrastTests}) +target_link_libraries(otbContrastTestDriver ${OTBContrast-Test_LIBRARIES}) +otb_module_target_label(otbContrastTestDriver) + +otb_add_test(NAME bfTuCLHistogramEqualizationFilterNew COMMAND otbContrastTestDriver + otbCLHistogramEqualizationFilterNew) + +otb_add_test(NAME bfTuApplyGainFilterNew COMMAND otbContrastTestDriver + otbApplyGainFilterNew) + +otb_add_test(NAME bfTuComputeGainLutFilterNew COMMAND otbContrastTestDriver + otbComputeGainLutFilterNew) + +otb_add_test(NAME bfTuComputeHistoFilterNew COMMAND otbContrastTestDriver + otbComputeHistoFilterNew) + +otb_add_test(NAME bfTvComputeHistoFilter COMMAND otbContrastTestDriver + --compare-image ${EPSILON_7} + ${BASELINE}/bfTvComputeHistoFilter.tif + ${TEMP}/bfTvComputeHistoFilter.tif + otbComputeHistoFilter + ${INPUTDATA}/QB_Suburb.png + ${TEMP}/bfTvComputeHistoFilter.tif + ) + +otb_add_test(NAME bfTvComputeGainLutFilter COMMAND otbContrastTestDriver + --compare-image ${EPSILON_7} + ${BASELINE}/bfTvComputeGainLutFilter.tif + ${TEMP}/bfTvComputeGainLutFilter.tif + otbComputeGainLutFilter + ${INPUTDATA}/QB_Suburb.png + ${INPUTDATA}/inputHisto_QB_Suburb.tif + ${TEMP}/bfTvComputeGainLutFilter.tif + ) + +otb_add_test(NAME bfTvApplyGainFilter COMMAND otbContrastTestDriver + --compare-image ${EPSILON_7} + ${BASELINE}/bfTvApplyGainFilter.tif + ${TEMP}/bfTvApplyGainFilter.tif + otbApplyGainFilter + ${INPUTDATA}/QB_Suburb.png + ${INPUTDATA}/inputLut_QB_Suburb.tif + ${TEMP}/bfTvApplyGainFilter.tif + ) + +otb_add_test(NAME bfTvCLHistoEqFilter COMMAND otbContrastTestDriver + --compare-image ${EPSILON_7} + ${BASELINE}/bfTvApplyGainFilter.tif + ${TEMP}/bfTvCLHistoEqFilter.tif + otbCLHistogramEqualizationFilter + ${INPUTDATA}/QB_Suburb.png + ${TEMP}/bfTvCLHistoEqFilter.tif + ) \ No newline at end of file diff --git a/Modules/Filtering/Contrast/test/otbApplyGainFilter.cxx b/Modules/Filtering/Contrast/test/otbApplyGainFilter.cxx new file mode 100644 index 0000000000000000000000000000000000000000..f16636f5fb378c666d110b1ff808f9901f6048c3 --- /dev/null +++ b/Modules/Filtering/Contrast/test/otbApplyGainFilter.cxx @@ -0,0 +1,66 @@ +/* + * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) + * + * This file is part of Orfeo Toolbox + * + * https://www.orfeo-toolbox.org/ + * + * 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 + * + * 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. + */ + +#include "otbImageFileWriter.h" +#include "otbImageFileReader.h" +#include "otbImage.h" +#include "otbVectorImage.h" +#include "otbApplyGainFilter.h" + +int otbApplyGainFilter(int itkNotUsed(argc), char * argv []) +{ + typedef int InputPixelType; + typedef double LutPixelType; + const unsigned int Dimension = 2; + + typedef otb::Image< InputPixelType , Dimension > InputImageType; + typedef otb::VectorImage< LutPixelType , Dimension > LutImageType; + typedef otb::ApplyGainFilter + < InputImageType , LutImageType , InputImageType > FilterType; + + + typedef otb::ImageFileReader< InputImageType > ReaderType; + typedef otb::ImageFileReader< LutImageType > ReaderLutType; + typedef otb::ImageFileWriter< InputImageType > WriterType; + ReaderType::Pointer reader ( ReaderType::New() ); + ReaderLutType::Pointer readerLut ( ReaderLutType::New() ); + WriterType::Pointer writer ( WriterType::New() ); + reader->SetFileName( argv[1] ); + readerLut->SetFileName( argv[2] ); + writer->SetFileName( argv[3] ); + reader->UpdateOutputInformation(); + readerLut->UpdateOutputInformation(); + + FilterType::Pointer appGain ( FilterType::New() ); + + appGain->SetInputLut( readerLut->GetOutput() ); + appGain->SetInputImage( reader->GetOutput() ); + appGain->SetMin(0); + appGain->SetMax(255); + appGain->ThumbSizeFromSpacingOn(); + + auto size = reader->GetOutput()->GetLargestPossibleRegion().GetSize(); + appGain->SetThumbSize(size); + + writer->SetInput( appGain->GetOutput() ); + writer->Update(); + + return EXIT_SUCCESS; +} \ No newline at end of file diff --git a/Modules/Filtering/Contrast/test/otbApplyGainFilterNew.cxx b/Modules/Filtering/Contrast/test/otbApplyGainFilterNew.cxx new file mode 100644 index 0000000000000000000000000000000000000000..47f5101513aae0d798aee4f7cccc66e44fab1fe1 --- /dev/null +++ b/Modules/Filtering/Contrast/test/otbApplyGainFilterNew.cxx @@ -0,0 +1,46 @@ +/* + * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) + * + * This file is part of Orfeo Toolbox + * + * https://www.orfeo-toolbox.org/ + * + * 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 + * + * 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. + */ + + + +#include "otbImage.h" +#include "otbVectorImage.h" +#include "otbApplyGainFilter.h" + +int otbApplyGainFilterNew(int itkNotUsed(argc), char * itkNotUsed(argv) []) +{ + typedef double InputPixelType; + typedef double OutputPixelType; + typedef float LutPixelType; + const unsigned int Dimension = 2; + + typedef otb::Image< InputPixelType , Dimension > InputImageType; + typedef otb::Image< OutputPixelType , Dimension > OutputImageType; + typedef otb::VectorImage< LutPixelType , Dimension > LutImageType; + typedef otb::ApplyGainFilter + < InputImageType , LutImageType , OutputImageType > FilterType; + + + FilterType::Pointer appGain ( FilterType::New() ); + + std::cout << appGain << std::endl; + + return EXIT_SUCCESS; +} \ No newline at end of file diff --git a/Modules/Filtering/Contrast/test/otbCLHistogramEqualizationFilter.cxx b/Modules/Filtering/Contrast/test/otbCLHistogramEqualizationFilter.cxx new file mode 100644 index 0000000000000000000000000000000000000000..7ddbe3e6fad80870d50a1f1fdf9edc3f63976e4f --- /dev/null +++ b/Modules/Filtering/Contrast/test/otbCLHistogramEqualizationFilter.cxx @@ -0,0 +1,57 @@ +/* + * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) + * + * This file is part of Orfeo Toolbox + * + * https://www.orfeo-toolbox.org/ + * + * 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 + * + * 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. + */ + +#include "otbImageFileReader.h" +#include "otbImageFileWriter.h" +#include "otbImage.h" +#include "otbCLHistogramEqualizationFilter.h" + +int otbCLHistogramEqualizationFilter(int itkNotUsed(argc), char * argv []) +{ + typedef int InputPixelType; + const unsigned int Dimension = 2; + + typedef otb::Image< InputPixelType , Dimension > InputImageType; + typedef otb::CLHistogramEqualizationFilter< InputImageType , InputImageType > + FilterType; + typedef otb::ImageFileReader< InputImageType > ReaderType; + typedef otb::ImageFileWriter< InputImageType > WriterType; + + ReaderType::Pointer reader ( ReaderType::New() ); + WriterType::Pointer writer ( WriterType::New() ); + reader->SetFileName( argv[1] ); + writer->SetFileName( argv[2] ); + reader->UpdateOutputInformation(); + + FilterType::Pointer histoEqualize ( FilterType::New() ); + + histoEqualize->SetInput( reader->GetOutput() ); + histoEqualize->SetMin(0); + histoEqualize->SetMax(255); + histoEqualize->SetNbBin(256); + auto size = reader->GetOutput()->GetLargestPossibleRegion().GetSize(); + size[0] /= 4; + size[1] /= 4; + histoEqualize->SetThumbSize( size ); + + writer->SetInput( histoEqualize->GetOutput() ); + writer->Update(); + return EXIT_SUCCESS; +} diff --git a/Modules/Filtering/Contrast/test/otbCLHistogramEqualizationFilterNew.cxx b/Modules/Filtering/Contrast/test/otbCLHistogramEqualizationFilterNew.cxx new file mode 100644 index 0000000000000000000000000000000000000000..4978121896f0d9480ea16427d590a3a69f55233c --- /dev/null +++ b/Modules/Filtering/Contrast/test/otbCLHistogramEqualizationFilterNew.cxx @@ -0,0 +1,41 @@ +/* + * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) + * + * This file is part of Orfeo Toolbox + * + * https://www.orfeo-toolbox.org/ + * + * 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 + * + * 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. + */ + + + +#include "otbImage.h" +#include "otbCLHistogramEqualizationFilter.h" + +int otbCLHistogramEqualizationFilterNew(int itkNotUsed(argc), char * itkNotUsed(argv) []) +{ + typedef double InputPixelType; + const unsigned int Dimension = 2; + + typedef otb::Image< InputPixelType , Dimension > InputImageType; + typedef otb::CLHistogramEqualizationFilter< InputImageType , InputImageType > + FilterType; + + + FilterType::Pointer histoEqualize ( FilterType::New() ); + + std::cout << histoEqualize << std::endl; + + return EXIT_SUCCESS; +} \ No newline at end of file diff --git a/Modules/Filtering/Contrast/test/otbComputeGainLutFilter.cxx b/Modules/Filtering/Contrast/test/otbComputeGainLutFilter.cxx new file mode 100644 index 0000000000000000000000000000000000000000..77a8c3fdeb108624b237ef71298aa7a89b0dc9ac --- /dev/null +++ b/Modules/Filtering/Contrast/test/otbComputeGainLutFilter.cxx @@ -0,0 +1,68 @@ +/* + * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) + * + * This file is part of Orfeo Toolbox + * + * https://www.orfeo-toolbox.org/ + * + * 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 + * + * 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. + */ + +#include "otbVectorImage.h" +#include "otbImage.h" +#include "otbImageFileWriter.h" +#include "otbImageFileReader.h" +#include "otbComputeGainLutFilter.h" + +int otbComputeGainLutFilter(int itkNotUsed(argc), char * argv []) +{ + typedef int InputPixelType; + typedef double OutputPixelType; + const unsigned int Dimension = 2; + + typedef otb::VectorImage< InputPixelType , Dimension > HistoImageType; + typedef otb::Image< InputPixelType , Dimension > InputImageType; + typedef otb::VectorImage< OutputPixelType , Dimension > LutImageType; + typedef otb::ComputeGainLutFilter< HistoImageType , LutImageType > FilterType; + + typedef otb::ImageFileReader< InputImageType > ReaderImageType; + typedef otb::ImageFileReader< HistoImageType > ReaderType; + typedef otb::ImageFileWriter< LutImageType > WriterType; + ReaderImageType::Pointer readerImage( ReaderImageType::New() ); + ReaderType::Pointer reader( ReaderType::New() ); + WriterType::Pointer writer ( WriterType::New() ); + readerImage->SetFileName( argv[1] ); + reader->SetFileName( argv[2] ); + writer->SetFileName( argv[3] ); + reader->UpdateOutputInformation(); + readerImage->UpdateOutputInformation(); + + FilterType::Pointer computeGainLut ( FilterType::New() ); + + computeGainLut->SetInput( reader->GetOutput() ); + computeGainLut->SetMin(0); + computeGainLut->SetMax(255); + auto size = readerImage->GetOutput()->GetLargestPossibleRegion().GetSize(); + size[0] /= 4; + size[1] /= 4; + auto nbPix = size[0]*size[1] ; + computeGainLut->SetNbPixel( nbPix ); + + writer->SetInput( computeGainLut->GetOutput() ); + writer->Update(); + + auto index = computeGainLut->GetOutput()->GetLargestPossibleRegion().GetIndex(); + std::cout<<computeGainLut->GetOutput()->GetPixel( index )<<std::endl; + + return EXIT_SUCCESS; +} diff --git a/Modules/Filtering/Contrast/test/otbComputeGainLutFilterNew.cxx b/Modules/Filtering/Contrast/test/otbComputeGainLutFilterNew.cxx new file mode 100644 index 0000000000000000000000000000000000000000..642e805fb0e5f8e83b5651bf393b42d1fc7e534f --- /dev/null +++ b/Modules/Filtering/Contrast/test/otbComputeGainLutFilterNew.cxx @@ -0,0 +1,40 @@ +/* + * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) + * + * This file is part of Orfeo Toolbox + * + * https://www.orfeo-toolbox.org/ + * + * 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 + * + * 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. + */ + +#include "otbVectorImage.h" +#include "otbComputeGainLutFilter.h" + +int otbComputeGainLutFilterNew(int itkNotUsed(argc), char * itkNotUsed(argv) []) +{ + typedef double InputPixelType; + typedef double OutputPixelType; + const unsigned int Dimension = 2; + + typedef otb::VectorImage< InputPixelType , Dimension > HistoImageType; + typedef otb::VectorImage< OutputPixelType , Dimension > LutImageType; + typedef otb::ComputeGainLutFilter< HistoImageType , LutImageType > FilterType; + + + FilterType::Pointer computeGainLut ( FilterType::New() ); + + std::cout << computeGainLut << std::endl; + + return EXIT_SUCCESS; +} \ No newline at end of file diff --git a/Modules/Filtering/Contrast/test/otbComputeHistoFilter.cxx b/Modules/Filtering/Contrast/test/otbComputeHistoFilter.cxx new file mode 100644 index 0000000000000000000000000000000000000000..f8407bee8d9123c51d012daa8240cdb45673bf04 --- /dev/null +++ b/Modules/Filtering/Contrast/test/otbComputeHistoFilter.cxx @@ -0,0 +1,59 @@ +/* + * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) + * + * This file is part of Orfeo Toolbox + * + * https://www.orfeo-toolbox.org/ + * + * 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 + * + * 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. + */ + +#include "otbImage.h" +#include "otbImageFileWriter.h" +#include "otbImageFileReader.h" +#include "otbVectorImage.h" +#include "otbComputeHistoFilter.h" + +int otbComputeHistoFilter(int itkNotUsed(argc), char * argv []) +{ + typedef int InputPixelType; + typedef int OutputPixelType; + const unsigned int Dimension = 2; + + typedef otb::Image< InputPixelType , Dimension > InputImageType; + typedef otb::VectorImage< OutputPixelType , Dimension > HistoImageType; + typedef otb::ComputeHistoFilter< InputImageType , HistoImageType > FilterType; + + typedef otb::ImageFileReader< InputImageType > ReaderType; + typedef otb::ImageFileWriter< HistoImageType > WriterType; + ReaderType::Pointer reader ( ReaderType::New() ); + WriterType::Pointer writer ( WriterType::New() ); + reader->SetFileName( argv[1] ); + writer->SetFileName( argv[2] ); + reader->UpdateOutputInformation(); + + FilterType::Pointer computeHisto ( FilterType::New() ); + + computeHisto->SetInput( reader->GetOutput() ); + computeHisto->SetMin(0); + computeHisto->SetMax(255); + computeHisto->SetNbBin(256); + auto size = reader->GetOutput()->GetLargestPossibleRegion().GetSize(); + size[0] /= 4; + size[1] /= 4; + computeHisto->SetThumbSize( size ); + + writer->SetInput( computeHisto->GetHistoOutput() ); + writer->Update(); + return EXIT_SUCCESS; +} \ No newline at end of file diff --git a/Modules/Filtering/Contrast/test/otbComputeHistoFilterNew.cxx b/Modules/Filtering/Contrast/test/otbComputeHistoFilterNew.cxx new file mode 100644 index 0000000000000000000000000000000000000000..4dd79d460e4b9747fa83dd22f1aadf1593ec0c7b --- /dev/null +++ b/Modules/Filtering/Contrast/test/otbComputeHistoFilterNew.cxx @@ -0,0 +1,43 @@ +/* + * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) + * + * This file is part of Orfeo Toolbox + * + * https://www.orfeo-toolbox.org/ + * + * 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 + * + * 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. + */ + + + +#include "otbImage.h" +#include "otbVectorImage.h" +#include "otbComputeHistoFilter.h" + +int otbComputeHistoFilterNew(int itkNotUsed(argc), char * itkNotUsed(argv) []) +{ + typedef double InputPixelType; + typedef double OutputPixelType; + const unsigned int Dimension = 2; + + typedef otb::Image< InputPixelType , Dimension > InputImageType; + typedef otb::VectorImage< OutputPixelType , Dimension > HistoImageType; + typedef otb::ComputeHistoFilter< InputImageType , HistoImageType > FilterType; + + + FilterType::Pointer computeHisto ( FilterType::New() ); + + std::cout << computeHisto << std::endl; + + return EXIT_SUCCESS; +} \ No newline at end of file diff --git a/Modules/Filtering/Contrast/test/otbContrastTestDriver.cxx b/Modules/Filtering/Contrast/test/otbContrastTestDriver.cxx new file mode 100644 index 0000000000000000000000000000000000000000..213bb726f31b703ac603c771dff16564b63ff856 --- /dev/null +++ b/Modules/Filtering/Contrast/test/otbContrastTestDriver.cxx @@ -0,0 +1,34 @@ +/* + * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) + * + * This file is part of Orfeo Toolbox + * + * https://www.orfeo-toolbox.org/ + * + * 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 + * + * 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. + */ + +#include "otbTestMain.h" + +void RegisterTests() +{ + REGISTER_TEST(otbComputeHistoFilterNew); + REGISTER_TEST(otbComputeHistoFilter); + REGISTER_TEST(otbComputeGainLutFilterNew); + REGISTER_TEST(otbComputeGainLutFilter); + REGISTER_TEST(otbApplyGainFilterNew); + REGISTER_TEST(otbApplyGainFilter); + REGISTER_TEST(otbCLHistogramEqualizationFilterNew); + REGISTER_TEST(otbCLHistogramEqualizationFilter); + REGISTER_TEST(otbHelperCLAHE); +} diff --git a/Modules/Filtering/Contrast/test/otbHelperCLAHE.cxx b/Modules/Filtering/Contrast/test/otbHelperCLAHE.cxx new file mode 100644 index 0000000000000000000000000000000000000000..b1fc8953f6322978d64f33492bbd41ca6c789dcd --- /dev/null +++ b/Modules/Filtering/Contrast/test/otbHelperCLAHE.cxx @@ -0,0 +1,92 @@ +/* + * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) + * + * This file is part of Orfeo Toolbox + * + * https://www.orfeo-toolbox.org/ + * + * 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 + * + * 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. + */ + +#include "otbImageFileReader.h" +#include "otbImageFileWriter.h" +#include "otbImage.h" +#include "otbVectorImage.h" +#include "itkImageRegionIterator.h" +#include "otbCLHistogramEqualizationFilter.h" + +int otbHelperCLAHE(int itkNotUsed(argc), char * argv []) +{ + typedef int InputPixelType; + const unsigned int Dimension = 2; + + typedef otb::Image< InputPixelType , Dimension > InputImageType; + typedef otb::ImageFileReader< InputImageType > ReaderType; + typedef otb::ImageFileWriter< InputImageType > WriterType; + ReaderType::Pointer reader ( ReaderType::New() ); + WriterType::Pointer writer ( WriterType::New() ); + reader->SetFileName( argv[1] ); + writer->SetFileName( argv[2] ); + reader->UpdateOutputInformation(); + + typedef otb::VectorImage< int , 2 > HistogramType; + typedef otb::VectorImage< double , 2 > LutType; + + typedef itk::StreamingImageFilter< LutType , LutType > + StreamingImageFilter; + + typedef otb::InPlacePassFilter < InputImageType > BufferFilter; + + typedef otb::ComputeHistoFilter< InputImageType , HistogramType > + HistoFilter; + + typedef otb::ComputeGainLutFilter< HistogramType , LutType > + GainLutFilter; + + typedef otb::ApplyGainFilter< InputImageType , LutType , InputImageType > + ApplyGainFilter; + + HistoFilter::Pointer histoFilter( HistoFilter::New() ); + GainLutFilter::Pointer lutFilter( GainLutFilter::New() ); + ApplyGainFilter::Pointer applyFilter( ApplyGainFilter::New() ); + BufferFilter::Pointer buffer( BufferFilter::New() ); + StreamingImageFilter::Pointer streamFilter( StreamingImageFilter::New() ); + InputImageType::SizeType size; + size = reader->GetOutput()->GetLargestPossibleRegion().GetSize(); + // histoEqualize->SetThreshold(100); + + histoFilter->SetInput( reader->GetOutput() ); + buffer->SetInput( reader->GetOutput() ); + lutFilter->SetInput( histoFilter->GetHistoOutput() ); + streamFilter->SetInput( lutFilter->GetOutput() ); + applyFilter->SetInputLut( streamFilter->GetOutput() ); + applyFilter->SetInputImage( buffer->GetOutput() ); + + histoFilter->SetMin(0); + histoFilter->SetMax(255); + lutFilter->SetMin(0); + lutFilter->SetMax(255); + applyFilter->SetMin(0); + applyFilter->SetMax(255); + + histoFilter->SetThumbSize(size); + applyFilter->SetThumbSize(size); + lutFilter->SetNbPixel(size[0]*size[1]); + + + + writer->SetInput(applyFilter->GetOutput()); + writer->Update(); + + return EXIT_SUCCESS; +} diff --git a/Modules/Filtering/ImageManipulation/include/otbInPlacePassFilter.h b/Modules/Filtering/ImageManipulation/include/otbInPlacePassFilter.h new file mode 100644 index 0000000000000000000000000000000000000000..22085dd00f9b5da1e3962550d4788e99d7855640 --- /dev/null +++ b/Modules/Filtering/ImageManipulation/include/otbInPlacePassFilter.h @@ -0,0 +1,99 @@ +/* + * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES) + * + * This file is part of Orfeo Toolbox + * + * https://www.orfeo-toolbox.org/ + * + * 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 + * + * 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 otbInPlacePassFilter_h +#define otbInPlacePassFilter_h + +#include "itkInPlaceImageFilter.h" +#include "otbImage.h" +#include "itkImageRegionIterator.h" + +namespace otb +{ + +/** \class InPlacePassFilter + * \brief This filter has the only purpose to recall regions + * + * This class is implemented to recall regions. Due to ITK implementation + * if the pipeline of the algorithm has branch (diamond) one might have + * an input with two requested regions : one from branch 1 (A) and one from + * branch 2 (B). Problem is after updating and generating data on branch 1 + * branch 2 will not propagate its region B again,and will use region A + * instead. By memorizing the region this buffer filter can be placed in + * in front of each branch so that the requested region will be saved.v + * + * \ingroup OTBImageManipulation + */ + +template < class TInputImage > +class ITK_EXPORT InPlacePassFilter : + public itk::InPlaceImageFilter < TInputImage , TInputImage > +{ +public: + /** typedef for standard classes. */ + typedef TInputImage InputImageType; + + typedef InPlacePassFilter Self; + typedef itk::InPlaceImageFilter< InputImageType, InputImageType > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self) + + /** Run-time type information (and related methods). */ + itkTypeMacro(InPlacePassFilter, InPlaceImageFilter) + +protected: + InPlacePassFilter() { + this->InPlaceOn(); + } + + ~InPlacePassFilter() override {} + + void ThreadedGenerateData( + const typename InputImageType::RegionType & + outputRegionForThread , + itk::ThreadIdType itkNotUsed(threadId) ) override + { + typename InputImageType::ConstPointer input ( this->GetInput() ); + typename InputImageType::Pointer output ( this->GetOutput() ); + itk::ImageRegionConstIterator < InputImageType > it ( input , + outputRegionForThread ); + itk::ImageRegionIterator < InputImageType > oit ( output , + outputRegionForThread ); + for ( oit.GoToBegin() , it.GoToBegin() ; !oit.IsAtEnd() || !it.IsAtEnd() ; + ++it , ++oit ) + { + oit.Set(it.Get()); + } + } +private: + InPlacePassFilter(const Self &) = delete ; + void operator =(const Self&) = delete ; + +}; + +} // End namespace otb + +#ifndef OTB_MANUAL_INSTANTIATION +#endif + +#endif