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