1//
2//  PRDFTFilter.m
3//  PRICE
4//
5//  Created by Riccardo Mottola on Tue Nov 18 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 "PRDFTFilter.h"
12#import "PRGrayscaleFilter.h"
13#import "FFT.h"
14#include <math.h>
15#import "PRProgressAction.h"
16
17extern double *cosinus;
18extern double *sinus;
19extern unsigned int *bitRevIndex;
20
21@implementation PRDFTFilter
22
23- (PRImage *)filterImage :(PRImage *)srcImage :(double **)filterMatRe :(double **)filterMatIm :(unsigned int) num :(BOOL)autoRange :(id)controller
24/* perform a convolution between the source image and the filter given by the two matrices */
25{
26  NSBitmapImageRep *srcImageRep;
27  PRImage *destImage;
28  NSBitmapImageRep *destImageRep;
29
30  unsigned int w, h;
31  unsigned int x, y; /* image scanning variables */
32  unsigned int i;    /* convolve matrix scanning */
33  unsigned char *srcData;
34  unsigned char *destData;
35  NSInteger destSamplesPerPixel;
36  NSInteger srcBytesPerRow;
37  NSInteger srcBytesPerPixel;
38
39  unsigned int leftBorder, topBorder;  /* borders to center the image */
40  unsigned int rightBorder, bottomBorder;
41
42  unsigned int bitNum;
43  double *Aa;
44  double *Ab;
45  double *ya;
46  double *yb;
47  double **Ma;   /* non quantized, float image */
48  double **Mb;
49
50  double min, max, scale;
51
52  if (controller)
53    {
54      [controller setActivity:@"initialize FFT"];
55      [controller advanceProgress];
56    }
57
58     /* get source image representation and associated information */
59    srcImageRep = [srcImage bitmapRep];
60
61    w = [srcImage width];
62    h = [srcImage height];
63
64    /* calculate back the bit number */
65    NSAssert(num > 0, @"Internal error. num <= 0");
66    bitNum = binaryLog(num);
67
68    sinus = (double*)calloc(num, sizeof(double));
69    cosinus = (double*)calloc(num, sizeof(double));
70    bitRevIndex = (unsigned int*)calloc(num, sizeof(unsigned int));
71    initTrigonometrics(num, bitNum);
72
73
74    if (![srcImage hasColor])
75      {
76        PRGrayscaleFilter *grayFilter;
77        printf("Color image\n");
78        grayFilter = [[PRGrayscaleFilter alloc] init];
79        srcImage = [grayFilter filterImage:srcImage :METHOD_AVERAGE];
80        [grayFilter release];
81        NSLog (@"done greyscale converting");
82        /* we reget previous information that is no longer valid */
83        /* get source image representation and associated information */
84        srcImageRep = [srcImage bitmapRep];
85      }
86    srcBytesPerRow = [srcImageRep bytesPerRow];
87    srcBytesPerPixel = [srcImageRep bitsPerPixel] / 8;
88
89    /* allocate 1-D FFT arrays */
90    Aa = (double*)calloc(num, sizeof(double));
91    Ab = (double*)calloc(num, sizeof(double));
92    ya = (double*)calloc(num, sizeof(double));
93    yb = (double*)calloc(num, sizeof(double));
94
95    /* allocate float destination matrix */
96
97    Ma = (double**)calloc(num, sizeof(double*));
98    for (i = 0; i < num; i++)
99    {
100        Ma[i] = (double*)calloc(num, sizeof(double));
101        if (!Ma[i])
102        {
103            printf("out of memory allocating float matrix?\n");
104            return srcImage;
105        }
106        memset(Ma[i], '\000', num);
107    }
108    Mb = (double**)calloc(num, sizeof(double*));
109    for (i = 0; i < num; i++)
110    {
111        Mb[i] = (double*)calloc(num, sizeof(double));
112        if (!Mb[i])
113        {
114            printf("out of memory allocating float matrix?\n");
115            return srcImage;
116        }
117        memset(Mb[i], '\000', num);
118    }
119
120    /* allocate destination image and its representation */
121    destSamplesPerPixel = 1;
122    destImage = [[PRImage alloc] initWithSize:NSMakeSize(w, h)];
123    destImageRep = [[NSBitmapImageRep alloc]
124                    initWithBitmapDataPlanes:NULL
125                    pixelsWide:w
126                    pixelsHigh:h
127                    bitsPerSample:8
128                    samplesPerPixel:destSamplesPerPixel
129                    hasAlpha:NO
130                    isPlanar:NO
131                    colorSpaceName:NSCalibratedWhiteColorSpace
132                    bytesPerRow:w*destSamplesPerPixel
133                    bitsPerPixel:0];
134
135    srcData = [srcImageRep bitmapData];
136    destData = [destImageRep bitmapData];
137
138    /* calculate borders */
139    leftBorder = (num - w) / 2;
140    topBorder = (num - h) / 2;
141    rightBorder =  (num - w) - leftBorder;
142    bottomBorder = (num - h) - topBorder;
143    printf("top %d left %d\n", topBorder, leftBorder);
144
145    /* copy the image to the float matrix */
146    for (y = 0; y < h; y++)
147        for (x = 0; x < w; x++)
148            Ma[y + topBorder][x + leftBorder] = (double)srcData[y * srcBytesPerRow + x*srcBytesPerPixel] / UCHAR_MAX; /* the image is greyscale, no bit-depth displace */
149
150
151    /* replicate borders */
152
153    /* top */
154    for (y = 0; y < topBorder; y++)
155        for (x = 0; x < w; x++)
156            Ma[y][x + leftBorder] = Ma[2*topBorder-y-1][x + leftBorder];
157
158    /* bottom */
159    for (y = 0; y < bottomBorder; y++)
160        for (x = 0; x < w; x++)
161            Ma[y+h+topBorder][x + leftBorder] = Ma[(h+topBorder)-y-1][x + leftBorder];
162
163    if (leftBorder > 0)
164    {
165        /* top left */
166        for (y = 0; y < topBorder; y++)
167            for (x = 0; x < leftBorder; x++)
168                Ma[y][x] = Ma[2*topBorder-y-1][2*leftBorder-x-1];
169        /* left */
170        for (y = 0; y < h; y++)
171            for (x = 0; x < leftBorder; x++)
172                Ma[y + topBorder][x] = Ma[y + topBorder][2*leftBorder-x-1];
173        /* bottom left */
174        for (y = 0; y < topBorder; y++)
175            for (x = 0; x < leftBorder; x++)
176                Ma[y+h+topBorder][x] = Ma[(h+topBorder)-y-1][2*leftBorder-x-1];
177    }
178
179    if (rightBorder > 0)
180    {
181        /* top right */
182        for (y = 0; y < topBorder; y++)
183            for (x = 0; x < rightBorder; x++)
184                Ma[y][x+w+leftBorder] = Ma[2*topBorder-y-1][(w+leftBorder)-x-1];
185        /* right */
186        for (y = 0; y < h; y++)
187            for (x = 0; x < rightBorder; x++)
188                Ma[y + topBorder][x+w+leftBorder] = Ma[y + topBorder][(w+rightBorder)-x-1];
189        /* bottom right */
190        for (y = 0; y < bottomBorder; y++)
191            for (x = 0; x < rightBorder; x++)
192                Ma[y+h+topBorder][x+w+leftBorder] = Ma[(h+topBorder)-y-1][(w+rightBorder)-x-1];
193    }
194
195
196
197    /* execute the actual filtering */
198
199    if (controller)
200    {
201        [controller setActivity:@"FFT forward pass"];
202        [controller advanceProgress];
203    }
204    /* forward pass */
205    /* horizontal pass */
206    for (y = 0; y < num; y++)
207    {
208        for (i = 0; i < num; i++)
209        {
210            Aa[i] = Ma[y][i];
211            Ab[i] = Mb[y][i];
212        }
213        fft(num, bitNum, Aa, Ab, ya, yb);
214        for (i = 0; i < num; i++)
215        {
216            Ma[y][i] = ya[i];
217            Mb[y][i] = yb[i];
218        }
219    }
220
221
222    /* vertical pass */
223    for (x = 0; x < num; x++)
224    {
225        for (i = 0; i < num; i++)
226        {
227            Aa[i] = Ma[i][x];
228            Ab[i] = Mb[i][x];
229        }
230        fft(num, bitNum, Aa, Ab, ya, yb);
231        for (i = 0; i < num; i++)
232        {
233            Ma[i][x] = ya[i];
234            Mb[i][x] = yb[i];
235        }
236    }
237
238//    printf("Fa %f %f\n", filterMatRe[0][0], filterMatRe[30][30]);
239
240    /* applying the transformation */
241    if (filterMatIm) /* we have a filter with I part */
242        for (x = 0; x < num; x++)
243            for (y = 0; y < num; y++)
244            {
245                Ma[x][y] = Ma[x][y]*filterMatRe[x][y] - Mb[x][y]*filterMatIm[x][y];
246                Mb[x][y] = Ma[x][y]*filterMatIm[x][y] + Mb[x][y]*filterMatRe[x][y];
247            }
248    else
249        for (x = 0; x < num; x++)
250            for (y = 0; y < num; y++)
251            {
252                Ma[x][y] = Ma[x][y]*filterMatRe[x][y];
253                Mb[x][y] = Mb[x][y]*filterMatRe[x][y];
254            }
255
256    if (controller)
257    {
258        [controller setActivity:@"FFT reverse pass"];
259        [controller advanceProgress];
260    }
261    /* reverse pass */
262    /* horizontal pass */
263    for (y = 0; y < num; y++)
264    {
265        for (i = 0; i < num; i++)
266        {
267            Aa[i] = Ma[y][i];
268            Ab[i] = Mb[y][i];
269        }
270        ifft(num, bitNum, Aa, Ab, ya, yb);
271        for (i = 0; i < num; i++)
272        {
273            Ma[y][i] = ya[i];
274            Mb[y][i] = yb[i];
275        }
276    }
277
278
279    /* vertical pass */
280    for (x = 0; x < num; x++)
281    {
282        for (i = 0; i < num; i++)
283        {
284            Aa[i] = Ma[i][x];
285            Ab[i] = Mb[i][x];
286        }
287        ifft(num, bitNum, Aa, Ab, ya, yb);
288        for (i = 0; i < num; i++)
289        {
290            Ma[i][x] = ya[i];
291            Mb[i][x] = yb[i];
292        }
293    }
294
295    if (autoRange)
296    {   /* AutoRange result */
297        /* now we find the range */
298        min = DBL_MAX;
299        max = DBL_MIN;
300        for (y = 0; y < h; y++)
301        {
302            for (x = 0; x < w; x++)
303            {
304                if (Ma[y+topBorder][x+leftBorder] > max)
305                    max = Ma[y+topBorder][x+leftBorder];
306                if (Ma[y+topBorder][x+leftBorder] < min)
307                    min = Ma[y+topBorder][x+leftBorder];
308            }
309        }
310        scale = fabs(max - min)/(double)UCHAR_MAX;
311        printf ("min = %f, max = %f, scale = %f\n", min, max, scale);
312        /* now we copy the result back, after scaling  */
313        for (y = 0; y < h; y++)
314            for (x = 0; x < w; x++)
315                destData[y*w + x] = (unsigned char) rint((Ma[y+topBorder][x+leftBorder]-min)/scale);
316    } else
317    { /* just scale back the result */
318        double sample;
319        /* now we copy the result back, after clipping  */
320        for (y = 0; y < h; y++)
321            for (x = 0; x < w; x++)
322            {
323                sample = rint(Ma[y+topBorder][x+leftBorder]*UCHAR_MAX);
324                if (sample > UCHAR_MAX)
325                    sample = UCHAR_MAX;
326                else if (sample < 0)
327                    sample = 0;
328                destData[y*w + x] = (unsigned char) sample;
329            }
330    }
331
332
333    [destImage setBitmapRep:destImageRep];
334    [destImageRep release];
335
336
337    /* free up the 1-D arrays */
338    free(Aa);
339    free(Ab);
340    free(ya);
341    free(yb);
342
343    /* now we free our float matrix */
344    for (i = 0; i < num; i++)
345        free(Ma[i]);
346    free(Ma);
347    for (i = 0; i < num; i++)
348        free(Mb[i]);
349    free(Mb);
350
351    /* free other arrays */
352    free(sinus);
353    free(cosinus);
354    free(bitRevIndex);
355
356    if (controller)
357    {
358        [controller setActivity:@"FFT done"];
359        [controller showProgress];
360    }
361
362    [destImage autorelease];
363    return destImage;
364}
365
366@end
367