Skip to content

GitLab

Commit e9c0d74d authored by Emmanuel Christophe's avatar Emmanuel Christophe
Browse files

STYLE: formatting indentation

parent fcf46b20
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 3709 additions and 3630 deletions
Code/BasicFilters/otbBSplineDecompositionImageFilter.h
View file @ e9c0d74d
......@@ -36,7 +36,7 @@ namespace otb
*/
template <class TInputImage, class TOutputImage>
class ITK_EXPORT BSplineDecompositionImageFilter :
public itk::ImageToImageFilter<TInputImage,TOutputImage>
public itk::ImageToImageFilter<TInputImage,TOutputImage>
{
public:
/** Standard class typedefs. */
......
Code/BasicFilters/otbBSplineDecompositionImageFilter.txx
View file @ e9c0d74d
......@@ -69,42 +69,42 @@ BSplineDecompositionImageFilter<TInputImage, TOutputImage>
double c0 = 1.0;
if (m_DataLength[m_IteratorDirection] == 1) //Required by mirror boundaries
{
{
return false;
}
}
// Compute overall gain
for (int k = 0; k < m_NumberOfPoles; k++)
{
{
// Note for cubic splines lambda = 6
c0 = c0 * (1.0 - m_SplinePoles[k]) * (1.0 - 1.0 / m_SplinePoles[k]);
}
}
// apply the gain
for (unsigned int n = 0; n < m_DataLength[m_IteratorDirection]; n++)
{
{
m_Scratch[n] *= c0;
}
}
// loop over all poles
for (int k = 0; k < m_NumberOfPoles; k++)
{
{
// causal initialization
this->SetInitialCausalCoefficient(m_SplinePoles[k]);
// causal recursion
for (unsigned int n = 1; n < m_DataLength[m_IteratorDirection]; n++)
{
{
m_Scratch[n] += m_SplinePoles[k] * m_Scratch[n - 1];
}
}
// anticausal initialization
this->SetInitialAntiCausalCoefficient(m_SplinePoles[k]);
// anticausal recursion
for ( int n = m_DataLength[m_IteratorDirection] - 2; 0 <= n; n--)
{
{
m_Scratch[n] = m_SplinePoles[k] * (m_Scratch[n + 1] - m_Scratch[n]);
}
}
}
return true;
}
......@@ -116,9 +116,9 @@ BSplineDecompositionImageFilter<TInputImage, TOutputImage>
::SetSplineOrder(unsigned int SplineOrder)
{
if (SplineOrder == m_SplineOrder)
{
{
return;
}
}
m_SplineOrder = SplineOrder;
this->SetPoles();
this->Modified();
......@@ -135,41 +135,41 @@ BSplineDecompositionImageFilter<TInputImage, TOutputImage>
// See also, Handbook of Medical Imaging, Processing and Analysis, Ed. Isaac N. Bankman,
// 2000, pg. 416.
switch (m_SplineOrder)
{
case 3:
m_NumberOfPoles = 1;
m_SplinePoles[0] = vcl_sqrt(3.0) - 2.0;
break;
case 0:
m_NumberOfPoles = 0;
break;
case 1:
m_NumberOfPoles = 0;
break;
case 2:
m_NumberOfPoles = 1;
m_SplinePoles[0] = vcl_sqrt(8.0) - 3.0;
break;
case 4:
m_NumberOfPoles = 2;
m_SplinePoles[0] = vcl_sqrt(664.0 - vcl_sqrt(438976.0)) + vcl_sqrt(304.0) - 19.0;
m_SplinePoles[1] = vcl_sqrt(664.0 + vcl_sqrt(438976.0)) - vcl_sqrt(304.0) - 19.0;
break;
case 5:
m_NumberOfPoles = 2;
m_SplinePoles[0] = vcl_sqrt(135.0 / 2.0 - vcl_sqrt(17745.0 / 4.0)) + vcl_sqrt(105.0 / 4.0)
- 13.0 / 2.0;
m_SplinePoles[1] = vcl_sqrt(135.0 / 2.0 + vcl_sqrt(17745.0 / 4.0)) - vcl_sqrt(105.0 / 4.0)
- 13.0 / 2.0;
break;
default:
// SplineOrder not implemented yet.
itk::ExceptionObject err(__FILE__, __LINE__);
err.SetLocation( ITK_LOCATION);
err.SetDescription( "SplineOrder must be between 0 and 5. Requested spline order has not been implemented yet." );
throw err;
break;
}
{
case 3:
m_NumberOfPoles = 1;
m_SplinePoles[0] = vcl_sqrt(3.0) - 2.0;
break;
case 0:
m_NumberOfPoles = 0;
break;
case 1:
m_NumberOfPoles = 0;
break;
case 2:
m_NumberOfPoles = 1;
m_SplinePoles[0] = vcl_sqrt(8.0) - 3.0;
break;
case 4:
m_NumberOfPoles = 2;
m_SplinePoles[0] = vcl_sqrt(664.0 - vcl_sqrt(438976.0)) + vcl_sqrt(304.0) - 19.0;
m_SplinePoles[1] = vcl_sqrt(664.0 + vcl_sqrt(438976.0)) - vcl_sqrt(304.0) - 19.0;
break;
case 5:
m_NumberOfPoles = 2;
m_SplinePoles[0] = vcl_sqrt(135.0 / 2.0 - vcl_sqrt(17745.0 / 4.0)) + vcl_sqrt(105.0 / 4.0)
- 13.0 / 2.0;
m_SplinePoles[1] = vcl_sqrt(135.0 / 2.0 + vcl_sqrt(17745.0 / 4.0)) - vcl_sqrt(105.0 / 4.0)
- 13.0 / 2.0;
break;
default:
// SplineOrder not implemented yet.
itk::ExceptionObject err(__FILE__, __LINE__);
err.SetLocation( ITK_LOCATION);
err.SetDescription( "SplineOrder must be between 0 and 5. Requested spline order has not been implemented yet." );
throw err;
break;
}
}
......@@ -188,33 +188,34 @@ BSplineDecompositionImageFilter<TInputImage, TOutputImage>
horizon = m_DataLength[m_IteratorDirection];
zn = z;
if (m_Tolerance > 0.0)
{
{
horizon = (long)vcl_ceil(log(m_Tolerance) / vcl_log(fabs(z)));
}
}
if (horizon < m_DataLength[m_IteratorDirection])
{
{
/* accelerated loop */
sum = m_Scratch[0]; // verify this
for (unsigned int n = 1; n < horizon; n++)
{
{
sum += zn * m_Scratch[n];
zn *= z;
}
m_Scratch[0] = sum;
}
else {
/* full loop */
iz = 1.0 / z;
z2n = vcl_pow(z, (double)(m_DataLength[m_IteratorDirection] - 1L));
sum = m_Scratch[0] + z2n * m_Scratch[m_DataLength[m_IteratorDirection] - 1L];
z2n *= z2n * iz;
for (unsigned int n = 1; n <= (m_DataLength[m_IteratorDirection] - 2); n++)
m_Scratch[0] = sum;
}
else
{
/* full loop */
iz = 1.0 / z;
z2n = vcl_pow(z, (double)(m_DataLength[m_IteratorDirection] - 1L));
sum = m_Scratch[0] + z2n * m_Scratch[m_DataLength[m_IteratorDirection] - 1L];
z2n *= z2n * iz;
for (unsigned int n = 1; n <= (m_DataLength[m_IteratorDirection] - 2); n++)
{
sum += (zn + z2n) * m_Scratch[n];
zn *= z;
z2n *= iz;
sum += (zn + z2n) * m_Scratch[n];
zn *= z;
z2n *= iz;
}
m_Scratch[0] = sum / (1.0 - zn * zn);
m_Scratch[0] = sum / (1.0 - zn * zn);
}
}
......@@ -250,7 +251,7 @@ BSplineDecompositionImageFilter<TInputImage, TOutputImage>
this->CopyImageToImage(); // Coefficients are initialized to the input data
for (unsigned int n=0; n < ImageDimension; n++)
{
{
m_IteratorDirection = n;
// Loop through each dimension
......@@ -259,7 +260,7 @@ BSplineDecompositionImageFilter<TInputImage, TOutputImage>
CIterator.SetDirection( m_IteratorDirection );
// For each data vector
while ( !CIterator.IsAtEnd() )
{
{
// Copy coefficients to scratch
this->CopyCoefficientsToScratch( CIterator );
......@@ -273,8 +274,8 @@ BSplineDecompositionImageFilter<TInputImage, TOutputImage>
this->CopyScratchToCoefficients( CIterator ); // m_Scratch = m_Image;
CIterator.NextLine();
progress.CompletedPixel();
}
}
}
}
......@@ -298,11 +299,11 @@ BSplineDecompositionImageFilter<TInputImage, TOutputImage>
outIt = outIt.Begin();
while ( !outIt.IsAtEnd() )
{
{
outIt.Set( static_cast<OutputPixelType>( inIt.Get() ) );
++inIt;
++outIt;
}
}
}
......@@ -318,11 +319,11 @@ BSplineDecompositionImageFilter<TInputImage, TOutputImage>
typedef typename TOutputImage::PixelType OutputPixelType;
unsigned long j = 0;
while ( !Iter.IsAtEndOfLine() )
{
{
Iter.Set( static_cast<OutputPixelType>( m_Scratch[j] ) );
++Iter;
++j;
}
}
}
......@@ -337,11 +338,11 @@ BSplineDecompositionImageFilter<TInputImage, TOutputImage>
{
unsigned long j = 0;
while ( !Iter.IsAtEndOfLine() )
{
{
m_Scratch[j] = static_cast<double>( Iter.Get() );
++Iter;
++j;
}
}
}
/**
......@@ -359,12 +360,12 @@ BSplineDecompositionImageFilter<TInputImage, TOutputImage>
unsigned long maxLength = 0;
for ( unsigned int n = 0; n < ImageDimension; n++ )
{
{
if ( m_DataLength[n] > maxLength )
{
{
maxLength = m_DataLength[n];
}
}
}
m_Scratch.resize( maxLength );
// Allocate memory for output image
......
Code/BasicFilters/otbBSplineInterpolateImageFunction.h
View file @ e9c0d74d
......@@ -41,11 +41,11 @@ namespace otb
* \ingroup ImageFunctions
*/
template <
class TImageType,
class TCoordRep = double,
class TCoefficientType = double >
class TImageType,
class TCoordRep = double,
class TCoefficientType = double >
class ITK_EXPORT BSplineInterpolateImageFunction :
public itk::InterpolateImageFunction<TImageType,TCoordRep>
public itk::InterpolateImageFunction<TImageType,TCoordRep>
{
public:
/** Standard class typedefs. */
......@@ -73,8 +73,8 @@ public:
/** Index typedef support. */
typedef typename Superclass::IndexType IndexType;
/** Region typedef support */
typedef typename InputImageType::RegionType RegionType;
/** Region typedef support */
typedef typename InputImageType::RegionType RegionType;
/** ContinuousIndex typedef support. */
typedef typename Superclass::ContinuousIndexType ContinuousIndexType;
......@@ -88,8 +88,8 @@ typedef typename InputImageType::RegionType RegionType;
/** Internal Coefficient typedef support */
typedef TCoefficientType CoefficientDataType;
typedef itk::Image<CoefficientDataType,
itkGetStaticConstMacro(ImageDimension)
> CoefficientImageType;
itkGetStaticConstMacro(ImageDimension)
> CoefficientImageType;
/** Define filter for calculating the BSpline coefficients */
typedef otb::BSplineDecompositionImageFilter<TImageType, CoefficientImageType>
......@@ -109,8 +109,8 @@ typedef typename InputImageType::RegionType RegionType;
/** Derivative typedef support */
typedef itk::CovariantVector<OutputType,
itkGetStaticConstMacro(ImageDimension)
> CovariantVectorType;
itkGetStaticConstMacro(ImageDimension)
> CovariantVectorType;
CovariantVectorType EvaluateDerivative( const PointType & point ) const
{
......@@ -135,7 +135,7 @@ typedef typename InputImageType::RegionType RegionType;
/** Update coefficients filter. Coefficient filter are computed over the buffered
region of the input image. */
virtual void UpdateCoefficientsFilter(void);
virtual void UpdateCoefficientsFilter(void);
protected:
BSplineInterpolateImageFunction();
......
Code/BasicFilters/otbBSplineInterpolateImageFunction.txx
View file @ e9c0d74d
......@@ -64,12 +64,12 @@ void
BSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
::UpdateCoefficientsFilter(void)
{
m_CoefficientFilter->GetOutput()->UpdateOutputInformation();
m_CoefficientFilter->GetOutput()->SetRequestedRegion(m_CoefficientFilter->GetInput()->GetBufferedRegion());
m_CoefficientFilter->GetOutput()->PropagateRequestedRegion();
m_CoefficientFilter->GetOutput()->UpdateOutputData();
m_Coefficients = m_CoefficientFilter->GetOutput();
m_CurrentBufferedRegion =m_CoefficientFilter->GetInput()->GetBufferedRegion();
m_CoefficientFilter->GetOutput()->UpdateOutputInformation();
m_CoefficientFilter->GetOutput()->SetRequestedRegion(m_CoefficientFilter->GetInput()->GetBufferedRegion());
m_CoefficientFilter->GetOutput()->PropagateRequestedRegion();
m_CoefficientFilter->GetOutput()->UpdateOutputData();
m_Coefficients = m_CoefficientFilter->GetOutput();
m_CurrentBufferedRegion =m_CoefficientFilter->GetInput()->GetBufferedRegion();
}
template <class TImageType, class TCoordRep, class TCoefficientType>
void
......@@ -77,8 +77,8 @@ BSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
::SetInputImage(const TImageType * inputData)
{
if ( inputData )
{
m_CoefficientFilter->SetInput(inputData);
{
m_CoefficientFilter->SetInput(inputData);
// the Coefficient Filter requires that the spline order and the input data be set.
// TODO: We need to ensure that this is only run once and only after both input and
......@@ -92,11 +92,11 @@ BSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
Superclass::SetInputImage(inputData);
m_DataLength = inputData->GetBufferedRegion().GetSize();
}
}
else
{
m_Coefficients = NULL;
}
{
m_Coefficients = NULL;
}
}
......@@ -106,18 +106,18 @@ BSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
::SetSplineOrder(unsigned int SplineOrder)
{
if (SplineOrder == m_SplineOrder)
{
{
return;
}
}
m_SplineOrder = SplineOrder;
m_CoefficientFilter->SetSplineOrder( SplineOrder );
//this->SetPoles();
m_MaxNumberInterpolationPoints = 1;
for (unsigned int n=0; n < ImageDimension; n++)
{
{
m_MaxNumberInterpolationPoints *= ( m_SplineOrder + 1);
}
}
this->GeneratePointsToIndex( );
}
......@@ -153,42 +153,42 @@ BSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
IndexType coefficientIndex;
// Step through eachpoint in the N-dimensional interpolation cube.
for (unsigned int p = 0; p < m_MaxNumberInterpolationPoints; p++)
{
{
// translate each step into the N-dimensional index.
// IndexType pointIndex = PointToIndex( p );
double w = 1.0;
for (unsigned int n = 0; n < ImageDimension; n++ )
{
{
w *= weights[n][ m_PointsToIndex[p][n] ];
coefficientIndex[n] = EvaluateIndex[n][m_PointsToIndex[p][n]]; // Build up ND index for coefficients.
//std::cout<<"From inside: "<<n<<" "<<p<<" "<<m_PointsToIndex[p][n]<<" "<< EvaluateIndex[n][m_PointsToIndex[p][n]]<<std::endl;
}
}
//std::cout<<"CoefficientIndex: "<<coefficientIndex<<std::endl;
// Convert our step p to the appropriate point in ND space in the
// m_Coefficients cube.
interpolated += w * m_Coefficients->GetPixel(coefficientIndex);
}
}
/* double interpolated = 0.0;
IndexType coefficientIndex;
// Step through eachpoint in the N-dimensional interpolation cube.
for (unsigned int sp = 0; sp <= m_SplineOrder; sp++)
{
for (unsigned int sp1=0; sp1 <= m_SplineOrder; sp1++)
/* double interpolated = 0.0;
IndexType coefficientIndex;
// Step through eachpoint in the N-dimensional interpolation cube.
for (unsigned int sp = 0; sp <= m_SplineOrder; sp++)
{
double w = 1.0;
for (unsigned int n1 = 0; n1 < ImageDimension; n1++ )
for (unsigned int sp1=0; sp1 <= m_SplineOrder; sp1++)
{
w *= weights[n1][ sp1 ];
coefficientIndex[n1] = EvaluateIndex[n1][sp]; // Build up ND index for coefficients.
}
interpolated += w * m_Coefficients->GetPixel(coefficientIndex);
double w = 1.0;
for (unsigned int n1 = 0; n1 < ImageDimension; n1++ )
{
w *= weights[n1][ sp1 ];
coefficientIndex[n1] = EvaluateIndex[n1][sp]; // Build up ND index for coefficients.
}
interpolated += w * m_Coefficients->GetPixel(coefficientIndex);
}
}
}
*/
*/
return(interpolated);
}
......@@ -223,30 +223,30 @@ BSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
double tempValue;
IndexType coefficientIndex;
for (unsigned int n = 0; n < ImageDimension; n++)
{
{
derivativeValue[n] = 0.0;
for (unsigned int p = 0; p < m_MaxNumberInterpolationPoints; p++)
{
{
tempValue = 1.0;
for (unsigned int n1 = 0; n1 < ImageDimension; n1++)
{
{
//coefficientIndex[n1] = EvaluateIndex[n1][sp];
coefficientIndex[n1] = EvaluateIndex[n1][m_PointsToIndex[p][n1]];
if (n1 == n)
{
{
//w *= weights[n][ m_PointsToIndex[p][n] ];
tempValue *= weightsDerivative[n1][ m_PointsToIndex[p][n1] ];
}
}
else
{
{
tempValue *= weights[n1][ m_PointsToIndex[p][n1] ];
}
}
derivativeValue[n] += m_Coefficients->GetPixel(coefficientIndex) * tempValue;
}
derivativeValue[n] /= this->GetInputImage()->GetSpacing()[n]; // take spacing into account
derivativeValue[n] += m_Coefficients->GetPixel(coefficientIndex) * tempValue;
}
derivativeValue[n] /= this->GetInputImage()->GetSpacing()[n]; // take spacing into account
}
return(derivativeValue);
......@@ -265,89 +265,89 @@ BSplineInterpolateImageFunction<TImageType,TCoordRep,TCoefficientType>
double w, w2, w4, t, t0, t1;
switch (splineOrder)
{
case 3:
for (unsigned int n = 0; n < ImageDimension; n++)
{
case 3:
for (unsigned int n = 0; n < ImageDimension; n++)
{
w = x[n] - (double) EvaluateIndex[n][1];
weights[n][3] = (1.0 / 6.0) * w * w * w;
weights[n][0] = (1.0 / 6.0) + 0.5 * w * (w - 1.0) - weights[n][3];
weights[n][2] = w + weights[n][0] - 2.0 * weights[n][3];
weights[n][1] = 1.0 - weights[n][0] - weights[n][2] - weights[n][3];
}
break;
case 0:
for (unsigned int n = 0; n < ImageDimension; n++)
{
weights[n][0] = 1; // implements nearest neighbor
}
break;
case 1:
for (unsigned int n = 0; n < ImageDimension; n++)
{
w = x[n] - (double) EvaluateIndex[n][0];
weights[n][1] = w;
weights[n][0] = 1.0 - w;
}
break;
case 2:
for (unsigned int n = 0; n < ImageDimension; n++)
{
/* x */
w = x[n] - (double)EvaluateIndex[n][1];
weights[n][1] = 0.75 - w * w;
weights[n][2] = 0.5 * (w - weights[n][1] + 1.0);
weights[n][0] = 1.0 - weights[n][1] - weights[n][2];
}
break;
case 4:
for (unsigned int n = 0; n < ImageDimension; n++)
{
/* x */
w = x[n] - (double)EvaluateIndex[n][2];
w2 = w * w;
t = (1.0 / 6.0) * w2;
weights[n][0] = 0.5 - w;
weights[n][0] *= weights[n][0];
weights[n][0] *= (1.0 / 24.0) * weights[n][0];
t0 = w * (t - 11.0 / 24.0);
t1 = 19.0 / 96.0 + w2 * (0.25 - t);
weights[n][1] = t1 + t0;
weights[n][3] = t1 - t0;
weights[n][4] = weights[n][0] + t0 + 0.5 * w;
weights[n][2] = 1.0 - weights[n][0] - weights[n][1] - weights[n][3] - weights[n][4];
}
break;
case 5:
for (unsigned int n = 0; n < ImageDimension; n++)
{
/* x */
w = x[n] - (double)EvaluateIndex[n][2];
w2 = w * w;
weights[n][5] = (1.0 / 120.0) * w * w2 * w2;
w2 -= w;
w4 = w2 * w2;
w -= 0.5;
t = w2 * (w2 - 3.0);
weights[n][0] = (1.0 / 24.0) * (1.0 / 5.0 + w2 + w4) - weights[n][5];
t0 = (1.0 / 24.0) * (w2 * (w2 - 5.0) + 46.0 / 5.0);
t1 = (-1.0 / 12.0) * w * (t + 4.0);
weights[n][2] = t0 + t1;
weights[n][3] = t0 - t1;
t0 = (1.0 / 16.0) * (9.0 / 5.0 - t);
t1 = (1.0 / 24.0) * w * (w4 - w2 - 5.0);
weights[n][1] = t0 + t1;
weights[n][4] = t0 - t1;
}
break;
default:
// SplineOrder not implemented yet.
itk::ExceptionObject err(__FILE__, __LINE__);
err.SetLocation( ITK_LOCATION );
err.SetDescription( "SplineOrder must be between 0 and 5. Requested spline order has not been implemented yet." );
throw err;