// ==========================================================
// Bitmap rotation using B-Splines
//
// Design and implementation by
// - Philippe Thévenaz (philippe.thevenaz@epfl.ch)
// Adaptation for FreeImage by
// - Hervé Drolon (drolon@infonie.fr)
//
// This file is part of FreeImage 3
//
// COVERED CODE IS PROVIDED UNDER THIS LICENSE ON AN "AS IS" BASIS, WITHOUT WARRANTY
// OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES
// THAT THE COVERED CODE IS FREE OF DEFECTS, MERCHANTABLE, FIT FOR A PARTICULAR PURPOSE
// OR NON-INFRINGING. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE COVERED
// CODE IS WITH YOU. SHOULD ANY COVERED CODE PROVE DEFECTIVE IN ANY RESPECT, YOU (NOT
// THE INITIAL DEVELOPER OR ANY OTHER CONTRIBUTOR) ASSUME THE COST OF ANY NECESSARY
// SERVICING, REPAIR OR CORRECTION. THIS DISCLAIMER OF WARRANTY CONSTITUTES AN ESSENTIAL
// PART OF THIS LICENSE. NO USE OF ANY COVERED CODE IS AUTHORIZED HEREUNDER EXCEPT UNDER
// THIS DISCLAIMER.
//
// Use at your own risk!
// ==========================================================
/* 
==========================================================
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
 @param Image Input image array
 @param y y coordinate of the selected line
 @param Line Output linear array
 @param Width Length of the line
*/
static void     
GetRow(double *Image, long y, double *Line, long Width) {
        long    x;
        Image = Image + (y * Width);
        for(x = 0L; x < Width; x++) {
                Line[x] = (double)*Image++;
        }
}
/**
 InitialAntiCausalCoefficient
 @param c Coefficients
 @param DataLength Number of samples or coefficients
 @param z Actual pole
 @return
*/
static double 
InitialAntiCausalCoefficient(double     *c, long DataLength, double     z) {
        // this initialization corresponds to mirror boundaries
        return((z / (z * z - 1.0)) * (z * c[DataLength - 2L] + c[DataLength - 1L]));
}
/**
 PutColumn
 @param Image Output image array
 @param Width Width of the image
 @param x x coordinate of the selected line
 @param Line Input linear array
 @param Height Length of the line and height of the image
*/
static void     
PutColumn(double *Image, long Width, long x, double *Line, long Height) {
        long    y;
        Image = Image + x;
        for(y = 0L; y < Height; y++) {
                *Image = (double)Line[y];
                Image += Width;
        }
}
/**
 PutRow
 @param Image Output image array
 @param y y coordinate of the selected line
 @param Line Input linear array
 @param Width length of the line and width of the image
*/
static void     
PutRow(double *Image, long y, double *Line, long Width) {
        long    x;
        Image = Image + (y * Width);
        for(x = 0L; x < Width; x++) {
                *Image++ = (double)Line[x];
        }
}
/**
 SamplesToCoefficients.<br>
 Implement the algorithm that converts the image samples into B-spline coefficients. 
 This efficient procedure essentially relies on the three papers cited above; 
 data are processed in-place. 
 Even though this algorithm is robust with respect to quantization, 
 we advocate the use of a floating-point format for the data. 
 @param Image Input / Output image (in-place processing)
 @param Width Width of the image
 @param Height Height of the image
 @param spline_degree Degree of the spline model
 @return Returns true if success, false otherwise
*/
static bool     
SamplesToCoefficients(double *Image, long Width, long Height, long spline_degree) {
        double  *Line;
        double  Pole[2];
        long    NbPoles;
        long    x, y;
        // recover the poles from a lookup table
        switch (spline_degree) {
                case 2L:
                        NbPoles = 1L;
                        Pole[0] = sqrt(8.0) - 3.0;
                        break;
                case 3L:
                        NbPoles = 1L;
                        Pole[0] = sqrt(3.0) - 2.0;
                        break;
                case 4L:
                        NbPoles = 2L;
                        Pole[0] = sqrt(664.0 - sqrt(438976.0)) + sqrt(304.0) - 19.0;
                        Pole[1] = sqrt(664.0 + sqrt(438976.0)) - sqrt(304.0) - 19.0;
                        break;
                case 5L:
                        NbPoles = 2L;
                        Pole[0] = sqrt(135.0 / 2.0 - sqrt(17745.0 / 4.0)) + sqrt(105.0 / 4.0)
                                - 13.0 / 2.0;
                        Pole[1] = sqrt(135.0 / 2.0 + sqrt(17745.0 / 4.0)) - sqrt(105.0 / 4.0)
                                - 13.0 / 2.0;
                        break;
                default:
                        // Invalid spline degree
                        return false;
        }
        // convert the image samples into interpolation coefficients 
        // in-place separable process, along x 
        Line = (double *)malloc(Width * sizeof(double));
        if (Line == NULL) {
                // Row allocation failed
                return false;
        }
        for (y = 0L; y < Height; y++) {
                GetRow(Image, y, Line, Width);
                ConvertToInterpolationCoefficients(Line, Width, Pole, NbPoles, DBL_EPSILON);
                PutRow(Image, y, Line, Width);
        }
        free(Line);
        // in-place separable process, along y 
        Line = (double *)malloc(Height * sizeof(double));
        if (Line == NULL) {
                // Column allocation failed
                return false;
        }
        for (x = 0L; x < Width; x++) {
                GetColumn(Image, Width, x, Line, Height);
                ConvertToInterpolationCoefficients(Line, Height, Pole, NbPoles, DBL_EPSILON);
                PutColumn(Image, Width, x, Line, Height);
        }
        free(Line);
        return true;
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Interpolation routines
/**
Perform the bidimensional interpolation of an image.
Given an array of spline coefficients, return the value of 
the underlying continuous spline model, sampled at the location (x, y). 
The model degree can be 2 (quadratic), 3 (cubic), 4 (quartic), or 5 (quintic).
@param Bcoeff Input B-spline array of coefficients
@param Width Width of the image
@param Height Height of the image
@param x x coordinate where to interpolate
@param y y coordinate where to interpolate
@param spline_degree Degree of the spline model
@return Returns the value of the underlying continuous spline model, 
sampled at the location (x, y)
*/
static double 
InterpolatedValue(double *Bcoeff, long Width, long Height, double x, double y, long spline_degree) {
        double  *p;
        double  xWeight[6], yWeight[6];
        double  interpolated;
        double  w, w2, w4, t, t0, t1;
        long    xIndex[6], yIndex[6];
        long    Width2 = 2L * Width - 2L, Height2 = 2L * Height - 2L;
        long    i, j, k;
        // compute the interpolation indexes
        if (spline_degree & 1L) {
                i = (long)floor(x) - spline_degree / 2L;
                j = (long)floor(y) - spline_degree / 2L;
                for(k = 0; k <= spline_degree; k++) {
                        xIndex[k] = i++;
                        yIndex[k] = j++;
                }
        }
        else {
                i = (long)floor(x + 0.5) - spline_degree / 2L;
                j = (long)floor(y + 0.5) - spline_degree / 2L;
                for (k = 0; k <= spline_degree; k++) {
                        xIndex[k] = i++;
                        yIndex[k] = j++;
                }
        }
        // compute the interpolation weights
        switch (spline_degree) {
                case 2L:
                        /* x */
                        w = x - (double)xIndex[1];
                        xWeight[1] = 3.0 / 4.0 - w * w;
                        xWeight[2] = (1.0 / 2.0) * (w - xWeight[1] + 1.0);
                        xWeight[0] = 1.0 - xWeight[1] - xWeight[2];
                        /* y */
                        w = y - (double)yIndex[1];
                        yWeight[1] = 3.0 / 4.0 - w * w;
                        yWeight[2] = (1.0 / 2.0) * (w - yWeight[1] + 1.0);
                        yWeight[0] = 1.0 - yWeight[1] - yWeight[2];
                        break;
                case 3L:
                        /* x */
                        w = x - (double)xIndex[1];
                        xWeight[3] = (1.0 / 6.0) * w * w * w;
                        xWeight[0] = (1.0 / 6.0) + (1.0 / 2.0) * w * (w - 1.0) - xWeight[3];
                        xWeight[2] = w + xWeight[0] - 2.0 * xWeight[3];
                        xWeight[1] = 1.0 - xWeight[0] - xWeight[2] - xWeight[3];
                        /* y */
                        w = y - (double)yIndex[1];
                        yWeight[3] = (1.0 / 6.0) * w * w * w;
                        yWeight[0] = (1.0 / 6.0) + (1.0 / 2.0) * w * (w - 1.0) - yWeight[3];
                        yWeight[2] = w + yWeight[0] - 2.0 * yWeight[3];
                        yWeight[1] = 1.0 - yWeight[0] - yWeight[2] - yWeight[3];
                        break;
                case 4L:
                        /* x */
                        w = x - (double)xIndex[2];
                        w2 = w * w;
                        t = (1.0 / 6.0) * w2;
                        xWeight[0] = 1.0 / 2.0 - w;
                        xWeight[0] *= xWeight[0];
                        xWeight[0] *= (1.0 / 24.0) * xWeight[0];
                        t0 = w * (t - 11.0 / 24.0);
                        t1 = 19.0 / 96.0 + w2 * (1.0 / 4.0 - t);
                        xWeight[1] = t1 + t0;
                        xWeight[3] = t1 - t0;
                        xWeight[4] = xWeight[0] + t0 + (1.0 / 2.0) * w;
                        xWeight[2] = 1.0 - xWeight[0] - xWeight[1] - xWeight[3] - xWeight[4];
                        /* y */
                        w = y - (double)yIndex[2];
                        w2 = w * w;
                        t = (1.0 / 6.0) * w2;
                        yWeight[0] = 1.0 / 2.0 - w;
                        yWeight[0] *= yWeight[0];
                        yWeight[0] *= (1.0 / 24.0) * yWeight[0];
                        t0 = w * (t - 11.0 / 24.0);
                        t1 = 19.0 / 96.0 + w2 * (1.0 / 4.0 - t);
                        yWeight[1] = t1 + t0;
                        yWeight[3] = t1 - t0;
                        yWeight[4] = yWeight[0] + t0 + (1.0 / 2.0) * w;
                        yWeight[2] = 1.0 - yWeight[0] - yWeight[1] - yWeight[3] - yWeight[4];
                        break;
                case 5L:
                        /* x */
                        w = x - (double)xIndex[2];
                        w2 = w * w;
                        xWeight[5] = (1.0 / 120.0) * w * w2 * w2;
                        w2 -= w;
                        w4 = w2 * w2;
                        w -= 1.0 / 2.0;
                        t = w2 * (w2 - 3.0);
                        xWeight[0] = (1.0 / 24.0) * (1.0 / 5.0 + w2 + w4) - xWeight[5];
                        t0 = (1.0 / 24.0) * (w2 * (w2 - 5.0) + 46.0 / 5.0);
                        t1 = (-1.0 / 12.0) * w * (t + 4.0);
                        xWeight[2] = t0 + t1;
                        xWeight[3] = t0 - t1;
                        t0 = (1.0 / 16.0) * (9.0 / 5.0 - t);
                        t1 = (1.0 / 24.0) * w * (w4 - w2 - 5.0);
                        xWeight[1] = t0 + t1;
                        xWeight[4] = t0 - t1;
                        /* y */
                        w = y - (double)yIndex[2];
                        w2 = w * w;
                        yWeight[5] = (1.0 / 120.0) * w * w2 * w2;
                        w2 -= w;
                        w4 = w2 * w2;
                        w -= 1.0 / 2.0;
                        t = w2 * (w2 - 3.0);
                        yWeight[0] = (1.0 / 24.0) * (1.0 / 5.0 + w2 + w4) - yWeight[5];
                        t0 = (1.0 / 24.0) * (w2 * (w2 - 5.0) + 46.0 / 5.0);
                        t1 = (-1.0 / 12.0) * w * (t + 4.0);
                        yWeight[2] = t0 + t1;
                        yWeight[3] = t0 - t1;
                        t0 = (1.0 / 16.0) * (9.0 / 5.0 - t);
                        t1 = (1.0 / 24.0) * w * (w4 - w2 - 5.0);
                        yWeight[1] = t0 + t1;
                        yWeight[4] = t0 - t1;
                        break;
                default:
                        // Invalid spline degree
                        return 0;
        }
        // apply the mirror boundary conditions
        for(k = 0; k <= spline_degree; k++) {
                xIndex[k] = (Width == 1L) ? (0L) : ((xIndex[k] < 0L) ?
                        (-xIndex[k] - Width2 * ((-xIndex[k]) / Width2))
                        : (xIndex[k] - Width2 * (xIndex[k] / Width2)));
                if (Width <= xIndex[k]) {
                        xIndex[k] = Width2 - xIndex[k];
                }
                yIndex[k] = (Height == 1L) ? (0L) : ((yIndex[k] < 0L) ?
                        (-yIndex[k] - Height2 * ((-yIndex[k]) / Height2))
                        : (yIndex[k] - Height2 * (yIndex[k] / Height2)));
                if (Height <= yIndex[k]) {
                        yIndex[k] = Height2 - yIndex[k];
                }
        }
        // perform interpolation
        interpolated = 0.0;
        for(j = 0; j <= spline_degree; j++) {
                p = Bcoeff + (yIndex[j] * Width);
                w = 0.0;
                for(i = 0; i <= spline_degree; i++) {
                        w += xWeight[i] * p[xIndex[i]];
                }
                interpolated += yWeight[j] * w;
        }
        return interpolated;
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
// FreeImage implementation
/** 
 Image translation and rotation using B-Splines.
 @param dib Input 8-bit greyscale image
 @param angle Output image rotation in degree
 @param x_shift Output image horizontal shift
 @param y_shift Output image vertical shift
 @param x_origin Output origin of the x-axis
 @param y_origin Output origin of the y-axis
 @param spline_degree Output degree of the B-spline model
 @param use_mask Whether or not to mask the image
 @return Returns the translated & rotated dib if successful, returns NULL otherwise
*/
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) {
        double  *ImageRasterArray;
        double  p;
        double  a11, a12, a21, a22;
        double  x0, y0, x1, y1;
        long    x, y;
        long    spline;
        bool    bResult;
        int bpp = FreeImage_GetBPP(dib);
        if(bpp != 8) {
                return NULL;
        }
         
        int width = FreeImage_GetWidth(dib);
        int height = FreeImage_GetHeight(dib);
        switch(spline_degree) {
                case ROTATE_QUADRATIC:
                        spline = 2L;    // Use splines of degree 2 (quadratic interpolation)
                        break;
                case ROTATE_CUBIC:
                        spline = 3L;    // Use splines of degree 3 (cubic interpolation)
                        break;
                case ROTATE_QUARTIC:
                        spline = 4L;    // Use splines of degree 4 (quartic interpolation)
                        break;
                case ROTATE_QUINTIC:
                        spline = 5L;    // Use splines of degree 5 (quintic interpolation)
                        break;
                default:
                        spline = 3L;
        }
        // allocate output image
        FIBITMAP *dst = FreeImage_Allocate(width, height, bpp);
        if(!dst)
                return NULL;
        // buid a grey scale palette
        RGBQUAD *pal = FreeImage_GetPalette(dst);
        for(int i = 0; i < 256; i++) {
                pal[i].rgbRed = pal[i].rgbGreen = pal[i].rgbBlue = (BYTE)i;
        }
        // allocate a temporary array
        ImageRasterArray = (double*)malloc(width * height * sizeof(double));
        if(!ImageRasterArray) {
                FreeImage_Unload(dst);
                return NULL;
        }
        // copy data samples
        for(y = 0; y < height; y++) {
                double *pImage = &ImageRasterArray[y*width];
                BYTE *src_bits = FreeImage_GetScanLine(dib, height-1-y);
                for(x = 0; x < width; x++) {
                        pImage[x] = (double)src_bits[x];
                }
        }
        // convert between a representation based on image samples
        // and a representation based on image B-spline coefficients
        bResult = SamplesToCoefficients(ImageRasterArray, width, height, spline);
        if(!bResult) {
                FreeImage_Unload(dst);
                free(ImageRasterArray);
                return NULL;
        }
        // prepare the geometry
        angle *= PI / 180.0;
        a11 = cos(angle);
        a12 = -sin(angle);
        a21 = sin(angle);
        a22 = cos(angle);
        x0 = a11 * (x_shift + x_origin) + a12 * (y_shift + y_origin);
        y0 = a21 * (x_shift + x_origin) + a22 * (y_shift + y_origin);
        x_shift = x_origin - x0;
        y_shift = y_origin - y0;
        // visit all pixels of the output image and assign their value
        for(y = 0; y < height; y++) {
                BYTE *dst_bits = FreeImage_GetScanLine(dst, height-1-y);
                 
                x0 = a12 * (double)y + x_shift;
                y0 = a22 * (double)y + y_shift;
                for(x = 0; x < width; x++) {
                        x1 = x0 + a11 * (double)x;
                        y1 = y0 + a21 * (double)x;
                        if(use_mask) {
                                if((x1 <= -0.5) || (((double)width - 0.5) <= x1) || (y1 <= -0.5) || (((double)height - 0.5) <= y1)) {
                                        p = 0;
                                }
                                else {
                                        p = (double)InterpolatedValue(ImageRasterArray, width, height, x1, y1, spline);
                                }
                        }
                        else {
                                p = (double)InterpolatedValue(ImageRasterArray, width, height, x1, y1, spline);
                        }
                        // clamp and convert to BYTE
                        dst_bits[x] = (BYTE)MIN(MAX((int)0, (int)(p + 0.5)), (int)255);
                }
        }
        // free working array and return
        free(ImageRasterArray);
        return dst;
}
/** 
 Image rotation using a 3rd order (cubic) B-Splines.
 @param dib Input dib (8, 24 or 32-bit)
 @param angle Output image rotation
 @param x_shift Output image horizontal shift
 @param y_shift Output image vertical shift
 @param x_origin Output origin of the x-axis
 @param y_origin Output origin of the y-axis
 @param use_mask Whether or not to mask the image
 @return Returns the translated & rotated dib if successful, returns NULL otherwise
*/
FIBITMAP * DLL_CALLCONV 
FreeImage_RotateEx(FIBITMAP *dib, double angle, double x_shift, double y_shift, double x_origin, double y_origin, BOOL use_mask) {
        int x, y, bpp;
        int channel, nb_channels;
        BYTE *src_bits, *dst_bits;
        FIBITMAP *src8 = NULL, *dst8 = NULL, *dst = NULL;
        if(!FreeImage_HasPixels(dib)) return NULL;
        try {
                bpp = FreeImage_GetBPP(dib);
                if(bpp == 8) {
                        FIBITMAP *dst_8 = Rotate8Bit(dib, angle, x_shift, y_shift, x_origin, y_origin, ROTATE_CUBIC, use_mask);
                        if(dst_8) {
                                // copy metadata from src to dst
                                FreeImage_CloneMetadata(dst_8, dib);
                        }
                        return dst_8;
                }
                if((bpp == 24) || (bpp == 32)) {
                        // allocate dst image
                        int width  = FreeImage_GetWidth(dib);
                        int height = FreeImage_GetHeight(dib);
                        if( bpp == 24 ) {
                                dst = FreeImage_Allocate(width, height, bpp, FI_RGBA_RED_MASK, FI_RGBA_GREEN_MASK, FI_RGBA_BLUE_MASK);
                        } else {
                                dst = FreeImage_Allocate(width, height, bpp, FI_RGBA_RED_MASK, FI_RGBA_GREEN_MASK, FI_RGBA_BLUE_MASK);
                        }
                        if(!dst) throw(1);
                        // allocate a temporary 8-bit dib (no need to build a palette)
                        src8 = FreeImage_Allocate(width, height, 8);
                        if(!src8) throw(1);
                        // process each channel separately
                        // -------------------------------
                        nb_channels = (bpp / 8);
                        for(channel = 0; channel < nb_channels; channel++) {
                                // extract channel from source dib
                                for(y = 0; y < height; y++) {
                                        src_bits = FreeImage_GetScanLine(dib, y);
                                        dst_bits = FreeImage_GetScanLine(src8, y);
                                        for(x = 0; x < width; x++) {
                                                dst_bits[x] = src_bits[channel];
                                                src_bits += nb_channels;
                                        }
                                }
                                // process channel
                                dst8 = Rotate8Bit(src8, angle, x_shift, y_shift, x_origin, y_origin, ROTATE_CUBIC, use_mask);
                                if(!dst8) throw(1);
                                // insert channel to destination dib
                                for(y = 0; y < height; y++) {
                                        src_bits = FreeImage_GetScanLine(dst8, y);
                                        dst_bits = FreeImage_GetScanLine(dst, y);
                                        for(x = 0; x < width; x++) {
                                                dst_bits[channel] = src_bits[x];
                                                dst_bits += nb_channels;
                                        }
                                }
                                FreeImage_Unload(dst8);
                        }
                        FreeImage_Unload(src8);
                        // copy metadata from src to dst
                        FreeImage_CloneMetadata(dst, dib);
                         
                        return dst;
                }
        } catch(int) {
                if(src8) FreeImage_Unload(src8);
                if(dst8) FreeImage_Unload(dst8);
                if(dst)  FreeImage_Unload(dst);
        }
        return NULL;
}