/*
 * interpolatedOps.c --
 *
 *      In this file are implemented those mage transformations 
 *      require interpolation. (scale and rotate)
 *
 * Copyright (c) 2016-2017 by Simon Bachmann <simonbachmann@bluewin.ch>
 * All rights reserved.
 *
 * Permission is hereby granted, without written agreement and without
 * license or royalty fees, to use, copy, modify, and distribute this
 * software and its documentation for any purpose, provided that the
 * above copyright notice and the following two paragraphs appear in
 * all copies of this software.
 * 
 */

#include <tcl.h>
#include <tk.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
#include <limits.h>

#ifdef _WINDOWS
   #include <windows.h>
#endif

//#include "config.h"
#include "imgtoolsDecl.h"
#include "catromcache.h"
#include "lanczos2cache.h"
#include "lanczos3cache.h"
#include "interpolatedOps.h"
#include "cputime.h"
#include "debugmsg.h"

#include <assert.h>

/* 
 * === Structures and Constants ===
 */

/*
 * Cache of weights and source coords for scale interpol function
 */ 
typedef struct {
        float *weights;     /* array of weights, pointsPerSample long */
                            /* may be in any order, but correspond to srcPos indices */
        float weightsSum;   /* sum of weights */
        int *srcPos;        /* indices of source pixels, size of array: pointsPerSample */
                            /* may be in any order, but correspond to weights */
        int pointsPerSample;/* number of points to use for interpolation */
        float posInSource;  /* position of dest pixel in source */
} Imgtools_ScaleCache;



/* 
 * === Prototypes ===
 */
static int Imgtools_ResizeImage(Tcl_Interp *interp, const Imgtools_Image *sourceData,
    const Imgtools_CommonOpts *options, Imgtools_Image *destData);
static int Imgtools_CalcScaleCache(Tcl_Interp *interp, int srcSize, int destSize,
    enum FilterType filter, Imgtools_ScaleCache *cachePtr);
static int Imgtools_FreeScaleCache(Tcl_Interp *interp, Imgtools_ScaleCache *cachePtr, int size);

static int Imgtools_RotateImage(Tcl_Interp *interp, const Imgtools_Image *sourceData, 
    double angle, const Imgtools_CommonOpts *options, const unsigned char outsidePixel[4], 
    Imgtools_Image *destData);
static int Imgtools_RotateImage90(Tcl_Interp *interp, const Imgtools_Image *sourceData, 
    double angle, const Imgtools_CommonOpts *options, const unsigned char outsidePixel[4], 
    Imgtools_Image *destData);

static int Imgtools_Interpolate(Tcl_Interp *interp, const Imgtools_Image *srourceData, 
    const Imgtools_CommonOpts *options, const unsigned char outsidePixel[4], 
    const Imgtools_FloatCoords *destPos, unsigned char *destPixelPtr);
static int Imgtools_NearestInterpol(Tcl_Interp *interp, const Imgtools_Image *srourceData, 
    const Imgtools_CommonOpts *options, const unsigned char outsidePixel[4], 
    const Imgtools_FloatCoords *destPos, unsigned char *destPixelPtr);

static float Imgtools_GeneralLanczosFilter (const float x, const int a);
static float Imgtools_Lanczos2Filter (const float x);
static float Imgtools_Lanczos3Filter (const float x);
static float Imgtools_CachedLanczos2Filter (const float x);
static float Imgtools_CachedLanczos3Filter (const float x);
static float Imgtools_CatromFilter (const float x);
static float Imgtools_LinearFilter (const float x);
static float Imgtools_CachedCatromFilter (const float x);
static float Imgtools_NearestFilter (const float x);


/*
 * Filters
 */

//filter width: delta of the region, where the filter has nonzero values
#ifdef IMGTOOLS_EXACT
const int filterWidth[FILTER_COUNT] = {2, 2, 4, 4, 6, 4, 4, 6};
#else
const int filterWidth[FILTER_COUNT] = {2, 2, 4, 4, 6};
#endif

#define MAX_FILTERWIDTH 6

/*
 * function pointer array for general interpol function
 */
int (*const Imgtools_InterpolFuncts[FILTER_COUNT])(Tcl_Interp *interp, const Imgtools_Image *, const Imgtools_CommonOpts *, const unsigned char *, const Imgtools_FloatCoords *, unsigned char *)= {
    Imgtools_NearestInterpol,
    Imgtools_Interpolate,   /* linear */
    Imgtools_Interpolate,   /* catrom */
    Imgtools_Interpolate,   /* lanczos2 */
    Imgtools_Interpolate    /* lanczos3 */
#ifdef IMGTOOLS_EXACT
    ,Imgtools_Interpolate,  /* catrom exact */
    Imgtools_Interpolate,   /* lanczos2 exact */
    Imgtools_Interpolate    /* lanczos3 exact */
#endif
};

/*
 * function pointer array for filter functions
 * (delta in, weight out)
 */
float (*const Imgtools_FilterFuncts[FILTER_COUNT])(const float) = {
    Imgtools_NearestFilter,
    Imgtools_LinearFilter,
    Imgtools_CachedCatromFilter,
    Imgtools_CachedLanczos2Filter,
    Imgtools_CachedLanczos3Filter
#ifdef IMGTOOLS_EXACT
    ,Imgtools_CatromFilter,
    Imgtools_Lanczos2Filter,
    Imgtools_Lanczos3Filter
#endif
};

/*
 * Filter names
 */
const char *const filterNames[] = {
    "nearest",
    "linear",
    "catrom",
    "lanczos2",
    "lanczos3",
#ifdef IMGTOOLS_EXACT
    "catromExact",
    "lanczos2Exact",
    "lanczos3Exact",
#endif
    NULL
};

/*
 * Settings for cputime
 */

enum Scopes {
    MAIN,
    CONVERT,
    ROTATE,
    ROTATE_LOOP,
    INTERPOL,
    SCALE,
    CACHE,
    SCALE_LOOP
};

const char *const scopeNameList[] = { 
    "Main function (option parsing, etc.)", 
    "Converting ops",
    "Rotate, pre-calculate and stuff",
    "Rotate main loop",
    "Interpolation",
    "Scale, pre-calculate and stuff",
    "Scale, create cache",
    "Scale, main loop",
    NULL 
};

/*
 * === forward Declarations ====
 */

/* 
 * Holds cached values for filters
 */
extern const float *const catromCache;
extern const float *const lanczos2Cache;
extern const float *const lanczos3Cache;

/*
 * === Scaling Images ===
 */

/*
 *--------------------------------------------------------------
 *
 * Imgtools_Scale --
 *
 *      Main Function for scaling images. Called by the user via
 *      the corresponding command (::imgtools::scale)
 *
 * Results:
 *      Standard Tcl result.
 *
 * Side effects:
 *      writes a modified image (according to options) to destImg. 
 *      destImg must exist, its contents will be partially or  
 *      completely overwritten. If the -width or -height options 
 *      for the dest image were set, these sizes are respected.
 *
 *--------------------------------------------------------------
 */
int Imgtools_Scale(clientData, interp, objc, objv)
    ClientData clientData;      /* not used */
    Tcl_Interp *interp;         /* current interpreter*/
    int objc;                   /* argument count */
    Tcl_Obj *const objv[];      /* arguments. See manpage for details */
{
    Cputime_SetScopes(interp, scopeNameList);
    Cputime_NewScope(interp, MAIN);
    
    int result = TCL_OK;
    
    Imgtools_Image srcData, destData;
    destData.memOwner = 0;  /* make sure we can safely call Imgtools_FreeImage() on these*/
    srcData.memOwner = 0;   /* before  they're initialized */
                                    
    /*
     * Parse arguments
     */
    Imgtools_CommonOpts commonOpts;
    int allowedOpts;
    
    Tk_PhotoHandle srcimg, destimg;
    Tcl_Obj *destImgName = NULL;
    
    if (objc < 4) {
        Tcl_WrongNumArgs(interp, 1, objv, "sourceImg sizespec ?options? destImg");
        Tcl_SetErrorCode(interp, "IMGTOOLS", "OPTIONS", "WRONG_NUM_ARGS", NULL);
        result = TCL_ERROR;
        goto tidyUpAndLeave;
    }
    if ((result = Imgtools_FindPhoto(interp, objv[1], &srcimg)) !=TCL_OK) {
        goto tidyUpAndLeave;
    }
    if ((result = Imgtools_FindPhoto(interp, objv[objc-1], &destimg)) != TCL_OK) {
        goto tidyUpAndLeave;
    }
    destImgName = objv[objc-1];

    /*
     * Parse the rest of options first, 
     * sizespec comes later
     */
    allowedOpts = OPT_COMPOSITE | OPT_SHRINK | OPT_FILTER | OPT_FROM | OPT_TO;
#ifdef IMGTOOLS_DEBUG
    allowedOpts |= OPT_FOLLOWPIXEL;
    int debugLevelSave = 0;
#endif

    if ((result = Imgtools_ParseCommonOpts(interp, objc-4, objv+3, allowedOpts, 
            &commonOpts)) != TCL_OK) {
        goto tidyUpAndLeave;
    }

    /*
     * Get the image data, adapt to -from option and calculate sizes
     */
    Cputime_NewScope(interp, CONVERT);
    if ((result = Imgtools_GetImage(interp, srcimg, &srcData)) != TCL_OK) {
        goto tidyUpAndLeave;
    }
    Cputime_EndScopeSecure(interp, CONVERT);
    
    if ((result = Imgtools_ApplySubregion(interp, commonOpts.fromCorners, &srcData)) != TCL_OK) 
    {
        goto tidyUpAndLeave;
    }
    
    Imgtools_IntCoords sourceSize, targetSize = {-1, -1};
    sourceSize.x = srcData.width;
    sourceSize.y = srcData.height;
    if (Imgtools_GetSizesFromSpec(interp, objv[2], sourceSize, &targetSize) != TCL_OK) {
        result = TCL_ERROR;
        goto tidyUpAndLeave;
    }
    
    if (targetSize.x < 0 || targetSize.y < 0) {
        /* 
         * This should never happen, as Imgtools_GetSizesFromSpec 
         * already checks for negative sizes
         */
        assert(0);
        targetSize.x = 0;
        targetSize.y = 0;
    }
    
#ifdef IMGTOOLS_DEBUG
    /*
     * If we got -followpixel, set appropriate debuglevel
     */
    debugLevelSave = Imgtools_GetDebugLevel(interp);
    if (commonOpts.options & OPT_FOLLOWPIXEL) {
        Imgtools_SetDebugLevel(interp, DEBUG_PARAMETERS | DEBUG_FOLLOW);
        debugmsg(interp, DEBUG_FOLLOW, "Data for dest pixel %d, %d\n", 
                 commonOpts.followPixel.x, commonOpts.followPixel.y);
        debugmsg(interp, DEBUG_FOLLOW, "=================================\n");
    }
#endif
    
    /*
     * We are ready to start the actual work!
     */
    Imgtools_Image *destDataPtr;
    destDataPtr = &destData;
    destData.width = targetSize.x;
    destData.height = targetSize.y;
    if (destData.width == 0 || destData.height == 0) {
        destData.pixelSize = 4;
        destData.pitch = destData.width * destData.pixelSize;
        destData.pixelPtr = NULL;
        destData.startPtr = NULL;
        destData.memOwner = 0;
    } else if (srcData.width == 0 || srcData.height == 0) {
        destDataPtr = &srcData;
    } else if (destData.width == srcData.width && destData.height == srcData.height) {
        destDataPtr = &srcData;
    } else {
        Imgtools_ResizeImage(interp, &srcData, &commonOpts, &destData);
    }
        
    /*
     * Put data to destImg and clean things up
     */
    Cputime_NewScope(interp, CONVERT);
    Imgtools_IntCoords start, end, size;
    size.x = destDataPtr->width;
    size.y = destDataPtr->height;
    Imgtools_GetRegion(interp, size, commonOpts.toCorners, ALLOW_LARGER,
        SIZE_IS_INCREMENT, &start, &end);
    Imgtools_GetSizeFromCorners(start, end, &size);
    if (Imgtools_PutImage(interp, commonOpts.composite, commonOpts.shrink, destDataPtr,
            start, size, Tcl_GetString(destImgName), destimg) != TCL_OK) {
        Tcl_SetObjResult(interp, Tcl_NewStringObj("Failed to put data to destination image.", -1));
        Tcl_SetErrorCode(interp, "IMGTOOLS", "DATA", "EXPORT_FAILED", NULL);
        result = TCL_ERROR;
        goto tidyUpAndLeave;
    }
    Cputime_EndScopeSecure(interp, CONVERT);
    
    Tcl_SetObjResult(interp, destImgName);
    result = TCL_OK;
    
#ifdef IMGTOOLS_DEBUG
    /*
     * Reset previous debuglevel settings
     */
    Imgtools_SetDebugLevel(interp, debugLevelSave);
#endif
    
    
tidyUpAndLeave:
    
    Imgtools_FreeImage(&destData);
    Imgtools_FreeImage(&srcData);
    
    Cputime_EndScopeSecure(interp, MAIN);
    Cputime_PrintSummary(interp);
    
    return result;
}

/*
 *--------------------------------------------------------------
 *
 * Imgtools_ResizeImage --
 *
 *      Resizes an image
 *      destData must have width and height set to desired values.
 *      the other fields in destData are ignored and overwritten with
 *      appropriate values.
 * 
 *      This function assumes that the filter is 0 at its endpoints 
 *      (e.g. -2 and 2 for catrom)
 *
 * Results:
 *      Standard Tcl Result. Writes the scaled Image to the destData struct.
 *
 * Side effects:
 *      Allocates needed Amount of memory for pixelPtr in dest Data.
 *      Remember to free.
 *
 *--------------------------------------------------------------
 */

static int Imgtools_ResizeImage(interp, srcData, options, destData)
Tcl_Interp *interp;                 /* used for debugmsg */
const Imgtools_Image *srcData;  /* The starting point*/
const Imgtools_CommonOpts *options; /* Filter, followpixel coords */
Imgtools_Image *destData;       /* Result is written here */
{
    /*
     * Initialize destData
     */
    destData->pixelSize = 4;
    destData->pitch = destData->width * 4;
    destData->pixelPtr = (unsigned char *)Tcl_Alloc(sizeof(unsigned char) 
        * destData->width * destData->height * 4 );
    destData->startPtr = destData->pixelPtr;
    destData->memOwner = 1;
            
    /*
     * Create a cache of all weights and weight sums for rows an columns
     */
    Cputime_NewScope(interp, CACHE);
    
    Imgtools_ScaleCache *xCache, *yCache;
    
    xCache = (Imgtools_ScaleCache *)Tcl_Alloc(sizeof(Imgtools_ScaleCache) * destData->width);
    yCache = (Imgtools_ScaleCache *)Tcl_Alloc(sizeof(Imgtools_ScaleCache) * destData->height);
    
    debugmsg(interp, DEBUG_PARAMETERS, "Parameters for horizontal (x) scaling:\n");
    Imgtools_CalcScaleCache(interp, srcData->width, destData->width, options->filter, xCache);
    debugmsg(interp, DEBUG_PARAMETERS, "Parameters for vertical (y) scaling:\n");
    Imgtools_CalcScaleCache(interp, srcData->height, destData->height, options->filter, yCache);
    
    Cputime_EndScopeSecure(interp, CACHE);
    Cputime_NewScope(interp, SCALE_LOOP);
        
    /*
     * Resample Rows
     */
    float *intermResults, *curIntermPtr;
    float curWeightAlpha;
    unsigned char *curPixelPtr;
    
    intermResults = (float*)Tcl_Alloc(sizeof(float) * destData->width * srcData->height * 4);
    memset(intermResults, 0, sizeof(float) * destData->width * srcData->height * 4);
    
#ifdef IMGTOOLS_DEBUG
    /*
     * Followpixel: calculate wich source rows to consider
     */
    int leftmostSrcCol, rightmostSrcCol;
    leftmostSrcCol = ceil(yCache[options->followPixel.y].posInSource 
                    - yCache[options->followPixel.y].pointsPerSample/2.0);
    rightmostSrcCol = leftmostSrcCol + yCache[options->followPixel.y].pointsPerSample -1;

#endif
    
    for (int y = srcData->height; y--; ) {
        for (int x = destData->width; x--; ) {
            curIntermPtr = intermResults + y*destData->pitch + x*destData->pixelSize;
#ifdef IMGTOOLS_DEBUG
            if (options->options & OPT_FOLLOWPIXEL && x == options->followPixel.x && 
                y >= leftmostSrcCol && y <= rightmostSrcCol
            ) {
                debugmsg(interp, DEBUG_FOLLOW, "Interpolating row %d\n", y);
                debugmsg(interp, DEBUG_FOLLOW, "Position in source (x): %f\n", xCache[x].posInSource); 
            }
#endif
            for (int p = xCache[x].pointsPerSample; p--; ) {
                if (xCache[x].srcPos[p] < 0 || xCache[x].srcPos[p] >= srcData->width) {
                    continue;
                }
                curPixelPtr = srcData->startPtr + y*srcData->pitch + xCache[x].srcPos[p]*srcData->pixelSize; 
                curWeightAlpha = xCache[x].weights[p] * curPixelPtr[ALPHA];
                curIntermPtr[RED] += curPixelPtr[RED] * curWeightAlpha;
                curIntermPtr[GREEN] += curPixelPtr[GREEN] * curWeightAlpha;
                curIntermPtr[BLUE] += curPixelPtr[BLUE] * curWeightAlpha;
                curIntermPtr[ALPHA] += curWeightAlpha;
#ifdef IMGTOOLS_DEBUG
                if (options->options & OPT_FOLLOWPIXEL && x == options->followPixel.x && 
                    y >= leftmostSrcCol && y <= rightmostSrcCol
                ) {
                    debugmsg(interp, DEBUG_FOLLOW, "\tSource pixel %d, %d\n", xCache[x].srcPos[p], y);
                    debugmsg(interp, DEBUG_FOLLOW, "\t\tValues: %d %d %d %d\n",
                             curPixelPtr[RED], curPixelPtr[GREEN], curPixelPtr[BLUE], curPixelPtr[ALPHA]);
                    debugmsg(interp, DEBUG_FOLLOW, "\t\tAssociated values: %d %d %d %d\n",
                             curPixelPtr[RED]*curPixelPtr[ALPHA], curPixelPtr[GREEN]*curPixelPtr[ALPHA],
                             curPixelPtr[BLUE]*curPixelPtr[ALPHA], curPixelPtr[ALPHA]);
                    debugmsg(interp, DEBUG_FOLLOW, "\t\tWeighted values: %f %f %f %f\n",
                             curPixelPtr[RED]*curWeightAlpha, curPixelPtr[GREEN]*curWeightAlpha,
                             curPixelPtr[BLUE]*curWeightAlpha, curWeightAlpha);
                    debugmsg(interp, DEBUG_FOLLOW, "\t\tWeight: %f\n", 
                             xCache[x].weights[p]);
                }
#endif
            }
#ifdef IMGTOOLS_DEBUG
            if (options->options & OPT_FOLLOWPIXEL && x == options->followPixel.x && 
                y >= leftmostSrcCol && y <= rightmostSrcCol
            ) {
                debugmsg(interp, DEBUG_FOLLOW, "\tWeighted sums: %f %f %f %f\n",
                         curIntermPtr[RED], curIntermPtr[GREEN], curIntermPtr[BLUE], curIntermPtr[ALPHA]);
            }
#endif
        }
    }
    
    /*
     * Resample Cols
     */
    float wSums[4], tmp;
    for (int x = destData->width; x--; ) {
        for (int y = destData->height; y--; ) {
            curPixelPtr = destData->startPtr + y*destData->pitch + x*destData->pixelSize;
            memset(wSums, 0, sizeof(float)*4);
#ifdef IMGTOOLS_DEBUG
            if (options->options & OPT_FOLLOWPIXEL && x == options->followPixel.x && 
                y == options->followPixel.y
            ) {
                debugmsg(interp, DEBUG_FOLLOW, "Interpolating dest pixel %d, %d\n", x, y);
                debugmsg(interp, DEBUG_FOLLOW, "Position in source (y): %f\n", yCache[y].posInSource);
            }
#endif
            
            for (int p = yCache[y].pointsPerSample; p--; ) {
                if (yCache[y].srcPos[p] < 0 || yCache[y].srcPos[p] >= srcData->height) {
                    continue;
                }
                curIntermPtr = intermResults + x*destData->pixelSize + yCache[y].srcPos[p]*destData->pitch;
                wSums[RED] += curIntermPtr[RED] * yCache[y].weights[p];
                wSums[GREEN] += curIntermPtr[GREEN] * yCache[y].weights[p];
                wSums[BLUE] += curIntermPtr[BLUE] * yCache[y].weights[p];
                wSums[ALPHA] += curIntermPtr[ALPHA] * yCache[y].weights[p];
#ifdef IMGTOOLS_DEBUG
                if (options->options & OPT_FOLLOWPIXEL 
                    && x == options->followPixel.x 
                    && y == options->followPixel.y
                ) {
                    debugmsg(interp, DEBUG_FOLLOW, "\tIntermediate pixel %d, %d\n", x, yCache[y].srcPos[p]);
                    debugmsg(interp, DEBUG_FOLLOW, "\t\tValues: %f %f %f %f\n", curIntermPtr[RED], 
                             curIntermPtr[GREEN], curIntermPtr[BLUE], curIntermPtr[ALPHA]);
                    debugmsg(interp, DEBUG_FOLLOW, "\t\tWeighted values: %f %f %f %f\n",
                             curIntermPtr[RED]*yCache[y].weights[p], curIntermPtr[GREEN]*yCache[y].weights[p], 
                             curIntermPtr[BLUE]*yCache[y].weights[p], curIntermPtr[ALPHA]*yCache[y].weights[p]);
                    debugmsg(interp, DEBUG_FOLLOW, "\t\tWeight: %f\n", 
                             yCache[x].weights[p]);
                }
#endif
            }
            /*
             * Divide by weights, round, clip and write to destData
             */
#ifdef IMGTOOLS_DEBUG
            if (options->options & OPT_FOLLOWPIXEL 
                && x == options->followPixel.x 
                && y == options->followPixel.y
            ) {
                debugmsg(interp, DEBUG_FOLLOW, "Raw weighted sums: %f %f %f %f\n", 
                         wSums[RED], wSums[GREEN], wSums[BLUE], wSums[ALPHA]);
            }
#endif
            if (wSums[ALPHA] == 0) {
                wSums[RED] = 0;
                wSums[GREEN] = 0;
                wSums[BLUE] = 0;
            } else {
                wSums[RED] /= wSums[ALPHA];
                wSums[GREEN] /= wSums[ALPHA];
                wSums[BLUE] /= wSums[ALPHA];
            }
            if ((yCache[y].weightsSum * xCache[x].weightsSum) == 0) {
                assert(0);
                wSums[ALPHA] = 0;
            } else {
                wSums[ALPHA] /= (yCache[y].weightsSum * xCache[x].weightsSum);
            }
            tmp = roundf(wSums[RED]);
            tmp = tmp < 0 ? 0 : tmp;
            tmp = tmp > UCHAR_MAX ? UCHAR_MAX : tmp;
            curPixelPtr[RED] = tmp;
            
            tmp = roundf(wSums[GREEN]);
            tmp = tmp < 0 ? 0 : tmp;
            tmp = tmp > UCHAR_MAX ? UCHAR_MAX : tmp;
            curPixelPtr[GREEN] = tmp;
            
            tmp = roundf(wSums[BLUE]);
            tmp = tmp < 0 ? 0 : tmp;
            tmp = tmp > UCHAR_MAX ? UCHAR_MAX : tmp;
            curPixelPtr[BLUE] = tmp;
            
            tmp = roundf(wSums[ALPHA]);
            tmp = tmp < 0 ? 0 : tmp;
            tmp = tmp > UCHAR_MAX ? UCHAR_MAX : tmp;
            curPixelPtr[ALPHA] = tmp;
#ifdef IMGTOOLS_DEBUG
            if (options->options & OPT_FOLLOWPIXEL 
                && x == options->followPixel.x 
                && y == options->followPixel.y
            ) {
                debugmsg(interp, DEBUG_FOLLOW, "Interpolated values (float): %f %f %f %f\n",
                         wSums[RED], wSums[GREEN], wSums[BLUE], wSums[ALPHA]);
                debugmsg(interp, DEBUG_FOLLOW, "Rounded and clipped values: %d %d %d %d\n\n",
                         curPixelPtr[RED], curPixelPtr[GREEN], curPixelPtr[BLUE], curPixelPtr[ALPHA]);
            }
#endif
        }
    }

    /*
     * Cleanup
     */
    Imgtools_FreeScaleCache(interp, xCache, destData->width);
    Imgtools_FreeScaleCache(interp, yCache, destData->height);
    Tcl_Free((char *)xCache);
    Tcl_Free((char *)yCache);
    Tcl_Free((char *)intermResults);
    
    Cputime_EndScopeSecure(interp, SCALE_LOOP);
    
    return TCL_OK;
}

/*
 *--------------------------------------------------------------
 *
 * Imgtools_CalcScaleCache --
 *
 *      Calculate weights for one direction for scaling interpolation
 *
 * Results:
 *      Standard Tcl result (always TCL_OK, actually).
 *
 * Side effects:
 *      Writes the weights and weightsSum to corresponding arrays
 *
 *--------------------------------------------------------------
 */
static int Imgtools_CalcScaleCache (interp, srcSize, destSize, filter, cachePtr)
Tcl_Interp *interp;             /* used for debugmsg */
int srcSize, destSize;          /* sizes (width or height) of src and dest images */
enum FilterType filter;         /* filter to use */
Imgtools_ScaleCache *cachePtr;  /* cache */
{
    /*
     * Calculate scale factors and points per sample
     */
    int pointsPerSample;        /* Numbers of source pixels to consider for one sample */
    double scaleFactor;         /* scale ratio */
    int leftmostSourcePx;       /* index of leftmost source pixel in sample */
    double destPosInSource;     /* position of dest pixel in source coodrs */
    double srcToFilterTransf;       /* coefficient to unse to transform a source distance to filter distance */
    double scaleAcoeff, scaleBcoeff;/*coefficients of the linear function used to calculate */
                                    /* relative pos of destination in source */
    
    scaleFactor = (double)destSize / srcSize;
    if (scaleFactor < 1) {
        pointsPerSample = ceil(filterWidth[filter] / scaleFactor);
        srcToFilterTransf = scaleFactor;
    } else {
        pointsPerSample = filterWidth[filter];
        srcToFilterTransf = 1;
    }
    scaleAcoeff = 1 / scaleFactor;
    scaleBcoeff = 1 / (2*scaleFactor) - 0.5;
    
    debugmsg(interp, DEBUG_PARAMETERS, "\tfilter: %s (%d), filterWidth: %d\n", filterNames[filter], filter, filterWidth[filter]);
    debugmsg(interp, DEBUG_PARAMETERS, "\tsrcSize: %d, destSize: %d, scaleFactor: %f\n", srcSize, destSize, scaleFactor);
    debugmsg(interp, DEBUG_PARAMETERS, "\tpointsPerSample: %d, \n", pointsPerSample);
    debugmsg(interp, DEBUG_PARAMETERS, "\tscale coefficients: a=%f, b=%f\n", scaleAcoeff, scaleBcoeff);
    
    /*
     * Alloc, initialize and calculate cache
     * To free the memory allocated here use cachePtr[0].weights/srcPos
     */
    float *weightsChunkPtr;
    int *srcPosChunkPtr;
    float distance;
    int smallestIndex = -1;
    float smallestDistance = -1; /* used for nearest interpol */
    weightsChunkPtr = (float *)Tcl_Alloc(sizeof(float) * destSize * pointsPerSample);
    srcPosChunkPtr = (int *)Tcl_Alloc(sizeof(int) * destSize * pointsPerSample);
    
    /*
     * Calculate weights and srcPos
     */    
    for (int i = 0; i < destSize; i++) {
        destPosInSource = scaleAcoeff * i + scaleBcoeff;
        leftmostSourcePx = ceil(destPosInSource - pointsPerSample/2.0);
        
        cachePtr[i].weights = weightsChunkPtr;
        weightsChunkPtr += pointsPerSample;
        cachePtr[i].srcPos = srcPosChunkPtr;
        srcPosChunkPtr += pointsPerSample;
        cachePtr[i].pointsPerSample = pointsPerSample;
        cachePtr[i].weightsSum = 0;
        cachePtr[i].posInSource = destPosInSource;
        smallestDistance = -1;
        smallestIndex = -1;
         
        debugmsg(interp, DEBUG_ALL, "i: %d\n", i);
        debugmsg(interp, DEBUG_ALL, "DestPosInSource: %f, leftmostSourcePx: %d\n", destPosInSource, leftmostSourcePx);
        debugmsg(interp, DEBUG_ALL, "distances (inverse order !):\n");
        
        /* 
         * The special condition needed for nearest neighbor "interpolation" 
         * here is kinda ugly. But could not come up with a better way...
         * The impact on performace should be negligible
         */
        for (int j = pointsPerSample; j--; ) {
            cachePtr[i].srcPos[j] = leftmostSourcePx + j;
            distance = (cachePtr[i].srcPos[j] - destPosInSource)*srcToFilterTransf;
            
            debugmsg(interp, DEBUG_ALL, "%6.3f ", distance);
            
            if (filter == NEAREST_NEIGHBOR) {
                cachePtr[i].weights[j] = 0;
                if (smallestDistance < 0 || fabsf(distance) < smallestDistance) {
                    smallestDistance = fabsf(distance);
                    smallestIndex = j;
                }
            } else if (leftmostSourcePx+j < 0 || leftmostSourcePx+j >= srcSize) {
                cachePtr[i].weights[j] = 0;
            } else {
                cachePtr[i].weights[j] = (*Imgtools_FilterFuncts[filter])(distance);
                cachePtr[i].weightsSum += cachePtr[i].weights[j];
            }
        }
        if (filter == NEAREST_NEIGHBOR) {
            assert(smallestIndex >= 0);
            cachePtr[i].weights[smallestIndex] = 1;
            cachePtr[i].weightsSum = 1;
        }

        debugmsg(interp, DEBUG_ALL, "\n\n");
    }
    
    /*
     * print a table with weights, distances and weight sums
     */
#ifdef IMGTOOLS_DEBUG
    debugmsg(interp, DEBUG_MORE, "weights and weightSums:\n");
    for(int i = 0; i < destSize; i++) {
        debugmsg(interp, DEBUG_MORE, "%2d: ", i);
        for (int j = 0; j < pointsPerSample; j++) {
            debugmsg(interp, DEBUG_MORE, "%6.3f ", cachePtr[i].weights[j]);
        }
        debugmsg(interp, DEBUG_MORE, " | %6.3f (sum)\n", cachePtr[i].weightsSum);
    }
    debugmsg(interp, DEBUG_MORE, "\n");
    debugmsg(interp, DEBUG_MORE, "source positions:\n");
    for(int i = 0; i < destSize; i++) {
        debugmsg(interp, DEBUG_MORE, "%d: ", i);
        for (int j = 0; j < pointsPerSample; j++) {
            debugmsg(interp, DEBUG_MORE, "%6d ", cachePtr[i].srcPos[j]);
        }
        debugmsg(interp, DEBUG_MORE, "\n");
    }
    debugmsg(interp, DEBUG_MORE, "\n");
#endif

    return TCL_OK;
}

/*
 *--------------------------------------------------------------
 *
 * Imgtools_FreeScaleCache --
 *
 *      Free the memory allocated for a scale cache
 *
 * Results:
 *      Standard Tcl result (always TCL_OK, actually).
 *
 * Side effects:
 *      Frees the memory allocated for scale cache
 *
 *--------------------------------------------------------------
 */
static int Imgtools_FreeScaleCache (interp, cachePtr, size)
Tcl_Interp *interp;             /* not used right now */
Imgtools_ScaleCache *cachePtr;  /* points to first element of cache */
int size;                       /* number of elements in cache */
{
    Tcl_Free((char *)cachePtr[0].weights);
    Tcl_Free((char *)cachePtr[0].srcPos);
    
    return TCL_OK;
}

/*
 * === Rotating Images ===
 */

/*
 *--------------------------------------------------------------
 *
 * Imgtools_Rotate --
 *
 *      Main Function for rotating images. Called by the user via
 *      the corresponding command (::imgtools::rotate)
 *
 * Results:
 *      TCL_OK or TCL_ERROR.
 *
 * Side effects:
 *      writes a modified image (according to options) to destImg. 
 *      destImg must exist, its contents will be partially or  
 *      completely overwritten. If the -width or -height options 
 *      for the dest image were set, these sizes are respected.
 *
 *--------------------------------------------------------------
 */
int Imgtools_Rotate(clientData, interp, objc, objv) 
ClientData clientData;  /* not used */
Tcl_Interp *interp;     /* current interpreter, used for error messages */
int objc;               /* argument count */
Tcl_Obj *const objv[];  /* argumetn values */
{
    
    Cputime_SetScopes(interp, scopeNameList);
    Cputime_NewScope(interp, MAIN);
    
    int result = TCL_OK;
    
    Imgtools_Image srcData, destData;
    srcData.memOwner = 0;   /* make sure we can safely call Imgtools_FreeImage() */
    destData.memOwner = 0;  /* on these before they're initialised (early return on rerros) */
    
    Imgtools_CommonOpts commonOpts;
    int allowedOpts;
    
    double angle;
    
    /*
     * Parse arguments
     */    
    Tk_PhotoHandle srcimg, destimg;
    Tcl_Obj *destImgName = NULL;
    
    if (objc < 4) {
        Tcl_WrongNumArgs(interp, 1, objv, "sourceImg angle ?options? destImg");
        Tcl_SetErrorCode(interp, "IMGTOOLS", "OPTIONS", "WRONG_NUM_ARGS", NULL);
        result = TCL_ERROR;
        goto cleanupAndReturn;
    }
    if ((result = Imgtools_FindPhoto(interp, objv[1], &srcimg)) !=TCL_OK) {
        goto cleanupAndReturn;
    }
    if ((result = Imgtools_FindPhoto(interp, objv[objc-1], &destimg)) != TCL_OK) {
        goto cleanupAndReturn;
    }
    destImgName = objv[objc-1];
    
    if ((result = Imgtools_GetAngle(interp, objv[2], &angle)) != TCL_OK) {
        Tcl_SetErrorCode(interp, "IMGTOOLS", "OPTIONS", "BAD_VALUE", NULL);
        goto cleanupAndReturn;
    }
    if ((result = Imgtools_ReduceAngle(interp, angle, &angle)) != TCL_OK) {
        Tcl_SetErrorCode(interp, "IMGTOOLS", "OPTIONS", "BAD_VALUE", NULL);
        goto cleanupAndReturn;
    }
    
    allowedOpts = OPT_ALPHA | OPT_COMPOSITE | OPT_SHRINK | OPT_FILTER 
        | OPT_BACKGROUND | OPT_CLIPPING | OPT_FROM | OPT_TO;
#ifdef IMGTOOLS_DEBUG
    allowedOpts |= OPT_FOLLOWPIXEL;
    int debugLevelSave = 0;
#endif
    if ((result = Imgtools_ParseCommonOpts(interp, objc-4, objv+3, allowedOpts, 
            &commonOpts)) != TCL_OK) {
        goto cleanupAndReturn;
    }
    
    /*
     * Get image data and adapt to 
     * -from option
     */
    Cputime_NewScope(interp, CONVERT);
    if ((result = Imgtools_GetImage(interp, srcimg, &srcData)) != TCL_OK) {
        goto cleanupAndReturn;
    }
    if ((result = Imgtools_ApplySubregion(interp, commonOpts.fromCorners, &srcData)) 
            != TCL_OK) {
        goto cleanupAndReturn;
    }
    Cputime_EndScopeSecure(interp, CONVERT);
    
    /*
     * with an empty source image, we don't have much to do
     */
    if (srcData.width == 0 || srcData.height == 0) {
        destData = srcData;
        destData.memOwner = 0;  /* Even if we own the memory for source data, */
                                /* this one has to be 0 to avoid freeing twice the same block */
        goto putDataAndReturn;
    }
    
    /*
     * With angle == 0, return source image
     */
    if (angle == 0) {
        destData = srcData;
        destData.memOwner = 0;
        goto putDataAndReturn;
    }
    
#ifdef IMGTOOLS_DEBUG
    /*
     * If we got -followpixel, save current debuglevel
     * the DEBUG_FOLLOW level is set only when we reach the target pixel in RotateImage
     */
    debugLevelSave = Imgtools_GetDebugLevel(interp);
    if (commonOpts.options & OPT_FOLLOWPIXEL) {
        Imgtools_SetDebugLevel(interp, DEBUG_PARAMETERS | DEBUG_FOLLOW);
        debugmsg(interp, DEBUG_FOLLOW, "Data for dest pixel %d, %d\n", 
                 commonOpts.followPixel.x, commonOpts.followPixel.y);
        debugmsg(interp, DEBUG_FOLLOW, "=================================\n");
        Imgtools_SetDebugLevel(interp, DEBUG_PARAMETERS);
    }
#endif

    /*
     * Time to interpolate!
     */
    unsigned char outsidePixelsColor[4];
    memcpy(outsidePixelsColor, commonOpts.background, 3 * sizeof(unsigned char));
    outsidePixelsColor[3] = roundf(commonOpts.alpha * UCHAR_MAX);
    
    Cputime_NewScope(interp, ROTATE);
    
    /*
     * If angle is a multiple of 90, we can use the faster rotate90 function
     */
    if (fmod(angle, 90) == 0) {
        if (Imgtools_RotateImage90(interp, &srcData, angle, &commonOpts, 
                outsidePixelsColor, &destData) != TCL_OK) {
            Tcl_SetObjResult(interp, Tcl_ObjPrintf("unexpected error rotating the image by %.0f°", angle));
            Tcl_SetErrorCode(interp, "IMGTOOLS", "ROTATE", "ROTATE_FAILED", NULL);
            result = TCL_ERROR;
            goto cleanupAndReturn;
        }
    } else {
        if (Imgtools_RotateImage(interp, &srcData, angle, &commonOpts, 
                outsidePixelsColor, &destData) != TCL_OK) {
            Tcl_SetObjResult(interp, Tcl_NewStringObj("something went wrong rotating the image. Sorry.", -1));
            Tcl_SetErrorCode(interp, "IMGTOOLS", "ROTATE", "ROTATE_FAILED", NULL);
            result = TCL_ERROR;
            goto cleanupAndReturn;
        }
    }

    Cputime_EndScopeSecure(interp, ROTATE);
    
#ifdef IMGTOOLS_DEBUG
    /*
     * Reset previous debuglevel settings
     */
    Imgtools_SetDebugLevel(interp, debugLevelSave);
#endif


    /*
     * Put data to destImg
     */
putDataAndReturn:
    Cputime_NewScope(interp, CONVERT);
    Imgtools_IntCoords start, end, size;
    size.x = destData.width;
    size.y = destData.height;
    Imgtools_GetRegion(interp, size, commonOpts.toCorners, ALLOW_LARGER,
        SIZE_IS_INCREMENT, &start, &end);
    Imgtools_GetSizeFromCorners(start, end, &size);
    if (Imgtools_PutImage(interp, commonOpts.composite, commonOpts.shrink, &destData, 
            start, size, Tcl_GetString(destImgName), destimg) != TCL_OK) {
        Tcl_SetObjResult(interp, Tcl_NewStringObj("failed to put data to destination image", -1));
        Tcl_SetErrorCode(interp, "IMGTOOLS", "DATA", "EXPORT_FAILED", NULL);
        result = TCL_ERROR;
        goto cleanupAndReturn;
    }

    Cputime_EndScopeSecure(interp, CONVERT);
    
    Tcl_SetObjResult(interp, destImgName);
    result = TCL_OK;

cleanupAndReturn:
    Imgtools_FreeImage(&srcData);
    Imgtools_FreeImage(&destData);
    
    Cputime_EndScopeSecure(interp, MAIN);
    Cputime_PrintSummary(interp);

    return result;
}

/*
 *--------------------------------------------------------------
 *
 * Imgtools_RotateImage --
 *
 *      Rotate the source image by a specified amount of degrees
 * 
 * Results:
 *      TCL_OK or TCL_ERROR.
 *
 * Side effects:
 *      Results are written to destData
 *      Allocates the needed memory for destData->pixelPtr
 *      The size of destImage might be bigger than source image,
 *      depending on the -clipping option.
 *      
 *--------------------------------------------------------------
 */
static int Imgtools_RotateImage(interp, sourceData, angle, options, outsidePixel, destData)
Tcl_Interp *interp;                 /* used for interp-global data / debugmsg */
const Imgtools_Image *sourceData;   /* image to rotate */
double angle;                       /* image is rotated by this amount of degrees */
const Imgtools_CommonOpts *options; /* filter, follow pixel coords, maybe more... */
const unsigned char outsidePixel[4];/* Values (RGBA) to use for Pixels that are outside the source image */
Imgtools_Image *destData;           /* rotated image is written here */
{
    // Calculate what we can outside the loops
    Imgtools_FloatCoords destSize = {0, 0};
    Imgtools_FloatCoords posCorrection, destPixelPos;
    Imgtools_IntCoords destCoords;
    double sinTheta, cosTheta, theta;
    double ySinTheta, yCosTheta;
    double r;
    enum FilterType filter = options->filter;
    
    /* 
     * the value we get for angle should already be ok
     */
    assert(angle <= 360 && angle >= -360);
    assert(sourceData->pixelSize == 4);
    
    theta = angle / 180 * acos(-1.0);
    sinTheta = sin(theta);
    cosTheta = cos(theta);
    
    /*
     * see http://tkimgtools.sourceforge.net/math.html for explanations about these formulas
     */
    switch (options->clipping) {
        case NONE:
            destSize.x = fabs(sinTheta) * sourceData->height + fabs(cosTheta) * sourceData->width;
            destSize.y = fabs(cosTheta) * sourceData->height + fabs(sinTheta) * sourceData->width;
            break;
        case KEEP:
            destSize.x = sourceData->width;
            destSize.y = sourceData->height;
            break;
        case ALL:
            r = (double)MIN(sourceData->width, sourceData->height) / MAX(sourceData->width, sourceData->height);
            if ( fabs(sin(2*theta)) < r) {
                destSize.x = (fabs(cosTheta) * sourceData->width - fabs(sinTheta) * sourceData->height) /
                    cos(2*theta);
                destSize.y = (fabs(cosTheta) * sourceData->height - fabs(sinTheta) * sourceData->width) /
                    cos(2*theta);
            }
            else if ( sourceData->width < sourceData->height ) {
                destSize.x = sourceData->width / (2 * fabs(cosTheta));
                destSize.y = sourceData->width / (2 * fabs(sinTheta));
            }
            else {
                destSize.x = sourceData->height / (2 * fabs(sinTheta));
                destSize.y = sourceData->height / (2 * fabs(cosTheta));
            }
            break;
        default:
            Tcl_Panic("RotateImage: unknown clipping %d\n", options->clipping);
    }
    
    /*
     * Coordinate translation and position correction are independent 
     * from the clipping method. We just need the dest size that
     * we already calculated
     */
    //posCorrection.x = -(destSize.x-1) / 2 * cosTheta - (destSize.y-1) / 2 * sinTheta
    //    + (sourceData->width-1) / 2;
    posCorrection.x = -(destSize.x-1) / 2 * cosTheta - (destSize.y-1) / 2 * sinTheta
        + ((double)sourceData->width-1) / 2;
    //posCorrection.x = ((sourceData->width - 1) - (destSize.x - 1)*cosTheta - (destSize.y-1)*sinTheta) / 2;
    posCorrection.y = (destSize.x-1) / 2 * sinTheta - (destSize.y-1) / 2 * cosTheta 
        + ((double)sourceData->height-1) / 2;
    
    debugmsg(interp, DEBUG_PARAMETERS, "Imgtools_RotateImage: theta: %f, sinTheta: %f, cosTheta: %f\n", theta, sinTheta, cosTheta);
    debugmsg(interp, DEBUG_PARAMETERS, "Imgtools_RotateImage: sourceData size: %dx%d\n", sourceData->width, sourceData->height);
    debugmsg(interp, DEBUG_PARAMETERS, "Imgtools_RotateImage: destSize: %fx%f\n", destSize.x, destSize.y);
    debugmsg(interp, DEBUG_PARAMETERS, "Imgtools_RotateImage: posCorrection: x=%.6f, y=%.6f\n", posCorrection.x, posCorrection.y);
    debugmsg(interp, DEBUG_PARAMETERS, "Imgtools_RotateImage: filter: %s (%d)\n", filterNames[filter], filter);
    
    /*
     * Initialize destData
     */
    destData->width = round(destSize.x);
    destData->height = round(destSize.y);
    destData->pixelSize = 4;
    destData->pitch = destData->pixelSize * destData->width;
    destData->pixelPtr = (unsigned char *)Tcl_Alloc(destData->width 
        * destData->height * destData->pixelSize * sizeof(unsigned char));
    destData->startPtr = destData->pixelPtr;
    destData->memOwner = 1;
    
    unsigned char *curPixelPtr;
    int halfFW = filterWidth[filter] / 2;
    
    Cputime_NewScope(interp, ROTATE_LOOP);

    for (destCoords.y = destData->height; destCoords.y--; ) {
        ySinTheta = destCoords.y * sinTheta;
        yCosTheta = destCoords.y * cosTheta;
        for (destCoords.x = destData->width; destCoords.x--; ) {
            destPixelPos.x = destCoords.x*cosTheta + ySinTheta + posCorrection.x;
            destPixelPos.y = yCosTheta - destCoords.x*sinTheta + posCorrection.y;
            curPixelPtr = destData->startPtr + destCoords.x*destData->pixelSize + destCoords.y*destData->pitch;

#ifdef IMGTOOLS_DEBUG_DISABLED
            // Causes severe hit on performance. Enable only if really needed...
            if (options->options & OPT_FOLLOWPIXEL
                && destCoords.x == options->followPixel.x 
                && destCoords.y == options->followPixel.y
            ) {
                Imgtools_SetDebugLevel(interp, DEBUG_FOLLOW);
            } else {
                Imgtools_SetDebugLevel(interp, DEBUG_NONE);
            }
#endif
            /*
             * permit destPixelPos to be up to half of filterWidth pixels outside source
             * this way we get interpolation with outsidePixel -> smoother edges
             */
            if (destPixelPos.x < -halfFW || destPixelPos.x >= sourceData->width + halfFW || destPixelPos.y < -halfFW || destPixelPos.y >= sourceData->height + halfFW) {
                curPixelPtr[RED] = outsidePixel[RED];
                curPixelPtr[GREEN] = outsidePixel[GREEN];
                curPixelPtr[BLUE] = outsidePixel[BLUE];
                curPixelPtr[ALPHA] = outsidePixel[ALPHA];
                /*debugmsg(interp, DEBUG_FOLLOW, "Assigned outsidePixel values to dest Pixel: %d %d %d %d\n\n",
                            curPixelPtr[RED], curPixelPtr[GREEN], curPixelPtr[BLUE], curPixelPtr[ALPHA]);*/
            } else {
                /*debugmsg(interp, DEBUG_FOLLOW, "Interpolating dest pixel %d, %d\n", destCoords.x, destCoords.y);
                debugmsg(interp, DEBUG_FOLLOW, "Position in source: %f, %f\n", destPixelPos.x, destPixelPos.y);*/
                // (*interpolFunc)(interp, sourceData, options, outsidePixel, &destPixelPos, curPixelPtr);
                (*Imgtools_InterpolFuncts[filter])
                    (interp, sourceData, options, outsidePixel, &destPixelPos, curPixelPtr);
            }
        }
    }
     
    Cputime_EndScopeSecure(interp, ROTATE_LOOP);

    return TCL_OK;    
}

/*
 *--------------------------------------------------------------
 *
 * Imgtools_RotateImage90 --
 *
 *      Rotate an image by a multiple of 90°.
 *      Does not support rotation by 0, 360, etc degrees! (no actual rotation cases)
 *      No need for interpolation here, just copy the pixels.
 *      
 * Results:
 *      Standard Tcl result.
 *
 * Side effects:
 *      The modified image is written to destImg.
 *      
 *--------------------------------------------------------------
 */

static int Imgtools_RotateImage90(interp, srcData, angle, options, outsidePixel, destData)
Tcl_Interp *interp;                 /* used for interp-global data / debugmsg */
const Imgtools_Image *srcData;      /* image to rotate */
double angle;                       /* image is rotated by this amount of degrees */
const Imgtools_CommonOpts *options; /* filter, follow pixel coords, maybe more... */
const unsigned char outsidePixel[4];/* Values (RGBA) to use for Pixels that are outside the source image */
Imgtools_Image *destData;           /* rotated image is written here */
{
    /*
     * Calculate what we can outside the loops
     */
    Imgtools_FloatCoords destSize = {0, 0};
    Imgtools_FloatCoords posCorrection;
    Imgtools_IntCoords destCoords, destPixelPos;
    int coordTranslation[4] = {0};  /* holds coefficients for coordinate translation (-1, 0 or 1) */
                                    /* element 0 and 1 are for x coords, 2 and 3 for y. */
    int intAngle;
    double theta;
    
    theta = angle / 180 * acos(-1.0);
    intAngle = angle;
    
    assert(intAngle < 360 && intAngle >= 0);
    assert(intAngle % 90 == 0);
    assert(intAngle != 0);
    
    switch (options->clipping) {
        case NONE:
        case ALL:
            if (intAngle == 180) {
                destSize.x = srcData->width;
                destSize.y = srcData->height;
            } else {
                destSize.x = srcData->height;
                destSize.y = srcData->width;
            }       
            break;
        case KEEP:
            destSize.x = srcData->width;
            destSize.y = srcData->height;
            break;
        default:
            Tcl_Panic("RotateImage90: unknown clipping %d\n", options->clipping);
    }
    
    /*
     * We need to round these, as they are gonna be used
     * in an int expression
     */
    posCorrection.x = round( -(destSize.x-1) / 2 * cos(theta) 
        - (destSize.y-1) / 2 * sin(theta) + ((double)srcData->width-1) / 2);
    posCorrection.y = round( (destSize.x-1) / 2 * sin(theta)
        - (destSize.y-1) / 2 * cos(theta) + ((double)srcData->height-1) / 2);
    
    /*
     * Coordinate translation is a little simpler here.
     */
    switch (intAngle) {
        case 90:
            coordTranslation[0] = 0;
            coordTranslation[1] = 1;
            coordTranslation[2] = -1;
            coordTranslation[3] = 0;
            break;
        case 180:
            coordTranslation[0] = -1;
            coordTranslation[1] = 0;
            coordTranslation[2] = 0;
            coordTranslation[3] = -1;
            break;
        case 270:
            coordTranslation[0] = 0;
            coordTranslation[1] = -1;
            coordTranslation[2] = 1;
            coordTranslation[3] = 0;
            break;
        default:
            Tcl_Panic("RotateImage90: got called with invalid angle");
    }
    
    /*
     * Initialize destData
     */
    destData->width = round(destSize.x);
    destData->height = round(destSize.y);
    destData->pixelSize = 4;
    destData->pitch = destData->pixelSize * destData->width;
    destData->pixelPtr = (unsigned char *)Tcl_Alloc(destData->width * destData->height * destData->pixelSize * sizeof(unsigned char));
    destData->startPtr = destData->pixelPtr;
    destData->memOwner = 1;
    
    unsigned char *destPixelPtr, *srcPixelPtr;
    
    for (destCoords.y = destData->height; destCoords.y--; ) {
        for (destCoords.x = destData->width; destCoords.x--; ) {
            destPixelPtr = destData->startPtr + destCoords.x*destData->pixelSize 
                + destCoords.y*destData->pitch;

            destPixelPos.x = destCoords.x * coordTranslation[0] 
                + destCoords.y * coordTranslation[1] + posCorrection.x;
            destPixelPos.y = destCoords.x * coordTranslation[2]
                + destCoords.y * coordTranslation[3] + posCorrection.y;
            srcPixelPtr = srcData->startPtr + destPixelPos.x * srcData->pixelSize
                + destPixelPos.y * srcData->pitch;
                
            if (destPixelPos.x < 0 || destPixelPos.x >= srcData->width 
                    || destPixelPos.y < 0 || destPixelPos.y >= srcData->height) {
                memcpy(destPixelPtr, outsidePixel, 4 * sizeof (unsigned char));
            } else {
                memcpy(destPixelPtr, srcPixelPtr, 4 * sizeof (unsigned char));
            }
        }
    }
     
    return TCL_OK;
    
}

/*
 * === Generic interpolation functions 
 */

/*
 *--------------------------------------------------------------
 *
 * Imgtools_Interpolate --
 *
 *      Calculate interpolation for one destination pixel.
 *      Assumes that the filter is scaled to source spacing (i.e. 
 *      the filter uses filterWidth[filter] source pixels).
 *      Thus, this function is not suited for downscaling images.
 *      
 * Results:
 *      Standard Tcl result
 *
 * Side effects:
 *      Results are written to destData
 *      
 *--------------------------------------------------------------
 */
inline static int Imgtools_Interpolate(interp, sourceData, options, outsidePixel, destPos, destPixelPtr)
Tcl_Interp *interp;
const Imgtools_Image *sourceData;       /* source Data. Must have pixelSize == 4 and */
                                        /* offset == {0, 1, 2, 3} */
const Imgtools_CommonOpts *options;     /* filter, followpixel, ... */
const unsigned char outsidePixel[4];    /* Source color to use for positions outside the bounds of source */
                                        /* RGBA values */
                                        /* if NULL, do not consider "outside pixels"*/
const Imgtools_FloatCoords *destPos;    /* position of destinaton pixel, relative to source data */
unsigned char *destPixelPtr;            /* adress of dest pixel - result is written there. */
                                        /* the offset for red, green, blue and alpha must be 0, 1, 2, and 3 */
{

    Cputime_NewScope(interp, INTERPOL);
    
    float xWeights[MAX_FILTERWIDTH], yWeights[MAX_FILTERWIDTH];
    float weightedSums[MAX_FILTERWIDTH] = {0};
    float weightsSum = 0;
    float weightedAlpha, curWeight;
    enum FilterType filter = options->filter;
            
    Imgtools_IntCoords startPos;
    startPos.x = floor(destPos->x) - filterWidth[filter]/2 + 1;
    startPos.y = floor(destPos->y) - filterWidth[filter]/2 + 1;
    
    /*
     * calculate weights for interpolation (same for each row / column) 
     */
    for (int i = filterWidth[filter]; i--; ) {
        xWeights[i] = (*Imgtools_FilterFuncts[filter])(startPos.x + i - destPos->x);        
        yWeights[i] = (*Imgtools_FilterFuncts[filter])(startPos.y + i - destPos->y);
    }
    
    /*
     * Calculate weighthed sums
     */
    const unsigned char *curPixelPtr;        
    for (int x = filterWidth[filter]; x--; ) {
        for (int y = filterWidth[filter]; y--; ) {
            if (startPos.x + x >= 0 && startPos.x + x < sourceData->width
                && startPos.y + y >= 0 && startPos.y + y < sourceData->height) {
                curPixelPtr = sourceData->startPtr + (startPos.x + x)*sourceData->pixelSize + (startPos.y + y)*sourceData->pitch;
            } else if (outsidePixel != NULL) {
                curPixelPtr = outsidePixel;
            } else {
                continue;
            }
            
            weightsSum += (curWeight = xWeights[x] * yWeights[y]);
            weightedAlpha = curPixelPtr[ALPHA] * curWeight;
            weightedSums[RED] += curPixelPtr[RED] * weightedAlpha;
            weightedSums[GREEN] += curPixelPtr[GREEN] * weightedAlpha;
            weightedSums[BLUE] += curPixelPtr[BLUE] * weightedAlpha;
            weightedSums[ALPHA] += weightedAlpha;

#ifdef IMGTOOLS_DEBUG_DISABLED
            /* 
             * The -followpixel feature is disabled in this function, as it has a
             * _huge_ impact on performance. And The need for it right now is not strong 
             * enough to implement something better.
             */
            
            /*
             * The decision, if we are to print followpixel messages is taken by the
             * caller (we don't know dest coords here), and controlled by setting
             * the debuglevel acordingly
             */
            debugmsg(interp, DEBUG_FOLLOW, "\tSource pixel %d, %d\n", x, y);
            debugmsg(interp, DEBUG_FOLLOW, "\t\tValues: %d %d %d %d\n", 
                     curPixelPtr[RED], curPixelPtr[GREEN], curPixelPtr[BLUE], curPixelPtr[ALPHA]);
            debugmsg(interp, DEBUG_FOLLOW, "\t\tAssociated alpha values: %d %d %d %d\n", 
                     curPixelPtr[RED]*curPixelPtr[ALPHA], curPixelPtr[GREEN]*curPixelPtr[ALPHA],
                     curPixelPtr[BLUE]*curPixelPtr[ALPHA], curPixelPtr[ALPHA]);
            debugmsg(interp, DEBUG_FOLLOW, "\t\txWeight: %f, yWeight: %f, combined %f\n",
                     xWeights[x], yWeights[y], xWeights[x]*yWeights[y]);
#endif
            
        }
    }
    
    /*
     * Divide by weights sum and disassociate alpha.
     */
#ifdef IMGTOOLS_DEBUG_DISABLED
    debugmsg(interp, DEBUG_FOLLOW, "Raw weighted sums: %f %f %f %f\n", 
             weightedSums[RED], weightedSums[GREEN], 
             weightedSums[BLUE], weightedSums[ALPHA]);
#endif
    
    if (weightedSums[ALPHA] == 0) {
        weightedSums[RED] = 0;
        weightedSums[GREEN] = 0;
        weightedSums[BLUE] = 0;
    } else {
        weightedSums[RED] /= weightedSums[ALPHA];
        weightedSums[GREEN] /= weightedSums[ALPHA];
        weightedSums[BLUE] /= weightedSums[ALPHA];
    }
    if (weightsSum == 0) {
        weightedSums[ALPHA] = 0;
    } else {
        weightedSums[ALPHA] /= weightsSum;
    }
    /*
     * round, clip and write to destData
     */
    float tmp;
    tmp = roundf(weightedSums[RED]);
    tmp = tmp < 0 ? 0 : tmp;
    tmp = tmp > UCHAR_MAX ? UCHAR_MAX : tmp;
    destPixelPtr[RED] = tmp;
    
    tmp = roundf(weightedSums[GREEN]);
    tmp = tmp < 0 ? 0 : tmp;
    tmp = tmp > UCHAR_MAX ? UCHAR_MAX : tmp;
    destPixelPtr[GREEN] = tmp;
    
    tmp = roundf(weightedSums[BLUE]);
    tmp = tmp < 0 ? 0 : tmp;
    tmp = tmp > UCHAR_MAX ? UCHAR_MAX : tmp;
    destPixelPtr[BLUE] = tmp;
    
    tmp = roundf(weightedSums[ALPHA]);
    tmp = tmp < 0 ? 0 : tmp;
    tmp = tmp > UCHAR_MAX ? UCHAR_MAX : tmp;
    destPixelPtr[ALPHA] = tmp;
    
#ifdef IMGTOOLS_DEBUG_DISABLED
    debugmsg(interp, DEBUG_FOLLOW, "Interpolated values (float): %f %f %f %f\n",
            weightedSums[RED], weightedSums[GREEN], 
            weightedSums[BLUE], weightedSums[ALPHA]);
    debugmsg(interp, DEBUG_FOLLOW, "Rounded and clipped values: %d %d %d %d\n\n",
            destPixelPtr[RED], destPixelPtr[GREEN], destPixelPtr[BLUE], destPixelPtr[ALPHA]);
#endif

    Cputime_EndScopeSecure(interp, INTERPOL);
            
    return TCL_OK;
}

/*
 *--------------------------------------------------------------
 *
 * Imgtools_NearestInterpol --
 *
 *      Calculate one destination pixel using nearest neighbor
 *      tecnique. 
 *      Assumes that the filter is scaled to source spacing (i.e. 
 *      there are only 4 source pixels to consider).
 *      Thus, this function is not suited for downsampling images.
 *      
 * Results:
 *      Standard Tcl result
 *
 * Side effects:
 *      Results are written to destData
 *      
 *--------------------------------------------------------------
 */
static int Imgtools_NearestInterpol(interp, sourceData, options, outsidePixel, destPos, destPixelPtr)
Tcl_Interp *interp;
const Imgtools_Image *sourceData;       /* source Data */
const Imgtools_CommonOpts * options;    /* used for followpixel debug option (not implemented here) */
const unsigned char outsidePixel[4];    /* Source color to use for positions outside the bounds of source */
                                        /* RGBA values */
                                        /* if NULL, do not consider "outside pixels"*/
const Imgtools_FloatCoords *destPos;    /* position of destinaton pixel, relative to source data */
unsigned char * destPixelPtr;           /* address dest pixel. The offsets for red, green, blue and alpha */
                                        /* must be 0, 1, 2 and 3 respectively */
{
    Cputime_NewScope(interp, INTERPOL);
    
    unsigned char fakeResult[4] = {0};  /* in case we're out of sourceData */
                                        /* and have no outside pixel, we return all 0 */
                                        /* just like the other interpol functions */
    const unsigned char *srcPixelPtr;
    Imgtools_IntCoords chosenPixel;
    
    chosenPixel.x = floor(destPos->x);
    chosenPixel.y = floor(destPos->y);
    
    if (destPos->x - chosenPixel.x > 0.5) {
        ++chosenPixel.x;
    }
    if (destPos->y - chosenPixel.y > 0.5) {
        ++chosenPixel.y;
    }
    
    if (chosenPixel.x < 0 || chosenPixel.x >= sourceData->width ||
        chosenPixel.y < 0 || chosenPixel.y >= sourceData->height) {
        
        srcPixelPtr = outsidePixel == NULL ? fakeResult : outsidePixel;
    } else {
        srcPixelPtr = sourceData->startPtr + chosenPixel.x*sourceData->pixelSize 
            + chosenPixel.y*sourceData->pitch;
    }
    
    memcpy(destPixelPtr, srcPixelPtr, sizeof(unsigned char) * 4);
        
    Cputime_EndScopeSecure(interp, INTERPOL);
    
    return TCL_OK;
}

/*
 *--------------------------------------------------------------
 *
 * Imgtools_GeneralLanczosFiler --
 *
 *      Calculate the value of LanczosN filter at the given point.
 *
 * Results:
 *      Returns the value of the filter at the given point
 *
 * Side effects:
 *      none
 *
 *--------------------------------------------------------------
 */
static float Imgtools_GeneralLanczosFilter (x, a)
const float x;   /* input value */
const int a;      /* half-width of the filter */
{
    float weight, xPi, xPiA;
    
    if (x <= -a || x >= a) {
        weight = 0;
    } else if (x == 0) {
        weight = 1;
    } else {
        xPi = M_PI * x;
        xPiA = xPi / a;
        weight = sinf(xPi) / xPi * sin(xPiA) / xPiA ;
    } 

    return weight;
}

/*
 *--------------------------------------------------------------
 *
 * Imgtools_Lanczos2Filer --
 *
 *      Calculate the value of Lanczos2 filter at the given point.
 *
 * Results:
 *      Returns the value of the filter at the given point
 *
 * Side effects:
 *      none
 *
 *--------------------------------------------------------------
 */
static float Imgtools_Lanczos2Filter (x)
float x;   /* input value */
{
    return Imgtools_GeneralLanczosFilter(x, 2);
}

/*
 *--------------------------------------------------------------
 *
 * Imgtools_CachedLanczos2Filer --
 *
 *      Get the value of Lanczos2 filter at the given point
 *      from the cache.
 *
 * Results:
 *      Returns the value of the filter at the given point
 *
 * Side effects:
 *      none
 *
 *--------------------------------------------------------------
 */
static float Imgtools_CachedLanczos2Filter (x)
float x;   /* input value */
{
    int idx;
    idx = roundf(x*LANCZOS2_CACHE_MULTIPLIER);
    idx = idx < LANCZOS2_CACHE_MININDEX ? LANCZOS2_CACHE_MININDEX : idx;
    idx = idx > LANCZOS2_CACHE_MAXINDEX ? LANCZOS2_CACHE_MAXINDEX : idx;
    return lanczos2Cache[idx];
}

/*
 *--------------------------------------------------------------
 *
 * Imgtools_Lanczos3Filer --
 *
 *      Calculate the value of Lanczos3 filter at the given point.
 *
 * Results:
 *      Returns the value of the filter at the given point
 *
 * Side effects:
 *      none
 *
 *--------------------------------------------------------------
 */

static float Imgtools_Lanczos3Filter (x)
float x;   /* input value */
{
    return Imgtools_GeneralLanczosFilter(x, 3);
}

/*
 *--------------------------------------------------------------
 *
 * Imgtools_CachedLanczos3Filer --
 *
 *      Get the value of Lanczos3 filter at the given point
 *      from the cached values.
 *
 * Results:
 *      Returns the value of the filter at the given point
 *
 * Side effects:
 *      none
 *
 *--------------------------------------------------------------
 */

static float Imgtools_CachedLanczos3Filter (x)
float x;   /* input value */
{
    int idx;
    idx = roundf(x*LANCZOS3_CACHE_MULTIPLIER);
    idx = idx < LANCZOS3_CACHE_MININDEX ? LANCZOS3_CACHE_MININDEX : idx;
    idx = idx > LANCZOS3_CACHE_MAXINDEX ? LANCZOS3_CACHE_MAXINDEX : idx;
    return lanczos3Cache[idx];
}

/*
 *--------------------------------------------------------------
 *
 * Imgtools_CatromFilter --
 *
 *      Calculate the value of Catmull-Rom filter at the given point.
 *
 * Results:
 *      Returns the value of the filter at the given point
 *
 * Side effects:
 *      none
 *
 *--------------------------------------------------------------
 */
static float Imgtools_CatromFilter (x)
float x;   /* input value */
{
    float weight;
    
    if (x < -2) {
        weight = 0;
    } else if (x < -1) {
        weight = 0.5*powf(x, 3) +2.5*powf(x, 2) +4*x +2;
    } else if (x < 0) {
        weight = -1.5*powf(x, 3) -2.5*powf(x, 2) +1;
    } else if (x < 1) {
        weight = 1.5*powf(x, 3) -2.5*powf(x, 2) +1;
    } else if (x < 2) {
        weight = -0.5*powf(x, 3) +2.5*powf(x, 2) -4*x +2;
    } else {
        weight = 0;
    }

    return weight;
}

/*
 *--------------------------------------------------------------
 *
 * Imgtools_CachedCatromFilter --
 *
 *      Returns the cached value of Catmull-Rom filter at the given point.
 *
 * Results:
 *      Returns the value of the filter at the given point
 *
 * Side effects:
 *      none
 *
 *--------------------------------------------------------------
 */
static float Imgtools_CachedCatromFilter (x)
const float x;  /* input value (distance) */
{
    int idx;
    idx = roundf(x*CATROM_CACHE_MULTIPLIER);
    idx = idx < CATROM_CACHE_MININDEX ? CATROM_CACHE_MININDEX : idx;
    idx = idx > CATROM_CACHE_MAXINDEX ? CATROM_CACHE_MAXINDEX : idx;
    return catromCache[idx];
}


/*
 *--------------------------------------------------------------
 *
 * Imgtools_LinearFilter --
 *
 *      Calculate the value of Linear filter at the given point.
 *
 * Results:
 *      Returns the value of the filter at the given point
 *
 * Side effects:
 *      none
 *
 *--------------------------------------------------------------
 */
static float Imgtools_LinearFilter (x)
float x;   /* input value */
{
    if (x < -1) {
        return 0;
    } else if (x < 0) {
        return x+1;
    } else if (x < 1) {
        return -x+1;
    } else {
        return 0;
    }
}

/*
 *--------------------------------------------------------------
 *
 * Imgtools_NearestFilter --
 *
 *      Calculate the value of nearest neighbor filter at the given point.
 *
 * Results:
 *      Returns the value of the filter at the given point
 *
 * Side effects:
 *      none
 *
 *--------------------------------------------------------------
 */
static float Imgtools_NearestFilter (x)
float x;   /* input value */
{
    return fabsf(x) < 0.5 ? 1 : 0;
}
