Skip to content
Snippets Groups Projects
Commit e868ff87 authored by Thomas Feuvrier's avatar Thomas Feuvrier
Browse files

nomsg

parent 4a469c41
No related branches found
No related tags found
No related merge requests found
/*=========================================================================
Programme : OTB (ORFEO ToolBox)
Auteurs : CS - P.Imbo
Language : C++
Date : 24 janvier 2006
Role : Filtre de débruitage de FROST sur une image
$Id$
=========================================================================*/
#ifndef __otbFrostImageFilter_txx
#define __otbFrostImageFilter_txx
#include "otbFrostImageFilter.h"
#include "itkDataObject.h"
#include "itkExceptionObject.h"
#include "itkConstNeighborhoodIterator.h"
#include "itkNeighborhoodInnerProduct.h"
#include "itkImageRegionIterator.h"
#include "itkNeighborhoodAlgorithm.h"
#include "itkZeroFluxNeumannBoundaryCondition.h"
#include "itkOffset.h"
#include "itkProgressReporter.h"
#include <math.h>
namespace otb
{
/**
*
*/
template <class TInputImage, class TOutputImage>
FrostImageFilter<TInputImage, TOutputImage>::FrostImageFilter()
{
m_Radius.Fill(1);
}
template <class TInputImage, class TOutputImage>
void FrostImageFilter<TInputImage, TOutputImage>::GenerateInputRequestedRegion() throw (itk::InvalidRequestedRegionError)
{
// call the superclass' implementation of this method
Superclass::GenerateInputRequestedRegion();
// get pointers to the input and output
typename Superclass::InputImagePointer inputPtr = const_cast< TInputImage * >( this->GetInput() );
typename Superclass::OutputImagePointer outputPtr = this->GetOutput();
if ( !inputPtr || !outputPtr )
{
return;
}
// get a copy of the input requested region (should equal the output
// requested region)
typename TInputImage::RegionType inputRequestedRegion;
inputRequestedRegion = inputPtr->GetRequestedRegion();
// pad the input requested region by the operator radius
inputRequestedRegion.PadByRadius( m_Radius );
// crop the input requested region at the input's largest possible region
if ( inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()) )
{
inputPtr->SetRequestedRegion( inputRequestedRegion );
return;
}
else
{
// Couldn't crop the region (requested region is outside the largest
// possible region). Throw an exception.
// store what we tried to request (prior to trying to crop)
inputPtr->SetRequestedRegion( inputRequestedRegion );
// build an exception
itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
itk::OStringStream msg;
msg << static_cast<const char *>(this->GetNameOfClass())
<< "::GenerateInputRequestedRegion()";
e.SetLocation(msg.str().c_str());
e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
e.SetDataObject(inputPtr);
throw e;
}
}
template< class TInputImage, class TOutputImage>
void FrostImageFilter< TInputImage, TOutputImage>::ThreadedGenerateData(
const OutputImageRegionType& outputRegionForThread,
int threadId
)
{
unsigned int i;
itk::ZeroFluxNeumannBoundaryCondition<InputImageType> nbc;
itk::ConstNeighborhoodIterator<InputImageType> bit;
typename itk::ConstNeighborhoodIterator<InputImageType>::OffsetType off;
itk::ImageRegionIterator<OutputImageType> it;
// Allocate output
typename OutputImageType::Pointer output = this->GetOutput();
typename InputImageType::ConstPointer input = this->GetInput();
// Find the data-set boundary "faces"
typename itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType>::FaceListType faceList;
typename itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType>::FaceListType::iterator fit;
itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType> bC;
faceList = bC(input, outputRegionForThread, m_Radius);
// support progress methods/callbacks
itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
// InputRealType pixel;
InputRealType sum;
InputRealType sum2;
double Moyenne, Variance;
double Alpha;
double NormFiltre;
double FrostFiltre;
double CoefFiltre;
double dPixel;
// Process each of the boundary faces. These are N-d regions which border
// the edge of the buffer.
for (fit=faceList.begin(); fit != faceList.end(); ++fit)
{
bit = itk::ConstNeighborhoodIterator<InputImageType>(m_Radius, input, *fit);
unsigned int neighborhoodSize = bit.Size();
unsigned int CenterPos = bit.GetCenterNeighborhoodIndex();
it = itk::ImageRegionIterator<OutputImageType>(output, *fit);
bit.OverrideBoundaryCondition(&nbc);
bit.GoToBegin();
while ( ! bit.IsAtEnd() )
{
sum = itk::NumericTraits<InputRealType>::Zero;
sum2 = itk::NumericTraits<InputRealType>::Zero;
for (i = 0; i < neighborhoodSize; ++i)
{
dPixel = static_cast<double>( bit.GetPixel(i) );
sum += dPixel;
sum2 += dPixel * dPixel;
}
Moyenne = sum / double(neighborhoodSize);
Variance = sum2 / double(neighborhoodSize) - Moyenne * Moyenne ;
Alpha = m_Deramp * Variance / (Moyenne * Moyenne) ;
NormFiltre = 0.0;
FrostFiltre = 0.0;
const double rad_x = static_cast<double>(m_Radius[0]);
const double rad_y = static_cast<double>(m_Radius[1]);
for (double x = -rad_x; x<= rad_x; x++)
{
for (double y = -rad_y; y <= rad_y; y++)
{
double Dist = double(sqrt(x*x+y*y));
off[0]= static_cast<int>(x);
off[1]= static_cast<int>(y);
// i = (unsigned int)((y+rad_y)*(2*rad_y+1)+(x+rad_x));
dPixel= static_cast<double>( bit.GetPixel(off));
// dPixel= static_cast<double>( bit.GetPixel(i));
CoefFiltre = Alpha * exp(-Alpha *Dist);
NormFiltre = NormFiltre + CoefFiltre;
FrostFiltre = FrostFiltre + (CoefFiltre * dPixel);
}
}
dPixel = (FrostFiltre/NormFiltre);
if (finite(dPixel)==0){
dPixel = 0.;
}
it.Set( static_cast<OutputPixelType>( dPixel ) );
++bit;
++it;
progress.CompletedPixel();
}
}
}
/**
* Standard "PrintSelf" method
*/
template <class TInputImage, class TOutput>
void
FrostImageFilter<TInputImage, TOutput>::PrintSelf(std::ostream& os, itk::Indent indent) const
{
Superclass::PrintSelf( os, indent );
os << indent << "Radius: " << m_Radius << std::endl;
}
} // end namespace otb
#endif
/*=========================================================================
Programme : OTB (ORFEO ToolBox)
Auteurs : CS - P.Imbo
Language : C++
Date : 23 janvier 2006
Role : Filtre de débruitage de LEE sur une image
$Id$
=========================================================================*/
#ifndef __otbLeeImageFilter_h
#define __otbLeeImageFilter_h
#include "itkImageToImageFilter.h"
#include "itkImage.h"
#include "itkNumericTraits.h"
namespace otb
{
/** \class LeeImageFilter
* \brief Applique un filtre de débruitage sur une image.
*
* Ce filtre correspond au filtre Anti-Speckle de LEE :
*
* R = E[I] + b(I-E[I]) avec b = C²r / ( C²r + C²v )
* Cv = 1 / sqrt(L) avec L le nombre de vues
* Cr = sqrt(Var(I)) / E[I] avec Var(I) = E[I²] - E[I]²
*
*/
template <class TInputImage, class TOutputImage>
class LeeImageFilter : public itk::ImageToImageFilter< TInputImage, TOutputImage >
{
public:
/** Extrait les dimensions aussi bien des images
d'entrée (Input) que de sortie (Output). */
itkStaticConstMacro( InputImageDimension,
unsigned int,
TInputImage::ImageDimension);
itkStaticConstMacro( OutputImageDimension,
unsigned int,
TOutputImage::ImageDimension);
/** "typedef" pour simplifier la définition et la déclaration de variables. */
typedef TInputImage InputImageType;
/** "typedef" pour simplifier la définition et la déclaration de variables. */
typedef TOutputImage OutputImageType;
/** "typedef" pour les classes standards. */
typedef LeeImageFilter Self;
typedef itk::ImageToImageFilter< InputImageType, OutputImageType> Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
/** Methode pour la gestion "object factory". */
itkNewMacro(Self);
/** Retourne le nom de la classe. */
itkTypeMacro(LeeImageFilter, ImageToImageFilter);
/** Définition des images supportées. */
typedef typename InputImageType::PixelType InputPixelType;
typedef typename OutputImageType::PixelType OutputPixelType;
/** "typedef" définissant un réel dans cette classe. */
typedef typename itk::NumericTraits<InputPixelType>::RealType InputRealType;
typedef typename InputImageType::RegionType InputImageRegionType;
typedef typename OutputImageType::RegionType OutputImageRegionType;
/** "typedef" définissant la taille d'une image. */
typedef typename InputImageType::SizeType SizeType;
/** Positionne le rayon définissant le voisinage utilisé pour le calcul du filtre. */
itkSetMacro(Radius, SizeType);
/** Récupère le rayon définissant le voisinage utilisé pour le calcul du filtre. */
itkGetConstReferenceMacro(Radius, SizeType);
/** Positionne le nombre de vues utilisé pour le calcul du filtre. */
itkSetMacro(NbVues, double);
/** Récupère le nombre de vues (référencé constant) utilisé pour le calcul du filtre. */
itkGetConstReferenceMacro(NbVues, double);
/** LeeImageFilter a besoin d'une zone de traitement plus large en entrée qu'en sortie
* pemettant une utilisation needs a larger input requested region than
* the output requested region. As such, LeeImageFilter needs
* to provide an implementation for GenerateInputRequestedRegion()
* in order to inform the pipeline execution model.
*
* \sa ImageToImageFilter::GenerateInputRequestedRegion() */
virtual void GenerateInputRequestedRegion() throw(itk::InvalidRequestedRegionError);
protected:
LeeImageFilter();
virtual ~LeeImageFilter() {};
void PrintSelf(std::ostream& os, itk::Indent indent) const;
/** LeeImageFilter peut etre implémentée pour un traitement de filtre multithreaded.
* Ainsi, cette implémentation fournit la méthode ThreadedGenerateData()
* qui est appelée pour chaque thread du process. Les données image sont allouées automatiquement
* par la classe "mère" en appelant la méthode ThreadedGenerateData(). ThreadedGenerateData peut seulement
* écrire la portion de l'image spécifiée par le paramètre "outputRegionForThread"
*
* Filtre de LEE :
* R = E[I] + b(I-E[I]) avec b = C²r / ( C²r + C²v )
* Cv = 1 / sqrt(L) avec L le nombre de vues
* Cr = sqrt(Var(I)) / E[I] avec Var(I) = E[I²] - E[I]²
*
* \sa ImageToImageFilter::ThreadedGenerateData(),
* ImageToImageFilter::GenerateData() */
void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
int threadId );
private:
LeeImageFilter(const Self&); //purposely not implemented
void operator=(const Self&); //purposely not implemented
/** Déclaration du rayon */
SizeType m_Radius;
/** Déclaration du nombre de vues du filtre */
double m_NbVues;
};
} // end namespace otb
#ifndef OTB_MANUAL_INSTANTIATION
#include "otbLeeImageFilter.txx"
#endif
#endif
/*=========================================================================
Programme : OTB (ORFEO ToolBox)
Auteurs : CS - P.Imbo
Language : C++
Date : 23 janvier 2006
Role : Filtre de débruitage de LEE sur une image
$Id$
=========================================================================*/
#ifndef __otbLeeImageFilter_txx
#define __otbLeeImageFilter_txx
#include "otbLeeImageFilter.h"
#include "itkDataObject.h"
#include "itkExceptionObject.h"
#include "itkConstNeighborhoodIterator.h"
#include "itkNeighborhoodInnerProduct.h"
#include "itkImageRegionIterator.h"
#include "itkNeighborhoodAlgorithm.h"
#include "itkZeroFluxNeumannBoundaryCondition.h"
#include "itkOffset.h"
#include "itkProgressReporter.h"
namespace otb
{
/**
*
*/
template <class TInputImage, class TOutputImage>
LeeImageFilter<TInputImage, TOutputImage>::LeeImageFilter()
{
m_Radius.Fill(1);
}
template <class TInputImage, class TOutputImage>
void LeeImageFilter<TInputImage, TOutputImage>::GenerateInputRequestedRegion() throw (itk::InvalidRequestedRegionError)
{
// call the superclass' implementation of this method
Superclass::GenerateInputRequestedRegion();
// get pointers to the input and output
typename Superclass::InputImagePointer inputPtr = const_cast< TInputImage * >( this->GetInput() );
typename Superclass::OutputImagePointer outputPtr = this->GetOutput();
if ( !inputPtr || !outputPtr )
{
return;
}
// get a copy of the input requested region (should equal the output
// requested region)
typename TInputImage::RegionType inputRequestedRegion;
inputRequestedRegion = inputPtr->GetRequestedRegion();
// pad the input requested region by the operator radius
inputRequestedRegion.PadByRadius( m_Radius );
// crop the input requested region at the input's largest possible region
if ( inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()) )
{
inputPtr->SetRequestedRegion( inputRequestedRegion );
return;
}
else
{
// Couldn't crop the region (requested region is outside the largest
// possible region). Throw an exception.
// store what we tried to request (prior to trying to crop)
inputPtr->SetRequestedRegion( inputRequestedRegion );
// build an exception
itk::InvalidRequestedRegionError e(__FILE__, __LINE__);
itk::OStringStream msg;
msg << static_cast<const char *>(this->GetNameOfClass())
<< "::GenerateInputRequestedRegion()";
e.SetLocation(msg.str().c_str());
e.SetDescription("Requested region is (at least partially) outside the largest possible region.");
e.SetDataObject(inputPtr);
throw e;
}
}
template< class TInputImage, class TOutputImage>
void LeeImageFilter< TInputImage, TOutputImage>::ThreadedGenerateData(
const OutputImageRegionType& outputRegionForThread,
int threadId
)
{
unsigned int i;
itk::ZeroFluxNeumannBoundaryCondition<InputImageType> nbc;
itk::ConstNeighborhoodIterator<InputImageType> bit;
itk::ImageRegionIterator<OutputImageType> it;
// Allocate output
typename OutputImageType::Pointer output = this->GetOutput();
typename InputImageType::ConstPointer input = this->GetInput();
// Find the data-set boundary "faces"
typename itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType>::FaceListType faceList;
typename itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType>::FaceListType::iterator fit;
itk::NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType> bC;
faceList = bC(input, outputRegionForThread, m_Radius);
// support progress methods/callbacks
itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
// InputRealType pixel;
InputRealType sum;
InputRealType sum2;
double Cr2, Cv2, E_I, I, Var_I, dPixel;
// dPixel = this->getNbVues();
dPixel = m_NbVues;
//Calcul du rapport
Cv2 = 1./(sqrt(dPixel));
Cv2 = Cv2*Cv2;
// Process each of the boundary faces. These are N-d regions which border
// the edge of the buffer.
for (fit=faceList.begin(); fit != faceList.end(); ++fit)
{
bit = itk::ConstNeighborhoodIterator<InputImageType>(m_Radius, input, *fit);
unsigned int neighborhoodSize = bit.Size();
it = itk::ImageRegionIterator<OutputImageType>(output, *fit);
bit.OverrideBoundaryCondition(&nbc);
bit.GoToBegin();
while ( ! bit.IsAtEnd() )
{
sum = itk::NumericTraits<InputRealType>::Zero;
sum2 = itk::NumericTraits<InputRealType>::Zero;
//Parcours du voisinage
for (i = 0; i < neighborhoodSize; ++i)
{
dPixel = static_cast<double>( bit.GetPixel(i) );
sum += dPixel;
sum2 += dPixel * dPixel;
}
E_I = sum / double(neighborhoodSize);
Var_I = sum2 / double(neighborhoodSize) - E_I*E_I ;
I = static_cast<double>( bit.GetCenterPixel() );
Cr2 = Var_I / (E_I*E_I);
dPixel = E_I + ((I-E_I)*(Cr2))/(Cr2 + Cv2);
// get the mean value
it.Set( static_cast<OutputPixelType>( dPixel ) );
++bit;
++it;
progress.CompletedPixel();
}
}
}
/**
* Standard "PrintSelf" method
*/
template <class TInputImage, class TOutput>
void
LeeImageFilter<TInputImage, TOutput>::PrintSelf(std::ostream& os, itk::Indent indent) const
{
Superclass::PrintSelf( os, indent );
os << indent << "Radius: " << m_Radius << std::endl;
}
} // end namespace otb
#endif
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment