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 "interpolation.h"
20 #include "palettes.h"
21 
adjust_percentage(double in)22 double adjust_percentage(double in) {
23 
24    if (in==0.0)
25       return(0.0);
26    else
27       return(pow(10.0, -log(1.0/in)/log(2)));
28 
29 }
30 
motion_funcs(int funcnum,double timeval)31 double motion_funcs(int funcnum, double timeval) {
32 
33    /* motion funcs should be cyclic, and equal to 0 at integral time values */
34    /* abs peak values should be not be greater than 1                       */
35    if (funcnum==MOTION_SIN) {
36       return (sin(2.0*M_PI*timeval));
37    } else if (funcnum==MOTION_TRIANGLE) {
38       double fr = fmod(timeval,1.0);
39 
40       if (fr<0) fr+= 1.0;
41 
42       if (fr<=.25)
43          fr = 4.0 * fr;
44       else if (fr<=.75)
45          fr = -4.0 * fr + 2.0;
46       else
47          fr = 4.0 * fr - 4.0;
48 
49       return(fr);
50    } else { //if (funcnum==MOTION_HILL) {
51       return( (1.0-cos(2.0*M_PI*timeval)) * 0.5);
52    }
53 
54 }
55 
det_matrix(double s[2][2])56 double det_matrix(double s[2][2]) {
57    return s[0][0] * s[1][1] - s[0][1] * s[1][0];
58 }
59 
id_matrix(double s[3][2])60 int id_matrix(double s[3][2]) {
61    return
62       (s[0][0] == 1.0) &&
63       (s[0][1] == 0.0) &&
64       (s[1][0] == 0.0) &&
65       (s[1][1] == 1.0) &&
66       (s[2][0] == 0.0) &&
67       (s[2][1] == 0.0);
68 }
69 
zero_matrix(double s[3][2])70 int zero_matrix(double s[3][2]) {
71    return
72       (s[0][0] == 0.0) &&
73       (s[0][1] == 0.0) &&
74       (s[1][0] == 0.0) &&
75       (s[1][1] == 0.0) &&
76       (s[2][0] == 0.0) &&
77       (s[2][1] == 0.0);
78 }
79 
copy_matrix(double to[3][2],double from[3][2])80 void copy_matrix(double to[3][2], double from[3][2]) {
81 
82    to[0][0] = from[0][0];
83    to[0][1] = from[0][1];
84    to[1][0] = from[1][0];
85    to[1][1] = from[1][1];
86    to[2][0] = from[2][0];
87    to[2][1] = from[2][1];
88 }
89 
90 
clear_matrix(double m[3][2])91 void clear_matrix(double m[3][2]) {
92    m[0][0] = 0.0;
93    m[0][1] = 0.0;
94    m[1][0] = 0.0;
95    m[1][1] = 0.0;
96    m[2][0] = 0.0;
97    m[2][1] = 0.0;
98 }
99 
sum_matrix(double s,double m1[3][2],double m2[3][2])100 void sum_matrix(double s, double m1[3][2], double m2[3][2]) {
101 
102    m2[0][0] += s * m1[0][0];
103    m2[0][1] += s * m1[0][1];
104    m2[1][0] += s * m1[1][0];
105    m2[1][1] += s * m1[1][1];
106    m2[2][0] += s * m1[2][0];
107    m2[2][1] += s * m1[2][1];
108 }
109 
mult_matrix(double s1[2][2],double s2[2][2],double d[2][2])110 void mult_matrix(double s1[2][2], double s2[2][2], double d[2][2]) {
111    d[0][0] = s1[0][0] * s2[0][0] + s1[1][0] * s2[0][1];
112    d[1][0] = s1[0][0] * s2[1][0] + s1[1][0] * s2[1][1];
113    d[0][1] = s1[0][1] * s2[0][0] + s1[1][1] * s2[0][1];
114    d[1][1] = s1[0][1] * s2[1][0] + s1[1][1] * s2[1][1];
115 }
116 
compare_xforms(const void * av,const void * bv)117 int compare_xforms(const void *av, const void *bv) {
118    flam3_xform *a = (flam3_xform *) av;
119    flam3_xform *b = (flam3_xform *) bv;
120    double aa[2][2];
121    double bb[2][2];
122    double ad, bd;
123 
124    aa[0][0] = a->c[0][0];
125    aa[0][1] = a->c[0][1];
126    aa[1][0] = a->c[1][0];
127    aa[1][1] = a->c[1][1];
128    bb[0][0] = b->c[0][0];
129    bb[0][1] = b->c[0][1];
130    bb[1][0] = b->c[1][0];
131    bb[1][1] = b->c[1][1];
132    ad = det_matrix(aa);
133    bd = det_matrix(bb);
134 
135    if (a->color_speed > b->color_speed) return 1;
136    if (a->color_speed < b->color_speed) return -1;
137    if (a->color_speed) {
138       if (ad < 0) return -1;
139       if (bd < 0) return 1;
140       ad = atan2(a->c[0][0], a->c[0][1]);
141       bd = atan2(b->c[0][0], b->c[0][1]);
142    }
143 
144    if (ad < bd) return -1;
145    if (ad > bd) return 1;
146    return 0;
147 }
148 
interpolate_cmap(flam3_palette cmap,double blend,int index0,double hue0,int index1,double hue1)149 void interpolate_cmap(flam3_palette cmap, double blend,
150                       int index0, double hue0, int index1, double hue1) {
151 
152    flam3_palette p0,p1;
153    int i, j, rcode;
154 
155    rcode = flam3_get_palette(index0, p0, hue0);
156    if (rcode<0)
157       fprintf(stderr,"unable to retrieve palette %d, setting to white\n", index0);
158    rcode = flam3_get_palette(index1, p1, hue1);
159    if (rcode<0)
160       fprintf(stderr,"unable to retrieve palette %d, setting to white\n", index1);
161 
162    for (i = 0; i < 256; i++) {
163       double t[5], s[5];
164 
165       rgb2hsv(p0[i].color, s);
166       rgb2hsv(p1[i].color, t);
167 
168       s[3] = p0[i].color[3];
169       t[3] = p1[i].color[3];
170 
171       s[4] = p0[i].index;
172       t[4] = p1[i].index;
173 
174       /* take the shorter way around the hue circle */
175 //      if (0 == i) {
176 //	fprintf(stderr, "xxx interpolating between hues, %g %g\n", s[0], t[0]);
177  //     }
178 
179       /* Correct the first hue to go the short way around */
180       if ((s[0] - t[0]) > 3.0)  /* first hue much bigger than second hue */
181         s[0] -= 6.0;
182       if ((s[0] - t[0]) < -3.0)  /* first hue much smaller than second hue */
183 	s[0] += 6.0;
184 
185       for (j = 0; j < 5; j++)
186          t[j] = ((1.0-blend) * s[j]) + (blend * t[j]);
187 
188       hsv2rgb(t, cmap[i].color);
189       cmap[i].color[3] = t[3];
190       cmap[i].index = t[4];
191    }
192 }
193 
interp_and_convert_back(double * c,int ncps,int xfi,double cxang[4][2],double cxmag[4][2],double cxtrn[4][2],double store_array[3][2])194 void interp_and_convert_back(double *c, int ncps, int xfi, double cxang[4][2],
195                              double cxmag[4][2], double cxtrn[4][2],double store_array[3][2]) {
196 
197    int i,col;
198 
199    double accang[2],accmag[2];
200    double expmag;
201    int accmode[2];
202 
203    accang[0] = 0.0;
204    accang[1] = 0.0;
205    accmag[0] = 0.0;
206    accmag[1] = 0.0;
207 
208    accmode[0]=accmode[1]=0;
209 
210    /* accumulation mode defaults to logarithmic, but in special */
211    /* cases we want to switch to linear accumulation            */
212    for (col=0; col<2; col++) {
213       for (i=0; i<ncps; i++) {
214          if (log(cxmag[i][col])<-10)
215             accmode[col]=1; // Mode set to linear interp
216       }
217    }
218 
219    for (i=0; i<ncps; i++) {
220       for (col=0; col<2; col++) {
221 
222          accang[col] += c[i] * cxang[i][col];
223 
224          if (accmode[col]==0)
225             accmag[col] += c[i] * log(cxmag[i][col]);
226          else
227             accmag[col] += c[i] * (cxmag[i][col]);
228 
229          /* translation is ready to go */
230          store_array[2][col] += c[i] * cxtrn[i][col];
231       }
232    }
233 
234    /* Convert the angle back to rectangular */
235    for (col=0;col<2;col++) {
236       if (accmode[col]==0)
237          expmag = exp(accmag[col]);
238       else
239          expmag = accmag[col];
240 
241       store_array[col][0] = expmag * cos(accang[col]);
242       store_array[col][1] = expmag * sin(accang[col]);
243    }
244 
245 }
246 
convert_linear_to_polar(flam3_genome * cp,int ncps,int xfi,int cflag,double cxang[4][2],double cxmag[4][2],double cxtrn[4][2])247 void convert_linear_to_polar(flam3_genome *cp, int ncps, int xfi, int cflag,
248                              double cxang[4][2], double cxmag[4][2], double cxtrn[4][2]) {
249 
250    double c1[2],d,t,refang;
251    int col,k;
252    int zlm[2];
253 
254    for (k=0; k<ncps;k++) {
255 
256       /* Establish the angles and magnitudes for each component */
257       /* Keep translation linear */
258       zlm[0]=zlm[1]=0;
259       for (col=0;col<2;col++) {
260 
261          if (cflag==0) {
262             c1[0] = cp[k].xform[xfi].c[col][0];
263             c1[1] = cp[k].xform[xfi].c[col][1];
264             t = cp[k].xform[xfi].c[2][col];
265          } else {
266             c1[0] = cp[k].xform[xfi].post[col][0];
267             c1[1] = cp[k].xform[xfi].post[col][1];
268             t = cp[k].xform[xfi].post[2][col];
269          }
270 
271          cxang[k][col] = atan2(c1[1],c1[0]);
272          cxmag[k][col] = sqrt(c1[0]*c1[0] + c1[1]*c1[1]);
273 
274          if (cxmag[k][col]== 0.0)
275             zlm[col]=1;
276 
277          cxtrn[k][col] = t;
278       }
279 
280       if (zlm[0]==1 && zlm[1]==0)
281          cxang[k][0] = cxang[k][1];
282       else if (zlm[0]==0 && zlm[1]==1)
283          cxang[k][1] = cxang[k][0];
284 
285    }
286 
287    /* Make sure the rotation is the shorter direction around the circle */
288    /* by adjusting each angle in succession, and rotate clockwise if 180 degrees */
289    for (col=0; col<2; col++) {
290       for (k=1;k<ncps;k++) {
291 
292          /* Adjust angles differently if we have an asymmetric case */
293          if (cp[k].xform[xfi].wind[col]>0 && cflag==0) {
294 
295             /* Adjust the angles to make sure that it's within wind:wind+2pi */
296             refang = cp[k].xform[xfi].wind[col] - 2*M_PI;
297 
298             /* Make sure both angles are within [refang refang+2*pi] */
299             while(cxang[k-1][col] < refang)
300                   cxang[k-1][col] += 2*M_PI;
301 
302             while(cxang[k-1][col] > refang + 2*M_PI)
303                   cxang[k-1][col] -= 2*M_PI;
304 
305             while(cxang[k][col] < refang)
306                   cxang[k][col] += 2*M_PI;
307 
308             while(cxang[k][col] > refang + 2*M_PI)
309                   cxang[k][col] -= 2*M_PI;
310 
311          } else {
312 
313             /* Normal way of adjusting angles */
314             d = cxang[k][col]-cxang[k-1][col];
315 
316             /* Adjust to avoid the -pi/pi discontinuity */
317             if (d > M_PI+EPS)
318                cxang[k][col] -= 2*M_PI;
319             else if (d < -(M_PI-EPS) ) /* Forces clockwise rotation at 180 */
320                cxang[k][col] += 2*M_PI;
321          }
322       }
323    }
324 }
325 
interpolate_catmull_rom(flam3_genome cps[],double t,flam3_genome * result)326 void interpolate_catmull_rom(flam3_genome cps[], double t, flam3_genome *result) {
327    double t2 = t * t;
328    double t3 = t2 * t;
329    double cmc[4];
330 
331    cmc[0] = (2*t2 - t - t3) / 2;
332    cmc[1] = (3*t3 - 5*t2 + 2) / 2;
333    cmc[2] = (4*t2 - 3*t3 + t) / 2;
334    cmc[3] = (t3 - t2) / 2;
335 
336    flam3_interpolate_n(result, 4, cps, cmc, 0);
337 }
338 
smoother(double t)339 double smoother(double t) {
340   return 3*t*t - 2*t*t*t;
341 }
342 
get_stagger_coef(double t,double stagger_prc,int num_xforms,int this_xform)343 double get_stagger_coef(double t, double stagger_prc, int num_xforms, int this_xform) {
344 
345    /* max_stag is the spacing between xform start times if stagger_prc = 1.0 */
346    double max_stag = (double)(num_xforms-1)/num_xforms;
347 
348    /* scale the spacing by stagger_prc */
349    double stag_scaled = stagger_prc * max_stag;
350 
351    /* t ranges from 1 to 0 (the contribution of cp[0] to the blend) */
352    /* the first line below makes the first xform interpolate first */
353    /* the second line makes the last xform interpolate first */
354    double st = stag_scaled * (num_xforms - 1 - this_xform) / (num_xforms-1);
355 //   double st = stag_scaled * (this_xform) / (num_xforms-1);
356    double et = st + (1-stag_scaled);
357 
358 //   printf("t=%f xf:%d st=%f et=%f : : %f\n",t,this_xform,st,et,smoother((t-st)/(1-stag_scaled)));
359 
360    if (t <= st)
361       return (0);
362    else if (t >= et)
363       return (1);
364    else
365       return ( smoother((t-st)/(1-stag_scaled)) );
366 
367 }
368 
369 
370 
371 /* all cpi and result must be aligned (have the same number of xforms,
372    and have final xform in the same slot) */
flam3_interpolate_n(flam3_genome * result,int ncp,flam3_genome * cpi,double * c,double stagger)373 void flam3_interpolate_n(flam3_genome *result, int ncp,
374           flam3_genome *cpi, double *c, double stagger) {
375    int i, j, k, l, numstd;
376 
377    if (flam3_palette_interpolation_sweep != cpi[0].palette_interpolation) {
378 
379       /* rgb, hsv or hsv_circular modes. */
380       double rgb_fraction = 0.0;  /* Assume that we are in plain hsv mode */
381       if (flam3_palette_interpolation_rgb == cpi[0].palette_interpolation)
382          rgb_fraction = 1.0;  /* All RGB output */
383       else if (flam3_palette_interpolation_hsv_circular == cpi[0].palette_interpolation)
384          rgb_fraction = cpi[0].hsv_rgb_palette_blend;
385 
386       for (i = 0; i < 256; i++) {
387          double col_rgb[3], col_hsv[3];
388          double new_rgb[3] = {0, 0, 0};
389          double new_hsv[3] = {0, 0, 0};
390          double new_count = 0, new_index = 0;
391          int alpha1 = 1;
392 
393          /* Loop over each control point's color at this index */
394          for (k = 0; k < ncp; k++) {
395 
396             /* Convert to hsv */
397             rgb2hsv(cpi[k].palette[i].color, col_hsv);
398 
399 	    /* Store the rgb */
400             for (l = 0; l < 3; l++)
401                col_rgb[l] = cpi[k].palette[i].color[l];
402 
403 	    if (2 == ncp && k == 0 && cpi[0].palette_interpolation == flam3_palette_interpolation_hsv_circular) {
404                /* only adjust the first coordinate based on the other control point's hue */
405                double second_color[3];
406                rgb2hsv(cpi[1].palette[i].color, second_color);
407 
408                /* Adjust the hue so that we go the shorter direction around the circle */
409                if ((second_color[0] - col_hsv[0]) > 3.0) {
410                   col_hsv[0] += 6.0;
411                } else if ((second_color[0] - col_hsv[0]) < -3.0) {
412                   col_hsv[0] -= 6.0;
413                }
414             }
415 
416             for (j = 0; j < 3; j++) {
417                new_rgb[j] += c[k] * col_rgb[j];
418                new_hsv[j] += c[k] * col_hsv[j];
419             }
420 
421             /* Compute the other two components of the color (count and index) */
422             new_count += c[k] * cpi[k].palette[i].color[3];
423             if (cpi[k].palette[i].color[3] != 1.0)
424                alpha1 = 0;
425             new_index += c[k] * cpi[k].palette[i].index;
426 
427          }
428 
429          if (alpha1 == 1)
430             new_count = 1.0;
431 
432          /* Convert the new hsv coord to back rgb */
433          double new_hsv_rgb[3];
434          hsv2rgb(new_hsv, new_hsv_rgb);
435 
436          /* Store the interpolated color in the new palette */
437          for (l = 0; l < 3; l++)
438             result->palette[i].color[l] = rgb_fraction * new_rgb[l] + (1.0-rgb_fraction) * new_hsv_rgb[l];
439 
440          result->palette[i].color[3] = new_count;
441          result->palette[i].index = new_index;
442 
443          /* Clip the new color appropriately */
444          for (j = 0; j < 4; j++) {
445             if (result->palette[i].color[j] < 0.0)
446                result->palette[i].color[j] = 0.0;
447             if (result->palette[i].color[j] > 1.0)
448                result->palette[i].color[j] = 1.0;
449          }
450 
451          if (result->palette[i].index < 0.0)
452             result->palette[i].index = 0.0;
453          if (result->palette[i].index > 255.0)
454             result->palette[i].index = 255.0;
455       }
456    } else {
457       /* Sweep - not the best option for float indices */
458       for (i = 0; i < 256; i++) {
459          j = (i < (256 * c[0])) ? 0 : 1;
460          result->palette[i] = cpi[j].palette[i];
461       }
462    }
463 
464    result->palette_index = flam3_palette_random;
465    result->symmetry = 0;
466    result->spatial_filter_select = cpi[0].spatial_filter_select;
467    result->temporal_filter_type = cpi[0].temporal_filter_type;
468    result->palette_mode = cpi[0].palette_mode;
469 
470    result->interpolation_type = cpi[0].interpolation_type;
471    result->palette_interpolation = cpi[0].palette_interpolation;
472    result->hsv_rgb_palette_blend = cpi[0].hsv_rgb_palette_blend;
473    INTERP(brightness);
474    INTERP(contrast);
475    INTERP(highlight_power);
476    INTERP(gamma);
477    INTERP(vibrancy);
478    INTERP(hue_rotation);
479    INTERI(width);
480    INTERI(height);
481    INTERI(spatial_oversample);
482    INTERP(center[0]);
483    INTERP(center[1]);
484    INTERP(rot_center[0]);
485    INTERP(rot_center[1]);
486    INTERP(background[0]);
487    INTERP(background[1]);
488    INTERP(background[2]);
489    INTERP(pixels_per_unit);
490    INTERP(spatial_filter_radius);
491    INTERP(temporal_filter_exp);
492    INTERP(temporal_filter_width);
493    INTERP(sample_density);
494    INTERP(zoom);
495    INTERP(rotate);
496    INTERI(nbatches);
497    INTERI(ntemporal_samples);
498    INTERP(estimator);
499    INTERP(estimator_minimum);
500    INTERP(estimator_curve);
501    INTERP(gam_lin_thresh);
502 
503    /* Interpolate the chaos array */
504    numstd = cpi[0].num_xforms - (cpi[0].final_xform_index >= 0);
505    for (i=0;i<numstd;i++) {
506       for (j=0;j<numstd;j++) {
507          INTERP(chaos[i][j]);
508          if (result->chaos[i][j]<0) result->chaos[i][j]=0;
509          //chaos can be > 1
510          //if (result->chaos[i][j]>1) result->chaos[i][j]=1.0;
511       }
512    }
513 
514    /* Interpolate each xform */
515    for (i = 0; i < cpi[0].num_xforms; i++) {
516 
517       double csave[2] = {0, 0};
518       double td;
519       int all_id;
520       int nx = cpi[0].num_xforms-(cpi[0].final_xform_index>=0);
521 
522       if (ncp==2 && stagger>0 && i!=cpi[0].final_xform_index) {
523          csave[0] = c[0];
524          csave[1] = c[1];
525          c[0] = get_stagger_coef(csave[0],stagger,nx,i);
526          c[1] = 1.0-c[0];
527       }
528 
529 
530       INTERP(xform[i].density);
531       td = result->xform[i].density;
532       result->xform[i].density = (td < 0.0) ? 0.0 : td;
533       INTERP(xform[i].color);
534       if (result->xform[i].color<0) result->xform[i].color=0;
535       if (result->xform[i].color>1) result->xform[i].color=1;
536 
537       INTERP(xform[i].color_speed);
538       if (result->xform[i].color_speed<0) result->xform[i].color_speed=0;
539       if (result->xform[i].color_speed>1) result->xform[i].color_speed=1;
540 
541       INTERP(xform[i].opacity);
542       INTERP(xform[i].animate);
543       INTERP(xform[i].blob_low);
544       INTERP(xform[i].blob_high);
545       INTERP(xform[i].blob_waves);
546       INTERP(xform[i].pdj_a);
547       INTERP(xform[i].pdj_b);
548       INTERP(xform[i].pdj_c);
549       INTERP(xform[i].pdj_d);
550       INTERP(xform[i].fan2_x);
551       INTERP(xform[i].fan2_y);
552       INTERP(xform[i].rings2_val);
553       INTERP(xform[i].perspective_angle);
554       INTERP(xform[i].perspective_dist);
555       INTERP(xform[i].julian_power);
556       INTERP(xform[i].julian_dist);
557       INTERP(xform[i].juliascope_power);
558       INTERP(xform[i].juliascope_dist);
559       INTERP(xform[i].radial_blur_angle);
560       INTERP(xform[i].pie_slices);
561       INTERP(xform[i].pie_rotation);
562       INTERP(xform[i].pie_thickness);
563       INTERP(xform[i].ngon_sides);
564       INTERP(xform[i].ngon_power);
565       INTERP(xform[i].ngon_circle);
566       INTERP(xform[i].ngon_corners);
567       INTERP(xform[i].curl_c1);
568       INTERP(xform[i].curl_c2);
569       INTERP(xform[i].rectangles_x);
570       INTERP(xform[i].rectangles_y);
571       INTERP(xform[i].amw_amp);
572       INTERP(xform[i].disc2_rot);
573       INTERP(xform[i].disc2_twist);
574       INTERP(xform[i].super_shape_rnd);
575       INTERP(xform[i].super_shape_m);
576       INTERP(xform[i].super_shape_n1);
577       INTERP(xform[i].super_shape_n2);
578       INTERP(xform[i].super_shape_n3);
579       INTERP(xform[i].super_shape_holes);
580       INTERP(xform[i].flower_petals);
581       INTERP(xform[i].flower_holes);
582       INTERP(xform[i].conic_eccentricity);
583       INTERP(xform[i].conic_holes);
584       INTERP(xform[i].parabola_height);
585       INTERP(xform[i].parabola_width);
586       INTERP(xform[i].bent2_x);
587       INTERP(xform[i].bent2_y);
588       INTERP(xform[i].bipolar_shift);
589       INTERP(xform[i].cell_size);
590       INTERP(xform[i].cpow_r);
591       INTERP(xform[i].cpow_i);
592       INTERP(xform[i].cpow_power);
593       INTERP(xform[i].curve_xamp);
594       INTERP(xform[i].curve_yamp);
595       INTERP(xform[i].curve_xlength);
596       INTERP(xform[i].curve_ylength);
597       INTERP(xform[i].escher_beta);
598       INTERP(xform[i].lazysusan_x);
599       INTERP(xform[i].lazysusan_y);
600       INTERP(xform[i].lazysusan_twist);
601       INTERP(xform[i].lazysusan_space);
602       INTERP(xform[i].lazysusan_spin);
603       INTERP(xform[i].modulus_x);
604       INTERP(xform[i].modulus_y);
605       INTERP(xform[i].oscope_separation);
606       INTERP(xform[i].oscope_frequency);
607       INTERP(xform[i].oscope_amplitude);
608       INTERP(xform[i].oscope_damping);
609       INTERP(xform[i].popcorn2_x);
610       INTERP(xform[i].popcorn2_y);
611       INTERP(xform[i].popcorn2_c);
612       INTERP(xform[i].separation_x);
613       INTERP(xform[i].separation_xinside);
614       INTERP(xform[i].separation_y);
615       INTERP(xform[i].separation_yinside);
616       INTERP(xform[i].split_xsize);
617       INTERP(xform[i].split_ysize);
618       INTERP(xform[i].splits_x);
619       INTERP(xform[i].splits_y);
620       INTERP(xform[i].stripes_space);
621       INTERP(xform[i].stripes_warp);
622       INTERP(xform[i].wedge_angle);
623       INTERP(xform[i].wedge_hole);
624       INTERP(xform[i].wedge_count);
625       INTERP(xform[i].wedge_swirl);
626       INTERP(xform[i].wedge_julia_angle);
627       INTERP(xform[i].wedge_julia_count);
628       INTERP(xform[i].wedge_julia_power);
629       INTERP(xform[i].wedge_julia_dist);
630       INTERP(xform[i].wedge_sph_angle);
631       INTERP(xform[i].wedge_sph_hole);
632       INTERP(xform[i].wedge_sph_count);
633       INTERP(xform[i].wedge_sph_swirl);
634       INTERP(xform[i].whorl_inside);
635       INTERP(xform[i].whorl_outside);
636       INTERP(xform[i].waves2_scalex);
637       INTERP(xform[i].waves2_scaley);
638       INTERP(xform[i].waves2_freqx);
639       INTERP(xform[i].waves2_freqy);
640       INTERP(xform[i].auger_sym);
641       INTERP(xform[i].auger_weight);
642       INTERP(xform[i].auger_freq);
643       INTERP(xform[i].auger_scale);
644       INTERP(xform[i].flux_spread);
645       INTERP(xform[i].mobius_re_a);
646       INTERP(xform[i].mobius_im_a);
647       INTERP(xform[i].mobius_re_b);
648       INTERP(xform[i].mobius_im_b);
649       INTERP(xform[i].mobius_re_c);
650       INTERP(xform[i].mobius_im_c);
651       INTERP(xform[i].mobius_re_d);
652       INTERP(xform[i].mobius_im_d);
653 
654       for (j = 0; j < flam3_nvariations; j++)
655          INTERP(xform[i].var[j]);
656 
657       if (flam3_inttype_log == cpi[0].interpolation_type) {
658          double cxmag[4][2];  // XXX why only 4? should be ncp
659          double cxang[4][2];
660          double cxtrn[4][2];
661 
662          /* affine part */
663          clear_matrix(result->xform[i].c);
664          convert_linear_to_polar(cpi,ncp,i,0,cxang,cxmag,cxtrn);
665          interp_and_convert_back(c, ncp, i, cxang, cxmag, cxtrn,result->xform[i].c);
666 
667          /* post part */
668          all_id = 1;
669          for (k=0; k<ncp; k++)
670             all_id &= id_matrix(cpi[k].xform[i].post);
671 
672          clear_matrix(result->xform[i].post);
673          if (all_id) {
674             result->xform[i].post[0][0] = 1.0;
675             result->xform[i].post[1][1] = 1.0;
676          } else {
677             convert_linear_to_polar(cpi,ncp,i,1,cxang,cxmag,cxtrn);
678             interp_and_convert_back(c, ncp, i, cxang, cxmag, cxtrn,result->xform[i].post);
679          }
680 
681       } else {
682 
683          /* Interpolate c matrix & post */
684          clear_matrix(result->xform[i].c);
685          clear_matrix(result->xform[i].post);
686          all_id = 1;
687          for (k = 0; k < ncp; k++) {
688             sum_matrix(c[k], cpi[k].xform[i].c, result->xform[i].c);
689             sum_matrix(c[k], cpi[k].xform[i].post, result->xform[i].post);
690 
691             all_id &= id_matrix(cpi[k].xform[i].post);
692 
693          }
694          if (all_id) {
695             clear_matrix(result->xform[i].post);
696             result->xform[i].post[0][0] = 1.0;
697             result->xform[i].post[1][1] = 1.0;
698          }
699       }
700 
701       if (ncp==2 && stagger>0 && i!=cpi[0].final_xform_index) {
702          c[0] = csave[0];
703          c[1] = csave[1];
704       }
705 
706    }
707 
708 }
709 
establish_asymmetric_refangles(flam3_genome * cp,int ncps)710 void establish_asymmetric_refangles(flam3_genome *cp, int ncps) {
711 
712    int k, xfi, col;
713 
714    double cxang[4][2],d,c1[2];
715 
716    for (xfi=0; xfi<cp[0].num_xforms; xfi++) {
717 
718       /* Final xforms don't rotate regardless of their symmetry */
719       if (cp[0].final_xform_enable==1 && xfi==cp[0].final_xform_index)
720          continue;
721 
722       for (k=0; k<ncps;k++) {
723 
724          /* Establish the angle for each component */
725          /* Should potentially functionalize */
726          for (col=0;col<2;col++) {
727 
728             c1[0] = cp[k].xform[xfi].c[col][0];
729             c1[1] = cp[k].xform[xfi].c[col][1];
730 
731             cxang[k][col] = atan2(c1[1],c1[0]);
732          }
733       }
734 
735       for (k=1; k<ncps; k++) {
736 
737          for (col=0;col<2;col++) {
738 
739             int sym0,sym1;
740             int padsymflag;
741 
742             d = cxang[k][col]-cxang[k-1][col];
743 
744             /* Adjust to avoid the -pi/pi discontinuity */
745             if (d > M_PI+EPS)
746             cxang[k][col] -= 2*M_PI;
747             else if (d < -(M_PI-EPS) )
748             cxang[k][col] += 2*M_PI;
749 
750             /* If this is an asymmetric case, store the NON-symmetric angle    */
751             /* Check them pairwise and store the reference angle in the second */
752             /* to avoid overwriting if asymmetric on both sides                */
753             padsymflag = 0;
754 
755             sym0 = (cp[k-1].xform[xfi].animate==0 || (cp[k-1].xform[xfi].padding==1 && padsymflag));
756             sym1 = (cp[k].xform[xfi].animate==0 || (cp[k].xform[xfi].padding==1 && padsymflag));
757 
758             if ( sym1 && !sym0 )
759                cp[k].xform[xfi].wind[col] = cxang[k-1][col] + 2*M_PI;
760             else if ( sym0 && !sym1 )
761                cp[k].xform[xfi].wind[col] = cxang[k][col] + 2*M_PI;
762 
763          }
764       }
765    }
766 }
767 
flam3_align(flam3_genome * dst,flam3_genome * src,int nsrc)768 void flam3_align(flam3_genome *dst, flam3_genome *src, int nsrc) {
769    int i, tfx, tnx, max_nx = 0, max_fx = 0;
770    int already_aligned=1;
771    int xf,j;
772    int ii,fnd;
773    double normed;
774    int usethisone;
775 
776    usethisone = (nsrc/2) - 1;
777 
778    max_nx = src[0].num_xforms - (src[0].final_xform_index >= 0);
779    max_fx = src[0].final_xform_enable;
780 
781    for (i = 1; i < nsrc; i++) {
782       tnx = src[i].num_xforms - (src[i].final_xform_index >= 0);
783       if (max_nx != tnx) {
784          already_aligned = 0;
785          if (tnx > max_nx) max_nx = tnx;
786       }
787 
788       tfx = src[i].final_xform_enable;
789       if (max_fx != tfx) {
790          already_aligned = 0;
791          max_fx |= tfx;
792       }
793    }
794 
795    /* Pad the cps to equal xforms */
796    for (i = 0; i < nsrc; i++) {
797       flam3_copyx(&dst[i], &src[i], max_nx, max_fx);
798    }
799 
800    /* Skip if this genome is compatibility mode */
801    if (dst[usethisone].interpolation_type == flam3_inttype_compat ||
802          dst[usethisone].interpolation_type == flam3_inttype_older)
803       return;
804 
805 
806    /* Check to see if there's a parametric variation present in one xform   */
807    /* but not in an aligned xform.  If this is the case, use the parameters */
808    /* from the xform with the variation as the defaults for the blank one.  */
809 
810    /* All genomes will have the same number of xforms at this point */
811    /* num = max_nx + max_fx */
812    for (i = 0; i<nsrc; i++) {
813 
814 
815       for (xf = 0; xf<max_nx+max_fx; xf++) {
816 
817          /* Loop over the variations to see which of them are set to 0 */
818          /* Note that there are no parametric variations < 23 */
819          for (j = 23; j < flam3_nvariations; j++) {
820 
821               if (dst[i].xform[xf].var[j]==0) {
822 
823                  if (i>0) {
824 
825                     /* Check to see if the prior genome's xform is populated */
826                     if (dst[i-1].xform[xf].var[j] != 0) {
827 
828                        /* Copy the prior genome's parameters and continue */
829                        flam3_copy_params(&(dst[i].xform[xf]), &(dst[i-1].xform[xf]), j);
830                        continue;
831                     }
832 
833                  } else if (i<nsrc-1) {
834 
835                     /* Check to see if the next genome's xform is populated */
836                     if (dst[i+1].xform[xf].var[j] != 0) {
837 
838                        /* Copy the next genome's parameters and continue */
839                        flam3_copy_params(&(dst[i].xform[xf]), &(dst[i+1].xform[xf]), j);
840                        continue;
841                     }
842                  }
843               }
844           } /* variations */
845 
846           if (dst[i].xform[xf].padding == 1 && !already_aligned) {
847 
848              /* This is a new xform.  Let's see if we can choose a better 'identity' xform. */
849              /* Check the neighbors to see if any of these variations are used: */
850              /* rings2, fan2, blob, perspective, julian, juliascope, ngon, curl, super_shape, split */
851              /* If so, we can use a better starting point for these */
852 
853              /* Remove linear from the list */
854              dst[i].xform[xf].var[0] = 0.0;
855 
856              /* Look through all of the 'companion' xforms to see if we get a match on any of these */
857              fnd=0;
858 
859              /* Only do the next substitution for log interpolation */
860              if ( (i==0 && dst[i].interpolation_type == flam3_inttype_log)
861                   || (i>0 && dst[i-1].interpolation_type==flam3_inttype_log) ) {
862 
863              for (ii=-1; ii<=1; ii+=2) {
864 
865                 /* Skip if out of bounds */
866                 if (i+ii<0 || i+ii>=nsrc)
867                    continue;
868 
869                 /* Skip if this is also padding */
870                 if (dst[i+ii].xform[xf].padding==1)
871                    continue;
872 
873                 /* Spherical / Ngon (trumps all others due to holes)       */
874                 /* Interpolate these against a 180 degree rotated identity */
875                 /* with weight -1.                                         */
876                 /* Added JULIAN/JULIASCOPE to get rid of black wedges      */
877                 if (dst[i+ii].xform[xf].var[VAR_SPHERICAL]>0 ||
878                       dst[i+ii].xform[xf].var[VAR_NGON]>0 ||
879                       dst[i+ii].xform[xf].var[VAR_JULIAN]>0 ||
880                       dst[i+ii].xform[xf].var[VAR_JULIASCOPE]>0 ||
881                       dst[i+ii].xform[xf].var[VAR_POLAR]>0 ||
882                       dst[i+ii].xform[xf].var[VAR_WEDGE_SPH]>0 ||
883                       dst[i+ii].xform[xf].var[VAR_WEDGE_JULIA]>0) {
884 
885                    dst[i].xform[xf].var[VAR_LINEAR] = -1.0;
886                    /* Set the coefs appropriately */
887                    dst[i].xform[xf].c[0][0] = -1.0;
888                    dst[i].xform[xf].c[0][1] = 0.0;
889                    dst[i].xform[xf].c[1][0] = 0.0;
890                    dst[i].xform[xf].c[1][1] = -1.0;
891                    dst[i].xform[xf].c[2][0] = 0.0;
892                    dst[i].xform[xf].c[2][1] = 0.0;
893                    fnd=-1;
894                 }
895              }
896 
897              }
898 
899              if (fnd==0) {
900 
901                 for (ii=-1; ii<=1; ii+=2) {
902 
903                    /* Skip if out of bounds */
904                    if (i+ii<0 || i+ii>=nsrc)
905                       continue;
906 
907                    /* Skip if also padding */
908                    if (dst[i+ii].xform[xf].padding==1)
909                       continue;
910 
911                    /* Rectangles */
912                    if (dst[i+ii].xform[xf].var[VAR_RECTANGLES]>0) {
913                       dst[i].xform[xf].var[VAR_RECTANGLES] = 1.0;
914                       dst[i].xform[xf].rectangles_x = 0.0;
915                       dst[i].xform[xf].rectangles_y = 0.0;
916                       fnd++;
917                    }
918 
919                    /* Rings 2 */
920                    if (dst[i+ii].xform[xf].var[VAR_RINGS2]>0) {
921                       dst[i].xform[xf].var[VAR_RINGS2] = 1.0;
922                       dst[i].xform[xf].rings2_val = 0.0;
923                       fnd++;
924                    }
925 
926                    /* Fan 2 */
927                    if (dst[i+ii].xform[xf].var[VAR_FAN2]>0) {
928                       dst[i].xform[xf].var[VAR_FAN2] = 1.0;
929                       dst[i].xform[xf].fan2_x = 0.0;
930                       dst[i].xform[xf].fan2_y = 0.0;
931                       fnd++;
932                    }
933 
934                    /* Blob */
935                    if (dst[i+ii].xform[xf].var[VAR_BLOB]>0) {
936                       dst[i].xform[xf].var[VAR_BLOB] = 1.0;
937                       dst[i].xform[xf].blob_low = 1.0;
938                       dst[i].xform[xf].blob_high = 1.0;
939                       dst[i].xform[xf].blob_waves = 1.0;
940                       fnd++;
941                    }
942 
943                    /* Perspective */
944                    if (dst[i+ii].xform[xf].var[VAR_PERSPECTIVE]>0) {
945                       dst[i].xform[xf].var[VAR_PERSPECTIVE] = 1.0;
946                       dst[i].xform[xf].perspective_angle = 0.0;
947                       /* Keep the perspective distance as-is */
948                       fnd++;
949                    }
950 
951                    /* Curl */
952                    if (dst[i+ii].xform[xf].var[VAR_CURL]>0) {
953                       dst[i].xform[xf].var[VAR_CURL] = 1.0;
954                       dst[i].xform[xf].curl_c1 = 0.0;
955                       dst[i].xform[xf].curl_c2 = 0.0;
956                       fnd++;
957                    }
958 
959                    /* Super-Shape */
960                    if (dst[i+ii].xform[xf].var[VAR_SUPER_SHAPE]>0) {
961                       dst[i].xform[xf].var[VAR_SUPER_SHAPE] = 1.0;
962                       /* Keep supershape_m the same */
963                       dst[i].xform[xf].super_shape_n1 = 2.0;
964                       dst[i].xform[xf].super_shape_n2 = 2.0;
965                       dst[i].xform[xf].super_shape_n3 = 2.0;
966                       dst[i].xform[xf].super_shape_rnd = 0.0;
967                       dst[i].xform[xf].super_shape_holes = 0.0;
968                       fnd++;
969                    }
970                 }
971              }
972 
973              /* If we didn't have any matches with those, */
974              /* try the affine ones, fan and rings        */
975              if (fnd==0) {
976 
977                 for (ii=-1; ii<=1; ii+=2) {
978 
979                    /* Skip if out of bounds */
980                    if (i+ii<0 || i+ii>=nsrc)
981                       continue;
982 
983                    /* Skip if also a padding xform */
984                    if (dst[i+ii].xform[xf].padding==1)
985                       continue;
986 
987                    /* Fan */
988                    if (dst[i+ii].xform[xf].var[VAR_FAN]>0) {
989                       dst[i].xform[xf].var[VAR_FAN] = 1.0;
990                       fnd++;
991                    }
992 
993                    /* Rings */
994                    if (dst[i+ii].xform[xf].var[VAR_RINGS]>0) {
995                       dst[i].xform[xf].var[VAR_RINGS] = 1.0;
996                       fnd++;
997                    }
998 
999                 }
1000 
1001                 if (fnd>0) {
1002                    /* Set the coefs appropriately */
1003                    dst[i].xform[xf].c[0][0] = 0.0;
1004                    dst[i].xform[xf].c[0][1] = 1.0;
1005                    dst[i].xform[xf].c[1][0] = 1.0;
1006                    dst[i].xform[xf].c[1][1] = 0.0;
1007                    dst[i].xform[xf].c[2][0] = 0.0;
1008                    dst[i].xform[xf].c[2][1] = 0.0;
1009                 }
1010              }
1011 
1012              /* If we still have no matches, switch back to linear */
1013              if (fnd==0)
1014 
1015                 dst[i].xform[xf].var[VAR_LINEAR] = 1.0;
1016 
1017              else if (fnd>0) {
1018 
1019                 /* Otherwise, go through and normalize the weights. */
1020                 normed = 0.0;
1021                 for (j = 0; j < flam3_nvariations; j++)
1022                    normed += dst[i].xform[xf].var[j];
1023 
1024                 for (j = 0; j < flam3_nvariations; j++)
1025                    dst[i].xform[xf].var[j] /= normed;
1026 
1027              }
1028           }
1029        } /* xforms */
1030    } /* genomes */
1031 
1032 }
1033 
1034