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