hide or rename any symbols that are likely to conflict with other
[imager.git] / filters.c
CommitLineData
02d1d628 1#include "image.h"
290bdf77 2#include "imagei.h"
02d1d628
AMH
3#include <stdlib.h>
4#include <math.h>
5
6
7/*
8=head1 NAME
9
10filters.c - implements filters that operate on images
11
12=head1 SYNOPSIS
13
14
15 i_contrast(im, 0.8);
16 i_hardinvert(im);
b6381851 17 i_unsharp_mask(im, 2.0, 1.0);
02d1d628
AMH
18 // and more
19
20=head1 DESCRIPTION
21
22filters.c implements basic filters for Imager. These filters
23should be accessible from the filter interface as defined in
24the pod for Imager.
25
26=head1 FUNCTION REFERENCE
27
28Some of these functions are internal.
29
b8c2033e 30=over
02d1d628
AMH
31
32=cut
33*/
34
35
36
37
b2778574
AMH
38/*
39=item saturate(in)
40
41Clamps the input value between 0 and 255. (internal)
42
43 in - input integer
44
45=cut
46*/
47
48static
49unsigned char
50saturate(int in) {
51 if (in>255) { return 255; }
52 else if (in>0) return in;
53 return 0;
54}
02d1d628
AMH
55
56
57
58/*
59=item i_contrast(im, intensity)
60
61Scales the pixel values by the amount specified.
62
63 im - image object
64 intensity - scalefactor
65
66=cut
67*/
68
69void
70i_contrast(i_img *im, float intensity) {
71 int x, y;
72 unsigned char ch;
73 unsigned int new_color;
74 i_color rcolor;
75
76 mm_log((1,"i_contrast(im %p, intensity %f)\n", im, intensity));
77
78 if(intensity < 0) return;
79
80 for(y = 0; y < im->ysize; y++) for(x = 0; x < im->xsize; x++) {
81 i_gpix(im, x, y, &rcolor);
82
83 for(ch = 0; ch < im->channels; ch++) {
84 new_color = (unsigned int) rcolor.channel[ch];
85 new_color *= intensity;
86
87 if(new_color > 255) {
88 new_color = 255;
89 }
90 rcolor.channel[ch] = (unsigned char) new_color;
91 }
92 i_ppix(im, x, y, &rcolor);
93 }
94}
95
96
97/*
98=item i_hardinvert(im)
99
100Inverts the pixel values of the input image.
101
102 im - image object
103
104=cut
105*/
106
107void
108i_hardinvert(i_img *im) {
109 int x, y;
110 unsigned char ch;
111
112 i_color rcolor;
113
114 mm_log((1,"i_hardinvert(im %p)\n", im));
115
116 for(y = 0; y < im->ysize; y++) {
117 for(x = 0; x < im->xsize; x++) {
118 i_gpix(im, x, y, &rcolor);
119
120 for(ch = 0; ch < im->channels; ch++) {
121 rcolor.channel[ch] = 255 - rcolor.channel[ch];
122 }
123
124 i_ppix(im, x, y, &rcolor);
125 }
126 }
127}
128
129
130
131/*
132=item i_noise(im, amount, type)
133
134Inverts the pixel values by the amount specified.
135
136 im - image object
137 amount - deviation in pixel values
138 type - noise individual for each channel if true
139
140=cut
141*/
142
143#ifdef _MSC_VER
144/* random() is non-ASCII, even if it is better than rand() */
145#define random() rand()
146#endif
147
148void
149i_noise(i_img *im, float amount, unsigned char type) {
150 int x, y;
151 unsigned char ch;
152 int new_color;
153 float damount = amount * 2;
154 i_color rcolor;
155 int color_inc = 0;
156
157 mm_log((1,"i_noise(im %p, intensity %.2f\n", im, amount));
158
159 if(amount < 0) return;
160
161 for(y = 0; y < im->ysize; y++) for(x = 0; x < im->xsize; x++) {
162 i_gpix(im, x, y, &rcolor);
163
164 if(type == 0) {
165 color_inc = (amount - (damount * ((float)random() / RAND_MAX)));
166 }
167
168 for(ch = 0; ch < im->channels; ch++) {
169 new_color = (int) rcolor.channel[ch];
170
171 if(type != 0) {
172 new_color += (amount - (damount * ((float)random() / RAND_MAX)));
173 } else {
174 new_color += color_inc;
175 }
176
177 if(new_color < 0) {
178 new_color = 0;
179 }
180 if(new_color > 255) {
181 new_color = 255;
182 }
183
184 rcolor.channel[ch] = (unsigned char) new_color;
185 }
186
187 i_ppix(im, x, y, &rcolor);
188 }
189}
190
191
192/*
193=item i_noise(im, amount, type)
194
195Inverts the pixel values by the amount specified.
196
197 im - image object
198 amount - deviation in pixel values
199 type - noise individual for each channel if true
200
201=cut
202*/
203
204
205/*
206=item i_applyimage(im, add_im, mode)
207
208Apply's an image to another image
209
210 im - target image
211 add_im - image that is applied to target
212 mode - what method is used in applying:
213
214 0 Normal
215 1 Multiply
216 2 Screen
217 3 Overlay
218 4 Soft Light
219 5 Hard Light
220 6 Color dodge
221 7 Color Burn
222 8 Darker
223 9 Lighter
224 10 Add
225 11 Subtract
226 12 Difference
227 13 Exclusion
228
229=cut
230*/
231
232void i_applyimage(i_img *im, i_img *add_im, unsigned char mode) {
233 int x, y;
234 int mx, my;
235
236 mm_log((1, "i_applyimage(im %p, add_im %p, mode %d", im, add_im, mode));
237
238 mx = (add_im->xsize <= im->xsize) ? add_im->xsize : add_im->xsize;
239 my = (add_im->ysize <= im->ysize) ? add_im->ysize : add_im->ysize;
240
241 for(x = 0; x < mx; x++) {
242 for(y = 0; y < my; y++) {
243 }
244 }
245}
246
247
248/*
249=item i_bumpmap(im, bump, channel, light_x, light_y, st)
250
251Makes a bumpmap on image im using the bump image as the elevation map.
252
253 im - target image
254 bump - image that contains the elevation info
255 channel - to take the elevation information from
256 light_x - x coordinate of light source
257 light_y - y coordinate of light source
258 st - length of shadow
259
260=cut
261*/
262
263void
264i_bumpmap(i_img *im, i_img *bump, int channel, int light_x, int light_y, int st) {
265 int x, y, ch;
266 int mx, my;
267 i_color x1_color, y1_color, x2_color, y2_color, dst_color;
268 double nX, nY;
269 double tX, tY, tZ;
270 double aX, aY, aL;
271 double fZ;
272 unsigned char px1, px2, py1, py2;
273
274 i_img new_im;
275
276 mm_log((1, "i_bumpmap(im %p, add_im %p, channel %d, light_x %d, light_y %d, st %d)\n",
277 im, bump, channel, light_x, light_y, st));
278
279
280 if(channel >= bump->channels) {
281 mm_log((1, "i_bumpmap: channel = %d while bump image only has %d channels\n", channel, bump->channels));
282 return;
283 }
b2778574 284
02d1d628
AMH
285 mx = (bump->xsize <= im->xsize) ? bump->xsize : im->xsize;
286 my = (bump->ysize <= im->ysize) ? bump->ysize : im->ysize;
287
288 i_img_empty_ch(&new_im, im->xsize, im->ysize, im->channels);
b2778574 289
02d1d628
AMH
290 aX = (light_x > (mx >> 1)) ? light_x : mx - light_x;
291 aY = (light_y > (my >> 1)) ? light_y : my - light_y;
292
293 aL = sqrt((aX * aX) + (aY * aY));
294
295 for(y = 1; y < my - 1; y++) {
296 for(x = 1; x < mx - 1; x++) {
297 i_gpix(bump, x + st, y, &x1_color);
298 i_gpix(bump, x, y + st, &y1_color);
299 i_gpix(bump, x - st, y, &x2_color);
300 i_gpix(bump, x, y - st, &y2_color);
301
302 i_gpix(im, x, y, &dst_color);
303
304 px1 = x1_color.channel[channel];
305 py1 = y1_color.channel[channel];
306 px2 = x2_color.channel[channel];
307 py2 = y2_color.channel[channel];
308
309 nX = px1 - px2;
310 nY = py1 - py2;
311
312 nX += 128;
313 nY += 128;
314
315 fZ = (sqrt((nX * nX) + (nY * nY)) / aL);
316
317 tX = abs(x - light_x) / aL;
318 tY = abs(y - light_y) / aL;
319
320 tZ = 1 - (sqrt((tX * tX) + (tY * tY)) * fZ);
321
322 if(tZ < 0) tZ = 0;
323 if(tZ > 2) tZ = 2;
324
325 for(ch = 0; ch < im->channels; ch++)
326 dst_color.channel[ch] = (unsigned char) (float)(dst_color.channel[ch] * tZ);
327
328 i_ppix(&new_im, x, y, &dst_color);
329 }
330 }
331
332 i_copyto(im, &new_im, 0, 0, (int)im->xsize, (int)im->ysize, 0, 0);
333
334 i_img_exorcise(&new_im);
335}
336
337
338
b2778574
AMH
339
340typedef struct {
341 float x,y,z;
342} fvec;
343
344
345static
346float
347dotp(fvec *a, fvec *b) {
348 return a->x*b->x+a->y*b->y+a->z*b->z;
349}
350
351static
352void
353normalize(fvec *a) {
354 double d = sqrt(dotp(a,a));
355 a->x /= d;
356 a->y /= d;
357 a->z /= d;
358}
359
360
361/*
362 positive directions:
363
364 x - right,
365 y - down
366 z - out of the plane
367
368 I = Ia + Ip*( cd*Scol(N.L) + cs*(R.V)^n )
369
370 Here, the variables are:
371
372 * Ia - ambient colour
373 * Ip - intensity of the point light source
374 * cd - diffuse coefficient
375 * Scol - surface colour
376 * cs - specular coefficient
377 * n - objects shinyness
378 * N - normal vector
379 * L - lighting vector
380 * R - reflection vector
381 * V - vision vector
382
383 static void fvec_dump(fvec *x) {
384 printf("(%.2f %.2f %.2f)", x->x, x->y, x->z);
385 }
386*/
387
388/* XXX: Should these return a code for success? */
389
390
391
392
393/*
394=item i_bumpmap_complex(im, bump, channel, tx, ty, Lx, Ly, Lz, Ip, cd, cs, n, Ia, Il, Is)
395
396Makes a bumpmap on image im using the bump image as the elevation map.
397
398 im - target image
399 bump - image that contains the elevation info
400 channel - to take the elevation information from
401 tx - shift in x direction of where to start applying bumpmap
402 ty - shift in y direction of where to start applying bumpmap
403 Lx - x position/direction of light
404 Ly - y position/direction of light
405 Lz - z position/direction of light
406 Ip - light intensity
407 cd - diffuse coefficient
408 cs - specular coefficient
409 n - surface shinyness
410 Ia - ambient colour
411 Il - light colour
412 Is - specular colour
413
414if z<0 then the L is taken to be the direction the light is shining in. Otherwise
415the L is taken to be the position of the Light, Relative to the image.
416
417=cut
418*/
419
420
421void
422i_bumpmap_complex(i_img *im,
423 i_img *bump,
424 int channel,
425 int tx,
426 int ty,
427 float Lx,
428 float Ly,
429 float Lz,
430 float cd,
431 float cs,
432 float n,
433 i_color *Ia,
434 i_color *Il,
435 i_color *Is) {
436 i_img new_im;
437
438 int inflight;
439 int x, y, ch;
440 int mx, Mx, my, My;
441
442 float cdc[MAXCHANNELS];
443 float csc[MAXCHANNELS];
444
445 i_color x1_color, y1_color, x2_color, y2_color;
446
447 i_color Scol; /* Surface colour */
448
449 fvec L; /* Light vector */
450 fvec N; /* surface normal */
451 fvec R; /* Reflection vector */
452 fvec V; /* Vision vector */
453
454 mm_log((1, "i_bumpmap_complex(im %p, bump %p, channel %d, tx %d, ty %d, Lx %.2f, Ly %.2f, Lz %.2f, cd %.2f, cs %.2f, n %.2f, Ia %p, Il %p, Is %p)\n",
455 im, bump, channel, tx, ty, Lx, Ly, Lz, cd, cs, n, Ia, Il, Is));
456
457 if (channel >= bump->channels) {
458 mm_log((1, "i_bumpmap_complex: channel = %d while bump image only has %d channels\n", channel, bump->channels));
459 return;
460 }
461
462 for(ch=0; ch<im->channels; ch++) {
463 cdc[ch] = (float)Il->channel[ch]*cd/255.f;
464 csc[ch] = (float)Is->channel[ch]*cs/255.f;
465 }
466
467 mx = 1;
468 my = 1;
469 Mx = bump->xsize-1;
470 My = bump->ysize-1;
471
472 V.x = 0;
473 V.y = 0;
474 V.z = 1;
475
476 if (Lz < 0) { /* Light specifies a direction vector, reverse it to get the vector from surface to light */
477 L.x = -Lx;
478 L.y = -Ly;
479 L.z = -Lz;
480 normalize(&L);
481 } else { /* Light is the position of the light source */
482 inflight = 0;
483 L.x = -0.2;
484 L.y = -0.4;
485 L.z = 1;
486 normalize(&L);
487 }
488
489 i_img_empty_ch(&new_im, im->xsize, im->ysize, im->channels);
490
491 for(y = 0; y < im->ysize; y++) {
492 for(x = 0; x < im->xsize; x++) {
493 double dp1, dp2;
494 double dx = 0, dy = 0;
495
496 /* Calculate surface normal */
497 if (mx<x && x<Mx && my<y && y<My) {
498 i_gpix(bump, x + 1, y, &x1_color);
499 i_gpix(bump, x - 1, y, &x2_color);
500 i_gpix(bump, x, y + 1, &y1_color);
501 i_gpix(bump, x, y - 1, &y2_color);
502 dx = x2_color.channel[channel] - x1_color.channel[channel];
503 dy = y2_color.channel[channel] - y1_color.channel[channel];
504 } else {
505 dx = 0;
506 dy = 0;
507 }
508 N.x = -dx * 0.015;
509 N.y = -dy * 0.015;
510 N.z = 1;
511 normalize(&N);
512
513 /* Calculate Light vector if needed */
514 if (Lz>=0) {
515 L.x = Lx - x;
516 L.y = Ly - y;
517 L.z = Lz;
518 normalize(&L);
519 }
520
521 dp1 = dotp(&L,&N);
522 R.x = -L.x + 2*dp1*N.x;
523 R.y = -L.y + 2*dp1*N.y;
524 R.z = -L.z + 2*dp1*N.z;
525
526 dp2 = dotp(&R,&V);
527
528 dp1 = dp1<0 ?0 : dp1;
529 dp2 = pow(dp2<0 ?0 : dp2,n);
530
531 i_gpix(im, x, y, &Scol);
532
533 for(ch = 0; ch < im->channels; ch++)
534 Scol.channel[ch] =
535 saturate( Ia->channel[ch] + cdc[ch]*Scol.channel[ch]*dp1 + csc[ch]*dp2 );
536
537 i_ppix(&new_im, x, y, &Scol);
538 }
539 }
540
541 i_copyto(im, &new_im, 0, 0, (int)im->xsize, (int)im->ysize, 0, 0);
542 i_img_exorcise(&new_im);
543}
544
545
02d1d628
AMH
546/*
547=item i_postlevels(im, levels)
548
549Quantizes Images to fewer levels.
550
551 im - target image
552 levels - number of levels
553
554=cut
555*/
556
557void
558i_postlevels(i_img *im, int levels) {
559 int x, y, ch;
560 float pv;
561 int rv;
562 float av;
563
564 i_color rcolor;
565
566 rv = (int) ((float)(256 / levels));
567 av = (float)levels;
568
569 for(y = 0; y < im->ysize; y++) for(x = 0; x < im->xsize; x++) {
570 i_gpix(im, x, y, &rcolor);
571
572 for(ch = 0; ch < im->channels; ch++) {
573 pv = (((float)rcolor.channel[ch] / 255)) * av;
574 pv = (int) ((int)pv * rv);
575
576 if(pv < 0) pv = 0;
577 else if(pv > 255) pv = 255;
578
579 rcolor.channel[ch] = (unsigned char) pv;
580 }
581 i_ppix(im, x, y, &rcolor);
582 }
583}
584
585
586/*
587=item i_mosaic(im, size)
588
589Makes an image looks like a mosaic with tilesize of size
590
591 im - target image
592 size - size of tiles
593
594=cut
595*/
596
597void
598i_mosaic(i_img *im, int size) {
599 int x, y, ch;
600 int lx, ly, z;
601 long sqrsize;
602
603 i_color rcolor;
604 long col[256];
605
606 sqrsize = size * size;
607
608 for(y = 0; y < im->ysize; y += size) for(x = 0; x < im->xsize; x += size) {
609 for(z = 0; z < 256; z++) col[z] = 0;
610
611 for(lx = 0; lx < size; lx++) {
612 for(ly = 0; ly < size; ly++) {
613 i_gpix(im, (x + lx), (y + ly), &rcolor);
614
615 for(ch = 0; ch < im->channels; ch++) {
616 col[ch] += rcolor.channel[ch];
617 }
618 }
619 }
620
621 for(ch = 0; ch < im->channels; ch++)
622 rcolor.channel[ch] = (int) ((float)col[ch] / sqrsize);
623
624
625 for(lx = 0; lx < size; lx++)
626 for(ly = 0; ly < size; ly++)
627 i_ppix(im, (x + lx), (y + ly), &rcolor);
628
629 }
630}
631
02d1d628
AMH
632
633/*
634=item i_watermark(im, wmark, tx, ty, pixdiff)
635
636Applies a watermark to the target image
637
638 im - target image
639 wmark - watermark image
640 tx - x coordinate of where watermark should be applied
641 ty - y coordinate of where watermark should be applied
642 pixdiff - the magnitude of the watermark, controls how visible it is
643
644=cut
645*/
646
647void
648i_watermark(i_img *im, i_img *wmark, int tx, int ty, int pixdiff) {
649 int vx, vy, ch;
650 i_color val, wval;
651
652 for(vx=0;vx<128;vx++) for(vy=0;vy<110;vy++) {
653
654 i_gpix(im, tx+vx, ty+vy,&val );
655 i_gpix(wmark, vx, vy, &wval);
656
657 for(ch=0;ch<im->channels;ch++)
658 val.channel[ch] = saturate( val.channel[ch] + (pixdiff* (wval.channel[0]-128) )/128 );
659
660 i_ppix(im,tx+vx,ty+vy,&val);
661 }
662}
663
664
665/*
666=item i_autolevels(im, lsat, usat, skew)
667
668Scales and translates each color such that it fills the range completely.
669Skew is not implemented yet - purpose is to control the color skew that can
670occur when changing the contrast.
671
672 im - target image
673 lsat - fraction of pixels that will be truncated at the lower end of the spectrum
674 usat - fraction of pixels that will be truncated at the higher end of the spectrum
675 skew - not used yet
676
677=cut
678*/
679
680void
681i_autolevels(i_img *im, float lsat, float usat, float skew) {
682 i_color val;
683 int i, x, y, rhist[256], ghist[256], bhist[256];
684 int rsum, rmin, rmax;
685 int gsum, gmin, gmax;
686 int bsum, bmin, bmax;
687 int rcl, rcu, gcl, gcu, bcl, bcu;
688
689 mm_log((1,"i_autolevels(im %p, lsat %f,usat %f,skew %f)\n", im, lsat,usat,skew));
690
691 rsum=gsum=bsum=0;
692 for(i=0;i<256;i++) rhist[i]=ghist[i]=bhist[i] = 0;
693 /* create histogram for each channel */
694 for(y = 0; y < im->ysize; y++) for(x = 0; x < im->xsize; x++) {
695 i_gpix(im, x, y, &val);
696 rhist[val.channel[0]]++;
697 ghist[val.channel[1]]++;
698 bhist[val.channel[2]]++;
699 }
700
701 for(i=0;i<256;i++) {
702 rsum+=rhist[i];
703 gsum+=ghist[i];
704 bsum+=bhist[i];
705 }
706
707 rmin = gmin = bmin = 0;
708 rmax = gmax = bmax = 255;
709
710 rcu = rcl = gcu = gcl = bcu = bcl = 0;
711
712 for(i=0; i<256; i++) {
713 rcl += rhist[i]; if ( (rcl<rsum*lsat) ) rmin=i;
714 rcu += rhist[255-i]; if ( (rcu<rsum*usat) ) rmax=255-i;
715
716 gcl += ghist[i]; if ( (gcl<gsum*lsat) ) gmin=i;
717 gcu += ghist[255-i]; if ( (gcu<gsum*usat) ) gmax=255-i;
718
719 bcl += bhist[i]; if ( (bcl<bsum*lsat) ) bmin=i;
720 bcu += bhist[255-i]; if ( (bcu<bsum*usat) ) bmax=255-i;
721 }
722
723 for(y = 0; y < im->ysize; y++) for(x = 0; x < im->xsize; x++) {
724 i_gpix(im, x, y, &val);
725 val.channel[0]=saturate((val.channel[0]-rmin)*255/(rmax-rmin));
726 val.channel[1]=saturate((val.channel[1]-gmin)*255/(gmax-gmin));
727 val.channel[2]=saturate((val.channel[2]-bmin)*255/(bmax-bmin));
728 i_ppix(im, x, y, &val);
729 }
730}
731
732/*
733=item Noise(x,y)
734
735Pseudo noise utility function used to generate perlin noise. (internal)
736
737 x - x coordinate
738 y - y coordinate
739
740=cut
741*/
742
743static
744float
745Noise(int x, int y) {
746 int n = x + y * 57;
747 n = (n<<13) ^ n;
748 return ( 1.0 - ( (n * (n * n * 15731 + 789221) + 1376312589) & 0x7fffffff) / 1073741824.0);
749}
750
751/*
752=item SmoothedNoise1(x,y)
753
754Pseudo noise utility function used to generate perlin noise. (internal)
755
756 x - x coordinate
757 y - y coordinate
758
759=cut
760*/
761
762static
763float
764SmoothedNoise1(float x, float y) {
765 float corners = ( Noise(x-1, y-1)+Noise(x+1, y-1)+Noise(x-1, y+1)+Noise(x+1, y+1) ) / 16;
766 float sides = ( Noise(x-1, y) +Noise(x+1, y) +Noise(x, y-1) +Noise(x, y+1) ) / 8;
767 float center = Noise(x, y) / 4;
768 return corners + sides + center;
769}
770
771
772/*
773=item G_Interpolate(a, b, x)
774
775Utility function used to generate perlin noise. (internal)
776
777=cut
778*/
779
780static
781float C_Interpolate(float a, float b, float x) {
782 /* float ft = x * 3.1415927; */
783 float ft = x * PI;
784 float f = (1 - cos(ft)) * .5;
785 return a*(1-f) + b*f;
786}
787
788
789/*
790=item InterpolatedNoise(x, y)
791
792Utility function used to generate perlin noise. (internal)
793
794=cut
795*/
796
797static
798float
799InterpolatedNoise(float x, float y) {
800
801 int integer_X = x;
802 float fractional_X = x - integer_X;
803 int integer_Y = y;
804 float fractional_Y = y - integer_Y;
805
806 float v1 = SmoothedNoise1(integer_X, integer_Y);
807 float v2 = SmoothedNoise1(integer_X + 1, integer_Y);
808 float v3 = SmoothedNoise1(integer_X, integer_Y + 1);
809 float v4 = SmoothedNoise1(integer_X + 1, integer_Y + 1);
810
811 float i1 = C_Interpolate(v1 , v2 , fractional_X);
812 float i2 = C_Interpolate(v3 , v4 , fractional_X);
813
814 return C_Interpolate(i1 , i2 , fractional_Y);
815}
816
817
818
819/*
820=item PerlinNoise_2D(x, y)
821
822Utility function used to generate perlin noise. (internal)
823
824=cut
825*/
826
827static
828float
829PerlinNoise_2D(float x, float y) {
830 int i,frequency;
831 float amplitude;
832 float total = 0;
833 int Number_Of_Octaves=6;
834 int n = Number_Of_Octaves - 1;
835
836 for(i=0;i<n;i++) {
837 frequency = 2*i;
838 amplitude = PI;
839 total = total + InterpolatedNoise(x * frequency, y * frequency) * amplitude;
840 }
841
842 return total;
843}
844
845
846/*
847=item i_radnoise(im, xo, yo, rscale, ascale)
848
849Perlin-like radial noise.
850
851 im - target image
852 xo - x coordinate of center
853 yo - y coordinate of center
854 rscale - radial scale
855 ascale - angular scale
856
857=cut
858*/
859
860void
861i_radnoise(i_img *im, int xo, int yo, float rscale, float ascale) {
862 int x, y, ch;
863 i_color val;
864 unsigned char v;
865 float xc, yc, r;
866 double a;
867
868 for(y = 0; y < im->ysize; y++) for(x = 0; x < im->xsize; x++) {
869 xc = (float)x-xo+0.5;
870 yc = (float)y-yo+0.5;
871 r = rscale*sqrt(xc*xc+yc*yc)+1.2;
872 a = (PI+atan2(yc,xc))*ascale;
873 v = saturate(128+100*(PerlinNoise_2D(a,r)));
874 /* v=saturate(120+12*PerlinNoise_2D(xo+(float)x/scale,yo+(float)y/scale)); Good soft marble */
875 for(ch=0; ch<im->channels; ch++) val.channel[ch]=v;
876 i_ppix(im, x, y, &val);
877 }
878}
879
880
881/*
882=item i_turbnoise(im, xo, yo, scale)
883
884Perlin-like 2d noise noise.
885
886 im - target image
887 xo - x coordinate translation
888 yo - y coordinate translation
889 scale - scale of noise
890
891=cut
892*/
893
894void
895i_turbnoise(i_img *im, float xo, float yo, float scale) {
896 int x,y,ch;
897 unsigned char v;
898 i_color val;
899
900 for(y = 0; y < im->ysize; y++) for(x = 0; x < im->xsize; x++) {
901 /* v=saturate(125*(1.0+PerlinNoise_2D(xo+(float)x/scale,yo+(float)y/scale))); */
902 v = saturate(120*(1.0+sin(xo+(float)x/scale+PerlinNoise_2D(xo+(float)x/scale,yo+(float)y/scale))));
903 for(ch=0; ch<im->channels; ch++) val.channel[ch] = v;
904 i_ppix(im, x, y, &val);
905 }
906}
907
908
909
910/*
911=item i_gradgen(im, num, xo, yo, ival, dmeasure)
912
913Gradient generating function.
914
915 im - target image
916 num - number of points given
917 xo - array of x coordinates
918 yo - array of y coordinates
919 ival - array of i_color objects
920 dmeasure - distance measure to be used.
921 0 = Euclidean
922 1 = Euclidean squared
923 2 = Manhattan distance
924
925=cut
926*/
927
928
929void
930i_gradgen(i_img *im, int num, int *xo, int *yo, i_color *ival, int dmeasure) {
931
932 i_color val;
933 int p, x, y, ch;
934 int channels = im->channels;
935 int xsize = im->xsize;
936 int ysize = im->ysize;
937
938 float *fdist;
939
940 mm_log((1,"i_gradgen(im %p, num %d, xo %p, yo %p, ival %p, dmeasure %d)\n", im, num, xo, yo, ival, dmeasure));
941
942 for(p = 0; p<num; p++) {
943 mm_log((1,"i_gradgen: (%d, %d)\n", xo[p], yo[p]));
944 ICL_info(&ival[p]);
945 }
946
947 fdist = mymalloc( sizeof(float) * num );
948
949 for(y = 0; y<ysize; y++) for(x = 0; x<xsize; x++) {
950 float cs = 0;
951 float csd = 0;
952 for(p = 0; p<num; p++) {
953 int xd = x-xo[p];
954 int yd = y-yo[p];
955 switch (dmeasure) {
956 case 0: /* euclidean */
957 fdist[p] = sqrt(xd*xd + yd*yd); /* euclidean distance */
958 break;
959 case 1: /* euclidean squared */
960 fdist[p] = xd*xd + yd*yd; /* euclidean distance */
961 break;
962 case 2: /* euclidean squared */
b33c08f8 963 fdist[p] = i_max(xd*xd, yd*yd); /* manhattan distance */
02d1d628
AMH
964 break;
965 default:
966 m_fatal(3,"i_gradgen: Unknown distance measure\n");
967 }
968 cs += fdist[p];
969 }
970
971 csd = 1/((num-1)*cs);
972
973 for(p = 0; p<num; p++) fdist[p] = (cs-fdist[p])*csd;
974
975 for(ch = 0; ch<channels; ch++) {
976 int tres = 0;
977 for(p = 0; p<num; p++) tres += ival[p].channel[ch] * fdist[p];
978 val.channel[ch] = saturate(tres);
979 }
980 i_ppix(im, x, y, &val);
981 }
a73aeb5f 982 myfree(fdist);
02d1d628
AMH
983
984}
985
02d1d628
AMH
986void
987i_nearest_color_foo(i_img *im, int num, int *xo, int *yo, i_color *ival, int dmeasure) {
988
a743c0a6 989 int p, x, y;
02d1d628
AMH
990 int xsize = im->xsize;
991 int ysize = im->ysize;
992
993 mm_log((1,"i_gradgen(im %p, num %d, xo %p, yo %p, ival %p, dmeasure %d)\n", im, num, xo, yo, ival, dmeasure));
994
995 for(p = 0; p<num; p++) {
996 mm_log((1,"i_gradgen: (%d, %d)\n", xo[p], yo[p]));
997 ICL_info(&ival[p]);
998 }
999
1000 for(y = 0; y<ysize; y++) for(x = 0; x<xsize; x++) {
1001 int midx = 0;
1002 float mindist = 0;
1003 float curdist = 0;
1004
1005 int xd = x-xo[0];
1006 int yd = y-yo[0];
1007
1008 switch (dmeasure) {
1009 case 0: /* euclidean */
1010 mindist = sqrt(xd*xd + yd*yd); /* euclidean distance */
1011 break;
1012 case 1: /* euclidean squared */
1013 mindist = xd*xd + yd*yd; /* euclidean distance */
1014 break;
1015 case 2: /* euclidean squared */
b33c08f8 1016 mindist = i_max(xd*xd, yd*yd); /* manhattan distance */
02d1d628
AMH
1017 break;
1018 default:
1019 m_fatal(3,"i_nearest_color: Unknown distance measure\n");
1020 }
1021
1022 for(p = 1; p<num; p++) {
1023 xd = x-xo[p];
1024 yd = y-yo[p];
1025 switch (dmeasure) {
1026 case 0: /* euclidean */
1027 curdist = sqrt(xd*xd + yd*yd); /* euclidean distance */
1028 break;
1029 case 1: /* euclidean squared */
1030 curdist = xd*xd + yd*yd; /* euclidean distance */
1031 break;
1032 case 2: /* euclidean squared */
b33c08f8 1033 curdist = i_max(xd*xd, yd*yd); /* manhattan distance */
02d1d628
AMH
1034 break;
1035 default:
1036 m_fatal(3,"i_nearest_color: Unknown distance measure\n");
1037 }
1038 if (curdist < mindist) {
1039 mindist = curdist;
1040 midx = p;
1041 }
1042 }
1043 i_ppix(im, x, y, &ival[midx]);
1044 }
1045}
1046
02d1d628
AMH
1047void
1048i_nearest_color(i_img *im, int num, int *xo, int *yo, i_color *oval, int dmeasure) {
1049 i_color *ival;
1050 float *tval;
1051 float c1, c2;
1052 i_color val;
1053 int p, x, y, ch;
02d1d628
AMH
1054 int xsize = im->xsize;
1055 int ysize = im->ysize;
1056 int *cmatch;
1057
d7f707d8 1058 mm_log((1,"i_nearest_color(im %p, num %d, xo %p, yo %p, ival %p, dmeasure %d)\n", im, num, xo, yo, oval, dmeasure));
02d1d628
AMH
1059
1060 tval = mymalloc( sizeof(float)*num*im->channels );
1061 ival = mymalloc( sizeof(i_color)*num );
1062 cmatch = mymalloc( sizeof(int)*num );
1063
1064 for(p = 0; p<num; p++) {
1065 for(ch = 0; ch<im->channels; ch++) tval[ p * im->channels + ch] = 0;
1066 cmatch[p] = 0;
1067 }
1068
1069
1070 for(y = 0; y<ysize; y++) for(x = 0; x<xsize; x++) {
1071 int midx = 0;
1072 float mindist = 0;
1073 float curdist = 0;
1074
1075 int xd = x-xo[0];
1076 int yd = y-yo[0];
1077
1078 switch (dmeasure) {
1079 case 0: /* euclidean */
1080 mindist = sqrt(xd*xd + yd*yd); /* euclidean distance */
1081 break;
1082 case 1: /* euclidean squared */
1083 mindist = xd*xd + yd*yd; /* euclidean distance */
1084 break;
1085 case 2: /* euclidean squared */
b33c08f8 1086 mindist = i_max(xd*xd, yd*yd); /* manhattan distance */
02d1d628
AMH
1087 break;
1088 default:
1089 m_fatal(3,"i_nearest_color: Unknown distance measure\n");
1090 }
1091
1092 for(p = 1; p<num; p++) {
1093 xd = x-xo[p];
1094 yd = y-yo[p];
1095 switch (dmeasure) {
1096 case 0: /* euclidean */
1097 curdist = sqrt(xd*xd + yd*yd); /* euclidean distance */
1098 break;
1099 case 1: /* euclidean squared */
1100 curdist = xd*xd + yd*yd; /* euclidean distance */
1101 break;
1102 case 2: /* euclidean squared */
b33c08f8 1103 curdist = i_max(xd*xd, yd*yd); /* manhattan distance */
02d1d628
AMH
1104 break;
1105 default:
1106 m_fatal(3,"i_nearest_color: Unknown distance measure\n");
1107 }
1108 if (curdist < mindist) {
1109 mindist = curdist;
1110 midx = p;
1111 }
1112 }
1113
1114 cmatch[midx]++;
1115 i_gpix(im, x, y, &val);
1116 c2 = 1.0/(float)(cmatch[midx]);
1117 c1 = 1.0-c2;
1118
02d1d628
AMH
1119 for(ch = 0; ch<im->channels; ch++)
1120 tval[midx*im->channels + ch] = c1*tval[midx*im->channels + ch] + c2 * (float) val.channel[ch];
3bb1c1f3 1121
02d1d628
AMH
1122
1123 }
1124
1125 for(p = 0; p<num; p++) for(ch = 0; ch<im->channels; ch++) ival[p].channel[ch] = tval[p*im->channels + ch];
1126
1127 i_nearest_color_foo(im, num, xo, yo, ival, dmeasure);
1128}
6607600c 1129
b6381851
TC
1130/*
1131=item i_unsharp_mask(im, stddev, scale)
1132
1133Perform an usharp mask, which is defined as subtracting the blurred
1134image from double the original.
1135
1136=cut
1137*/
1138void i_unsharp_mask(i_img *im, double stddev, double scale) {
1139 i_img copy;
1140 int x, y, ch;
1141
1142 if (scale < 0)
1143 return;
1144 /* it really shouldn't ever be more than 1.0, but maybe ... */
1145 if (scale > 100)
1146 scale = 100;
1147
1148 i_copy(&copy, im);
1149 i_gaussian(&copy, stddev);
1150 if (im->bits == i_8_bits) {
1151 i_color *blur = mymalloc(im->xsize * sizeof(i_color) * 2);
1152 i_color *out = blur + im->xsize;
1153
1154 for (y = 0; y < im->ysize; ++y) {
1155 i_glin(&copy, 0, copy.xsize, y, blur);
1156 i_glin(im, 0, im->xsize, y, out);
1157 for (x = 0; x < im->xsize; ++x) {
1158 for (ch = 0; ch < im->channels; ++ch) {
1159 /*int temp = out[x].channel[ch] +
1160 scale * (out[x].channel[ch] - blur[x].channel[ch]);*/
1161 int temp = out[x].channel[ch] * 2 - blur[x].channel[ch];
1162 if (temp < 0)
1163 temp = 0;
1164 else if (temp > 255)
1165 temp = 255;
1166 out[x].channel[ch] = temp;
1167 }
1168 }
1169 i_plin(im, 0, im->xsize, y, out);
1170 }
1171
1172 myfree(blur);
1173 }
1174 else {
1175 i_fcolor *blur = mymalloc(im->xsize * sizeof(i_fcolor) * 2);
1176 i_fcolor *out = blur + im->xsize;
1177
1178 for (y = 0; y < im->ysize; ++y) {
1179 i_glinf(&copy, 0, copy.xsize, y, blur);
1180 i_glinf(im, 0, im->xsize, y, out);
1181 for (x = 0; x < im->xsize; ++x) {
1182 for (ch = 0; ch < im->channels; ++ch) {
1183 double temp = out[x].channel[ch] +
1184 scale * (out[x].channel[ch] - blur[x].channel[ch]);
1185 if (temp < 0)
1186 temp = 0;
1187 else if (temp > 1.0)
1188 temp = 1.0;
1189 out[x].channel[ch] = temp;
1190 }
1191 }
1192 i_plinf(im, 0, im->xsize, y, out);
1193 }
1194
1195 myfree(blur);
1196 }
1197 i_img_exorcise(&copy);
1198}
1199
dff75dee
TC
1200/*
1201=item i_diff_image(im1, im2, mindiff)
1202
1203Creates a new image that is transparent, except where the pixel in im2
1204is different from im1, where it is the pixel from im2.
1205
1206The samples must differ by at least mindiff to be considered different.
1207
1208=cut
1209*/
1210
1211i_img *
1212i_diff_image(i_img *im1, i_img *im2, int mindiff) {
1213 i_img *out;
1214 int outchans, diffchans;
1215 int xsize, ysize;
1216 i_img temp;
1217
1218 i_clear_error();
1219 if (im1->channels != im2->channels) {
1220 i_push_error(0, "different number of channels");
1221 return NULL;
1222 }
1223
1224 outchans = diffchans = im1->channels;
1225 if (outchans == 1 || outchans == 3)
1226 ++outchans;
1227
b33c08f8
TC
1228 xsize = i_min(im1->xsize, im2->xsize);
1229 ysize = i_min(im1->ysize, im2->ysize);
dff75dee
TC
1230
1231 out = i_sametype_chans(im1, xsize, ysize, outchans);
1232
1233 if (im1->bits == i_8_bits && im2->bits == i_8_bits) {
1234 i_color *line1 = mymalloc(2 * xsize * sizeof(*line1));
1235 i_color *line2 = line1 + xsize;
1236 i_color empty;
1237 int x, y, ch;
1238
1239 for (ch = 0; ch < MAXCHANNELS; ++ch)
1240 empty.channel[ch] = 0;
1241
1242 for (y = 0; y < ysize; ++y) {
1243 i_glin(im1, 0, xsize, y, line1);
1244 i_glin(im2, 0, xsize, y, line2);
1245 for (x = 0; x < xsize; ++x) {
1246 int diff = 0;
1247 for (ch = 0; ch < diffchans; ++ch) {
1248 if (line1[x].channel[ch] != line2[x].channel[ch]
1249 && abs(line1[x].channel[ch] - line2[x].channel[ch]) > mindiff) {
1250 diff = 1;
1251 break;
1252 }
1253 }
1254 if (!diff)
1255 line2[x] = empty;
1256 }
1257 i_plin(out, 0, xsize, y, line2);
1258 }
1259 myfree(line1);
1260 }
1261 else {
1262 i_fcolor *line1 = mymalloc(2 * xsize * sizeof(*line1));
1263 i_fcolor *line2 = line1 + xsize;
1264 i_fcolor empty;
1265 int x, y, ch;
1266 double dist = mindiff / 255;
1267
1268 for (ch = 0; ch < MAXCHANNELS; ++ch)
1269 empty.channel[ch] = 0;
1270
1271 for (y = 0; y < ysize; ++y) {
1272 i_glinf(im1, 0, xsize, y, line1);
1273 i_glinf(im2, 0, xsize, y, line2);
1274 for (x = 0; x < xsize; ++x) {
1275 int diff = 0;
1276 for (ch = 0; ch < diffchans; ++ch) {
1277 if (line1[x].channel[ch] != line2[x].channel[ch]
1278 && abs(line1[x].channel[ch] - line2[x].channel[ch]) > dist) {
1279 diff = 1;
1280 break;
1281 }
1282 }
1283 if (!diff)
1284 line2[x] = empty;
1285 }
1286 i_plinf(out, 0, xsize, y, line2);
1287 }
1288 myfree(line1);
1289 }
1290
1291 return out;
1292}
1293
f1ac5027 1294struct fount_state;
6607600c
TC
1295static double linear_fount_f(double x, double y, struct fount_state *state);
1296static double bilinear_fount_f(double x, double y, struct fount_state *state);
1297static double radial_fount_f(double x, double y, struct fount_state *state);
1298static double square_fount_f(double x, double y, struct fount_state *state);
1299static double revolution_fount_f(double x, double y,
1300 struct fount_state *state);
1301static double conical_fount_f(double x, double y, struct fount_state *state);
1302
1303typedef double (*fount_func)(double, double, struct fount_state *);
1304static fount_func fount_funcs[] =
1305{
1306 linear_fount_f,
1307 bilinear_fount_f,
1308 radial_fount_f,
1309 square_fount_f,
1310 revolution_fount_f,
1311 conical_fount_f,
1312};
1313
1314static double linear_interp(double pos, i_fountain_seg *seg);
1315static double sine_interp(double pos, i_fountain_seg *seg);
1316static double sphereup_interp(double pos, i_fountain_seg *seg);
1317static double spheredown_interp(double pos, i_fountain_seg *seg);
1318typedef double (*fount_interp)(double pos, i_fountain_seg *seg);
1319static fount_interp fount_interps[] =
1320{
1321 linear_interp,
1322 linear_interp,
1323 sine_interp,
1324 sphereup_interp,
1325 spheredown_interp,
1326};
1327
1328static void direct_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg);
1329static void hue_up_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg);
1330static void hue_down_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg);
1331typedef void (*fount_cinterp)(i_fcolor *out, double pos, i_fountain_seg *seg);
1332static fount_cinterp fount_cinterps[] =
1333{
1334 direct_cinterp,
1335 hue_up_cinterp,
1336 hue_down_cinterp,
1337};
1338
1339typedef double (*fount_repeat)(double v);
1340static double fount_r_none(double v);
1341static double fount_r_sawtooth(double v);
1342static double fount_r_triangle(double v);
1343static double fount_r_saw_both(double v);
1344static double fount_r_tri_both(double v);
1345static fount_repeat fount_repeats[] =
1346{
1347 fount_r_none,
1348 fount_r_sawtooth,
1349 fount_r_triangle,
1350 fount_r_saw_both,
1351 fount_r_tri_both,
1352};
1353
f1ac5027
TC
1354static int simple_ssample(i_fcolor *out, double x, double y,
1355 struct fount_state *state);
1356static int random_ssample(i_fcolor *out, double x, double y,
1357 struct fount_state *state);
1358static int circle_ssample(i_fcolor *out, double x, double y,
1359 struct fount_state *state);
1360typedef int (*fount_ssample)(i_fcolor *out, double x, double y,
1361 struct fount_state *state);
6607600c
TC
1362static fount_ssample fount_ssamples[] =
1363{
1364 NULL,
1365 simple_ssample,
1366 random_ssample,
1367 circle_ssample,
1368};
1369
1370static int
f1ac5027
TC
1371fount_getat(i_fcolor *out, double x, double y, struct fount_state *state);
1372
1373/*
1374 Keep state information used by each type of fountain fill
1375*/
1376struct fount_state {
1377 /* precalculated for the equation of the line perpendicular to the line AB */
1378 double lA, lB, lC;
1379 double AB;
1380 double sqrtA2B2;
1381 double mult;
1382 double cos;
1383 double sin;
1384 double theta;
1385 int xa, ya;
1386 void *ssample_data;
1387 fount_func ffunc;
1388 fount_repeat rpfunc;
1389 fount_ssample ssfunc;
1390 double parm;
1391 i_fountain_seg *segs;
1392 int count;
1393};
1394
1395static void
1396fount_init_state(struct fount_state *state, double xa, double ya,
1397 double xb, double yb, i_fountain_type type,
1398 i_fountain_repeat repeat, int combine, int super_sample,
1399 double ssample_param, int count, i_fountain_seg *segs);
1400
1401static void
1402fount_finish_state(struct fount_state *state);
6607600c
TC
1403
1404#define EPSILON (1e-6)
1405
1406/*
1407=item i_fountain(im, xa, ya, xb, yb, type, repeat, combine, super_sample, ssample_param, count, segs)
1408
1409Draws a fountain fill using A(xa, ya) and B(xb, yb) as reference points.
1410
1411I<type> controls how the reference points are used:
1412
1413=over
1414
1415=item i_ft_linear
1416
1417linear, where A is 0 and B is 1.
1418
1419=item i_ft_bilinear
1420
1421linear in both directions from A.
1422
1423=item i_ft_radial
1424
1425circular, where A is the centre of the fill, and B is a point
1426on the radius.
1427
1428=item i_ft_radial_square
1429
1430where A is the centre of the fill and B is the centre of
1431one side of the square.
1432
1433=item i_ft_revolution
1434
1435where A is the centre of the fill and B defines the 0/1.0
1436angle of the fill.
1437
1438=item i_ft_conical
1439
1440similar to i_ft_revolution, except that the revolution goes in both
1441directions
1442
1443=back
1444
1445I<repeat> can be one of:
1446
1447=over
1448
1449=item i_fr_none
1450
1451values < 0 are treated as zero, values > 1 are treated as 1.
1452
1453=item i_fr_sawtooth
1454
1455negative values are treated as 0, positive values are modulo 1.0
1456
1457=item i_fr_triangle
1458
1459negative values are treated as zero, if (int)value is odd then the value is treated as 1-(value
1460mod 1.0), otherwise the same as for sawtooth.
1461
1462=item i_fr_saw_both
1463
1464like i_fr_sawtooth, except that the sawtooth pattern repeats into
1465negative values.
1466
1467=item i_fr_tri_both
1468
1469Like i_fr_triangle, except that negative values are handled as their
1470absolute values.
1471
1472=back
1473
1474If combine is non-zero then non-opaque values are combined with the
1475underlying color.
1476
1477I<super_sample> controls super sampling, if any. At some point I'll
1478probably add a adaptive super-sampler. Current possible values are:
1479
1480=over
1481
1482=item i_fts_none
1483
1484No super-sampling is done.
1485
1486=item i_fts_grid
1487
1488A square grid of points withing the pixel are sampled.
1489
1490=item i_fts_random
1491
1492Random points within the pixel are sampled.
1493
1494=item i_fts_circle
1495
1496Points on the radius of a circle are sampled. This produces fairly
1497good results, but is fairly slow since sin() and cos() are evaluated
1498for each point.
1499
1500=back
1501
1502I<ssample_param> is intended to be roughly the number of points
1503sampled within the pixel.
1504
1505I<count> and I<segs> define the segments of the fill.
1506
1507=cut
1508
1509*/
1510
1511void
1512i_fountain(i_img *im, double xa, double ya, double xb, double yb,
1513 i_fountain_type type, i_fountain_repeat repeat,
1514 int combine, int super_sample, double ssample_param,
1515 int count, i_fountain_seg *segs) {
1516 struct fount_state state;
6607600c
TC
1517 int x, y;
1518 i_fcolor *line = mymalloc(sizeof(i_fcolor) * im->xsize);
efdc2568 1519 i_fcolor *work = NULL;
290bdf77 1520
f1ac5027 1521 i_fountain_seg *my_segs;
efdc2568
TC
1522 i_fill_combine_f combine_func = NULL;
1523 i_fill_combinef_f combinef_func = NULL;
1524
1525 i_get_combine(combine, &combine_func, &combinef_func);
1526 if (combinef_func)
1527 work = mymalloc(sizeof(i_fcolor) * im->xsize);
f1ac5027
TC
1528
1529 fount_init_state(&state, xa, ya, xb, yb, type, repeat, combine,
1530 super_sample, ssample_param, count, segs);
1531 my_segs = state.segs;
1532
1533 for (y = 0; y < im->ysize; ++y) {
1534 i_glinf(im, 0, im->xsize, y, line);
1535 for (x = 0; x < im->xsize; ++x) {
1536 i_fcolor c;
1537 int got_one;
f1ac5027
TC
1538 if (super_sample == i_fts_none)
1539 got_one = fount_getat(&c, x, y, &state);
1540 else
1541 got_one = state.ssfunc(&c, x, y, &state);
1542 if (got_one) {
efdc2568
TC
1543 if (combine)
1544 work[x] = c;
f1ac5027
TC
1545 else
1546 line[x] = c;
1547 }
1548 }
efdc2568
TC
1549 if (combine)
1550 combinef_func(line, work, im->channels, im->xsize);
f1ac5027
TC
1551 i_plinf(im, 0, im->xsize, y, line);
1552 }
1553 fount_finish_state(&state);
a73aeb5f 1554 if (work) myfree(work);
f1ac5027
TC
1555 myfree(line);
1556}
1557
1558typedef struct {
1559 i_fill_t base;
1560 struct fount_state state;
1561} i_fill_fountain_t;
1562
1563static void
1564fill_fountf(i_fill_t *fill, int x, int y, int width, int channels,
43c5dacb 1565 i_fcolor *data);
f1ac5027
TC
1566static void
1567fount_fill_destroy(i_fill_t *fill);
1568
1569/*
1570=item i_new_fount(xa, ya, xb, yb, type, repeat, combine, super_sample, ssample_param, count, segs)
1571
773bc121
TC
1572Creates a new general fill which fills with a fountain fill.
1573
f1ac5027
TC
1574=cut
1575*/
1576
1577i_fill_t *
1578i_new_fill_fount(double xa, double ya, double xb, double yb,
1579 i_fountain_type type, i_fountain_repeat repeat,
1580 int combine, int super_sample, double ssample_param,
1581 int count, i_fountain_seg *segs) {
1582 i_fill_fountain_t *fill = mymalloc(sizeof(i_fill_fountain_t));
1583
1584 fill->base.fill_with_color = NULL;
1585 fill->base.fill_with_fcolor = fill_fountf;
1586 fill->base.destroy = fount_fill_destroy;
efdc2568
TC
1587 if (combine)
1588 i_get_combine(combine, &fill->base.combine, &fill->base.combinef);
1589 else {
1590 fill->base.combine = NULL;
1591 fill->base.combinef = NULL;
1592 }
f1ac5027
TC
1593 fount_init_state(&fill->state, xa, ya, xb, yb, type, repeat, combine,
1594 super_sample, ssample_param, count, segs);
1595
1596 return &fill->base;
1597}
1598
1599/*
1600=back
1601
1602=head1 INTERNAL FUNCTIONS
1603
1604=over
1605
1606=item fount_init_state(...)
1607
1608Used by both the fountain fill filter and the fountain fill.
1609
1610=cut
1611*/
1612
1613static void
1614fount_init_state(struct fount_state *state, double xa, double ya,
1615 double xb, double yb, i_fountain_type type,
1616 i_fountain_repeat repeat, int combine, int super_sample,
1617 double ssample_param, int count, i_fountain_seg *segs) {
6607600c
TC
1618 int i, j;
1619 i_fountain_seg *my_segs = mymalloc(sizeof(i_fountain_seg) * count);
f1ac5027 1620 /*int have_alpha = im->channels == 2 || im->channels == 4;*/
f1ac5027
TC
1621
1622 memset(state, 0, sizeof(*state));
6607600c
TC
1623 /* we keep a local copy that we can adjust for speed */
1624 for (i = 0; i < count; ++i) {
1625 i_fountain_seg *seg = my_segs + i;
1626
1627 *seg = segs[i];
4c033fd4
TC
1628 if (seg->type < 0 || seg->type >= i_fst_end)
1629 seg->type = i_fst_linear;
6607600c
TC
1630 if (seg->color < 0 || seg->color >= i_fc_end)
1631 seg->color = i_fc_direct;
1632 if (seg->color == i_fc_hue_up || seg->color == i_fc_hue_down) {
1633 /* so we don't have to translate to HSV on each request, do it here */
1634 for (j = 0; j < 2; ++j) {
1635 i_rgb_to_hsvf(seg->c+j);
1636 }
1637 if (seg->color == i_fc_hue_up) {
1638 if (seg->c[1].channel[0] <= seg->c[0].channel[0])
1639 seg->c[1].channel[0] += 1.0;
1640 }
1641 else {
1642 if (seg->c[0].channel[0] <= seg->c[0].channel[1])
1643 seg->c[0].channel[0] += 1.0;
1644 }
1645 }
1646 /*printf("start %g mid %g end %g c0(%g,%g,%g,%g) c1(%g,%g,%g,%g) type %d color %d\n",
1647 seg->start, seg->middle, seg->end, seg->c[0].channel[0],
1648 seg->c[0].channel[1], seg->c[0].channel[2], seg->c[0].channel[3],
1649 seg->c[1].channel[0], seg->c[1].channel[1], seg->c[1].channel[2],
1650 seg->c[1].channel[3], seg->type, seg->color);*/
1651
1652 }
1653
1654 /* initialize each engine */
1655 /* these are so common ... */
f1ac5027
TC
1656 state->lA = xb - xa;
1657 state->lB = yb - ya;
1658 state->AB = sqrt(state->lA * state->lA + state->lB * state->lB);
1659 state->xa = xa;
1660 state->ya = ya;
6607600c
TC
1661 switch (type) {
1662 default:
1663 type = i_ft_linear; /* make the invalid value valid */
1664 case i_ft_linear:
1665 case i_ft_bilinear:
f1ac5027
TC
1666 state->lC = ya * ya - ya * yb + xa * xa - xa * xb;
1667 state->mult = 1;
1668 state->mult = 1/linear_fount_f(xb, yb, state);
6607600c
TC
1669 break;
1670
1671 case i_ft_radial:
f1ac5027
TC
1672 state->mult = 1.0 / sqrt((double)(xb-xa)*(xb-xa)
1673 + (double)(yb-ya)*(yb-ya));
6607600c
TC
1674 break;
1675
1676 case i_ft_radial_square:
f1ac5027
TC
1677 state->cos = state->lA / state->AB;
1678 state->sin = state->lB / state->AB;
1679 state->mult = 1.0 / state->AB;
6607600c
TC
1680 break;
1681
1682 case i_ft_revolution:
f1ac5027
TC
1683 state->theta = atan2(yb-ya, xb-xa);
1684 state->mult = 1.0 / (PI * 2);
6607600c
TC
1685 break;
1686
1687 case i_ft_conical:
f1ac5027
TC
1688 state->theta = atan2(yb-ya, xb-xa);
1689 state->mult = 1.0 / PI;
6607600c
TC
1690 break;
1691 }
f1ac5027 1692 state->ffunc = fount_funcs[type];
6607600c 1693 if (super_sample < 0
290bdf77 1694 || super_sample >= (int)(sizeof(fount_ssamples)/sizeof(*fount_ssamples))) {
6607600c
TC
1695 super_sample = 0;
1696 }
f1ac5027 1697 state->ssample_data = NULL;
6607600c
TC
1698 switch (super_sample) {
1699 case i_fts_grid:
1700 ssample_param = floor(0.5 + sqrt(ssample_param));
f1ac5027 1701 state->ssample_data = mymalloc(sizeof(i_fcolor) * ssample_param * ssample_param);
6607600c
TC
1702 break;
1703
1704 case i_fts_random:
1705 case i_fts_circle:
1706 ssample_param = floor(0.5+ssample_param);
f1ac5027 1707 state->ssample_data = mymalloc(sizeof(i_fcolor) * ssample_param);
6607600c
TC
1708 break;
1709 }
f1ac5027
TC
1710 state->parm = ssample_param;
1711 state->ssfunc = fount_ssamples[super_sample];
6607600c
TC
1712 if (repeat < 0 || repeat >= (sizeof(fount_repeats)/sizeof(*fount_repeats)))
1713 repeat = 0;
f1ac5027
TC
1714 state->rpfunc = fount_repeats[repeat];
1715 state->segs = my_segs;
1716 state->count = count;
6607600c
TC
1717}
1718
f1ac5027
TC
1719static void
1720fount_finish_state(struct fount_state *state) {
1721 if (state->ssample_data)
1722 myfree(state->ssample_data);
1723 myfree(state->segs);
1724}
6607600c 1725
6607600c 1726
f1ac5027 1727/*
6607600c
TC
1728=item fount_getat(out, x, y, ffunc, rpfunc, state, segs, count)
1729
1730Evaluates the fountain fill at the given point.
1731
1732This is called by both the non-super-sampling and super-sampling code.
1733
1734You might think that it would make sense to sample the fill parameter
1735instead, and combine those, but this breaks badly.
1736
1737=cut
1738*/
1739
1740static int
f1ac5027
TC
1741fount_getat(i_fcolor *out, double x, double y, struct fount_state *state) {
1742 double v = (state->rpfunc)((state->ffunc)(x, y, state));
6607600c
TC
1743 int i;
1744
1745 i = 0;
f1ac5027
TC
1746 while (i < state->count
1747 && (v < state->segs[i].start || v > state->segs[i].end)) {
6607600c
TC
1748 ++i;
1749 }
f1ac5027
TC
1750 if (i < state->count) {
1751 v = (fount_interps[state->segs[i].type])(v, state->segs+i);
1752 (fount_cinterps[state->segs[i].color])(out, v, state->segs+i);
6607600c
TC
1753 return 1;
1754 }
1755 else
1756 return 0;
1757}
1758
1759/*
1760=item linear_fount_f(x, y, state)
1761
1762Calculate the fill parameter for a linear fountain fill.
1763
1764Uses the point to line distance function, with some precalculation
1765done in i_fountain().
1766
1767=cut
1768*/
1769static double
1770linear_fount_f(double x, double y, struct fount_state *state) {
1771 return (state->lA * x + state->lB * y + state->lC) / state->AB * state->mult;
1772}
1773
1774/*
1775=item bilinear_fount_f(x, y, state)
1776
1777Calculate the fill parameter for a bi-linear fountain fill.
1778
1779=cut
1780*/
1781static double
1782bilinear_fount_f(double x, double y, struct fount_state *state) {
1783 return fabs((state->lA * x + state->lB * y + state->lC) / state->AB * state->mult);
1784}
1785
1786/*
1787=item radial_fount_f(x, y, state)
1788
1789Calculate the fill parameter for a radial fountain fill.
1790
1791Simply uses the distance function.
1792
1793=cut
1794 */
1795static double
1796radial_fount_f(double x, double y, struct fount_state *state) {
1797 return sqrt((double)(state->xa-x)*(state->xa-x)
1798 + (double)(state->ya-y)*(state->ya-y)) * state->mult;
1799}
1800
1801/*
1802=item square_fount_f(x, y, state)
1803
1804Calculate the fill parameter for a square fountain fill.
1805
1806Works by rotating the reference co-ordinate around the centre of the
1807square.
1808
1809=cut
1810*/
1811static double
1812square_fount_f(double x, double y, struct fount_state *state) {
1813 int xc, yc; /* centred on A */
1814 double xt, yt; /* rotated by theta */
1815 xc = x - state->xa;
1816 yc = y - state->ya;
1817 xt = fabs(xc * state->cos + yc * state->sin);
1818 yt = fabs(-xc * state->sin + yc * state->cos);
1819 return (xt > yt ? xt : yt) * state->mult;
1820}
1821
1822/*
1823=item revolution_fount_f(x, y, state)
1824
1825Calculates the fill parameter for the revolution fountain fill.
1826
1827=cut
1828*/
1829static double
1830revolution_fount_f(double x, double y, struct fount_state *state) {
1831 double angle = atan2(y - state->ya, x - state->xa);
1832
1833 angle -= state->theta;
1834 if (angle < 0) {
1835 angle = fmod(angle+ PI * 4, PI*2);
1836 }
1837
1838 return angle * state->mult;
1839}
1840
1841/*
1842=item conical_fount_f(x, y, state)
1843
1844Calculates the fill parameter for the conical fountain fill.
1845
1846=cut
1847*/
1848static double
1849conical_fount_f(double x, double y, struct fount_state *state) {
1850 double angle = atan2(y - state->ya, x - state->xa);
1851
1852 angle -= state->theta;
1853 if (angle < -PI)
1854 angle += PI * 2;
1855 else if (angle > PI)
1856 angle -= PI * 2;
1857
1858 return fabs(angle) * state->mult;
1859}
1860
1861/*
1862=item linear_interp(pos, seg)
1863
1864Calculates linear interpolation on the fill parameter. Breaks the
1865segment into 2 regions based in the I<middle> value.
1866
1867=cut
1868*/
1869static double
1870linear_interp(double pos, i_fountain_seg *seg) {
1871 if (pos < seg->middle) {
1872 double len = seg->middle - seg->start;
1873 if (len < EPSILON)
1874 return 0.0;
1875 else
1876 return (pos - seg->start) / len / 2;
1877 }
1878 else {
1879 double len = seg->end - seg->middle;
1880 if (len < EPSILON)
1881 return 1.0;
1882 else
1883 return 0.5 + (pos - seg->middle) / len / 2;
1884 }
1885}
1886
1887/*
1888=item sine_interp(pos, seg)
1889
1890Calculates sine function interpolation on the fill parameter.
1891
1892=cut
1893*/
1894static double
1895sine_interp(double pos, i_fountain_seg *seg) {
1896 /* I wonder if there's a simple way to smooth the transition for this */
1897 double work = linear_interp(pos, seg);
1898
1899 return (1-cos(work * PI))/2;
1900}
1901
1902/*
1903=item sphereup_interp(pos, seg)
1904
1905Calculates spherical interpolation on the fill parameter, with the cusp
1906at the low-end.
1907
1908=cut
1909*/
1910static double
1911sphereup_interp(double pos, i_fountain_seg *seg) {
1912 double work = linear_interp(pos, seg);
1913
1914 return sqrt(1.0 - (1-work) * (1-work));
1915}
1916
1917/*
1918=item spheredown_interp(pos, seg)
1919
1920Calculates spherical interpolation on the fill parameter, with the cusp
1921at the high-end.
1922
1923=cut
1924*/
1925static double
1926spheredown_interp(double pos, i_fountain_seg *seg) {
1927 double work = linear_interp(pos, seg);
1928
1929 return 1-sqrt(1.0 - work * work);
1930}
1931
1932/*
1933=item direct_cinterp(out, pos, seg)
1934
1935Calculates the fountain color based on direct scaling of the channels
1936of the color channels.
1937
1938=cut
1939*/
1940static void
1941direct_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
1942 int ch;
1943 for (ch = 0; ch < MAXCHANNELS; ++ch) {
1944 out->channel[ch] = seg->c[0].channel[ch] * (1 - pos)
1945 + seg->c[1].channel[ch] * pos;
1946 }
1947}
1948
1949/*
1950=item hue_up_cinterp(put, pos, seg)
1951
1952Calculates the fountain color based on scaling a HSV value. The hue
1953increases as the fill parameter increases.
1954
1955=cut
1956*/
1957static void
1958hue_up_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
1959 int ch;
1960 for (ch = 0; ch < MAXCHANNELS; ++ch) {
1961 out->channel[ch] = seg->c[0].channel[ch] * (1 - pos)
1962 + seg->c[1].channel[ch] * pos;
1963 }
1964 i_hsv_to_rgbf(out);
1965}
1966
1967/*
1968=item hue_down_cinterp(put, pos, seg)
1969
1970Calculates the fountain color based on scaling a HSV value. The hue
1971decreases as the fill parameter increases.
1972
1973=cut
1974*/
1975static void
1976hue_down_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
1977 int ch;
1978 for (ch = 0; ch < MAXCHANNELS; ++ch) {
1979 out->channel[ch] = seg->c[0].channel[ch] * (1 - pos)
1980 + seg->c[1].channel[ch] * pos;
1981 }
1982 i_hsv_to_rgbf(out);
1983}
1984
1985/*
1986=item simple_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
1987
1988Simple grid-based super-sampling.
1989
1990=cut
1991*/
1992static int
f1ac5027 1993simple_ssample(i_fcolor *out, double x, double y, struct fount_state *state) {
6607600c
TC
1994 i_fcolor *work = state->ssample_data;
1995 int dx, dy;
f1ac5027 1996 int grid = state->parm;
6607600c
TC
1997 double base = -0.5 + 0.5 / grid;
1998 double step = 1.0 / grid;
1999 int ch, i;
2000 int samp_count = 0;
2001
2002 for (dx = 0; dx < grid; ++dx) {
2003 for (dy = 0; dy < grid; ++dy) {
2004 if (fount_getat(work+samp_count, x + base + step * dx,
f1ac5027 2005 y + base + step * dy, state)) {
6607600c
TC
2006 ++samp_count;
2007 }
2008 }
2009 }
2010 for (ch = 0; ch < MAXCHANNELS; ++ch) {
2011 out->channel[ch] = 0;
2012 for (i = 0; i < samp_count; ++i) {
2013 out->channel[ch] += work[i].channel[ch];
2014 }
2015 /* we divide by 4 rather than samp_count since if there's only one valid
2016 sample it should be mostly transparent */
2017 out->channel[ch] /= grid * grid;
2018 }
2019 return samp_count;
2020}
2021
2022/*
2023=item random_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
2024
2025Random super-sampling.
2026
2027=cut
2028*/
2029static int
f1ac5027
TC
2030random_ssample(i_fcolor *out, double x, double y,
2031 struct fount_state *state) {
6607600c
TC
2032 i_fcolor *work = state->ssample_data;
2033 int i, ch;
f1ac5027 2034 int maxsamples = state->parm;
6607600c
TC
2035 double rand_scale = 1.0 / RAND_MAX;
2036 int samp_count = 0;
2037 for (i = 0; i < maxsamples; ++i) {
2038 if (fount_getat(work+samp_count, x - 0.5 + rand() * rand_scale,
f1ac5027 2039 y - 0.5 + rand() * rand_scale, state)) {
6607600c
TC
2040 ++samp_count;
2041 }
2042 }
2043 for (ch = 0; ch < MAXCHANNELS; ++ch) {
2044 out->channel[ch] = 0;
2045 for (i = 0; i < samp_count; ++i) {
2046 out->channel[ch] += work[i].channel[ch];
2047 }
2048 /* we divide by maxsamples rather than samp_count since if there's
2049 only one valid sample it should be mostly transparent */
2050 out->channel[ch] /= maxsamples;
2051 }
2052 return samp_count;
2053}
2054
2055/*
2056=item circle_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
2057
2058Super-sampling around the circumference of a circle.
2059
2060I considered saving the sin()/cos() values and transforming step-size
2061around the circle, but that's inaccurate, though it may not matter
2062much.
2063
2064=cut
2065 */
2066static int
f1ac5027
TC
2067circle_ssample(i_fcolor *out, double x, double y,
2068 struct fount_state *state) {
6607600c
TC
2069 i_fcolor *work = state->ssample_data;
2070 int i, ch;
f1ac5027 2071 int maxsamples = state->parm;
6607600c
TC
2072 double angle = 2 * PI / maxsamples;
2073 double radius = 0.3; /* semi-random */
2074 int samp_count = 0;
2075 for (i = 0; i < maxsamples; ++i) {
2076 if (fount_getat(work+samp_count, x + radius * cos(angle * i),
f1ac5027 2077 y + radius * sin(angle * i), state)) {
6607600c
TC
2078 ++samp_count;
2079 }
2080 }
2081 for (ch = 0; ch < MAXCHANNELS; ++ch) {
2082 out->channel[ch] = 0;
2083 for (i = 0; i < samp_count; ++i) {
2084 out->channel[ch] += work[i].channel[ch];
2085 }
2086 /* we divide by maxsamples rather than samp_count since if there's
2087 only one valid sample it should be mostly transparent */
2088 out->channel[ch] /= maxsamples;
2089 }
2090 return samp_count;
2091}
2092
2093/*
2094=item fount_r_none(v)
2095
2096Implements no repeats. Simply clamps the fill value.
2097
2098=cut
2099*/
2100static double
2101fount_r_none(double v) {
2102 return v < 0 ? 0 : v > 1 ? 1 : v;
2103}
2104
2105/*
2106=item fount_r_sawtooth(v)
2107
2108Implements sawtooth repeats. Clamps negative values and uses fmod()
2109on others.
2110
2111=cut
2112*/
2113static double
2114fount_r_sawtooth(double v) {
2115 return v < 0 ? 0 : fmod(v, 1.0);
2116}
2117
2118/*
2119=item fount_r_triangle(v)
2120
2121Implements triangle repeats. Clamps negative values, uses fmod to get
2122a range 0 through 2 and then adjusts values > 1.
2123
2124=cut
2125*/
2126static double
2127fount_r_triangle(double v) {
2128 if (v < 0)
2129 return 0;
2130 else {
2131 v = fmod(v, 2.0);
2132 return v > 1.0 ? 2.0 - v : v;
2133 }
2134}
2135
2136/*
2137=item fount_r_saw_both(v)
2138
2139Implements sawtooth repeats in the both postive and negative directions.
2140
2141Adjusts the value to be postive and then just uses fmod().
2142
2143=cut
2144*/
2145static double
2146fount_r_saw_both(double v) {
2147 if (v < 0)
2148 v += 1+(int)(-v);
2149 return fmod(v, 1.0);
2150}
2151
2152/*
2153=item fount_r_tri_both(v)
2154
2155Implements triangle repeats in the both postive and negative directions.
2156
2157Uses fmod on the absolute value, and then adjusts values > 1.
2158
2159=cut
2160*/
2161static double
2162fount_r_tri_both(double v) {
2163 v = fmod(fabs(v), 2.0);
2164 return v > 1.0 ? 2.0 - v : v;
2165}
2166
f1ac5027
TC
2167/*
2168=item fill_fountf(fill, x, y, width, channels, data)
2169
2170The fill function for fountain fills.
2171
2172=cut
2173*/
2174static void
2175fill_fountf(i_fill_t *fill, int x, int y, int width, int channels,
43c5dacb 2176 i_fcolor *data) {
f1ac5027 2177 i_fill_fountain_t *f = (i_fill_fountain_t *)fill;
efdc2568 2178
43c5dacb
TC
2179 while (width--) {
2180 i_fcolor c;
2181 int got_one;
2182
2183 if (f->state.ssfunc)
2184 got_one = f->state.ssfunc(&c, x, y, &f->state);
2185 else
2186 got_one = fount_getat(&c, x, y, &f->state);
2187
2188 *data++ = c;
2189
2190 ++x;
f1ac5027
TC
2191 }
2192}
2193
2194/*
2195=item fount_fill_destroy(fill)
2196
2197=cut
2198*/
2199static void
2200fount_fill_destroy(i_fill_t *fill) {
2201 i_fill_fountain_t *f = (i_fill_fountain_t *)fill;
2202 fount_finish_state(&f->state);
2203}
2204
6607600c
TC
2205/*
2206=back
2207
2208=head1 AUTHOR
2209
2210Arnar M. Hrafnkelsson <addi@umich.edu>
2211
2212Tony Cook <tony@develop-help.com> (i_fountain())
2213
2214=head1 SEE ALSO
2215
2216Imager(3)
2217
2218=cut
2219*/