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