view release on metacpan or  search on metacpan

src/Source/FreeImageToolkit/BSplineRotate.cpp  view on Meta::CPAN

/* 
==========================================================
This code was taken and adapted from the following reference : 

[1] Philippe Thévenaz, Spline interpolation, a C source code 
implementation. http://bigwww.epfl.ch/thevenaz/

It implements ideas described in the following papers : 

[2] Unser M., Splines: A Perfect Fit for Signal and Image Processing. 
IEEE Signal Processing Magazine, vol. 16, no. 6, pp. 22-38, November 1999. 

[3] Unser M., Aldroubi A., Eden M., B-Spline Signal Processing: Part I--Theory.
IEEE Transactions on Signal Processing, vol. 41, no. 2, pp. 821-832, February 1993. 

[4] Unser M., Aldroubi A., Eden M., B-Spline Signal Processing: Part II--Efficient Design and Applications.
IEEE Transactions on Signal Processing, vol. 41, no. 2, pp. 834-848, February 1993.

========================================================== 
*/


#include <float.h>
#include "FreeImage.h"
#include "Utilities.h"

#define PI	((double)3.14159265358979323846264338327950288419716939937510)

#define ROTATE_QUADRATIC 2L	// Use B-splines of degree 2 (quadratic interpolation)
#define ROTATE_CUBIC     3L	// Use B-splines of degree 3 (cubic interpolation)
#define ROTATE_QUARTIC   4L	// Use B-splines of degree 4 (quartic interpolation)
#define ROTATE_QUINTIC   5L	// Use B-splines of degree 5 (quintic interpolation)


/////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Prototypes definition

static void ConvertToInterpolationCoefficients(double *c, long DataLength, double *z, long NbPoles,	double Tolerance);
static double InitialCausalCoefficient(double *c, long DataLength, double z, double Tolerance);
static void GetColumn(double *Image, long Width, long x, double *Line, long Height);
static void	GetRow(double *Image, long y, double *Line, long Width);
static double InitialAntiCausalCoefficient(double *c, long DataLength, double z);
static void	PutColumn(double *Image, long Width, long x, double *Line, long Height);
static void	PutRow(double *Image, long y, double *Line, long Width);
static bool SamplesToCoefficients(double *Image, long Width, long Height, long spline_degree);
static double InterpolatedValue(double *Bcoeff, long Width, long Height, double x, double y, long spline_degree);

static FIBITMAP * Rotate8Bit(FIBITMAP *dib, double angle, double x_shift, double y_shift, double x_origin, double y_origin, long spline_degree, BOOL use_mask);

/////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Coefficients routines

/**
 ConvertToInterpolationCoefficients

 @param c Input samples --> output coefficients
 @param DataLength Number of samples or coefficients
 @param z Poles
 @param NbPoles Number of poles
 @param Tolerance Admissible relative error
*/
static void 
ConvertToInterpolationCoefficients(double *c, long DataLength, double *z, long NbPoles,	double Tolerance) {
	double	Lambda = 1;
	long	n, k;

	// special case required by mirror boundaries
	if(DataLength == 1L) {
		return;
	}
	// compute the overall gain
	for(k = 0L; k < NbPoles; k++) {
		Lambda = Lambda * (1.0 - z[k]) * (1.0 - 1.0 / z[k]);
	}
	// apply the gain 
	for (n = 0L; n < DataLength; n++) {
		c[n] *= Lambda;
	}
	// loop over all poles 
	for (k = 0L; k < NbPoles; k++) {
		// causal initialization 
		c[0] = InitialCausalCoefficient(c, DataLength, z[k], Tolerance);
		// causal recursion 
		for (n = 1L; n < DataLength; n++) {
			c[n] += z[k] * c[n - 1L];
		}
		// anticausal initialization 
		c[DataLength - 1L] = InitialAntiCausalCoefficient(c, DataLength, z[k]);
		// anticausal recursion 
		for (n = DataLength - 2L; 0 <= n; n--) {
			c[n] = z[k] * (c[n + 1L] - c[n]);
		}
	}
} 

/**
 InitialCausalCoefficient

 @param c Coefficients
 @param DataLength Number of coefficients
 @param z Actual pole
 @param Tolerance Admissible relative error
 @return
*/
static double 
InitialCausalCoefficient(double	*c, long DataLength, double	z, double Tolerance) {
	double	Sum, zn, z2n, iz;
	long	n, Horizon;

	// this initialization corresponds to mirror boundaries 
	Horizon = DataLength;
	if(Tolerance > 0) {
		Horizon = (long)ceil(log(Tolerance) / log(fabs(z)));
	}
	if(Horizon < DataLength) {
		// accelerated loop
		zn = z;
		Sum = c[0];
		for (n = 1L; n < Horizon; n++) {
			Sum += zn * c[n];
			zn *= z;
		}
		return(Sum);
	}
	else {
		// full loop 
		zn = z;
		iz = 1.0 / z;
		z2n = pow(z, (double)(DataLength - 1L));
		Sum = c[0] + z2n * c[DataLength - 1L];
		z2n *= z2n * iz;
		for (n = 1L; n <= DataLength - 2L; n++) {
			Sum += (zn + z2n) * c[n];
			zn *= z;
			z2n *= iz;
		}
		return(Sum / (1.0 - zn * zn));
	}
}

/**
 GetColumn

 @param Image Input image array
 @param Width Width of the image
 @param x x coordinate of the selected line
 @param Line Output linear array
 @param Height Length of the line
*/
static void 
GetColumn(double *Image, long Width, long x, double *Line, long Height) {
	long y;

	Image = Image + x;
	for(y = 0L; y < Height; y++) {
		Line[y] = (double)*Image;
		Image += Width;
	}
}

/**
 GetRow