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