1//
2//  PRFourier.m
3//  PRICE
4//
5//  Created by Riccardo Mottola on Fri Jan 24 2003.
6//  Copyright (c) 2003-2014 Carduus. All rights reserved.
7//
8// This application is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.
9// 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.
10
11#import "PRFourier.h"
12#import "PRGrayscaleFilter.h"
13#include "FFT.h"
14
15extern double *cosinus;
16extern double *sinus;
17extern unsigned int *bitRevIndex;
18
19@implementation PRFourier
20
21- (PRImage *)filterImage:(PRImage *)image with:(NSArray *)parameters progressPanel:(PRCProgress *)progressPanel
22{
23
24    /* interpret the parameters */
25
26    return [self transformImage:image];
27}
28
29- (NSString *)actionName
30{
31    return @"FFT";
32}
33
34- (PRImage *)transformImage:(PRImage *)srcImage
35{
36    NSBitmapImageRep *srcImageRep;
37    PRImage *destImage;
38    NSBitmapImageRep *destImageRep;
39    unsigned int w, h;
40    unsigned int x, y; /* image scanning variables */
41    unsigned int halfW, halfH; /* the middle */
42    unsigned int i;    /* convolve matrix scanning */
43    unsigned char *srcData;
44    unsigned char *destData;
45    unsigned char *p1, *p2;
46    int destSamplesPerPixel;
47    NSInteger srcBytesPerRow;
48    NSInteger srcBytesPerPixel;
49
50    int num;
51    unsigned int bitNum;
52    double *Aa;
53    double *Ab;
54    double *ya;
55    double *yb;
56    double **Ma;   /* non quantized output */
57    double **Mb;
58
59    double min, max, scale;
60
61
62    /* get source image representation and associated information */
63    srcImageRep = [srcImage bitmapRep];
64
65    w = [srcImage width];
66    h = [srcImage height];
67
68    if (w > h)
69        bitNum = binaryLog(w);
70    else
71        bitNum = binaryLog(h);
72    num = binpow(bitNum);
73    NSAssert(num > 0, @"PRFourier. Internal error. num <= 0");
74
75    /* find the center, maybe a better way is needed */
76    halfW = (unsigned int) rint(w / 2);
77    halfH = (unsigned int) rint(h / 2);
78
79    sinus = (double*)calloc(num, sizeof(double));
80    cosinus = (double*)calloc(num, sizeof(double));
81    bitRevIndex = (unsigned int*)calloc(num, sizeof(unsigned int));
82    initTrigonometrics(num, bitNum);
83
84    Aa = (double*)calloc(num, sizeof(double));
85    Ab = (double*)calloc(num, sizeof(double));
86    ya = (double*)calloc(num, sizeof(double));
87    yb = (double*)calloc(num, sizeof(double));
88
89
90    if ([srcImage hasColor])
91      {
92        PRGrayscaleFilter *grayFilter;
93        printf("Color image\n");
94        grayFilter = [[PRGrayscaleFilter alloc] init];
95        srcImage = [grayFilter filterImage:srcImage :METHOD_AVERAGE];
96        [grayFilter release];
97
98        srcImageRep = [srcImage bitmapRep];
99      }
100    srcBytesPerRow = [srcImageRep bytesPerRow];
101    srcBytesPerPixel = [srcImageRep bitsPerPixel] / 8;
102    destSamplesPerPixel = 1;
103
104    /* allocate float destination matrix */
105    Ma = (double**)calloc(num, sizeof(double*));
106    for (i = 0; i < num; i++)
107    {
108        Ma[i] = (double*)calloc(num, sizeof(double));
109        if (!Ma[i])
110        {
111            printf("out of memory allocating float matrix?\n");
112            return srcImage;
113        }
114        memset(Ma[i], '\000', num);
115    }
116    Mb = (double**)calloc(num, sizeof(double*));
117    for (i = 0; i < num; i++)
118    {
119        Mb[i] = (double*)calloc(num, sizeof(double));
120        if (!Mb[i])
121        {
122            printf("out of memory allocating float matrix?\n");
123            return srcImage;
124        }
125        memset(Mb[i], '\000', num);
126    }
127
128    /* allocate destination image and its representation */
129    destImage = [[PRImage alloc] initWithSize:NSMakeSize(w, h)];
130    destImageRep = [[NSBitmapImageRep alloc]
131                    initWithBitmapDataPlanes:NULL
132                    pixelsWide:w
133                    pixelsHigh:h
134                    bitsPerSample:8
135                    samplesPerPixel:destSamplesPerPixel
136                    hasAlpha:NO
137                    isPlanar:NO
138                    colorSpaceName:NSCalibratedWhiteColorSpace
139                    bytesPerRow:w*destSamplesPerPixel
140                    bitsPerPixel:0] ;
141
142    srcData = [srcImageRep bitmapData];
143    destData = [destImageRep bitmapData];
144
145
146    /* copy the image to the float matrix */
147    for (y = 0; y < h; y++)
148    {
149        for (x = 0; x < w; x++)
150        {
151            p1 = srcData +  (y * srcBytesPerRow + x*srcBytesPerPixel);
152            Ma[y][x] = (double) *p1 / UCHAR_MAX;
153        }
154    }
155    /* replicate borders */
156    if (w < num)
157    {
158        for (y = 0; y < h; y++)
159        {
160            for (x = w; x < num; x++)
161            {
162                p1 = srcData + (y * srcBytesPerRow + (srcBytesPerRow - (x*srcBytesPerPixel-srcBytesPerRow+1)));
163                Ma[y][x] = (double) *p1 / UCHAR_MAX;
164            }
165        }
166    }
167    if (h < num)
168    {
169        for (y = h; y < num; y++)
170        {
171            for (x = 0; x < w; x++)
172            {
173                p1 = srcData + ((h - (y-h+1)) * srcBytesPerRow + x*srcBytesPerPixel);
174                Ma[y][x] = (double) *p1 / UCHAR_MAX;
175            }
176        }
177    }
178    /* execute the actual filtering */
179    /* horizontal pass */
180    for (y = 0; y < num; y++)
181    {
182        for (i = 0; i < num; i++)
183        {
184            Aa[i] = Ma[y][i];
185            Ab[i] = Mb[y][i];
186        }
187        fft(num, bitNum, Aa, Ab, ya, yb);
188        for (i = 0; i < num; i++)
189        {
190            Ma[y][i] = ya[i];
191            Mb[y][i] = yb[i];
192        }
193    }
194
195    /* vertical pass */
196    for (x = 0; x < num; x++)
197    {
198        for (i = 0; i < num; i++)
199        {
200            Aa[i] = Ma[i][x];
201            Ab[i] = Mb[i][x];
202        }
203        fft(num, bitNum, Aa, Ab, ya, yb);
204        for (i = 0; i < num; i++)
205        {
206            Ma[i][x] = ya[i];
207            Mb[i][x] = yb[i];
208        }
209    }
210
211    /* log of rised square modulus */
212    for (y = 0; y < num; y++)
213    {
214        for (x = 0; x < num; x++)
215        {
216            double temp;
217            temp = Ma[y][x]*Ma[y][x] + Mb[y][x]*Mb[y][x] + 1;
218            if (temp <= 1)
219                Ma[y][x] = 0;
220            else
221                Ma[y][x] = log(temp);
222        }
223    }
224
225    /* now we find the range */
226    min = DBL_MAX;
227    max = DBL_MIN;
228    y = 0;
229    for (x = 1; x < w; x++)
230    {
231        if (Ma[y][x] > max)
232            max = Ma[y][x];
233        if (Ma[y][x] < min)
234            min = Ma[y][x];
235    }
236    for (y = 1; y < h; y++)
237    {
238        for (x = 0; x < w; x++)
239        {
240            if (Ma[y][x] > max)
241                max = Ma[y][x];
242            if (Ma[y][x] < min)
243                min = Ma[y][x];
244        }
245    }
246
247    scale = fabs(max - min)/(double)UCHAR_MAX;
248    printf ("min = %f, max = %f, scale = %f\n", min, max, scale);
249
250    /* now we copy the result back, after scaling  */
251
252    /* we recenter the FFT */
253    /* top left */
254    for (y = 0; y < halfH; y++)
255    {
256        p2 = destData + (y * w);
257        for (x = 0; x < halfW; x++)
258        {
259            p2[x] = (unsigned char) rint((Ma[num-1-halfH+y][num-1-halfW+x]-min)/scale);
260        }
261    }
262    /* top right */
263    for (y = 0; y < halfH; y++)
264    {
265        p2 = destData + (y * w);
266        for (x = halfW; x < w; x++)
267        {
268            p2[x] = (unsigned char) rint((Ma[num-1-halfH+y][x-halfW]-min)/scale);
269        }
270    }
271    /* bottom left */
272    for (y = halfH; y < h; y++)
273    {
274        p2 = destData + (y * w);
275        for (x = 0; x < halfW; x++)
276        {
277            p2[x] = (unsigned char) rint((Ma[y-halfH][num-1-halfW+x]-min)/scale);
278        }
279    }
280    /* bottom right */
281    for (y = halfH; y < h; y++)
282    {
283        p2 = destData + (y * w);
284        for (x = halfW; x < w; x++)
285        {
286            p2[x] = (unsigned char) rint((Ma[y-halfH][x-halfW]-min)/scale);
287        }
288    }
289
290    [destImage setBitmapRep:destImageRep];
291    [destImageRep release];
292
293    /* now we free our float matrix */
294    for (i = 0; i < num; i++)
295        free(Ma[i]);
296    free(Ma);
297    for (i = 0; i < num; i++)
298        free(Mb[i]);
299    free(Mb);
300
301    [destImage autorelease];
302    return destImage;
303}
304
305- (BOOL)displayProgress
306{
307    return NO;
308}
309
310
311
312@end
313