//
//  PRFourier.m
//  PRICE
//
//  Created by Riccardo Mottola on Fri Jan 24 2003.
//  Copyright (c) 2003-2005 Carduus. All rights reserved.
//
// This program is free software; you can redistribute it and/or modify it under the terms of the version 2 of the GNU General Public License as published by the Free Software Foundation.
// This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

#import "PRFourier.h"
#import "PRGrayscaleFilter.h"
#include "FFT.h"

extern double *cosinus;
extern double *sinus;
extern unsigned int *bitRevIndex;

@implementation PRFourier

- (PRImage *)transformImage:(PRImage *)srcImage
{
    NSBitmapImageRep *srcImageRep;
    PRImage *destImage;
    NSBitmapImageRep *destImageRep;
    unsigned int w, h;
    unsigned int x, y; /* image scanning variables */
    unsigned int halfW, halfH; /* the middle */
    unsigned int i;    /* convolve matrix scanning */
    unsigned char *srcData;
    unsigned char *destData;
    unsigned char *p1, *p2;
    int bytesPerPixel;
    
    int num;
    unsigned int bitNum;
    double *Aa;
    double *Ab;
    double *ya;
    double *yb;
    double **Ma;   /* non quantized output */
    double **Mb;
    
    double min, max, scale;
    
    char tempStr[256];
    
    /* get source image representation and associated information */
    srcImageRep = [srcImage tiffRep];

    w = [srcImage width];
    h = [srcImage height];
    bytesPerPixel = [srcImage bitsPP] / 8;

    if (w > h)
        bitNum = binaryLog(w);
    else
        bitNum = binaryLog(h);
    num = binpow(bitNum);
    
    /* find the center, maybe a better way is needed */
    halfW = (unsigned int) rint(w / 2);
    halfH = (unsigned int) rint(h / 2);
    
    sinus = (double*)calloc(num, sizeof(double));
    cosinus = (double*)calloc(num, sizeof(double));
    bitRevIndex = (unsigned int*)calloc(num, sizeof(unsigned int));
    initTrigonometrics(num, bitNum);
    
    Aa = (double*)calloc(num, sizeof(double));
    Ab = (double*)calloc(num, sizeof(double));
    ya = (double*)calloc(num, sizeof(double));
    yb = (double*)calloc(num, sizeof(double));
    
    
    if ([srcImageRep hasAlpha])
    {
        if ([srcImageRep samplesPerPixel] == 2)
        {
            printf("Grayscale image + alpha\n");
        } else
        {
            printf("Color image + alpha\n");
            return srcImage;
        }
    } else
    {
        if ([srcImageRep samplesPerPixel] == 1)
        {
            printf("Grayscale image\n");
        } else
        {
            PRGrayscaleFilter *grayFilter;
            printf("Color image\n");
            grayFilter = [[PRGrayscaleFilter alloc] init];
            srcImage = [grayFilter filterImage:srcImage];
            [grayFilter release];
            NSLog (@"done greyscale converting");
            /* we reget previous information that is no longer valid */
            /* get source image representation and associated information */
            srcImageRep = [NSBitmapImageRep imageRepWithData:[srcImage TIFFRepresentation]];
            bytesPerPixel = [srcImageRep bitsPerPixel] /8;
        }
    }
    
    /* allocate float destination matrix */
    Ma = (double**)calloc(num, sizeof(double*));
    for (i = 0; i < num; i++)
    {
        Ma[i] = (double*)calloc(num, sizeof(double));
        if (!Ma[i])
        {
            printf("out of memory allocating float matrix?\n");
            return srcImage;
        }
        memset(Ma[i], '\000', num);
    }
    Mb = (double**)calloc(num, sizeof(double*));
    for (i = 0; i < num; i++)
    {
        Mb[i] = (double*)calloc(num, sizeof(double));
        if (!Mb[i])
        {
            printf("out of memory allocating float matrix?\n");
            return srcImage;
        }
        memset(Mb[i], '\000', num);
    }
    
    /* allocate destination image and its representation */
    destImage = [[PRImage alloc] initWithSize:NSMakeSize(w, h)];
    destImageRep = [[NSBitmapImageRep alloc] \
                    initWithBitmapDataPlanes:NULL \
                    pixelsWide:w \
                    pixelsHigh:h \
                    bitsPerSample:8 \
                    samplesPerPixel:1 \
                    hasAlpha:NO \
                    isPlanar:NO \
                    colorSpaceName:NSCalibratedWhiteColorSpace \
                    bytesPerRow:0 \
                    bitsPerPixel:0] ;
    
    srcData = [srcImageRep bitmapData];
    destData = [destImageRep bitmapData];
    
    [[srcImageRep colorSpaceName] getCString:tempStr];
//    printf("image space: %s\n", tempStr);


    /* copy the image to the float matrix */
    for (y = 0; y < h; y++)
    {
        for (x = 0; x < w; x++)
        {
            p1 = srcData +  (y * w + x);  /* the image is greyscale, no bit-depth displace */
            Ma[y][x] = (double) *p1 / UCHAR_MAX;
        }
    }
    /* replicate borders */
    if (w < num)
    {
        for (y = 0; y < h; y++)
        {
            for (x = w; x < num; x++)
            {
                p1 = srcData + (y * w + (w - (x-w+1)));
                Ma[y][x] = (double) *p1 / UCHAR_MAX;
            }
        }
    }
    if (h < num)
    {
        for (y = h; y < num; y++)
        {
            for (x = 0; x < w; x++)
            {
                p1 = srcData + ((h - (y-h+1)) * w + x);
                Ma[y][x] = (double) *p1 / UCHAR_MAX;
            }
        }
    }
    /* execute the actual filtering */
    /* horizontal pass */
    for (y = 0; y < num; y++)
    {
        for (i = 0; i < num; i++)
        {
            Aa[i] = Ma[y][i];
            Ab[i] = Mb[y][i];
        }
        fft(num, bitNum, Aa, Ab, ya, yb);
        for (i = 0; i < num; i++)
        {
            Ma[y][i] = ya[i];
            Mb[y][i] = yb[i];
        }
    }

    /* vertical pass */
    for (x = 0; x < num; x++)
    {
        for (i = 0; i < num; i++)
        {
            Aa[i] = Ma[i][x];
            Ab[i] = Mb[i][x];
        }
        fft(num, bitNum, Aa, Ab, ya, yb);
        for (i = 0; i < num; i++)
        {
            Ma[i][x] = ya[i];
            Mb[i][x] = yb[i];
        }
    }
    
    /* log of rised square modulus */
    for (y = 0; y < num; y++)
    {
        for (x = 0; x < num; x++)
        {
            double temp;
            temp = Ma[y][x]*Ma[y][x] + Mb[y][x]*Mb[y][x] + 1;
            if (temp <= 1)
                Ma[y][x] = 0;
            else
                Ma[y][x] = log(temp);
        }
    }

    /* now we find the range */
    min = DBL_MAX;
    max = DBL_MIN;
    y = 0;
    for (x = 1; x < w; x++)
    {
        if (Ma[y][x] > max)
            max = Ma[y][x];
        if (Ma[y][x] < min)
            min = Ma[y][x];
    }
    for (y = 1; y < h; y++)
    {
        for (x = 0; x < w; x++)
        {
            if (Ma[y][x] > max)
                max = Ma[y][x];
            if (Ma[y][x] < min)
                min = Ma[y][x];
        }
    }
    
    scale = fabs(max - min)/(double)UCHAR_MAX;
    printf ("min = %f, max = %f, scale = %f\n", min, max, scale);

    /* now we copy the result back, after scaling  */
    
    /* we recenter the FFT */
    /* top left */
    for (y = 0; y < halfH; y++)
    {
        p2 = destData + (y * w);
        for (x = 0; x < halfW; x++)
        {
            p2[x] = (unsigned char) rint((Ma[num-1-halfH+y][num-1-halfW+x]-min)/scale);
        }
    }
    /* top right */
    for (y = 0; y < halfH; y++)
    {
        p2 = destData + (y * w);
        for (x = halfW; x < w; x++)
        {
            p2[x] = (unsigned char) rint((Ma[num-1-halfH+y][x-halfW]-min)/scale);
        }
    }
    /* bottom left */
    for (y = halfH; y < h; y++)
    {
        p2 = destData + (y * w);
        for (x = 0; x < halfW; x++)
        {
            p2[x] = (unsigned char) rint((Ma[y-halfH][num-1-halfW+x]-min)/scale);
        }
    }
    /* bottom right */
    for (y = halfH; y < h; y++)
    {
        p2 = destData + (y * w);
        for (x = halfW; x < w; x++)
        {
            p2[x] = (unsigned char) rint((Ma[y-halfH][x-halfW]-min)/scale);
        }
    }

    [destImage addRepresentation:destImageRep];
    [destImageRep release];

    /* now we free our float matrix */
    for (i = 0; i < num; i++)
        free(Ma[i]);
    free(Ma);
    for (i = 0; i < num; i++)
        free(Mb[i]);
    free(Mb);

    [destImage autorelease];
    return destImage;
}


@end