1 // ==========================================================
2 // Bitmap rotation using B-Splines
3 //
4 // Design and implementation by
5 // - Philippe Th�venaz (philippe.thevenaz@epfl.ch)
6 // Adaptation for FreeImage by
7 // - Herv� Drolon (drolon@infonie.fr)
8 //
9 // This file is part of FreeImage 3
10 //
11 // COVERED CODE IS PROVIDED UNDER THIS LICENSE ON AN "AS IS" BASIS, WITHOUT WARRANTY
12 // OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES
13 // THAT THE COVERED CODE IS FREE OF DEFECTS, MERCHANTABLE, FIT FOR A PARTICULAR PURPOSE
14 // OR NON-INFRINGING. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE COVERED
15 // CODE IS WITH YOU. SHOULD ANY COVERED CODE PROVE DEFECTIVE IN ANY RESPECT, YOU (NOT
16 // THE INITIAL DEVELOPER OR ANY OTHER CONTRIBUTOR) ASSUME THE COST OF ANY NECESSARY
17 // SERVICING, REPAIR OR CORRECTION. THIS DISCLAIMER OF WARRANTY CONSTITUTES AN ESSENTIAL
18 // PART OF THIS LICENSE. NO USE OF ANY COVERED CODE IS AUTHORIZED HEREUNDER EXCEPT UNDER
19 // THIS DISCLAIMER.
20 //
21 // Use at your own risk!
22 // ==========================================================
23 
24 /*
25 ==========================================================
26 This code was taken and adapted from the following reference :
27 
28 [1] Philippe Th�venaz, Spline interpolation, a C source code
29 implementation. http://bigwww.epfl.ch/thevenaz/
30 
31 It implements ideas described in the following papers :
32 
33 [2] Unser M., Splines: A Perfect Fit for Signal and Image Processing.
34 IEEE Signal Processing Magazine, vol. 16, no. 6, pp. 22-38, November 1999.
35 
36 [3] Unser M., Aldroubi A., Eden M., B-Spline Signal Processing: Part I--Theory.
37 IEEE Transactions on Signal Processing, vol. 41, no. 2, pp. 821-832, February 1993.
38 
39 [4] Unser M., Aldroubi A., Eden M., B-Spline Signal Processing: Part II--Efficient Design and Applications.
40 IEEE Transactions on Signal Processing, vol. 41, no. 2, pp. 834-848, February 1993.
41 
42 ==========================================================
43 */
44 
45 
46 #include <float.h>
47 #include "FreeImage.h"
48 #include "Utilities.h"
49 
50 #define PI	((double)3.14159265358979323846264338327950288419716939937510)
51 
52 #define ROTATE_QUADRATIC 2L	// Use B-splines of degree 2 (quadratic interpolation)
53 #define ROTATE_CUBIC     3L	// Use B-splines of degree 3 (cubic interpolation)
54 #define ROTATE_QUARTIC   4L	// Use B-splines of degree 4 (quartic interpolation)
55 #define ROTATE_QUINTIC   5L	// Use B-splines of degree 5 (quintic interpolation)
56 
57 
58 /////////////////////////////////////////////////////////////////////////////////////////////////////////////
59 // Prototypes definition
60 
61 static void ConvertToInterpolationCoefficients(double *c, long DataLength, double *z, long NbPoles,	double Tolerance);
62 static double InitialCausalCoefficient(double *c, long DataLength, double z, double Tolerance);
63 static void GetColumn(double *Image, long Width, long x, double *Line, long Height);
64 static void	GetRow(double *Image, long y, double *Line, long Width);
65 static double InitialAntiCausalCoefficient(double *c, long DataLength, double z);
66 static void	PutColumn(double *Image, long Width, long x, double *Line, long Height);
67 static void	PutRow(double *Image, long y, double *Line, long Width);
68 static bool SamplesToCoefficients(double *Image, long Width, long Height, long spline_degree);
69 static double InterpolatedValue(double *Bcoeff, long Width, long Height, double x, double y, long spline_degree);
70 
71 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);
72 
73 /////////////////////////////////////////////////////////////////////////////////////////////////////////////
74 // Coefficients routines
75 
76 /**
77  ConvertToInterpolationCoefficients
78 
79  @param c Input samples --> output coefficients
80  @param DataLength Number of samples or coefficients
81  @param z Poles
82  @param NbPoles Number of poles
83  @param Tolerance Admissible relative error
84 */
85 static void
ConvertToInterpolationCoefficients(double * c,long DataLength,double * z,long NbPoles,double Tolerance)86 ConvertToInterpolationCoefficients(double *c, long DataLength, double *z, long NbPoles,	double Tolerance) {
87 	double	Lambda = 1;
88 	long	n, k;
89 
90 	// special case required by mirror boundaries
91 	if(DataLength == 1L) {
92 		return;
93 	}
94 	// compute the overall gain
95 	for(k = 0L; k < NbPoles; k++) {
96 		Lambda = Lambda * (1.0 - z[k]) * (1.0 - 1.0 / z[k]);
97 	}
98 	// apply the gain
99 	for (n = 0L; n < DataLength; n++) {
100 		c[n] *= Lambda;
101 	}
102 	// loop over all poles
103 	for (k = 0L; k < NbPoles; k++) {
104 		// causal initialization
105 		c[0] = InitialCausalCoefficient(c, DataLength, z[k], Tolerance);
106 		// causal recursion
107 		for (n = 1L; n < DataLength; n++) {
108 			c[n] += z[k] * c[n - 1L];
109 		}
110 		// anticausal initialization
111 		c[DataLength - 1L] = InitialAntiCausalCoefficient(c, DataLength, z[k]);
112 		// anticausal recursion
113 		for (n = DataLength - 2L; 0 <= n; n--) {
114 			c[n] = z[k] * (c[n + 1L] - c[n]);
115 		}
116 	}
117 }
118 
119 /**
120  InitialCausalCoefficient
121 
122  @param c Coefficients
123  @param DataLength Number of coefficients
124  @param z Actual pole
125  @param Tolerance Admissible relative error
126  @return
127 */
128 static double
InitialCausalCoefficient(double * c,long DataLength,double z,double Tolerance)129 InitialCausalCoefficient(double	*c, long DataLength, double	z, double Tolerance) {
130 	double	Sum, zn, z2n, iz;
131 	long	n, Horizon;
132 
133 	// this initialization corresponds to mirror boundaries
134 	Horizon = DataLength;
135 	if(Tolerance > 0) {
136 		Horizon = (long)ceil(log(Tolerance) / log(fabs(z)));
137 	}
138 	if(Horizon < DataLength) {
139 		// accelerated loop
140 		zn = z;
141 		Sum = c[0];
142 		for (n = 1L; n < Horizon; n++) {
143 			Sum += zn * c[n];
144 			zn *= z;
145 		}
146 		return(Sum);
147 	}
148 	else {
149 		// full loop
150 		zn = z;
151 		iz = 1.0 / z;
152 		z2n = pow(z, (double)(DataLength - 1L));
153 		Sum = c[0] + z2n * c[DataLength - 1L];
154 		z2n *= z2n * iz;
155 		for (n = 1L; n <= DataLength - 2L; n++) {
156 			Sum += (zn + z2n) * c[n];
157 			zn *= z;
158 			z2n *= iz;
159 		}
160 		return(Sum / (1.0 - zn * zn));
161 	}
162 }
163 
164 /**
165  GetColumn
166 
167  @param Image Input image array
168  @param Width Width of the image
169  @param x x coordinate of the selected line
170  @param Line Output linear array
171  @param Height Length of the line
172 */
173 static void
GetColumn(double * Image,long Width,long x,double * Line,long Height)174 GetColumn(double *Image, long Width, long x, double *Line, long Height) {
175 	long y;
176 
177 	Image = Image + x;
178 	for(y = 0L; y < Height; y++) {
179 		Line[y] = (double)*Image;
180 		Image += Width;
181 	}
182 }
183 
184 /**
185  GetRow
186 
187  @param Image Input image array
188  @param y y coordinate of the selected line
189  @param Line Output linear array
190  @param Width Length of the line
191 */
192 static void
GetRow(double * Image,long y,double * Line,long Width)193 GetRow(double *Image, long y, double *Line, long Width) {
194 	long	x;
195 
196 	Image = Image + (y * Width);
197 	for(x = 0L; x < Width; x++) {
198 		Line[x] = (double)*Image++;
199 	}
200 }
201 
202 /**
203  InitialAntiCausalCoefficient
204 
205  @param c Coefficients
206  @param DataLength Number of samples or coefficients
207  @param z Actual pole
208  @return
209 */
210 static double
InitialAntiCausalCoefficient(double * c,long DataLength,double z)211 InitialAntiCausalCoefficient(double	*c, long DataLength, double	z) {
212 	// this initialization corresponds to mirror boundaries
213 	return((z / (z * z - 1.0)) * (z * c[DataLength - 2L] + c[DataLength - 1L]));
214 }
215 
216 /**
217  PutColumn
218 
219  @param Image Output image array
220  @param Width Width of the image
221  @param x x coordinate of the selected line
222  @param Line Input linear array
223  @param Height Length of the line and height of the image
224 */
225 static void
PutColumn(double * Image,long Width,long x,double * Line,long Height)226 PutColumn(double *Image, long Width, long x, double *Line, long Height) {
227 	long	y;
228 
229 	Image = Image + x;
230 	for(y = 0L; y < Height; y++) {
231 		*Image = (double)Line[y];
232 		Image += Width;
233 	}
234 }
235 
236 /**
237  PutRow
238 
239  @param Image Output image array
240  @param y y coordinate of the selected line
241  @param Line Input linear array
242  @param Width length of the line and width of the image
243 */
244 static void
PutRow(double * Image,long y,double * Line,long Width)245 PutRow(double *Image, long y, double *Line, long Width) {
246 	long	x;
247 
248 	Image = Image + (y * Width);
249 	for(x = 0L; x < Width; x++) {
250 		*Image++ = (double)Line[x];
251 	}
252 }
253 
254 /**
255  SamplesToCoefficients.<br>
256  Implement the algorithm that converts the image samples into B-spline coefficients.
257  This efficient procedure essentially relies on the three papers cited above;
258  data are processed in-place.
259  Even though this algorithm is robust with respect to quantization,
260  we advocate the use of a floating-point format for the data.
261 
262  @param Image Input / Output image (in-place processing)
263  @param Width Width of the image
264  @param Height Height of the image
265  @param spline_degree Degree of the spline model
266  @return Returns true if success, false otherwise
267 */
268 static bool
SamplesToCoefficients(double * Image,long Width,long Height,long spline_degree)269 SamplesToCoefficients(double *Image, long Width, long Height, long spline_degree) {
270 	double	*Line;
271 	double	Pole[2];
272 	long	NbPoles;
273 	long	x, y;
274 
275 	// recover the poles from a lookup table
276 	switch (spline_degree) {
277 		case 2L:
278 			NbPoles = 1L;
279 			Pole[0] = sqrt(8.0) - 3.0;
280 			break;
281 		case 3L:
282 			NbPoles = 1L;
283 			Pole[0] = sqrt(3.0) - 2.0;
284 			break;
285 		case 4L:
286 			NbPoles = 2L;
287 			Pole[0] = sqrt(664.0 - sqrt(438976.0)) + sqrt(304.0) - 19.0;
288 			Pole[1] = sqrt(664.0 + sqrt(438976.0)) - sqrt(304.0) - 19.0;
289 			break;
290 		case 5L:
291 			NbPoles = 2L;
292 			Pole[0] = sqrt(135.0 / 2.0 - sqrt(17745.0 / 4.0)) + sqrt(105.0 / 4.0)
293 				- 13.0 / 2.0;
294 			Pole[1] = sqrt(135.0 / 2.0 + sqrt(17745.0 / 4.0)) - sqrt(105.0 / 4.0)
295 				- 13.0 / 2.0;
296 			break;
297 		default:
298 			// Invalid spline degree
299 			return false;
300 	}
301 
302 	// convert the image samples into interpolation coefficients
303 
304 	// in-place separable process, along x
305 	Line = (double *)malloc(Width * sizeof(double));
306 	if (Line == NULL) {
307 		// Row allocation failed
308 		return false;
309 	}
310 	for (y = 0L; y < Height; y++) {
311 		GetRow(Image, y, Line, Width);
312 		ConvertToInterpolationCoefficients(Line, Width, Pole, NbPoles, DBL_EPSILON);
313 		PutRow(Image, y, Line, Width);
314 	}
315 	free(Line);
316 
317 	// in-place separable process, along y
318 	Line = (double *)malloc(Height * sizeof(double));
319 	if (Line == NULL) {
320 		// Column allocation failed
321 		return false;
322 	}
323 	for (x = 0L; x < Width; x++) {
324 		GetColumn(Image, Width, x, Line, Height);
325 		ConvertToInterpolationCoefficients(Line, Height, Pole, NbPoles, DBL_EPSILON);
326 		PutColumn(Image, Width, x, Line, Height);
327 	}
328 	free(Line);
329 
330 	return true;
331 }
332 
333 /////////////////////////////////////////////////////////////////////////////////////////////////////////////
334 // Interpolation routines
335 
336 /**
337 Perform the bidimensional interpolation of an image.
338 Given an array of spline coefficients, return the value of
339 the underlying continuous spline model, sampled at the location (x, y).
340 The model degree can be 2 (quadratic), 3 (cubic), 4 (quartic), or 5 (quintic).
341 
342 @param Bcoeff Input B-spline array of coefficients
343 @param Width Width of the image
344 @param Height Height of the image
345 @param x x coordinate where to interpolate
346 @param y y coordinate where to interpolate
347 @param spline_degree Degree of the spline model
348 @return Returns the value of the underlying continuous spline model,
349 sampled at the location (x, y)
350 */
351 static double
InterpolatedValue(double * Bcoeff,long Width,long Height,double x,double y,long spline_degree)352 InterpolatedValue(double *Bcoeff, long Width, long Height, double x, double y, long spline_degree) {
353 	double	*p;
354 	double	xWeight[6], yWeight[6];
355 	double	interpolated;
356 	double	w, w2, w4, t, t0, t1;
357 	long	xIndex[6], yIndex[6];
358 	long	Width2 = 2L * Width - 2L, Height2 = 2L * Height - 2L;
359 	long	i, j, k;
360 
361 	// compute the interpolation indexes
362 	if (spline_degree & 1L) {
363 		i = (long)floor(x) - spline_degree / 2L;
364 		j = (long)floor(y) - spline_degree / 2L;
365 		for(k = 0; k <= spline_degree; k++) {
366 			xIndex[k] = i++;
367 			yIndex[k] = j++;
368 		}
369 	}
370 	else {
371 		i = (long)floor(x + 0.5) - spline_degree / 2L;
372 		j = (long)floor(y + 0.5) - spline_degree / 2L;
373 		for (k = 0; k <= spline_degree; k++) {
374 			xIndex[k] = i++;
375 			yIndex[k] = j++;
376 		}
377 	}
378 
379 	// compute the interpolation weights
380 	switch (spline_degree) {
381 		case 2L:
382 			/* x */
383 			w = x - (double)xIndex[1];
384 			xWeight[1] = 3.0 / 4.0 - w * w;
385 			xWeight[2] = (1.0 / 2.0) * (w - xWeight[1] + 1.0);
386 			xWeight[0] = 1.0 - xWeight[1] - xWeight[2];
387 			/* y */
388 			w = y - (double)yIndex[1];
389 			yWeight[1] = 3.0 / 4.0 - w * w;
390 			yWeight[2] = (1.0 / 2.0) * (w - yWeight[1] + 1.0);
391 			yWeight[0] = 1.0 - yWeight[1] - yWeight[2];
392 			break;
393 		case 3L:
394 			/* x */
395 			w = x - (double)xIndex[1];
396 			xWeight[3] = (1.0 / 6.0) * w * w * w;
397 			xWeight[0] = (1.0 / 6.0) + (1.0 / 2.0) * w * (w - 1.0) - xWeight[3];
398 			xWeight[2] = w + xWeight[0] - 2.0 * xWeight[3];
399 			xWeight[1] = 1.0 - xWeight[0] - xWeight[2] - xWeight[3];
400 			/* y */
401 			w = y - (double)yIndex[1];
402 			yWeight[3] = (1.0 / 6.0) * w * w * w;
403 			yWeight[0] = (1.0 / 6.0) + (1.0 / 2.0) * w * (w - 1.0) - yWeight[3];
404 			yWeight[2] = w + yWeight[0] - 2.0 * yWeight[3];
405 			yWeight[1] = 1.0 - yWeight[0] - yWeight[2] - yWeight[3];
406 			break;
407 		case 4L:
408 			/* x */
409 			w = x - (double)xIndex[2];
410 			w2 = w * w;
411 			t = (1.0 / 6.0) * w2;
412 			xWeight[0] = 1.0 / 2.0 - w;
413 			xWeight[0] *= xWeight[0];
414 			xWeight[0] *= (1.0 / 24.0) * xWeight[0];
415 			t0 = w * (t - 11.0 / 24.0);
416 			t1 = 19.0 / 96.0 + w2 * (1.0 / 4.0 - t);
417 			xWeight[1] = t1 + t0;
418 			xWeight[3] = t1 - t0;
419 			xWeight[4] = xWeight[0] + t0 + (1.0 / 2.0) * w;
420 			xWeight[2] = 1.0 - xWeight[0] - xWeight[1] - xWeight[3] - xWeight[4];
421 			/* y */
422 			w = y - (double)yIndex[2];
423 			w2 = w * w;
424 			t = (1.0 / 6.0) * w2;
425 			yWeight[0] = 1.0 / 2.0 - w;
426 			yWeight[0] *= yWeight[0];
427 			yWeight[0] *= (1.0 / 24.0) * yWeight[0];
428 			t0 = w * (t - 11.0 / 24.0);
429 			t1 = 19.0 / 96.0 + w2 * (1.0 / 4.0 - t);
430 			yWeight[1] = t1 + t0;
431 			yWeight[3] = t1 - t0;
432 			yWeight[4] = yWeight[0] + t0 + (1.0 / 2.0) * w;
433 			yWeight[2] = 1.0 - yWeight[0] - yWeight[1] - yWeight[3] - yWeight[4];
434 			break;
435 		case 5L:
436 			/* x */
437 			w = x - (double)xIndex[2];
438 			w2 = w * w;
439 			xWeight[5] = (1.0 / 120.0) * w * w2 * w2;
440 			w2 -= w;
441 			w4 = w2 * w2;
442 			w -= 1.0 / 2.0;
443 			t = w2 * (w2 - 3.0);
444 			xWeight[0] = (1.0 / 24.0) * (1.0 / 5.0 + w2 + w4) - xWeight[5];
445 			t0 = (1.0 / 24.0) * (w2 * (w2 - 5.0) + 46.0 / 5.0);
446 			t1 = (-1.0 / 12.0) * w * (t + 4.0);
447 			xWeight[2] = t0 + t1;
448 			xWeight[3] = t0 - t1;
449 			t0 = (1.0 / 16.0) * (9.0 / 5.0 - t);
450 			t1 = (1.0 / 24.0) * w * (w4 - w2 - 5.0);
451 			xWeight[1] = t0 + t1;
452 			xWeight[4] = t0 - t1;
453 			/* y */
454 			w = y - (double)yIndex[2];
455 			w2 = w * w;
456 			yWeight[5] = (1.0 / 120.0) * w * w2 * w2;
457 			w2 -= w;
458 			w4 = w2 * w2;
459 			w -= 1.0 / 2.0;
460 			t = w2 * (w2 - 3.0);
461 			yWeight[0] = (1.0 / 24.0) * (1.0 / 5.0 + w2 + w4) - yWeight[5];
462 			t0 = (1.0 / 24.0) * (w2 * (w2 - 5.0) + 46.0 / 5.0);
463 			t1 = (-1.0 / 12.0) * w * (t + 4.0);
464 			yWeight[2] = t0 + t1;
465 			yWeight[3] = t0 - t1;
466 			t0 = (1.0 / 16.0) * (9.0 / 5.0 - t);
467 			t1 = (1.0 / 24.0) * w * (w4 - w2 - 5.0);
468 			yWeight[1] = t0 + t1;
469 			yWeight[4] = t0 - t1;
470 			break;
471 		default:
472 			// Invalid spline degree
473 			return 0;
474 	}
475 
476 	// apply the mirror boundary conditions
477 	for(k = 0; k <= spline_degree; k++) {
478 		xIndex[k] = (Width == 1L) ? (0L) : ((xIndex[k] < 0L) ?
479 			(-xIndex[k] - Width2 * ((-xIndex[k]) / Width2))
480 			: (xIndex[k] - Width2 * (xIndex[k] / Width2)));
481 		if (Width <= xIndex[k]) {
482 			xIndex[k] = Width2 - xIndex[k];
483 		}
484 		yIndex[k] = (Height == 1L) ? (0L) : ((yIndex[k] < 0L) ?
485 			(-yIndex[k] - Height2 * ((-yIndex[k]) / Height2))
486 			: (yIndex[k] - Height2 * (yIndex[k] / Height2)));
487 		if (Height <= yIndex[k]) {
488 			yIndex[k] = Height2 - yIndex[k];
489 		}
490 	}
491 
492 	// perform interpolation
493 	interpolated = 0.0;
494 	for(j = 0; j <= spline_degree; j++) {
495 		p = Bcoeff + (yIndex[j] * Width);
496 		w = 0.0;
497 		for(i = 0; i <= spline_degree; i++) {
498 			w += xWeight[i] * p[xIndex[i]];
499 		}
500 		interpolated += yWeight[j] * w;
501 	}
502 
503 	return interpolated;
504 }
505 
506 /////////////////////////////////////////////////////////////////////////////////////////////////////////////
507 // FreeImage implementation
508 
509 
510 /**
511  Image translation and rotation using B-Splines.
512 
513  @param dib Input 8-bit greyscale image
514  @param angle Output image rotation in degree
515  @param x_shift Output image horizontal shift
516  @param y_shift Output image vertical shift
517  @param x_origin Output origin of the x-axis
518  @param y_origin Output origin of the y-axis
519  @param spline_degree Output degree of the B-spline model
520  @param use_mask Whether or not to mask the image
521  @return Returns the translated & rotated dib if successful, returns NULL otherwise
522 */
523 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)524 Rotate8Bit(FIBITMAP *dib, double angle, double x_shift, double y_shift, double x_origin, double y_origin, long spline_degree, BOOL use_mask) {
525 	double	*ImageRasterArray;
526 	double	p;
527 	double	a11, a12, a21, a22;
528 	double	x0, y0, x1, y1;
529 	long	x, y;
530 	long	spline;
531 	bool	bResult;
532 
533 	int bpp = FreeImage_GetBPP(dib);
534 	if(bpp != 8) {
535 		return NULL;
536 	}
537 
538 	int width = FreeImage_GetWidth(dib);
539 	int height = FreeImage_GetHeight(dib);
540 	switch(spline_degree) {
541 		case ROTATE_QUADRATIC:
542 			spline = 2L;	// Use splines of degree 2 (quadratic interpolation)
543 			break;
544 		case ROTATE_CUBIC:
545 			spline = 3L;	// Use splines of degree 3 (cubic interpolation)
546 			break;
547 		case ROTATE_QUARTIC:
548 			spline = 4L;	// Use splines of degree 4 (quartic interpolation)
549 			break;
550 		case ROTATE_QUINTIC:
551 			spline = 5L;	// Use splines of degree 5 (quintic interpolation)
552 			break;
553 		default:
554 			spline = 3L;
555 	}
556 
557 	// allocate output image
558 	FIBITMAP *dst = FreeImage_Allocate(width, height, bpp);
559 	if(!dst)
560 		return NULL;
561 	// buid a grey scale palette
562 	RGBQUAD *pal = FreeImage_GetPalette(dst);
563 	for(int i = 0; i < 256; i++) {
564 		pal[i].rgbRed = pal[i].rgbGreen = pal[i].rgbBlue = (BYTE)i;
565 	}
566 
567 	// allocate a temporary array
568 	ImageRasterArray = (double*)malloc(width * height * sizeof(double));
569 	if(!ImageRasterArray) {
570 		FreeImage_Unload(dst);
571 		return NULL;
572 	}
573 	// copy data samples
574 	for(y = 0; y < height; y++) {
575 		double *pImage = &ImageRasterArray[y*width];
576 		BYTE *src_bits = FreeImage_GetScanLine(dib, height-1-y);
577 
578 		for(x = 0; x < width; x++) {
579 			pImage[x] = (double)src_bits[x];
580 		}
581 	}
582 
583 	// convert between a representation based on image samples
584 	// and a representation based on image B-spline coefficients
585 	bResult = SamplesToCoefficients(ImageRasterArray, width, height, spline);
586 	if(!bResult) {
587 		FreeImage_Unload(dst);
588 		free(ImageRasterArray);
589 		return NULL;
590 	}
591 
592 	// prepare the geometry
593 	angle *= PI / 180.0;
594 	a11 = cos(angle);
595 	a12 = -sin(angle);
596 	a21 = sin(angle);
597 	a22 = cos(angle);
598 	x0 = a11 * (x_shift + x_origin) + a12 * (y_shift + y_origin);
599 	y0 = a21 * (x_shift + x_origin) + a22 * (y_shift + y_origin);
600 	x_shift = x_origin - x0;
601 	y_shift = y_origin - y0;
602 
603 	// visit all pixels of the output image and assign their value
604 	for(y = 0; y < height; y++) {
605 		BYTE *dst_bits = FreeImage_GetScanLine(dst, height-1-y);
606 
607 		x0 = a12 * (double)y + x_shift;
608 		y0 = a22 * (double)y + y_shift;
609 
610 		for(x = 0; x < width; x++) {
611 			x1 = x0 + a11 * (double)x;
612 			y1 = y0 + a21 * (double)x;
613 			if(use_mask) {
614 				if((x1 <= -0.5) || (((double)width - 0.5) <= x1) || (y1 <= -0.5) || (((double)height - 0.5) <= y1)) {
615 					p = 0;
616 				}
617 				else {
618 					p = (double)InterpolatedValue(ImageRasterArray, width, height, x1, y1, spline);
619 				}
620 			}
621 			else {
622 				p = (double)InterpolatedValue(ImageRasterArray, width, height, x1, y1, spline);
623 			}
624 			// clamp and convert to BYTE
625 			dst_bits[x] = (BYTE)MIN(MAX((int)0, (int)(p + 0.5)), (int)255);
626 		}
627 	}
628 
629 	// free working array and return
630 	free(ImageRasterArray);
631 
632 	return dst;
633 }
634 
635 /**
636  Image rotation using a 3rd order (cubic) B-Splines.
637 
638  @param dib Input dib (8, 24 or 32-bit)
639  @param angle Output image rotation
640  @param x_shift Output image horizontal shift
641  @param y_shift Output image vertical shift
642  @param x_origin Output origin of the x-axis
643  @param y_origin Output origin of the y-axis
644  @param use_mask Whether or not to mask the image
645  @return Returns the translated & rotated dib if successful, returns NULL otherwise
646 */
647 FIBITMAP * DLL_CALLCONV
FreeImage_RotateEx(FIBITMAP * dib,double angle,double x_shift,double y_shift,double x_origin,double y_origin,BOOL use_mask)648 FreeImage_RotateEx(FIBITMAP *dib, double angle, double x_shift, double y_shift, double x_origin, double y_origin, BOOL use_mask) {
649 
650 	int x, y, bpp;
651 	int channel, nb_channels;
652 	BYTE *src_bits, *dst_bits;
653 	FIBITMAP *src8 = NULL, *dst8 = NULL, *dst = NULL;
654 
655 	if(!FreeImage_HasPixels(dib)) return NULL;
656 
657 	try {
658 
659 		bpp = FreeImage_GetBPP(dib);
660 
661 		if(bpp == 8) {
662 			FIBITMAP *dst_8 = Rotate8Bit(dib, angle, x_shift, y_shift, x_origin, y_origin, ROTATE_CUBIC, use_mask);
663 			if(dst_8) {
664 				// copy metadata from src to dst
665 				FreeImage_CloneMetadata(dst_8, dib);
666 			}
667 			return dst_8;
668 		}
669 		if((bpp == 24) || (bpp == 32)) {
670 			// allocate dst image
671 			int width  = FreeImage_GetWidth(dib);
672 			int height = FreeImage_GetHeight(dib);
673 			if( bpp == 24 ) {
674 				dst = FreeImage_Allocate(width, height, bpp, FI_RGBA_RED_MASK, FI_RGBA_GREEN_MASK, FI_RGBA_BLUE_MASK);
675 			} else {
676 				dst = FreeImage_Allocate(width, height, bpp, FI_RGBA_RED_MASK, FI_RGBA_GREEN_MASK, FI_RGBA_BLUE_MASK);
677 			}
678 			if(!dst) throw(1);
679 
680 			// allocate a temporary 8-bit dib (no need to build a palette)
681 			src8 = FreeImage_Allocate(width, height, 8);
682 			if(!src8) throw(1);
683 
684 			// process each channel separately
685 			// -------------------------------
686 			nb_channels = (bpp / 8);
687 
688 			for(channel = 0; channel < nb_channels; channel++) {
689 				// extract channel from source dib
690 				for(y = 0; y < height; y++) {
691 					src_bits = FreeImage_GetScanLine(dib, y);
692 					dst_bits = FreeImage_GetScanLine(src8, y);
693 					for(x = 0; x < width; x++) {
694 						dst_bits[x] = src_bits[channel];
695 						src_bits += nb_channels;
696 					}
697 				}
698 
699 				// process channel
700 				dst8 = Rotate8Bit(src8, angle, x_shift, y_shift, x_origin, y_origin, ROTATE_CUBIC, use_mask);
701 				if(!dst8) throw(1);
702 
703 				// insert channel to destination dib
704 				for(y = 0; y < height; y++) {
705 					src_bits = FreeImage_GetScanLine(dst8, y);
706 					dst_bits = FreeImage_GetScanLine(dst, y);
707 					for(x = 0; x < width; x++) {
708 						dst_bits[channel] = src_bits[x];
709 						dst_bits += nb_channels;
710 					}
711 				}
712 
713 				FreeImage_Unload(dst8);
714 			}
715 
716 			FreeImage_Unload(src8);
717 
718 			// copy metadata from src to dst
719 			FreeImage_CloneMetadata(dst, dib);
720 
721 			return dst;
722 		}
723 	} catch(int) {
724 		if(src8) FreeImage_Unload(src8);
725 		if(dst8) FreeImage_Unload(dst8);
726 		if(dst)  FreeImage_Unload(dst);
727 	}
728 
729 	return NULL;
730 }
731