Fixed i_transform2() so malloc(0) doesn't happen. Also corrected pod errors and
[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 */
963 fdist[p] = max(xd*xd, yd*yd); /* manhattan distance */
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 }
982
983}
984
02d1d628
AMH
985void
986i_nearest_color_foo(i_img *im, int num, int *xo, int *yo, i_color *ival, int dmeasure) {
987
a743c0a6 988 int p, x, y;
02d1d628
AMH
989 int xsize = im->xsize;
990 int ysize = im->ysize;
991
992 mm_log((1,"i_gradgen(im %p, num %d, xo %p, yo %p, ival %p, dmeasure %d)\n", im, num, xo, yo, ival, dmeasure));
993
994 for(p = 0; p<num; p++) {
995 mm_log((1,"i_gradgen: (%d, %d)\n", xo[p], yo[p]));
996 ICL_info(&ival[p]);
997 }
998
999 for(y = 0; y<ysize; y++) for(x = 0; x<xsize; x++) {
1000 int midx = 0;
1001 float mindist = 0;
1002 float curdist = 0;
1003
1004 int xd = x-xo[0];
1005 int yd = y-yo[0];
1006
1007 switch (dmeasure) {
1008 case 0: /* euclidean */
1009 mindist = sqrt(xd*xd + yd*yd); /* euclidean distance */
1010 break;
1011 case 1: /* euclidean squared */
1012 mindist = xd*xd + yd*yd; /* euclidean distance */
1013 break;
1014 case 2: /* euclidean squared */
1015 mindist = max(xd*xd, yd*yd); /* manhattan distance */
1016 break;
1017 default:
1018 m_fatal(3,"i_nearest_color: Unknown distance measure\n");
1019 }
1020
1021 for(p = 1; p<num; p++) {
1022 xd = x-xo[p];
1023 yd = y-yo[p];
1024 switch (dmeasure) {
1025 case 0: /* euclidean */
1026 curdist = sqrt(xd*xd + yd*yd); /* euclidean distance */
1027 break;
1028 case 1: /* euclidean squared */
1029 curdist = xd*xd + yd*yd; /* euclidean distance */
1030 break;
1031 case 2: /* euclidean squared */
1032 curdist = max(xd*xd, yd*yd); /* manhattan distance */
1033 break;
1034 default:
1035 m_fatal(3,"i_nearest_color: Unknown distance measure\n");
1036 }
1037 if (curdist < mindist) {
1038 mindist = curdist;
1039 midx = p;
1040 }
1041 }
1042 i_ppix(im, x, y, &ival[midx]);
1043 }
1044}
1045
02d1d628
AMH
1046void
1047i_nearest_color(i_img *im, int num, int *xo, int *yo, i_color *oval, int dmeasure) {
1048 i_color *ival;
1049 float *tval;
1050 float c1, c2;
1051 i_color val;
1052 int p, x, y, ch;
02d1d628
AMH
1053 int xsize = im->xsize;
1054 int ysize = im->ysize;
1055 int *cmatch;
1056
d7f707d8 1057 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
1058
1059 tval = mymalloc( sizeof(float)*num*im->channels );
1060 ival = mymalloc( sizeof(i_color)*num );
1061 cmatch = mymalloc( sizeof(int)*num );
1062
1063 for(p = 0; p<num; p++) {
1064 for(ch = 0; ch<im->channels; ch++) tval[ p * im->channels + ch] = 0;
1065 cmatch[p] = 0;
1066 }
1067
1068
1069 for(y = 0; y<ysize; y++) for(x = 0; x<xsize; x++) {
1070 int midx = 0;
1071 float mindist = 0;
1072 float curdist = 0;
1073
1074 int xd = x-xo[0];
1075 int yd = y-yo[0];
1076
1077 switch (dmeasure) {
1078 case 0: /* euclidean */
1079 mindist = sqrt(xd*xd + yd*yd); /* euclidean distance */
1080 break;
1081 case 1: /* euclidean squared */
1082 mindist = xd*xd + yd*yd; /* euclidean distance */
1083 break;
1084 case 2: /* euclidean squared */
1085 mindist = max(xd*xd, yd*yd); /* manhattan distance */
1086 break;
1087 default:
1088 m_fatal(3,"i_nearest_color: Unknown distance measure\n");
1089 }
1090
1091 for(p = 1; p<num; p++) {
1092 xd = x-xo[p];
1093 yd = y-yo[p];
1094 switch (dmeasure) {
1095 case 0: /* euclidean */
1096 curdist = sqrt(xd*xd + yd*yd); /* euclidean distance */
1097 break;
1098 case 1: /* euclidean squared */
1099 curdist = xd*xd + yd*yd; /* euclidean distance */
1100 break;
1101 case 2: /* euclidean squared */
1102 curdist = max(xd*xd, yd*yd); /* manhattan distance */
1103 break;
1104 default:
1105 m_fatal(3,"i_nearest_color: Unknown distance measure\n");
1106 }
1107 if (curdist < mindist) {
1108 mindist = curdist;
1109 midx = p;
1110 }
1111 }
1112
1113 cmatch[midx]++;
1114 i_gpix(im, x, y, &val);
1115 c2 = 1.0/(float)(cmatch[midx]);
1116 c1 = 1.0-c2;
1117
02d1d628
AMH
1118 for(ch = 0; ch<im->channels; ch++)
1119 tval[midx*im->channels + ch] = c1*tval[midx*im->channels + ch] + c2 * (float) val.channel[ch];
3bb1c1f3 1120
02d1d628
AMH
1121
1122 }
1123
1124 for(p = 0; p<num; p++) for(ch = 0; ch<im->channels; ch++) ival[p].channel[ch] = tval[p*im->channels + ch];
1125
1126 i_nearest_color_foo(im, num, xo, yo, ival, dmeasure);
1127}
6607600c 1128
b6381851
TC
1129/*
1130=item i_unsharp_mask(im, stddev, scale)
1131
1132Perform an usharp mask, which is defined as subtracting the blurred
1133image from double the original.
1134
1135=cut
1136*/
1137void i_unsharp_mask(i_img *im, double stddev, double scale) {
1138 i_img copy;
1139 int x, y, ch;
1140
1141 if (scale < 0)
1142 return;
1143 /* it really shouldn't ever be more than 1.0, but maybe ... */
1144 if (scale > 100)
1145 scale = 100;
1146
1147 i_copy(&copy, im);
1148 i_gaussian(&copy, stddev);
1149 if (im->bits == i_8_bits) {
1150 i_color *blur = mymalloc(im->xsize * sizeof(i_color) * 2);
1151 i_color *out = blur + im->xsize;
1152
1153 for (y = 0; y < im->ysize; ++y) {
1154 i_glin(&copy, 0, copy.xsize, y, blur);
1155 i_glin(im, 0, im->xsize, y, out);
1156 for (x = 0; x < im->xsize; ++x) {
1157 for (ch = 0; ch < im->channels; ++ch) {
1158 /*int temp = out[x].channel[ch] +
1159 scale * (out[x].channel[ch] - blur[x].channel[ch]);*/
1160 int temp = out[x].channel[ch] * 2 - blur[x].channel[ch];
1161 if (temp < 0)
1162 temp = 0;
1163 else if (temp > 255)
1164 temp = 255;
1165 out[x].channel[ch] = temp;
1166 }
1167 }
1168 i_plin(im, 0, im->xsize, y, out);
1169 }
1170
1171 myfree(blur);
1172 }
1173 else {
1174 i_fcolor *blur = mymalloc(im->xsize * sizeof(i_fcolor) * 2);
1175 i_fcolor *out = blur + im->xsize;
1176
1177 for (y = 0; y < im->ysize; ++y) {
1178 i_glinf(&copy, 0, copy.xsize, y, blur);
1179 i_glinf(im, 0, im->xsize, y, out);
1180 for (x = 0; x < im->xsize; ++x) {
1181 for (ch = 0; ch < im->channels; ++ch) {
1182 double temp = out[x].channel[ch] +
1183 scale * (out[x].channel[ch] - blur[x].channel[ch]);
1184 if (temp < 0)
1185 temp = 0;
1186 else if (temp > 1.0)
1187 temp = 1.0;
1188 out[x].channel[ch] = temp;
1189 }
1190 }
1191 i_plinf(im, 0, im->xsize, y, out);
1192 }
1193
1194 myfree(blur);
1195 }
1196 i_img_exorcise(&copy);
1197}
1198
f1ac5027 1199struct fount_state;
6607600c
TC
1200static double linear_fount_f(double x, double y, struct fount_state *state);
1201static double bilinear_fount_f(double x, double y, struct fount_state *state);
1202static double radial_fount_f(double x, double y, struct fount_state *state);
1203static double square_fount_f(double x, double y, struct fount_state *state);
1204static double revolution_fount_f(double x, double y,
1205 struct fount_state *state);
1206static double conical_fount_f(double x, double y, struct fount_state *state);
1207
1208typedef double (*fount_func)(double, double, struct fount_state *);
1209static fount_func fount_funcs[] =
1210{
1211 linear_fount_f,
1212 bilinear_fount_f,
1213 radial_fount_f,
1214 square_fount_f,
1215 revolution_fount_f,
1216 conical_fount_f,
1217};
1218
1219static double linear_interp(double pos, i_fountain_seg *seg);
1220static double sine_interp(double pos, i_fountain_seg *seg);
1221static double sphereup_interp(double pos, i_fountain_seg *seg);
1222static double spheredown_interp(double pos, i_fountain_seg *seg);
1223typedef double (*fount_interp)(double pos, i_fountain_seg *seg);
1224static fount_interp fount_interps[] =
1225{
1226 linear_interp,
1227 linear_interp,
1228 sine_interp,
1229 sphereup_interp,
1230 spheredown_interp,
1231};
1232
1233static void direct_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg);
1234static void hue_up_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg);
1235static void hue_down_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg);
1236typedef void (*fount_cinterp)(i_fcolor *out, double pos, i_fountain_seg *seg);
1237static fount_cinterp fount_cinterps[] =
1238{
1239 direct_cinterp,
1240 hue_up_cinterp,
1241 hue_down_cinterp,
1242};
1243
1244typedef double (*fount_repeat)(double v);
1245static double fount_r_none(double v);
1246static double fount_r_sawtooth(double v);
1247static double fount_r_triangle(double v);
1248static double fount_r_saw_both(double v);
1249static double fount_r_tri_both(double v);
1250static fount_repeat fount_repeats[] =
1251{
1252 fount_r_none,
1253 fount_r_sawtooth,
1254 fount_r_triangle,
1255 fount_r_saw_both,
1256 fount_r_tri_both,
1257};
1258
f1ac5027
TC
1259static int simple_ssample(i_fcolor *out, double x, double y,
1260 struct fount_state *state);
1261static int random_ssample(i_fcolor *out, double x, double y,
1262 struct fount_state *state);
1263static int circle_ssample(i_fcolor *out, double x, double y,
1264 struct fount_state *state);
1265typedef int (*fount_ssample)(i_fcolor *out, double x, double y,
1266 struct fount_state *state);
6607600c
TC
1267static fount_ssample fount_ssamples[] =
1268{
1269 NULL,
1270 simple_ssample,
1271 random_ssample,
1272 circle_ssample,
1273};
1274
1275static int
f1ac5027
TC
1276fount_getat(i_fcolor *out, double x, double y, struct fount_state *state);
1277
1278/*
1279 Keep state information used by each type of fountain fill
1280*/
1281struct fount_state {
1282 /* precalculated for the equation of the line perpendicular to the line AB */
1283 double lA, lB, lC;
1284 double AB;
1285 double sqrtA2B2;
1286 double mult;
1287 double cos;
1288 double sin;
1289 double theta;
1290 int xa, ya;
1291 void *ssample_data;
1292 fount_func ffunc;
1293 fount_repeat rpfunc;
1294 fount_ssample ssfunc;
1295 double parm;
1296 i_fountain_seg *segs;
1297 int count;
1298};
1299
1300static void
1301fount_init_state(struct fount_state *state, double xa, double ya,
1302 double xb, double yb, i_fountain_type type,
1303 i_fountain_repeat repeat, int combine, int super_sample,
1304 double ssample_param, int count, i_fountain_seg *segs);
1305
1306static void
1307fount_finish_state(struct fount_state *state);
6607600c
TC
1308
1309#define EPSILON (1e-6)
1310
1311/*
1312=item i_fountain(im, xa, ya, xb, yb, type, repeat, combine, super_sample, ssample_param, count, segs)
1313
1314Draws a fountain fill using A(xa, ya) and B(xb, yb) as reference points.
1315
1316I<type> controls how the reference points are used:
1317
1318=over
1319
1320=item i_ft_linear
1321
1322linear, where A is 0 and B is 1.
1323
1324=item i_ft_bilinear
1325
1326linear in both directions from A.
1327
1328=item i_ft_radial
1329
1330circular, where A is the centre of the fill, and B is a point
1331on the radius.
1332
1333=item i_ft_radial_square
1334
1335where A is the centre of the fill and B is the centre of
1336one side of the square.
1337
1338=item i_ft_revolution
1339
1340where A is the centre of the fill and B defines the 0/1.0
1341angle of the fill.
1342
1343=item i_ft_conical
1344
1345similar to i_ft_revolution, except that the revolution goes in both
1346directions
1347
1348=back
1349
1350I<repeat> can be one of:
1351
1352=over
1353
1354=item i_fr_none
1355
1356values < 0 are treated as zero, values > 1 are treated as 1.
1357
1358=item i_fr_sawtooth
1359
1360negative values are treated as 0, positive values are modulo 1.0
1361
1362=item i_fr_triangle
1363
1364negative values are treated as zero, if (int)value is odd then the value is treated as 1-(value
1365mod 1.0), otherwise the same as for sawtooth.
1366
1367=item i_fr_saw_both
1368
1369like i_fr_sawtooth, except that the sawtooth pattern repeats into
1370negative values.
1371
1372=item i_fr_tri_both
1373
1374Like i_fr_triangle, except that negative values are handled as their
1375absolute values.
1376
1377=back
1378
1379If combine is non-zero then non-opaque values are combined with the
1380underlying color.
1381
1382I<super_sample> controls super sampling, if any. At some point I'll
1383probably add a adaptive super-sampler. Current possible values are:
1384
1385=over
1386
1387=item i_fts_none
1388
1389No super-sampling is done.
1390
1391=item i_fts_grid
1392
1393A square grid of points withing the pixel are sampled.
1394
1395=item i_fts_random
1396
1397Random points within the pixel are sampled.
1398
1399=item i_fts_circle
1400
1401Points on the radius of a circle are sampled. This produces fairly
1402good results, but is fairly slow since sin() and cos() are evaluated
1403for each point.
1404
1405=back
1406
1407I<ssample_param> is intended to be roughly the number of points
1408sampled within the pixel.
1409
1410I<count> and I<segs> define the segments of the fill.
1411
1412=cut
1413
1414*/
1415
1416void
1417i_fountain(i_img *im, double xa, double ya, double xb, double yb,
1418 i_fountain_type type, i_fountain_repeat repeat,
1419 int combine, int super_sample, double ssample_param,
1420 int count, i_fountain_seg *segs) {
1421 struct fount_state state;
6607600c
TC
1422 int x, y;
1423 i_fcolor *line = mymalloc(sizeof(i_fcolor) * im->xsize);
efdc2568 1424 i_fcolor *work = NULL;
290bdf77 1425
f1ac5027 1426 i_fountain_seg *my_segs;
efdc2568
TC
1427 i_fill_combine_f combine_func = NULL;
1428 i_fill_combinef_f combinef_func = NULL;
1429
1430 i_get_combine(combine, &combine_func, &combinef_func);
1431 if (combinef_func)
1432 work = mymalloc(sizeof(i_fcolor) * im->xsize);
f1ac5027
TC
1433
1434 fount_init_state(&state, xa, ya, xb, yb, type, repeat, combine,
1435 super_sample, ssample_param, count, segs);
1436 my_segs = state.segs;
1437
1438 for (y = 0; y < im->ysize; ++y) {
1439 i_glinf(im, 0, im->xsize, y, line);
1440 for (x = 0; x < im->xsize; ++x) {
1441 i_fcolor c;
1442 int got_one;
f1ac5027
TC
1443 if (super_sample == i_fts_none)
1444 got_one = fount_getat(&c, x, y, &state);
1445 else
1446 got_one = state.ssfunc(&c, x, y, &state);
1447 if (got_one) {
efdc2568
TC
1448 if (combine)
1449 work[x] = c;
f1ac5027
TC
1450 else
1451 line[x] = c;
1452 }
1453 }
efdc2568
TC
1454 if (combine)
1455 combinef_func(line, work, im->channels, im->xsize);
f1ac5027
TC
1456 i_plinf(im, 0, im->xsize, y, line);
1457 }
1458 fount_finish_state(&state);
1459 myfree(line);
1460}
1461
1462typedef struct {
1463 i_fill_t base;
1464 struct fount_state state;
1465} i_fill_fountain_t;
1466
1467static void
1468fill_fountf(i_fill_t *fill, int x, int y, int width, int channels,
efdc2568 1469 i_fcolor *data, i_fcolor *work);
f1ac5027
TC
1470static void
1471fount_fill_destroy(i_fill_t *fill);
1472
1473/*
1474=item i_new_fount(xa, ya, xb, yb, type, repeat, combine, super_sample, ssample_param, count, segs)
1475
773bc121
TC
1476Creates a new general fill which fills with a fountain fill.
1477
f1ac5027
TC
1478=cut
1479*/
1480
1481i_fill_t *
1482i_new_fill_fount(double xa, double ya, double xb, double yb,
1483 i_fountain_type type, i_fountain_repeat repeat,
1484 int combine, int super_sample, double ssample_param,
1485 int count, i_fountain_seg *segs) {
1486 i_fill_fountain_t *fill = mymalloc(sizeof(i_fill_fountain_t));
1487
1488 fill->base.fill_with_color = NULL;
1489 fill->base.fill_with_fcolor = fill_fountf;
1490 fill->base.destroy = fount_fill_destroy;
efdc2568
TC
1491 if (combine)
1492 i_get_combine(combine, &fill->base.combine, &fill->base.combinef);
1493 else {
1494 fill->base.combine = NULL;
1495 fill->base.combinef = NULL;
1496 }
f1ac5027
TC
1497 fount_init_state(&fill->state, xa, ya, xb, yb, type, repeat, combine,
1498 super_sample, ssample_param, count, segs);
1499
1500 return &fill->base;
1501}
1502
1503/*
1504=back
1505
1506=head1 INTERNAL FUNCTIONS
1507
1508=over
1509
1510=item fount_init_state(...)
1511
1512Used by both the fountain fill filter and the fountain fill.
1513
1514=cut
1515*/
1516
1517static void
1518fount_init_state(struct fount_state *state, double xa, double ya,
1519 double xb, double yb, i_fountain_type type,
1520 i_fountain_repeat repeat, int combine, int super_sample,
1521 double ssample_param, int count, i_fountain_seg *segs) {
6607600c
TC
1522 int i, j;
1523 i_fountain_seg *my_segs = mymalloc(sizeof(i_fountain_seg) * count);
f1ac5027 1524 /*int have_alpha = im->channels == 2 || im->channels == 4;*/
f1ac5027
TC
1525
1526 memset(state, 0, sizeof(*state));
6607600c
TC
1527 /* we keep a local copy that we can adjust for speed */
1528 for (i = 0; i < count; ++i) {
1529 i_fountain_seg *seg = my_segs + i;
1530
1531 *seg = segs[i];
4c033fd4
TC
1532 if (seg->type < 0 || seg->type >= i_fst_end)
1533 seg->type = i_fst_linear;
6607600c
TC
1534 if (seg->color < 0 || seg->color >= i_fc_end)
1535 seg->color = i_fc_direct;
1536 if (seg->color == i_fc_hue_up || seg->color == i_fc_hue_down) {
1537 /* so we don't have to translate to HSV on each request, do it here */
1538 for (j = 0; j < 2; ++j) {
1539 i_rgb_to_hsvf(seg->c+j);
1540 }
1541 if (seg->color == i_fc_hue_up) {
1542 if (seg->c[1].channel[0] <= seg->c[0].channel[0])
1543 seg->c[1].channel[0] += 1.0;
1544 }
1545 else {
1546 if (seg->c[0].channel[0] <= seg->c[0].channel[1])
1547 seg->c[0].channel[0] += 1.0;
1548 }
1549 }
1550 /*printf("start %g mid %g end %g c0(%g,%g,%g,%g) c1(%g,%g,%g,%g) type %d color %d\n",
1551 seg->start, seg->middle, seg->end, seg->c[0].channel[0],
1552 seg->c[0].channel[1], seg->c[0].channel[2], seg->c[0].channel[3],
1553 seg->c[1].channel[0], seg->c[1].channel[1], seg->c[1].channel[2],
1554 seg->c[1].channel[3], seg->type, seg->color);*/
1555
1556 }
1557
1558 /* initialize each engine */
1559 /* these are so common ... */
f1ac5027
TC
1560 state->lA = xb - xa;
1561 state->lB = yb - ya;
1562 state->AB = sqrt(state->lA * state->lA + state->lB * state->lB);
1563 state->xa = xa;
1564 state->ya = ya;
6607600c
TC
1565 switch (type) {
1566 default:
1567 type = i_ft_linear; /* make the invalid value valid */
1568 case i_ft_linear:
1569 case i_ft_bilinear:
f1ac5027
TC
1570 state->lC = ya * ya - ya * yb + xa * xa - xa * xb;
1571 state->mult = 1;
1572 state->mult = 1/linear_fount_f(xb, yb, state);
6607600c
TC
1573 break;
1574
1575 case i_ft_radial:
f1ac5027
TC
1576 state->mult = 1.0 / sqrt((double)(xb-xa)*(xb-xa)
1577 + (double)(yb-ya)*(yb-ya));
6607600c
TC
1578 break;
1579
1580 case i_ft_radial_square:
f1ac5027
TC
1581 state->cos = state->lA / state->AB;
1582 state->sin = state->lB / state->AB;
1583 state->mult = 1.0 / state->AB;
6607600c
TC
1584 break;
1585
1586 case i_ft_revolution:
f1ac5027
TC
1587 state->theta = atan2(yb-ya, xb-xa);
1588 state->mult = 1.0 / (PI * 2);
6607600c
TC
1589 break;
1590
1591 case i_ft_conical:
f1ac5027
TC
1592 state->theta = atan2(yb-ya, xb-xa);
1593 state->mult = 1.0 / PI;
6607600c
TC
1594 break;
1595 }
f1ac5027 1596 state->ffunc = fount_funcs[type];
6607600c 1597 if (super_sample < 0
290bdf77 1598 || super_sample >= (int)(sizeof(fount_ssamples)/sizeof(*fount_ssamples))) {
6607600c
TC
1599 super_sample = 0;
1600 }
f1ac5027 1601 state->ssample_data = NULL;
6607600c
TC
1602 switch (super_sample) {
1603 case i_fts_grid:
1604 ssample_param = floor(0.5 + sqrt(ssample_param));
f1ac5027 1605 state->ssample_data = mymalloc(sizeof(i_fcolor) * ssample_param * ssample_param);
6607600c
TC
1606 break;
1607
1608 case i_fts_random:
1609 case i_fts_circle:
1610 ssample_param = floor(0.5+ssample_param);
f1ac5027 1611 state->ssample_data = mymalloc(sizeof(i_fcolor) * ssample_param);
6607600c
TC
1612 break;
1613 }
f1ac5027
TC
1614 state->parm = ssample_param;
1615 state->ssfunc = fount_ssamples[super_sample];
6607600c
TC
1616 if (repeat < 0 || repeat >= (sizeof(fount_repeats)/sizeof(*fount_repeats)))
1617 repeat = 0;
f1ac5027
TC
1618 state->rpfunc = fount_repeats[repeat];
1619 state->segs = my_segs;
1620 state->count = count;
6607600c
TC
1621}
1622
f1ac5027
TC
1623static void
1624fount_finish_state(struct fount_state *state) {
1625 if (state->ssample_data)
1626 myfree(state->ssample_data);
1627 myfree(state->segs);
1628}
6607600c 1629
6607600c 1630
f1ac5027 1631/*
6607600c
TC
1632=item fount_getat(out, x, y, ffunc, rpfunc, state, segs, count)
1633
1634Evaluates the fountain fill at the given point.
1635
1636This is called by both the non-super-sampling and super-sampling code.
1637
1638You might think that it would make sense to sample the fill parameter
1639instead, and combine those, but this breaks badly.
1640
1641=cut
1642*/
1643
1644static int
f1ac5027
TC
1645fount_getat(i_fcolor *out, double x, double y, struct fount_state *state) {
1646 double v = (state->rpfunc)((state->ffunc)(x, y, state));
6607600c
TC
1647 int i;
1648
1649 i = 0;
f1ac5027
TC
1650 while (i < state->count
1651 && (v < state->segs[i].start || v > state->segs[i].end)) {
6607600c
TC
1652 ++i;
1653 }
f1ac5027
TC
1654 if (i < state->count) {
1655 v = (fount_interps[state->segs[i].type])(v, state->segs+i);
1656 (fount_cinterps[state->segs[i].color])(out, v, state->segs+i);
6607600c
TC
1657 return 1;
1658 }
1659 else
1660 return 0;
1661}
1662
1663/*
1664=item linear_fount_f(x, y, state)
1665
1666Calculate the fill parameter for a linear fountain fill.
1667
1668Uses the point to line distance function, with some precalculation
1669done in i_fountain().
1670
1671=cut
1672*/
1673static double
1674linear_fount_f(double x, double y, struct fount_state *state) {
1675 return (state->lA * x + state->lB * y + state->lC) / state->AB * state->mult;
1676}
1677
1678/*
1679=item bilinear_fount_f(x, y, state)
1680
1681Calculate the fill parameter for a bi-linear fountain fill.
1682
1683=cut
1684*/
1685static double
1686bilinear_fount_f(double x, double y, struct fount_state *state) {
1687 return fabs((state->lA * x + state->lB * y + state->lC) / state->AB * state->mult);
1688}
1689
1690/*
1691=item radial_fount_f(x, y, state)
1692
1693Calculate the fill parameter for a radial fountain fill.
1694
1695Simply uses the distance function.
1696
1697=cut
1698 */
1699static double
1700radial_fount_f(double x, double y, struct fount_state *state) {
1701 return sqrt((double)(state->xa-x)*(state->xa-x)
1702 + (double)(state->ya-y)*(state->ya-y)) * state->mult;
1703}
1704
1705/*
1706=item square_fount_f(x, y, state)
1707
1708Calculate the fill parameter for a square fountain fill.
1709
1710Works by rotating the reference co-ordinate around the centre of the
1711square.
1712
1713=cut
1714*/
1715static double
1716square_fount_f(double x, double y, struct fount_state *state) {
1717 int xc, yc; /* centred on A */
1718 double xt, yt; /* rotated by theta */
1719 xc = x - state->xa;
1720 yc = y - state->ya;
1721 xt = fabs(xc * state->cos + yc * state->sin);
1722 yt = fabs(-xc * state->sin + yc * state->cos);
1723 return (xt > yt ? xt : yt) * state->mult;
1724}
1725
1726/*
1727=item revolution_fount_f(x, y, state)
1728
1729Calculates the fill parameter for the revolution fountain fill.
1730
1731=cut
1732*/
1733static double
1734revolution_fount_f(double x, double y, struct fount_state *state) {
1735 double angle = atan2(y - state->ya, x - state->xa);
1736
1737 angle -= state->theta;
1738 if (angle < 0) {
1739 angle = fmod(angle+ PI * 4, PI*2);
1740 }
1741
1742 return angle * state->mult;
1743}
1744
1745/*
1746=item conical_fount_f(x, y, state)
1747
1748Calculates the fill parameter for the conical fountain fill.
1749
1750=cut
1751*/
1752static double
1753conical_fount_f(double x, double y, struct fount_state *state) {
1754 double angle = atan2(y - state->ya, x - state->xa);
1755
1756 angle -= state->theta;
1757 if (angle < -PI)
1758 angle += PI * 2;
1759 else if (angle > PI)
1760 angle -= PI * 2;
1761
1762 return fabs(angle) * state->mult;
1763}
1764
1765/*
1766=item linear_interp(pos, seg)
1767
1768Calculates linear interpolation on the fill parameter. Breaks the
1769segment into 2 regions based in the I<middle> value.
1770
1771=cut
1772*/
1773static double
1774linear_interp(double pos, i_fountain_seg *seg) {
1775 if (pos < seg->middle) {
1776 double len = seg->middle - seg->start;
1777 if (len < EPSILON)
1778 return 0.0;
1779 else
1780 return (pos - seg->start) / len / 2;
1781 }
1782 else {
1783 double len = seg->end - seg->middle;
1784 if (len < EPSILON)
1785 return 1.0;
1786 else
1787 return 0.5 + (pos - seg->middle) / len / 2;
1788 }
1789}
1790
1791/*
1792=item sine_interp(pos, seg)
1793
1794Calculates sine function interpolation on the fill parameter.
1795
1796=cut
1797*/
1798static double
1799sine_interp(double pos, i_fountain_seg *seg) {
1800 /* I wonder if there's a simple way to smooth the transition for this */
1801 double work = linear_interp(pos, seg);
1802
1803 return (1-cos(work * PI))/2;
1804}
1805
1806/*
1807=item sphereup_interp(pos, seg)
1808
1809Calculates spherical interpolation on the fill parameter, with the cusp
1810at the low-end.
1811
1812=cut
1813*/
1814static double
1815sphereup_interp(double pos, i_fountain_seg *seg) {
1816 double work = linear_interp(pos, seg);
1817
1818 return sqrt(1.0 - (1-work) * (1-work));
1819}
1820
1821/*
1822=item spheredown_interp(pos, seg)
1823
1824Calculates spherical interpolation on the fill parameter, with the cusp
1825at the high-end.
1826
1827=cut
1828*/
1829static double
1830spheredown_interp(double pos, i_fountain_seg *seg) {
1831 double work = linear_interp(pos, seg);
1832
1833 return 1-sqrt(1.0 - work * work);
1834}
1835
1836/*
1837=item direct_cinterp(out, pos, seg)
1838
1839Calculates the fountain color based on direct scaling of the channels
1840of the color channels.
1841
1842=cut
1843*/
1844static void
1845direct_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
1846 int ch;
1847 for (ch = 0; ch < MAXCHANNELS; ++ch) {
1848 out->channel[ch] = seg->c[0].channel[ch] * (1 - pos)
1849 + seg->c[1].channel[ch] * pos;
1850 }
1851}
1852
1853/*
1854=item hue_up_cinterp(put, pos, seg)
1855
1856Calculates the fountain color based on scaling a HSV value. The hue
1857increases as the fill parameter increases.
1858
1859=cut
1860*/
1861static void
1862hue_up_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
1863 int ch;
1864 for (ch = 0; ch < MAXCHANNELS; ++ch) {
1865 out->channel[ch] = seg->c[0].channel[ch] * (1 - pos)
1866 + seg->c[1].channel[ch] * pos;
1867 }
1868 i_hsv_to_rgbf(out);
1869}
1870
1871/*
1872=item hue_down_cinterp(put, pos, seg)
1873
1874Calculates the fountain color based on scaling a HSV value. The hue
1875decreases as the fill parameter increases.
1876
1877=cut
1878*/
1879static void
1880hue_down_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
1881 int ch;
1882 for (ch = 0; ch < MAXCHANNELS; ++ch) {
1883 out->channel[ch] = seg->c[0].channel[ch] * (1 - pos)
1884 + seg->c[1].channel[ch] * pos;
1885 }
1886 i_hsv_to_rgbf(out);
1887}
1888
1889/*
1890=item simple_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
1891
1892Simple grid-based super-sampling.
1893
1894=cut
1895*/
1896static int
f1ac5027 1897simple_ssample(i_fcolor *out, double x, double y, struct fount_state *state) {
6607600c
TC
1898 i_fcolor *work = state->ssample_data;
1899 int dx, dy;
f1ac5027 1900 int grid = state->parm;
6607600c
TC
1901 double base = -0.5 + 0.5 / grid;
1902 double step = 1.0 / grid;
1903 int ch, i;
1904 int samp_count = 0;
1905
1906 for (dx = 0; dx < grid; ++dx) {
1907 for (dy = 0; dy < grid; ++dy) {
1908 if (fount_getat(work+samp_count, x + base + step * dx,
f1ac5027 1909 y + base + step * dy, state)) {
6607600c
TC
1910 ++samp_count;
1911 }
1912 }
1913 }
1914 for (ch = 0; ch < MAXCHANNELS; ++ch) {
1915 out->channel[ch] = 0;
1916 for (i = 0; i < samp_count; ++i) {
1917 out->channel[ch] += work[i].channel[ch];
1918 }
1919 /* we divide by 4 rather than samp_count since if there's only one valid
1920 sample it should be mostly transparent */
1921 out->channel[ch] /= grid * grid;
1922 }
1923 return samp_count;
1924}
1925
1926/*
1927=item random_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
1928
1929Random super-sampling.
1930
1931=cut
1932*/
1933static int
f1ac5027
TC
1934random_ssample(i_fcolor *out, double x, double y,
1935 struct fount_state *state) {
6607600c
TC
1936 i_fcolor *work = state->ssample_data;
1937 int i, ch;
f1ac5027 1938 int maxsamples = state->parm;
6607600c
TC
1939 double rand_scale = 1.0 / RAND_MAX;
1940 int samp_count = 0;
1941 for (i = 0; i < maxsamples; ++i) {
1942 if (fount_getat(work+samp_count, x - 0.5 + rand() * rand_scale,
f1ac5027 1943 y - 0.5 + rand() * rand_scale, state)) {
6607600c
TC
1944 ++samp_count;
1945 }
1946 }
1947 for (ch = 0; ch < MAXCHANNELS; ++ch) {
1948 out->channel[ch] = 0;
1949 for (i = 0; i < samp_count; ++i) {
1950 out->channel[ch] += work[i].channel[ch];
1951 }
1952 /* we divide by maxsamples rather than samp_count since if there's
1953 only one valid sample it should be mostly transparent */
1954 out->channel[ch] /= maxsamples;
1955 }
1956 return samp_count;
1957}
1958
1959/*
1960=item circle_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
1961
1962Super-sampling around the circumference of a circle.
1963
1964I considered saving the sin()/cos() values and transforming step-size
1965around the circle, but that's inaccurate, though it may not matter
1966much.
1967
1968=cut
1969 */
1970static int
f1ac5027
TC
1971circle_ssample(i_fcolor *out, double x, double y,
1972 struct fount_state *state) {
6607600c
TC
1973 i_fcolor *work = state->ssample_data;
1974 int i, ch;
f1ac5027 1975 int maxsamples = state->parm;
6607600c
TC
1976 double angle = 2 * PI / maxsamples;
1977 double radius = 0.3; /* semi-random */
1978 int samp_count = 0;
1979 for (i = 0; i < maxsamples; ++i) {
1980 if (fount_getat(work+samp_count, x + radius * cos(angle * i),
f1ac5027 1981 y + radius * sin(angle * i), state)) {
6607600c
TC
1982 ++samp_count;
1983 }
1984 }
1985 for (ch = 0; ch < MAXCHANNELS; ++ch) {
1986 out->channel[ch] = 0;
1987 for (i = 0; i < samp_count; ++i) {
1988 out->channel[ch] += work[i].channel[ch];
1989 }
1990 /* we divide by maxsamples rather than samp_count since if there's
1991 only one valid sample it should be mostly transparent */
1992 out->channel[ch] /= maxsamples;
1993 }
1994 return samp_count;
1995}
1996
1997/*
1998=item fount_r_none(v)
1999
2000Implements no repeats. Simply clamps the fill value.
2001
2002=cut
2003*/
2004static double
2005fount_r_none(double v) {
2006 return v < 0 ? 0 : v > 1 ? 1 : v;
2007}
2008
2009/*
2010=item fount_r_sawtooth(v)
2011
2012Implements sawtooth repeats. Clamps negative values and uses fmod()
2013on others.
2014
2015=cut
2016*/
2017static double
2018fount_r_sawtooth(double v) {
2019 return v < 0 ? 0 : fmod(v, 1.0);
2020}
2021
2022/*
2023=item fount_r_triangle(v)
2024
2025Implements triangle repeats. Clamps negative values, uses fmod to get
2026a range 0 through 2 and then adjusts values > 1.
2027
2028=cut
2029*/
2030static double
2031fount_r_triangle(double v) {
2032 if (v < 0)
2033 return 0;
2034 else {
2035 v = fmod(v, 2.0);
2036 return v > 1.0 ? 2.0 - v : v;
2037 }
2038}
2039
2040/*
2041=item fount_r_saw_both(v)
2042
2043Implements sawtooth repeats in the both postive and negative directions.
2044
2045Adjusts the value to be postive and then just uses fmod().
2046
2047=cut
2048*/
2049static double
2050fount_r_saw_both(double v) {
2051 if (v < 0)
2052 v += 1+(int)(-v);
2053 return fmod(v, 1.0);
2054}
2055
2056/*
2057=item fount_r_tri_both(v)
2058
2059Implements triangle repeats in the both postive and negative directions.
2060
2061Uses fmod on the absolute value, and then adjusts values > 1.
2062
2063=cut
2064*/
2065static double
2066fount_r_tri_both(double v) {
2067 v = fmod(fabs(v), 2.0);
2068 return v > 1.0 ? 2.0 - v : v;
2069}
2070
f1ac5027
TC
2071/*
2072=item fill_fountf(fill, x, y, width, channels, data)
2073
2074The fill function for fountain fills.
2075
2076=cut
2077*/
2078static void
2079fill_fountf(i_fill_t *fill, int x, int y, int width, int channels,
efdc2568 2080 i_fcolor *data, i_fcolor *work) {
f1ac5027 2081 i_fill_fountain_t *f = (i_fill_fountain_t *)fill;
efdc2568
TC
2082
2083 if (fill->combinef) {
2084 i_fcolor *wstart = work;
2085 int count = width;
f1ac5027 2086
efdc2568
TC
2087 while (width--) {
2088 i_fcolor c;
2089 int got_one;
290bdf77 2090
efdc2568
TC
2091 if (f->state.ssfunc)
2092 got_one = f->state.ssfunc(&c, x, y, &f->state);
2093 else
2094 got_one = fount_getat(&c, x, y, &f->state);
2095
2096 *work++ = c;
2097
2098 ++x;
2099 }
2100 (fill->combinef)(data, wstart, channels, count);
2101 }
2102 else {
2103 while (width--) {
2104 i_fcolor c;
2105 int got_one;
290bdf77 2106
efdc2568
TC
2107 if (f->state.ssfunc)
2108 got_one = f->state.ssfunc(&c, x, y, &f->state);
2109 else
2110 got_one = fount_getat(&c, x, y, &f->state);
2111
2112 *data++ = c;
2113
2114 ++x;
f1ac5027 2115 }
f1ac5027
TC
2116 }
2117}
2118
2119/*
2120=item fount_fill_destroy(fill)
2121
2122=cut
2123*/
2124static void
2125fount_fill_destroy(i_fill_t *fill) {
2126 i_fill_fountain_t *f = (i_fill_fountain_t *)fill;
2127 fount_finish_state(&f->state);
2128}
2129
6607600c
TC
2130/*
2131=back
2132
2133=head1 AUTHOR
2134
2135Arnar M. Hrafnkelsson <addi@umich.edu>
2136
2137Tony Cook <tony@develop-help.com> (i_fountain())
2138
2139=head1 SEE ALSO
2140
2141Imager(3)
2142
2143=cut
2144*/