most numeric parameters to the XS implementation now throw an exception if supplied...
[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);
1450
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
f1ac5027 1527struct fount_state;
6607600c
TC
1528static double linear_fount_f(double x, double y, struct fount_state *state);
1529static double bilinear_fount_f(double x, double y, struct fount_state *state);
1530static double radial_fount_f(double x, double y, struct fount_state *state);
1531static double square_fount_f(double x, double y, struct fount_state *state);
1532static double revolution_fount_f(double x, double y,
1533 struct fount_state *state);
1534static double conical_fount_f(double x, double y, struct fount_state *state);
1535
1536typedef double (*fount_func)(double, double, struct fount_state *);
1537static fount_func fount_funcs[] =
1538{
1539 linear_fount_f,
1540 bilinear_fount_f,
1541 radial_fount_f,
1542 square_fount_f,
1543 revolution_fount_f,
1544 conical_fount_f,
1545};
1546
1547static double linear_interp(double pos, i_fountain_seg *seg);
1548static double sine_interp(double pos, i_fountain_seg *seg);
1549static double sphereup_interp(double pos, i_fountain_seg *seg);
1550static double spheredown_interp(double pos, i_fountain_seg *seg);
1551typedef double (*fount_interp)(double pos, i_fountain_seg *seg);
1552static fount_interp fount_interps[] =
1553{
1554 linear_interp,
1555 linear_interp,
1556 sine_interp,
1557 sphereup_interp,
1558 spheredown_interp,
1559};
1560
1561static void direct_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg);
1562static void hue_up_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg);
1563static void hue_down_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg);
1564typedef void (*fount_cinterp)(i_fcolor *out, double pos, i_fountain_seg *seg);
1565static fount_cinterp fount_cinterps[] =
1566{
1567 direct_cinterp,
1568 hue_up_cinterp,
1569 hue_down_cinterp,
1570};
1571
1572typedef double (*fount_repeat)(double v);
1573static double fount_r_none(double v);
1574static double fount_r_sawtooth(double v);
1575static double fount_r_triangle(double v);
1576static double fount_r_saw_both(double v);
1577static double fount_r_tri_both(double v);
1578static fount_repeat fount_repeats[] =
1579{
1580 fount_r_none,
1581 fount_r_sawtooth,
1582 fount_r_triangle,
1583 fount_r_saw_both,
1584 fount_r_tri_both,
1585};
1586
f1ac5027
TC
1587static int simple_ssample(i_fcolor *out, double x, double y,
1588 struct fount_state *state);
1589static int random_ssample(i_fcolor *out, double x, double y,
1590 struct fount_state *state);
1591static int circle_ssample(i_fcolor *out, double x, double y,
1592 struct fount_state *state);
1593typedef int (*fount_ssample)(i_fcolor *out, double x, double y,
1594 struct fount_state *state);
6607600c
TC
1595static fount_ssample fount_ssamples[] =
1596{
1597 NULL,
1598 simple_ssample,
1599 random_ssample,
1600 circle_ssample,
1601};
1602
1603static int
f1ac5027
TC
1604fount_getat(i_fcolor *out, double x, double y, struct fount_state *state);
1605
1606/*
1607 Keep state information used by each type of fountain fill
1608*/
1609struct fount_state {
1610 /* precalculated for the equation of the line perpendicular to the line AB */
1611 double lA, lB, lC;
1612 double AB;
1613 double sqrtA2B2;
1614 double mult;
1615 double cos;
1616 double sin;
1617 double theta;
8d14daab 1618 i_img_dim xa, ya;
f1ac5027
TC
1619 void *ssample_data;
1620 fount_func ffunc;
1621 fount_repeat rpfunc;
1622 fount_ssample ssfunc;
1623 double parm;
1624 i_fountain_seg *segs;
1625 int count;
1626};
1627
1628static void
1629fount_init_state(struct fount_state *state, double xa, double ya,
1630 double xb, double yb, i_fountain_type type,
1631 i_fountain_repeat repeat, int combine, int super_sample,
1632 double ssample_param, int count, i_fountain_seg *segs);
1633
1634static void
1635fount_finish_state(struct fount_state *state);
6607600c
TC
1636
1637#define EPSILON (1e-6)
1638
1639/*
1640=item i_fountain(im, xa, ya, xb, yb, type, repeat, combine, super_sample, ssample_param, count, segs)
1641
1642Draws a fountain fill using A(xa, ya) and B(xb, yb) as reference points.
1643
1644I<type> controls how the reference points are used:
1645
1646=over
1647
1648=item i_ft_linear
1649
1650linear, where A is 0 and B is 1.
1651
1652=item i_ft_bilinear
1653
1654linear in both directions from A.
1655
1656=item i_ft_radial
1657
1658circular, where A is the centre of the fill, and B is a point
1659on the radius.
1660
1661=item i_ft_radial_square
1662
1663where A is the centre of the fill and B is the centre of
1664one side of the square.
1665
1666=item i_ft_revolution
1667
1668where A is the centre of the fill and B defines the 0/1.0
1669angle of the fill.
1670
1671=item i_ft_conical
1672
1673similar to i_ft_revolution, except that the revolution goes in both
1674directions
1675
1676=back
1677
1678I<repeat> can be one of:
1679
1680=over
1681
1682=item i_fr_none
1683
1684values < 0 are treated as zero, values > 1 are treated as 1.
1685
1686=item i_fr_sawtooth
1687
1688negative values are treated as 0, positive values are modulo 1.0
1689
1690=item i_fr_triangle
1691
1692negative values are treated as zero, if (int)value is odd then the value is treated as 1-(value
1693mod 1.0), otherwise the same as for sawtooth.
1694
1695=item i_fr_saw_both
1696
1697like i_fr_sawtooth, except that the sawtooth pattern repeats into
1698negative values.
1699
1700=item i_fr_tri_both
1701
1702Like i_fr_triangle, except that negative values are handled as their
1703absolute values.
1704
1705=back
1706
1707If combine is non-zero then non-opaque values are combined with the
1708underlying color.
1709
1710I<super_sample> controls super sampling, if any. At some point I'll
1711probably add a adaptive super-sampler. Current possible values are:
1712
1713=over
1714
1715=item i_fts_none
1716
1717No super-sampling is done.
1718
1719=item i_fts_grid
1720
1721A square grid of points withing the pixel are sampled.
1722
1723=item i_fts_random
1724
1725Random points within the pixel are sampled.
1726
1727=item i_fts_circle
1728
1729Points on the radius of a circle are sampled. This produces fairly
1730good results, but is fairly slow since sin() and cos() are evaluated
1731for each point.
1732
1733=back
1734
1735I<ssample_param> is intended to be roughly the number of points
1736sampled within the pixel.
1737
1738I<count> and I<segs> define the segments of the fill.
1739
1740=cut
1741
1742*/
1743
e310e5f9 1744int
6607600c
TC
1745i_fountain(i_img *im, double xa, double ya, double xb, double yb,
1746 i_fountain_type type, i_fountain_repeat repeat,
1747 int combine, int super_sample, double ssample_param,
1748 int count, i_fountain_seg *segs) {
1749 struct fount_state state;
8d14daab 1750 i_img_dim x, y;
e310e5f9 1751 i_fcolor *line = NULL;
efdc2568 1752 i_fcolor *work = NULL;
8d14daab 1753 size_t line_bytes;
efdc2568
TC
1754 i_fill_combine_f combine_func = NULL;
1755 i_fill_combinef_f combinef_func = NULL;
d17ee717 1756 dIMCTXim(im);
efdc2568 1757
e310e5f9
TC
1758 i_clear_error();
1759
1760 /* i_fountain() allocates floating colors even for 8-bit images,
1761 so we need to do this check */
1762 line_bytes = sizeof(i_fcolor) * im->xsize;
1763 if (line_bytes / sizeof(i_fcolor) != im->xsize) {
1764 i_push_error(0, "integer overflow calculating memory allocation");
1765 return 0;
1766 }
1767
1768 line = mymalloc(line_bytes); /* checked 17feb2005 tonyc */
1769
efdc2568 1770 i_get_combine(combine, &combine_func, &combinef_func);
56ce7446 1771 if (combinef_func) {
e310e5f9 1772 work = mymalloc(line_bytes); /* checked 17feb2005 tonyc */
56ce7446 1773 }
f1ac5027
TC
1774
1775 fount_init_state(&state, xa, ya, xb, yb, type, repeat, combine,
e4330a7b 1776 super_sample, ssample_param, count, segs);
f1ac5027
TC
1777
1778 for (y = 0; y < im->ysize; ++y) {
1779 i_glinf(im, 0, im->xsize, y, line);
1780 for (x = 0; x < im->xsize; ++x) {
1781 i_fcolor c;
1782 int got_one;
f1ac5027
TC
1783 if (super_sample == i_fts_none)
1784 got_one = fount_getat(&c, x, y, &state);
1785 else
1786 got_one = state.ssfunc(&c, x, y, &state);
1787 if (got_one) {
56ce7446 1788 if (combinef_func)
efdc2568 1789 work[x] = c;
f1ac5027
TC
1790 else
1791 line[x] = c;
1792 }
1793 }
56ce7446 1794 if (combinef_func)
efdc2568 1795 combinef_func(line, work, im->channels, im->xsize);
f1ac5027
TC
1796 i_plinf(im, 0, im->xsize, y, line);
1797 }
1798 fount_finish_state(&state);
56ce7446 1799 myfree(work);
f1ac5027 1800 myfree(line);
e310e5f9
TC
1801
1802 return 1;
f1ac5027
TC
1803}
1804
1805typedef struct {
1806 i_fill_t base;
1807 struct fount_state state;
1808} i_fill_fountain_t;
1809
1810static void
50c75381 1811fill_fountf(i_fill_t *fill, i_img_dim x, i_img_dim y, i_img_dim width, int channels,
43c5dacb 1812 i_fcolor *data);
f1ac5027
TC
1813static void
1814fount_fill_destroy(i_fill_t *fill);
1815
9b1ec2b8
TC
1816static i_fill_fountain_t
1817fount_fill_proto =
1818 {
1819 {
1820 NULL,
1821 fill_fountf,
1822 fount_fill_destroy
1823 }
1824 };
1825
1826
f1ac5027 1827/*
5715f7c3 1828=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
1829
1830=category Fills
1831=synopsis fill = i_new_fill_fount(0, 0, 100, 100, i_ft_linear, i_ft_linear,
1832=synopsis i_fr_triangle, 0, i_fts_grid, 9, 1, segs);
1833
f1ac5027 1834
773bc121
TC
1835Creates a new general fill which fills with a fountain fill.
1836
f1ac5027
TC
1837=cut
1838*/
1839
1840i_fill_t *
1841i_new_fill_fount(double xa, double ya, double xb, double yb,
1842 i_fountain_type type, i_fountain_repeat repeat,
1843 int combine, int super_sample, double ssample_param,
1844 int count, i_fountain_seg *segs) {
1845 i_fill_fountain_t *fill = mymalloc(sizeof(i_fill_fountain_t));
1846
9b1ec2b8 1847 *fill = fount_fill_proto;
efdc2568
TC
1848 if (combine)
1849 i_get_combine(combine, &fill->base.combine, &fill->base.combinef);
1850 else {
1851 fill->base.combine = NULL;
1852 fill->base.combinef = NULL;
1853 }
f1ac5027
TC
1854 fount_init_state(&fill->state, xa, ya, xb, yb, type, repeat, combine,
1855 super_sample, ssample_param, count, segs);
1856
1857 return &fill->base;
1858}
1859
1860/*
1861=back
1862
1863=head1 INTERNAL FUNCTIONS
1864
1865=over
1866
1867=item fount_init_state(...)
1868
1869Used by both the fountain fill filter and the fountain fill.
1870
1871=cut
1872*/
1873
1874static void
1875fount_init_state(struct fount_state *state, double xa, double ya,
1876 double xb, double yb, i_fountain_type type,
1877 i_fountain_repeat repeat, int combine, int super_sample,
1878 double ssample_param, int count, i_fountain_seg *segs) {
6607600c 1879 int i, j;
8d14daab 1880 size_t bytes;
d3d01aeb 1881 i_fountain_seg *my_segs = mymalloc(sizeof(i_fountain_seg) * count); /* checked 2jul06 - duplicating original */
f1ac5027 1882 /*int have_alpha = im->channels == 2 || im->channels == 4;*/
f1ac5027
TC
1883
1884 memset(state, 0, sizeof(*state));
6607600c
TC
1885 /* we keep a local copy that we can adjust for speed */
1886 for (i = 0; i < count; ++i) {
1887 i_fountain_seg *seg = my_segs + i;
1888
1889 *seg = segs[i];
4c033fd4
TC
1890 if (seg->type < 0 || seg->type >= i_fst_end)
1891 seg->type = i_fst_linear;
6607600c
TC
1892 if (seg->color < 0 || seg->color >= i_fc_end)
1893 seg->color = i_fc_direct;
1894 if (seg->color == i_fc_hue_up || seg->color == i_fc_hue_down) {
1895 /* so we don't have to translate to HSV on each request, do it here */
1896 for (j = 0; j < 2; ++j) {
1897 i_rgb_to_hsvf(seg->c+j);
1898 }
1899 if (seg->color == i_fc_hue_up) {
1900 if (seg->c[1].channel[0] <= seg->c[0].channel[0])
1901 seg->c[1].channel[0] += 1.0;
1902 }
1903 else {
1904 if (seg->c[0].channel[0] <= seg->c[0].channel[1])
1905 seg->c[0].channel[0] += 1.0;
1906 }
1907 }
1908 /*printf("start %g mid %g end %g c0(%g,%g,%g,%g) c1(%g,%g,%g,%g) type %d color %d\n",
1909 seg->start, seg->middle, seg->end, seg->c[0].channel[0],
1910 seg->c[0].channel[1], seg->c[0].channel[2], seg->c[0].channel[3],
1911 seg->c[1].channel[0], seg->c[1].channel[1], seg->c[1].channel[2],
1912 seg->c[1].channel[3], seg->type, seg->color);*/
1913
1914 }
1915
1916 /* initialize each engine */
1917 /* these are so common ... */
f1ac5027
TC
1918 state->lA = xb - xa;
1919 state->lB = yb - ya;
1920 state->AB = sqrt(state->lA * state->lA + state->lB * state->lB);
1921 state->xa = xa;
1922 state->ya = ya;
6607600c
TC
1923 switch (type) {
1924 default:
1925 type = i_ft_linear; /* make the invalid value valid */
1926 case i_ft_linear:
1927 case i_ft_bilinear:
f1ac5027
TC
1928 state->lC = ya * ya - ya * yb + xa * xa - xa * xb;
1929 state->mult = 1;
1930 state->mult = 1/linear_fount_f(xb, yb, state);
6607600c
TC
1931 break;
1932
1933 case i_ft_radial:
f1ac5027
TC
1934 state->mult = 1.0 / sqrt((double)(xb-xa)*(xb-xa)
1935 + (double)(yb-ya)*(yb-ya));
6607600c
TC
1936 break;
1937
1938 case i_ft_radial_square:
f1ac5027
TC
1939 state->cos = state->lA / state->AB;
1940 state->sin = state->lB / state->AB;
1941 state->mult = 1.0 / state->AB;
6607600c
TC
1942 break;
1943
1944 case i_ft_revolution:
f1ac5027
TC
1945 state->theta = atan2(yb-ya, xb-xa);
1946 state->mult = 1.0 / (PI * 2);
6607600c
TC
1947 break;
1948
1949 case i_ft_conical:
f1ac5027
TC
1950 state->theta = atan2(yb-ya, xb-xa);
1951 state->mult = 1.0 / PI;
6607600c
TC
1952 break;
1953 }
f1ac5027 1954 state->ffunc = fount_funcs[type];
6607600c 1955 if (super_sample < 0
290bdf77 1956 || super_sample >= (int)(sizeof(fount_ssamples)/sizeof(*fount_ssamples))) {
6607600c
TC
1957 super_sample = 0;
1958 }
f1ac5027 1959 state->ssample_data = NULL;
6607600c
TC
1960 switch (super_sample) {
1961 case i_fts_grid:
1962 ssample_param = floor(0.5 + sqrt(ssample_param));
d3d01aeb
TC
1963 bytes = ssample_param * ssample_param * sizeof(i_fcolor);
1964 if (bytes / sizeof(i_fcolor) == ssample_param * ssample_param) {
1965 state->ssample_data = mymalloc(sizeof(i_fcolor) * ssample_param * ssample_param); /* checked 1jul06 tonyc */
1966 }
1967 else {
1968 super_sample = i_fts_none;
1969 }
6607600c
TC
1970 break;
1971
1972 case i_fts_random:
1973 case i_fts_circle:
1974 ssample_param = floor(0.5+ssample_param);
a2c2d400
TC
1975 if (im_size_t_max / sizeof(i_fcolor) > ssample_param) {
1976 bytes = sizeof(i_fcolor) * ssample_param;
1977 state->ssample_data = mymalloc(bytes);
d3d01aeb
TC
1978 }
1979 else {
a2c2d400
TC
1980 dIMCTX;
1981 im_log((aIMCTX, 1,"size_t overflow calculating super-sample array size for random or circl"));
d3d01aeb
TC
1982 super_sample = i_fts_none;
1983 }
6607600c
TC
1984 break;
1985 }
f1ac5027
TC
1986 state->parm = ssample_param;
1987 state->ssfunc = fount_ssamples[super_sample];
6607600c
TC
1988 if (repeat < 0 || repeat >= (sizeof(fount_repeats)/sizeof(*fount_repeats)))
1989 repeat = 0;
f1ac5027
TC
1990 state->rpfunc = fount_repeats[repeat];
1991 state->segs = my_segs;
1992 state->count = count;
6607600c
TC
1993}
1994
f1ac5027
TC
1995static void
1996fount_finish_state(struct fount_state *state) {
1997 if (state->ssample_data)
1998 myfree(state->ssample_data);
1999 myfree(state->segs);
2000}
6607600c 2001
6607600c 2002
f1ac5027 2003/*
6607600c
TC
2004=item fount_getat(out, x, y, ffunc, rpfunc, state, segs, count)
2005
2006Evaluates the fountain fill at the given point.
2007
2008This is called by both the non-super-sampling and super-sampling code.
2009
2010You might think that it would make sense to sample the fill parameter
2011instead, and combine those, but this breaks badly.
2012
2013=cut
2014*/
2015
2016static int
f1ac5027
TC
2017fount_getat(i_fcolor *out, double x, double y, struct fount_state *state) {
2018 double v = (state->rpfunc)((state->ffunc)(x, y, state));
6607600c
TC
2019 int i;
2020
2021 i = 0;
f1ac5027
TC
2022 while (i < state->count
2023 && (v < state->segs[i].start || v > state->segs[i].end)) {
6607600c
TC
2024 ++i;
2025 }
f1ac5027
TC
2026 if (i < state->count) {
2027 v = (fount_interps[state->segs[i].type])(v, state->segs+i);
2028 (fount_cinterps[state->segs[i].color])(out, v, state->segs+i);
6607600c
TC
2029 return 1;
2030 }
2031 else
2032 return 0;
2033}
2034
2035/*
2036=item linear_fount_f(x, y, state)
2037
2038Calculate the fill parameter for a linear fountain fill.
2039
2040Uses the point to line distance function, with some precalculation
2041done in i_fountain().
2042
2043=cut
2044*/
2045static double
2046linear_fount_f(double x, double y, struct fount_state *state) {
2047 return (state->lA * x + state->lB * y + state->lC) / state->AB * state->mult;
2048}
2049
2050/*
2051=item bilinear_fount_f(x, y, state)
2052
2053Calculate the fill parameter for a bi-linear fountain fill.
2054
2055=cut
2056*/
2057static double
2058bilinear_fount_f(double x, double y, struct fount_state *state) {
2059 return fabs((state->lA * x + state->lB * y + state->lC) / state->AB * state->mult);
2060}
2061
2062/*
2063=item radial_fount_f(x, y, state)
2064
2065Calculate the fill parameter for a radial fountain fill.
2066
2067Simply uses the distance function.
2068
2069=cut
2070 */
2071static double
2072radial_fount_f(double x, double y, struct fount_state *state) {
2073 return sqrt((double)(state->xa-x)*(state->xa-x)
2074 + (double)(state->ya-y)*(state->ya-y)) * state->mult;
2075}
2076
2077/*
2078=item square_fount_f(x, y, state)
2079
2080Calculate the fill parameter for a square fountain fill.
2081
2082Works by rotating the reference co-ordinate around the centre of the
2083square.
2084
2085=cut
2086*/
2087static double
2088square_fount_f(double x, double y, struct fount_state *state) {
8d14daab 2089 i_img_dim xc, yc; /* centred on A */
6607600c
TC
2090 double xt, yt; /* rotated by theta */
2091 xc = x - state->xa;
2092 yc = y - state->ya;
2093 xt = fabs(xc * state->cos + yc * state->sin);
2094 yt = fabs(-xc * state->sin + yc * state->cos);
2095 return (xt > yt ? xt : yt) * state->mult;
2096}
2097
2098/*
2099=item revolution_fount_f(x, y, state)
2100
2101Calculates the fill parameter for the revolution fountain fill.
2102
2103=cut
2104*/
2105static double
2106revolution_fount_f(double x, double y, struct fount_state *state) {
2107 double angle = atan2(y - state->ya, x - state->xa);
2108
2109 angle -= state->theta;
2110 if (angle < 0) {
2111 angle = fmod(angle+ PI * 4, PI*2);
2112 }
2113
2114 return angle * state->mult;
2115}
2116
2117/*
2118=item conical_fount_f(x, y, state)
2119
2120Calculates the fill parameter for the conical fountain fill.
2121
2122=cut
2123*/
2124static double
2125conical_fount_f(double x, double y, struct fount_state *state) {
2126 double angle = atan2(y - state->ya, x - state->xa);
2127
2128 angle -= state->theta;
2129 if (angle < -PI)
2130 angle += PI * 2;
2131 else if (angle > PI)
2132 angle -= PI * 2;
2133
2134 return fabs(angle) * state->mult;
2135}
2136
2137/*
2138=item linear_interp(pos, seg)
2139
2140Calculates linear interpolation on the fill parameter. Breaks the
2141segment into 2 regions based in the I<middle> value.
2142
2143=cut
2144*/
2145static double
2146linear_interp(double pos, i_fountain_seg *seg) {
2147 if (pos < seg->middle) {
2148 double len = seg->middle - seg->start;
2149 if (len < EPSILON)
2150 return 0.0;
2151 else
2152 return (pos - seg->start) / len / 2;
2153 }
2154 else {
2155 double len = seg->end - seg->middle;
2156 if (len < EPSILON)
2157 return 1.0;
2158 else
2159 return 0.5 + (pos - seg->middle) / len / 2;
2160 }
2161}
2162
2163/*
2164=item sine_interp(pos, seg)
2165
2166Calculates sine function interpolation on the fill parameter.
2167
2168=cut
2169*/
2170static double
2171sine_interp(double pos, i_fountain_seg *seg) {
2172 /* I wonder if there's a simple way to smooth the transition for this */
2173 double work = linear_interp(pos, seg);
2174
2175 return (1-cos(work * PI))/2;
2176}
2177
2178/*
2179=item sphereup_interp(pos, seg)
2180
2181Calculates spherical interpolation on the fill parameter, with the cusp
2182at the low-end.
2183
2184=cut
2185*/
2186static double
2187sphereup_interp(double pos, i_fountain_seg *seg) {
2188 double work = linear_interp(pos, seg);
2189
2190 return sqrt(1.0 - (1-work) * (1-work));
2191}
2192
2193/*
2194=item spheredown_interp(pos, seg)
2195
2196Calculates spherical interpolation on the fill parameter, with the cusp
2197at the high-end.
2198
2199=cut
2200*/
2201static double
2202spheredown_interp(double pos, i_fountain_seg *seg) {
2203 double work = linear_interp(pos, seg);
2204
2205 return 1-sqrt(1.0 - work * work);
2206}
2207
2208/*
2209=item direct_cinterp(out, pos, seg)
2210
2211Calculates the fountain color based on direct scaling of the channels
2212of the color channels.
2213
2214=cut
2215*/
2216static void
2217direct_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
2218 int ch;
2219 for (ch = 0; ch < MAXCHANNELS; ++ch) {
2220 out->channel[ch] = seg->c[0].channel[ch] * (1 - pos)
2221 + seg->c[1].channel[ch] * pos;
2222 }
2223}
2224
2225/*
2226=item hue_up_cinterp(put, pos, seg)
2227
2228Calculates the fountain color based on scaling a HSV value. The hue
2229increases as the fill parameter increases.
2230
2231=cut
2232*/
2233static void
2234hue_up_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
2235 int ch;
2236 for (ch = 0; ch < MAXCHANNELS; ++ch) {
2237 out->channel[ch] = seg->c[0].channel[ch] * (1 - pos)
2238 + seg->c[1].channel[ch] * pos;
2239 }
2240 i_hsv_to_rgbf(out);
2241}
2242
2243/*
2244=item hue_down_cinterp(put, pos, seg)
2245
2246Calculates the fountain color based on scaling a HSV value. The hue
2247decreases as the fill parameter increases.
2248
2249=cut
2250*/
2251static void
2252hue_down_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
2253 int ch;
2254 for (ch = 0; ch < MAXCHANNELS; ++ch) {
2255 out->channel[ch] = seg->c[0].channel[ch] * (1 - pos)
2256 + seg->c[1].channel[ch] * pos;
2257 }
2258 i_hsv_to_rgbf(out);
2259}
2260
2261/*
2262=item simple_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
2263
2264Simple grid-based super-sampling.
2265
2266=cut
2267*/
2268static int
f1ac5027 2269simple_ssample(i_fcolor *out, double x, double y, struct fount_state *state) {
6607600c 2270 i_fcolor *work = state->ssample_data;
8d14daab 2271 i_img_dim dx, dy;
f1ac5027 2272 int grid = state->parm;
6607600c
TC
2273 double base = -0.5 + 0.5 / grid;
2274 double step = 1.0 / grid;
2275 int ch, i;
2276 int samp_count = 0;
2277
2278 for (dx = 0; dx < grid; ++dx) {
2279 for (dy = 0; dy < grid; ++dy) {
2280 if (fount_getat(work+samp_count, x + base + step * dx,
f1ac5027 2281 y + base + step * dy, state)) {
6607600c
TC
2282 ++samp_count;
2283 }
2284 }
2285 }
2286 for (ch = 0; ch < MAXCHANNELS; ++ch) {
2287 out->channel[ch] = 0;
2288 for (i = 0; i < samp_count; ++i) {
2289 out->channel[ch] += work[i].channel[ch];
2290 }
2291 /* we divide by 4 rather than samp_count since if there's only one valid
2292 sample it should be mostly transparent */
2293 out->channel[ch] /= grid * grid;
2294 }
2295 return samp_count;
2296}
2297
2298/*
2299=item random_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
2300
2301Random super-sampling.
2302
2303=cut
2304*/
2305static int
f1ac5027
TC
2306random_ssample(i_fcolor *out, double x, double y,
2307 struct fount_state *state) {
6607600c
TC
2308 i_fcolor *work = state->ssample_data;
2309 int i, ch;
f1ac5027 2310 int maxsamples = state->parm;
6607600c
TC
2311 double rand_scale = 1.0 / RAND_MAX;
2312 int samp_count = 0;
2313 for (i = 0; i < maxsamples; ++i) {
2314 if (fount_getat(work+samp_count, x - 0.5 + rand() * rand_scale,
f1ac5027 2315 y - 0.5 + rand() * rand_scale, state)) {
6607600c
TC
2316 ++samp_count;
2317 }
2318 }
2319 for (ch = 0; ch < MAXCHANNELS; ++ch) {
2320 out->channel[ch] = 0;
2321 for (i = 0; i < samp_count; ++i) {
2322 out->channel[ch] += work[i].channel[ch];
2323 }
2324 /* we divide by maxsamples rather than samp_count since if there's
2325 only one valid sample it should be mostly transparent */
2326 out->channel[ch] /= maxsamples;
2327 }
2328 return samp_count;
2329}
2330
2331/*
2332=item circle_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
2333
2334Super-sampling around the circumference of a circle.
2335
2336I considered saving the sin()/cos() values and transforming step-size
2337around the circle, but that's inaccurate, though it may not matter
2338much.
2339
2340=cut
2341 */
2342static int
f1ac5027
TC
2343circle_ssample(i_fcolor *out, double x, double y,
2344 struct fount_state *state) {
6607600c
TC
2345 i_fcolor *work = state->ssample_data;
2346 int i, ch;
f1ac5027 2347 int maxsamples = state->parm;
6607600c
TC
2348 double angle = 2 * PI / maxsamples;
2349 double radius = 0.3; /* semi-random */
2350 int samp_count = 0;
2351 for (i = 0; i < maxsamples; ++i) {
2352 if (fount_getat(work+samp_count, x + radius * cos(angle * i),
f1ac5027 2353 y + radius * sin(angle * i), state)) {
6607600c
TC
2354 ++samp_count;
2355 }
2356 }
2357 for (ch = 0; ch < MAXCHANNELS; ++ch) {
2358 out->channel[ch] = 0;
2359 for (i = 0; i < samp_count; ++i) {
2360 out->channel[ch] += work[i].channel[ch];
2361 }
2362 /* we divide by maxsamples rather than samp_count since if there's
2363 only one valid sample it should be mostly transparent */
2364 out->channel[ch] /= maxsamples;
2365 }
2366 return samp_count;
2367}
2368
2369/*
2370=item fount_r_none(v)
2371
2372Implements no repeats. Simply clamps the fill value.
2373
2374=cut
2375*/
2376static double
2377fount_r_none(double v) {
2378 return v < 0 ? 0 : v > 1 ? 1 : v;
2379}
2380
2381/*
2382=item fount_r_sawtooth(v)
2383
2384Implements sawtooth repeats. Clamps negative values and uses fmod()
2385on others.
2386
2387=cut
2388*/
2389static double
2390fount_r_sawtooth(double v) {
2391 return v < 0 ? 0 : fmod(v, 1.0);
2392}
2393
2394/*
2395=item fount_r_triangle(v)
2396
2397Implements triangle repeats. Clamps negative values, uses fmod to get
2398a range 0 through 2 and then adjusts values > 1.
2399
2400=cut
2401*/
2402static double
2403fount_r_triangle(double v) {
2404 if (v < 0)
2405 return 0;
2406 else {
2407 v = fmod(v, 2.0);
2408 return v > 1.0 ? 2.0 - v : v;
2409 }
2410}
2411
2412/*
2413=item fount_r_saw_both(v)
2414
2415Implements sawtooth repeats in the both postive and negative directions.
2416
2417Adjusts the value to be postive and then just uses fmod().
2418
2419=cut
2420*/
2421static double
2422fount_r_saw_both(double v) {
2423 if (v < 0)
2424 v += 1+(int)(-v);
2425 return fmod(v, 1.0);
2426}
2427
2428/*
2429=item fount_r_tri_both(v)
2430
2431Implements triangle repeats in the both postive and negative directions.
2432
2433Uses fmod on the absolute value, and then adjusts values > 1.
2434
2435=cut
2436*/
2437static double
2438fount_r_tri_both(double v) {
2439 v = fmod(fabs(v), 2.0);
2440 return v > 1.0 ? 2.0 - v : v;
2441}
2442
f1ac5027
TC
2443/*
2444=item fill_fountf(fill, x, y, width, channels, data)
2445
2446The fill function for fountain fills.
2447
2448=cut
2449*/
2450static void
50c75381
TC
2451fill_fountf(i_fill_t *fill, i_img_dim x, i_img_dim y, i_img_dim width,
2452 int channels, i_fcolor *data) {
f1ac5027 2453 i_fill_fountain_t *f = (i_fill_fountain_t *)fill;
efdc2568 2454
43c5dacb
TC
2455 while (width--) {
2456 i_fcolor c;
2457 int got_one;
2458
2459 if (f->state.ssfunc)
2460 got_one = f->state.ssfunc(&c, x, y, &f->state);
2461 else
2462 got_one = fount_getat(&c, x, y, &f->state);
98747309
TC
2463
2464 if (got_one)
2465 *data++ = c;
43c5dacb
TC
2466
2467 ++x;
f1ac5027
TC
2468 }
2469}
2470
2471/*
2472=item fount_fill_destroy(fill)
2473
2474=cut
2475*/
2476static void
2477fount_fill_destroy(i_fill_t *fill) {
2478 i_fill_fountain_t *f = (i_fill_fountain_t *)fill;
2479 fount_finish_state(&f->state);
2480}
2481
6607600c
TC
2482/*
2483=back
2484
2485=head1 AUTHOR
2486
2487Arnar M. Hrafnkelsson <addi@umich.edu>
2488
2489Tony Cook <tony@develop-help.com> (i_fountain())
2490
2491=head1 SEE ALSO
2492
2493Imager(3)
2494
2495=cut
2496*/