1 /*
2 FLAM3 - cosmic recursive fractal flames
3 Copyright (C) 1992-2009 Spotworks LLC
4
5 This program is free software; you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation; either version 3 of the License, or
8 (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19 #include "filters.h"
20
21
22 /*
23 * filter function definitions
24 * from Graphics Gems III code
25 * and ImageMagick resize.c
26 */
27
28
29 double flam3_spatial_support[flam3_num_spatialfilters] = {
30
31 1.5, /* gaussian */
32 1.0, /* hermite */
33 0.5, /* box */
34 1.0, /* triangle */
35 1.5, /* bell */
36 2.0, /* b spline */
37 2.0, /* mitchell */
38 1.0, /* blackman */
39 2.0, /* catrom */
40 1.0, /* hanning */
41 1.0, /* hamming */
42 3.0, /* lanczos3 */
43 2.0, /* lanczos2 */
44 1.5 /* quadratic */
45 };
46
flam3_hermite_filter(double t)47 double flam3_hermite_filter(double t) {
48 /* f(t) = 2|t|^3 - 3|t|^2 + 1, -1 <= t <= 1 */
49 if(t < 0.0) t = -t;
50 if(t < 1.0) return((2.0 * t - 3.0) * t * t + 1.0);
51 return(0.0);
52 }
53
flam3_box_filter(double t)54 double flam3_box_filter(double t) {
55 if((t > -0.5) && (t <= 0.5)) return(1.0);
56 return(0.0);
57 }
58
flam3_triangle_filter(double t)59 double flam3_triangle_filter(double t) {
60 if(t < 0.0) t = -t;
61 if(t < 1.0) return(1.0 - t);
62 return(0.0);
63 }
64
flam3_bell_filter(double t)65 double flam3_bell_filter(double t) {
66 /* box (*) box (*) box */
67 if(t < 0) t = -t;
68 if(t < .5) return(.75 - (t * t));
69 if(t < 1.5) {
70 t = (t - 1.5);
71 return(.5 * (t * t));
72 }
73 return(0.0);
74 }
75
flam3_b_spline_filter(double t)76 double flam3_b_spline_filter(double t) {
77
78 /* box (*) box (*) box (*) box */
79 double tt;
80
81 if(t < 0) t = -t;
82 if(t < 1) {
83 tt = t * t;
84 return((.5 * tt * t) - tt + (2.0 / 3.0));
85 } else if(t < 2) {
86 t = 2 - t;
87 return((1.0 / 6.0) * (t * t * t));
88 }
89 return(0.0);
90 }
91
flam3_sinc(double x)92 double flam3_sinc(double x) {
93 x *= M_PI;
94 if(x != 0) return(sin(x) / x);
95 return(1.0);
96 }
97
flam3_blackman_filter(double x)98 double flam3_blackman_filter(double x) {
99 return(0.42+0.5*cos(M_PI*x)+0.08*cos(2*M_PI*x));
100 }
101
flam3_catrom_filter(double x)102 double flam3_catrom_filter(double x) {
103 if (x < -2.0)
104 return(0.0);
105 if (x < -1.0)
106 return(0.5*(4.0+x*(8.0+x*(5.0+x))));
107 if (x < 0.0)
108 return(0.5*(2.0+x*x*(-5.0-3.0*x)));
109 if (x < 1.0)
110 return(0.5*(2.0+x*x*(-5.0+3.0*x)));
111 if (x < 2.0)
112 return(0.5*(4.0+x*(-8.0+x*(5.0-x))));
113 return(0.0);
114 }
115
flam3_mitchell_filter(double t)116 double flam3_mitchell_filter(double t) {
117 double tt;
118
119 tt = t * t;
120 if(t < 0) t = -t;
121 if(t < 1.0) {
122 t = (((12.0 - 9.0 * flam3_mitchell_b - 6.0 * flam3_mitchell_c) * (t * tt))
123 + ((-18.0 + 12.0 * flam3_mitchell_b + 6.0 * flam3_mitchell_c) * tt)
124 + (6.0 - 2 * flam3_mitchell_b));
125 return(t / 6.0);
126 } else if(t < 2.0) {
127 t = (((-1.0 * flam3_mitchell_b - 6.0 * flam3_mitchell_c) * (t * tt))
128 + ((6.0 * flam3_mitchell_b + 30.0 * flam3_mitchell_c) * tt)
129 + ((-12.0 * flam3_mitchell_b - 48.0 * flam3_mitchell_c) * t)
130 + (8.0 * flam3_mitchell_b + 24 * flam3_mitchell_c));
131 return(t / 6.0);
132 }
133 return(0.0);
134 }
135
flam3_hanning_filter(double x)136 double flam3_hanning_filter(double x) {
137 return(0.5+0.5*cos(M_PI*x));
138 }
139
flam3_hamming_filter(double x)140 double flam3_hamming_filter(double x) {
141 return(0.54+0.46*cos(M_PI*x));
142 }
143
flam3_lanczos3_filter(double t)144 double flam3_lanczos3_filter(double t) {
145 if(t < 0) t = -t;
146 if(t < 3.0) return(flam3_sinc(t) * flam3_sinc(t/3.0));
147 return(0.0);
148 }
149
flam3_lanczos2_filter(double t)150 double flam3_lanczos2_filter(double t) {
151 if(t < 0) t = -t;
152 if(t < 2.0) return(flam3_sinc(t) * flam3_sinc(t/2.0));
153 return(0.0);
154 }
155
flam3_gaussian_filter(double x)156 double flam3_gaussian_filter(double x) {
157 return(exp((-2.0*x*x))*sqrt(2.0/M_PI));
158 }
159
flam3_quadratic_filter(double x)160 double flam3_quadratic_filter(double x) {
161 if (x < -1.5)
162 return(0.0);
163 if (x < -0.5)
164 return(0.5*(x+1.5)*(x+1.5));
165 if (x < 0.5)
166 return(0.75-x*x);
167 if (x < 1.5)
168 return(0.5*(x-1.5)*(x-1.5));
169 return(0.0);
170 }
171
flam3_spatial_filter(int knum,double x)172 double flam3_spatial_filter(int knum, double x) {
173
174 if (knum==0)
175 return flam3_gaussian_filter(x);
176 else if (knum==1)
177 return flam3_hermite_filter(x);
178 else if (knum==2)
179 return flam3_box_filter(x);
180 else if (knum==3)
181 return flam3_triangle_filter(x);
182 else if (knum==4)
183 return flam3_bell_filter(x);
184 else if (knum==5)
185 return flam3_b_spline_filter(x);
186 else if (knum==6)
187 return flam3_mitchell_filter(x);
188 else if (knum==7)
189 return flam3_sinc(x)*flam3_blackman_filter(x);
190 else if (knum==8)
191 return flam3_catrom_filter(x);
192 else if (knum==9)
193 return flam3_sinc(x)*flam3_hanning_filter(x);
194 else if (knum==10)
195 return flam3_sinc(x)*flam3_hamming_filter(x);
196 else if (knum==11)
197 return flam3_lanczos3_filter(x)*flam3_sinc(x/3.0);
198 else if (knum==12)
199 return flam3_lanczos2_filter(x)*flam3_sinc(x/2.0);
200 else // if (knum==13)
201 return flam3_quadratic_filter(x);
202 }
203
normalize_vector(double * v,int n)204 int normalize_vector(double *v, int n) {
205 double t = 0.0;
206 int i;
207 for (i = 0; i < n; i++)
208 t += v[i];
209 if (0.0 == t) return 1;
210 t = 1.0 / t;
211 for (i = 0; i < n; i++)
212 v[i] *= t;
213 return 0;
214 }
215
216
flam3_create_spatial_filter(flam3_frame * spec,int field,double ** filter)217 int flam3_create_spatial_filter(flam3_frame *spec, int field, double **filter) {
218
219 int sf_kernel = spec->genomes[0].spatial_filter_select;
220 int supersample = spec->genomes[0].spatial_oversample;
221 double sf_radius = spec->genomes[0].spatial_filter_radius;
222 double aspect_ratio = spec->pixel_aspect_ratio;
223 double sf_supp = flam3_spatial_support[sf_kernel];
224
225 double fw = 2.0 * sf_supp * supersample * sf_radius / aspect_ratio;
226 double adjust, ii, jj;
227
228 int fwidth = ((int) fw) + 1;
229 int i,j;
230
231
232 /* Make sure the filter kernel has same parity as oversample */
233 if ((fwidth ^ supersample) & 1)
234 fwidth++;
235
236 /* Calculate the coordinate scaling factor for the kernel values */
237 if (fw > 0.0)
238 adjust = sf_supp * fwidth / fw;
239 else
240 adjust = 1.0;
241
242 /* Calling function MUST FREE THE RETURNED KERNEL, lest ye leak memory */
243 (*filter) = (double *)calloc(fwidth * fwidth,sizeof(double));
244
245 /* fill in the coefs */
246 for (i = 0; i < fwidth; i++)
247 for (j = 0; j < fwidth; j++) {
248
249 /* Calculate the function inputs for the kernel function */
250 ii = ((2.0 * i + 1.0) / (double)fwidth - 1.0)*adjust;
251 jj = ((2.0 * j + 1.0) / (double)fwidth - 1.0)*adjust;
252
253 /* Scale for scanlines */
254 if (field) jj *= 2.0;
255
256 /* Adjust for aspect ratio */
257 jj /= aspect_ratio;
258
259 (*filter)[i + j * fwidth] =
260 flam3_spatial_filter(sf_kernel,ii) * flam3_spatial_filter(sf_kernel,jj);
261 }
262
263
264 if (normalize_vector((*filter), fwidth * fwidth)) {
265 fprintf(stderr, "Spatial filter value is too small: %g. Terminating.\n",sf_radius);
266 return(-1);
267 }
268
269 return (fwidth);
270 }
271
flam3_create_de_filters(double max_rad,double min_rad,double curve,int ss)272 flam3_de_helper flam3_create_de_filters(double max_rad, double min_rad, double curve, int ss) {
273
274 flam3_de_helper de;
275 double comp_max_radius, comp_min_radius;
276 double num_de_filters_d;
277 int num_de_filters,de_max_ind;
278 int de_row_size, de_half_size;
279 int filtloop;
280 int keep_thresh=100;
281
282 de.kernel_size=-1;
283
284 if (curve <= 0.0) {
285 fprintf(stderr,"estimator curve must be > 0\n");
286 return(de);
287 }
288
289 if (max_rad < min_rad) {
290 fprintf(stderr,"estimator must be larger than estimator_minimum.\n");
291 fprintf(stderr,"(%f > %f) ? \n",max_rad,min_rad);
292 return(de);
293 }
294
295 /* We should scale the filter width by the oversample */
296 /* The '+1' comes from the assumed distance to the first pixel */
297 comp_max_radius = max_rad*ss + 1;
298 comp_min_radius = min_rad*ss + 1;
299
300 /* Calculate how many filter kernels we need based on the decay function */
301 /* */
302 /* num filters = (de_max_width / de_min_width)^(1/estimator_curve) */
303 /* */
304 num_de_filters_d = pow( comp_max_radius/comp_min_radius, 1.0/curve );
305 if (num_de_filters_d>1e7) {
306 fprintf(stderr,"too many filters required in this configuration (%g)\n",num_de_filters_d);
307 return(de);
308 }
309 num_de_filters = (int)ceil(num_de_filters_d);
310
311 /* Condense the smaller kernels to save space */
312 if (num_de_filters>keep_thresh) {
313 de_max_ind = (int)ceil(DE_THRESH + pow(num_de_filters-DE_THRESH,curve))+1;
314 de.max_filtered_counts = (int)pow( (double)(de_max_ind-DE_THRESH), 1.0/curve) + DE_THRESH;
315 } else {
316 de_max_ind = num_de_filters;
317 de.max_filtered_counts = de_max_ind;
318 }
319
320 /* Allocate the memory for these filters */
321 /* and the hit/width lookup vector */
322 de_row_size = (int)(2*ceil(comp_max_radius)-1);
323 de_half_size = (de_row_size-1)/2;
324 de.kernel_size = (de_half_size+1)*(2+de_half_size)/2;
325
326 de.filter_coefs = (double *) calloc (de_max_ind * de.kernel_size,sizeof(double));
327 de.filter_widths = (double *) calloc (de_max_ind,sizeof(double));
328
329 /* Generate the filter coefficients */
330 de.max_filter_index = 0;
331 for (filtloop=0;filtloop<de_max_ind;filtloop++) {
332
333 double de_filt_sum=0.0, de_filt_d;
334 double de_filt_h;
335 int dej,dek;
336 double adjloop;
337 int filter_coef_idx;
338
339 /* Calculate the filter width for this number of hits in a bin */
340 if (filtloop<keep_thresh)
341 de_filt_h = (comp_max_radius / pow(filtloop+1,curve));
342 else {
343 adjloop = pow(filtloop-keep_thresh,(1.0/curve)) + keep_thresh;
344 de_filt_h = (comp_max_radius / pow(adjloop+1,curve));
345 }
346
347 /* Once we've reached the min radius, don't populate any more */
348 if (de_filt_h <= comp_min_radius) {
349 de_filt_h = comp_min_radius;
350 de.max_filter_index = filtloop;
351 }
352
353 de.filter_widths[filtloop] = de_filt_h;
354
355 /* Calculate norm of kernel separately (easier) */
356 for (dej=-de_half_size; dej<=de_half_size; dej++) {
357 for (dek=-de_half_size; dek<=de_half_size; dek++) {
358
359 de_filt_d = sqrt( (double)(dej*dej+dek*dek) ) / de_filt_h;
360
361 /* Only populate the coefs within this radius */
362 if (de_filt_d <= 1.0) {
363
364 /* Gaussian */
365 de_filt_sum += flam3_spatial_filter(flam3_gaussian_kernel,
366 flam3_spatial_support[flam3_gaussian_kernel]*de_filt_d);
367
368 /* Epanichnikov */
369 // de_filt_sum += (1.0 - (de_filt_d * de_filt_d));
370 }
371 }
372 }
373
374 filter_coef_idx = filtloop*de.kernel_size;
375
376 /* Calculate the unique entries of the kernel */
377 for (dej=0; dej<=de_half_size; dej++) {
378 for (dek=0; dek<=dej; dek++) {
379 de_filt_d = sqrt( (double)(dej*dej+dek*dek) ) / de_filt_h;
380
381 /* Only populate the coefs within this radius */
382 if (de_filt_d>1.0)
383 de.filter_coefs[filter_coef_idx] = 0.0;
384 else {
385
386 /* Gaussian */
387 de.filter_coefs[filter_coef_idx] = flam3_spatial_filter(flam3_gaussian_kernel,
388 flam3_spatial_support[flam3_gaussian_kernel]*de_filt_d)/de_filt_sum;
389
390 /* Epanichnikov */
391 // de_filter_coefs[filter_coef_idx] = (1.0 - (de_filt_d * de_filt_d))/de_filt_sum;
392 }
393
394 filter_coef_idx ++;
395 }
396 }
397
398 if (de.max_filter_index>0)
399 break;
400 }
401
402 if (de.max_filter_index==0)
403 de.max_filter_index = de_max_ind-1;
404
405
406 return(de);
407 }
408
flam3_create_temporal_filter(int numsteps,int filter_type,double filter_exp,double filter_width,double ** temporal_filter,double ** temporal_deltas)409 double flam3_create_temporal_filter(int numsteps, int filter_type, double filter_exp, double filter_width,
410 double **temporal_filter, double **temporal_deltas) {
411
412 double maxfilt = 0.0;
413 double sumfilt = 0.0;
414 double slpx,halfsteps;
415 double *deltas, *filter;
416
417 int i;
418
419 /* Allocate memory - this must be freed in the calling routine! */
420 deltas = (double *)malloc(numsteps*sizeof(double));
421 filter = (double *)malloc(numsteps*sizeof(double));
422
423 /* Deal with only one step */
424 if (numsteps==1) {
425 deltas[0] = 0;
426 filter[0] = 1.0;
427 *temporal_deltas = deltas;
428 *temporal_filter = filter;
429 return(1.0);
430 }
431
432 /* Define the temporal deltas */
433 for (i = 0; i < numsteps; i++)
434 deltas[i] = ((double)i /(double)(numsteps - 1) - 0.5)*filter_width;
435
436 /* Define the filter coefs */
437 if (flam3_temporal_exp == filter_type) {
438
439 for (i=0; i < numsteps; i++) {
440
441 if (filter_exp>=0)
442 slpx = ((double)i+1.0)/numsteps;
443 else
444 slpx = (double)(numsteps - i)/numsteps;
445
446 /* Scale the color based on these values */
447 filter[i] = pow(slpx,fabs(filter_exp));
448
449 /* Keep the max */
450 if (filter[i]>maxfilt)
451 maxfilt = filter[i];
452 }
453
454 } else if (flam3_temporal_gaussian == filter_type) {
455
456 halfsteps = numsteps/2.0;
457 for (i=0; i < numsteps; i++) {
458
459 /* Gaussian */
460 filter[i] = flam3_spatial_filter(flam3_gaussian_kernel,
461 flam3_spatial_support[flam3_gaussian_kernel]*fabs(i - halfsteps)/halfsteps);
462 /* Keep the max */
463 if (filter[i]>maxfilt)
464 maxfilt = filter[i];
465 }
466
467 } else { // (flam3_temporal_box)
468
469 for (i=0; i < numsteps; i++)
470 filter[i] = 1.0;
471
472 maxfilt = 1.0;
473
474 }
475
476 /* Adjust the filter so that the max is 1.0, and */
477 /* calculate the K2 scaling factor */
478 for (i=0;i<numsteps;i++) {
479 filter[i] /= maxfilt;
480 sumfilt += filter[i];
481 }
482
483 sumfilt /= numsteps;
484
485 *temporal_deltas = deltas;
486 *temporal_filter = filter;
487
488 return(sumfilt);
489 }
490
491