avoid an unneeded check in the FT1 has_chars() method implementation
[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
TC
1770 i_get_combine(combine, &combine_func, &combinef_func);
1771 if (combinef_func)
e310e5f9 1772 work = mymalloc(line_bytes); /* checked 17feb2005 tonyc */
f1ac5027
TC
1773
1774 fount_init_state(&state, xa, ya, xb, yb, type, repeat, combine,
e4330a7b 1775 super_sample, ssample_param, count, segs);
f1ac5027
TC
1776
1777 for (y = 0; y < im->ysize; ++y) {
1778 i_glinf(im, 0, im->xsize, y, line);
1779 for (x = 0; x < im->xsize; ++x) {
1780 i_fcolor c;
1781 int got_one;
f1ac5027
TC
1782 if (super_sample == i_fts_none)
1783 got_one = fount_getat(&c, x, y, &state);
1784 else
1785 got_one = state.ssfunc(&c, x, y, &state);
1786 if (got_one) {
efdc2568
TC
1787 if (combine)
1788 work[x] = c;
f1ac5027
TC
1789 else
1790 line[x] = c;
1791 }
1792 }
efdc2568
TC
1793 if (combine)
1794 combinef_func(line, work, im->channels, im->xsize);
f1ac5027
TC
1795 i_plinf(im, 0, im->xsize, y, line);
1796 }
1797 fount_finish_state(&state);
a73aeb5f 1798 if (work) myfree(work);
f1ac5027 1799 myfree(line);
e310e5f9
TC
1800
1801 return 1;
f1ac5027
TC
1802}
1803
1804typedef struct {
1805 i_fill_t base;
1806 struct fount_state state;
1807} i_fill_fountain_t;
1808
1809static void
50c75381 1810fill_fountf(i_fill_t *fill, i_img_dim x, i_img_dim y, i_img_dim width, int channels,
43c5dacb 1811 i_fcolor *data);
f1ac5027
TC
1812static void
1813fount_fill_destroy(i_fill_t *fill);
1814
9b1ec2b8
TC
1815static i_fill_fountain_t
1816fount_fill_proto =
1817 {
1818 {
1819 NULL,
1820 fill_fountf,
1821 fount_fill_destroy
1822 }
1823 };
1824
1825
f1ac5027 1826/*
5715f7c3 1827=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
1828
1829=category Fills
1830=synopsis fill = i_new_fill_fount(0, 0, 100, 100, i_ft_linear, i_ft_linear,
1831=synopsis i_fr_triangle, 0, i_fts_grid, 9, 1, segs);
1832
f1ac5027 1833
773bc121
TC
1834Creates a new general fill which fills with a fountain fill.
1835
f1ac5027
TC
1836=cut
1837*/
1838
1839i_fill_t *
1840i_new_fill_fount(double xa, double ya, double xb, double yb,
1841 i_fountain_type type, i_fountain_repeat repeat,
1842 int combine, int super_sample, double ssample_param,
1843 int count, i_fountain_seg *segs) {
1844 i_fill_fountain_t *fill = mymalloc(sizeof(i_fill_fountain_t));
1845
9b1ec2b8 1846 *fill = fount_fill_proto;
efdc2568
TC
1847 if (combine)
1848 i_get_combine(combine, &fill->base.combine, &fill->base.combinef);
1849 else {
1850 fill->base.combine = NULL;
1851 fill->base.combinef = NULL;
1852 }
f1ac5027
TC
1853 fount_init_state(&fill->state, xa, ya, xb, yb, type, repeat, combine,
1854 super_sample, ssample_param, count, segs);
1855
1856 return &fill->base;
1857}
1858
1859/*
1860=back
1861
1862=head1 INTERNAL FUNCTIONS
1863
1864=over
1865
1866=item fount_init_state(...)
1867
1868Used by both the fountain fill filter and the fountain fill.
1869
1870=cut
1871*/
1872
1873static void
1874fount_init_state(struct fount_state *state, double xa, double ya,
1875 double xb, double yb, i_fountain_type type,
1876 i_fountain_repeat repeat, int combine, int super_sample,
1877 double ssample_param, int count, i_fountain_seg *segs) {
6607600c 1878 int i, j;
8d14daab 1879 size_t bytes;
d3d01aeb 1880 i_fountain_seg *my_segs = mymalloc(sizeof(i_fountain_seg) * count); /* checked 2jul06 - duplicating original */
f1ac5027 1881 /*int have_alpha = im->channels == 2 || im->channels == 4;*/
f1ac5027
TC
1882
1883 memset(state, 0, sizeof(*state));
6607600c
TC
1884 /* we keep a local copy that we can adjust for speed */
1885 for (i = 0; i < count; ++i) {
1886 i_fountain_seg *seg = my_segs + i;
1887
1888 *seg = segs[i];
4c033fd4
TC
1889 if (seg->type < 0 || seg->type >= i_fst_end)
1890 seg->type = i_fst_linear;
6607600c
TC
1891 if (seg->color < 0 || seg->color >= i_fc_end)
1892 seg->color = i_fc_direct;
1893 if (seg->color == i_fc_hue_up || seg->color == i_fc_hue_down) {
1894 /* so we don't have to translate to HSV on each request, do it here */
1895 for (j = 0; j < 2; ++j) {
1896 i_rgb_to_hsvf(seg->c+j);
1897 }
1898 if (seg->color == i_fc_hue_up) {
1899 if (seg->c[1].channel[0] <= seg->c[0].channel[0])
1900 seg->c[1].channel[0] += 1.0;
1901 }
1902 else {
1903 if (seg->c[0].channel[0] <= seg->c[0].channel[1])
1904 seg->c[0].channel[0] += 1.0;
1905 }
1906 }
1907 /*printf("start %g mid %g end %g c0(%g,%g,%g,%g) c1(%g,%g,%g,%g) type %d color %d\n",
1908 seg->start, seg->middle, seg->end, seg->c[0].channel[0],
1909 seg->c[0].channel[1], seg->c[0].channel[2], seg->c[0].channel[3],
1910 seg->c[1].channel[0], seg->c[1].channel[1], seg->c[1].channel[2],
1911 seg->c[1].channel[3], seg->type, seg->color);*/
1912
1913 }
1914
1915 /* initialize each engine */
1916 /* these are so common ... */
f1ac5027
TC
1917 state->lA = xb - xa;
1918 state->lB = yb - ya;
1919 state->AB = sqrt(state->lA * state->lA + state->lB * state->lB);
1920 state->xa = xa;
1921 state->ya = ya;
6607600c
TC
1922 switch (type) {
1923 default:
1924 type = i_ft_linear; /* make the invalid value valid */
1925 case i_ft_linear:
1926 case i_ft_bilinear:
f1ac5027
TC
1927 state->lC = ya * ya - ya * yb + xa * xa - xa * xb;
1928 state->mult = 1;
1929 state->mult = 1/linear_fount_f(xb, yb, state);
6607600c
TC
1930 break;
1931
1932 case i_ft_radial:
f1ac5027
TC
1933 state->mult = 1.0 / sqrt((double)(xb-xa)*(xb-xa)
1934 + (double)(yb-ya)*(yb-ya));
6607600c
TC
1935 break;
1936
1937 case i_ft_radial_square:
f1ac5027
TC
1938 state->cos = state->lA / state->AB;
1939 state->sin = state->lB / state->AB;
1940 state->mult = 1.0 / state->AB;
6607600c
TC
1941 break;
1942
1943 case i_ft_revolution:
f1ac5027
TC
1944 state->theta = atan2(yb-ya, xb-xa);
1945 state->mult = 1.0 / (PI * 2);
6607600c
TC
1946 break;
1947
1948 case i_ft_conical:
f1ac5027
TC
1949 state->theta = atan2(yb-ya, xb-xa);
1950 state->mult = 1.0 / PI;
6607600c
TC
1951 break;
1952 }
f1ac5027 1953 state->ffunc = fount_funcs[type];
6607600c 1954 if (super_sample < 0
290bdf77 1955 || super_sample >= (int)(sizeof(fount_ssamples)/sizeof(*fount_ssamples))) {
6607600c
TC
1956 super_sample = 0;
1957 }
f1ac5027 1958 state->ssample_data = NULL;
6607600c
TC
1959 switch (super_sample) {
1960 case i_fts_grid:
1961 ssample_param = floor(0.5 + sqrt(ssample_param));
d3d01aeb
TC
1962 bytes = ssample_param * ssample_param * sizeof(i_fcolor);
1963 if (bytes / sizeof(i_fcolor) == ssample_param * ssample_param) {
1964 state->ssample_data = mymalloc(sizeof(i_fcolor) * ssample_param * ssample_param); /* checked 1jul06 tonyc */
1965 }
1966 else {
1967 super_sample = i_fts_none;
1968 }
6607600c
TC
1969 break;
1970
1971 case i_fts_random:
1972 case i_fts_circle:
1973 ssample_param = floor(0.5+ssample_param);
d3d01aeb
TC
1974 bytes = sizeof(i_fcolor) * ssample_param;
1975 if (bytes / sizeof(i_fcolor) == ssample_param) {
1976 state->ssample_data = mymalloc(sizeof(i_fcolor) * ssample_param);
1977 }
1978 else {
1979 super_sample = i_fts_none;
1980 }
6607600c
TC
1981 break;
1982 }
f1ac5027
TC
1983 state->parm = ssample_param;
1984 state->ssfunc = fount_ssamples[super_sample];
6607600c
TC
1985 if (repeat < 0 || repeat >= (sizeof(fount_repeats)/sizeof(*fount_repeats)))
1986 repeat = 0;
f1ac5027
TC
1987 state->rpfunc = fount_repeats[repeat];
1988 state->segs = my_segs;
1989 state->count = count;
6607600c
TC
1990}
1991
f1ac5027
TC
1992static void
1993fount_finish_state(struct fount_state *state) {
1994 if (state->ssample_data)
1995 myfree(state->ssample_data);
1996 myfree(state->segs);
1997}
6607600c 1998
6607600c 1999
f1ac5027 2000/*
6607600c
TC
2001=item fount_getat(out, x, y, ffunc, rpfunc, state, segs, count)
2002
2003Evaluates the fountain fill at the given point.
2004
2005This is called by both the non-super-sampling and super-sampling code.
2006
2007You might think that it would make sense to sample the fill parameter
2008instead, and combine those, but this breaks badly.
2009
2010=cut
2011*/
2012
2013static int
f1ac5027
TC
2014fount_getat(i_fcolor *out, double x, double y, struct fount_state *state) {
2015 double v = (state->rpfunc)((state->ffunc)(x, y, state));
6607600c
TC
2016 int i;
2017
2018 i = 0;
f1ac5027
TC
2019 while (i < state->count
2020 && (v < state->segs[i].start || v > state->segs[i].end)) {
6607600c
TC
2021 ++i;
2022 }
f1ac5027
TC
2023 if (i < state->count) {
2024 v = (fount_interps[state->segs[i].type])(v, state->segs+i);
2025 (fount_cinterps[state->segs[i].color])(out, v, state->segs+i);
6607600c
TC
2026 return 1;
2027 }
2028 else
2029 return 0;
2030}
2031
2032/*
2033=item linear_fount_f(x, y, state)
2034
2035Calculate the fill parameter for a linear fountain fill.
2036
2037Uses the point to line distance function, with some precalculation
2038done in i_fountain().
2039
2040=cut
2041*/
2042static double
2043linear_fount_f(double x, double y, struct fount_state *state) {
2044 return (state->lA * x + state->lB * y + state->lC) / state->AB * state->mult;
2045}
2046
2047/*
2048=item bilinear_fount_f(x, y, state)
2049
2050Calculate the fill parameter for a bi-linear fountain fill.
2051
2052=cut
2053*/
2054static double
2055bilinear_fount_f(double x, double y, struct fount_state *state) {
2056 return fabs((state->lA * x + state->lB * y + state->lC) / state->AB * state->mult);
2057}
2058
2059/*
2060=item radial_fount_f(x, y, state)
2061
2062Calculate the fill parameter for a radial fountain fill.
2063
2064Simply uses the distance function.
2065
2066=cut
2067 */
2068static double
2069radial_fount_f(double x, double y, struct fount_state *state) {
2070 return sqrt((double)(state->xa-x)*(state->xa-x)
2071 + (double)(state->ya-y)*(state->ya-y)) * state->mult;
2072}
2073
2074/*
2075=item square_fount_f(x, y, state)
2076
2077Calculate the fill parameter for a square fountain fill.
2078
2079Works by rotating the reference co-ordinate around the centre of the
2080square.
2081
2082=cut
2083*/
2084static double
2085square_fount_f(double x, double y, struct fount_state *state) {
8d14daab 2086 i_img_dim xc, yc; /* centred on A */
6607600c
TC
2087 double xt, yt; /* rotated by theta */
2088 xc = x - state->xa;
2089 yc = y - state->ya;
2090 xt = fabs(xc * state->cos + yc * state->sin);
2091 yt = fabs(-xc * state->sin + yc * state->cos);
2092 return (xt > yt ? xt : yt) * state->mult;
2093}
2094
2095/*
2096=item revolution_fount_f(x, y, state)
2097
2098Calculates the fill parameter for the revolution fountain fill.
2099
2100=cut
2101*/
2102static double
2103revolution_fount_f(double x, double y, struct fount_state *state) {
2104 double angle = atan2(y - state->ya, x - state->xa);
2105
2106 angle -= state->theta;
2107 if (angle < 0) {
2108 angle = fmod(angle+ PI * 4, PI*2);
2109 }
2110
2111 return angle * state->mult;
2112}
2113
2114/*
2115=item conical_fount_f(x, y, state)
2116
2117Calculates the fill parameter for the conical fountain fill.
2118
2119=cut
2120*/
2121static double
2122conical_fount_f(double x, double y, struct fount_state *state) {
2123 double angle = atan2(y - state->ya, x - state->xa);
2124
2125 angle -= state->theta;
2126 if (angle < -PI)
2127 angle += PI * 2;
2128 else if (angle > PI)
2129 angle -= PI * 2;
2130
2131 return fabs(angle) * state->mult;
2132}
2133
2134/*
2135=item linear_interp(pos, seg)
2136
2137Calculates linear interpolation on the fill parameter. Breaks the
2138segment into 2 regions based in the I<middle> value.
2139
2140=cut
2141*/
2142static double
2143linear_interp(double pos, i_fountain_seg *seg) {
2144 if (pos < seg->middle) {
2145 double len = seg->middle - seg->start;
2146 if (len < EPSILON)
2147 return 0.0;
2148 else
2149 return (pos - seg->start) / len / 2;
2150 }
2151 else {
2152 double len = seg->end - seg->middle;
2153 if (len < EPSILON)
2154 return 1.0;
2155 else
2156 return 0.5 + (pos - seg->middle) / len / 2;
2157 }
2158}
2159
2160/*
2161=item sine_interp(pos, seg)
2162
2163Calculates sine function interpolation on the fill parameter.
2164
2165=cut
2166*/
2167static double
2168sine_interp(double pos, i_fountain_seg *seg) {
2169 /* I wonder if there's a simple way to smooth the transition for this */
2170 double work = linear_interp(pos, seg);
2171
2172 return (1-cos(work * PI))/2;
2173}
2174
2175/*
2176=item sphereup_interp(pos, seg)
2177
2178Calculates spherical interpolation on the fill parameter, with the cusp
2179at the low-end.
2180
2181=cut
2182*/
2183static double
2184sphereup_interp(double pos, i_fountain_seg *seg) {
2185 double work = linear_interp(pos, seg);
2186
2187 return sqrt(1.0 - (1-work) * (1-work));
2188}
2189
2190/*
2191=item spheredown_interp(pos, seg)
2192
2193Calculates spherical interpolation on the fill parameter, with the cusp
2194at the high-end.
2195
2196=cut
2197*/
2198static double
2199spheredown_interp(double pos, i_fountain_seg *seg) {
2200 double work = linear_interp(pos, seg);
2201
2202 return 1-sqrt(1.0 - work * work);
2203}
2204
2205/*
2206=item direct_cinterp(out, pos, seg)
2207
2208Calculates the fountain color based on direct scaling of the channels
2209of the color channels.
2210
2211=cut
2212*/
2213static void
2214direct_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
2215 int ch;
2216 for (ch = 0; ch < MAXCHANNELS; ++ch) {
2217 out->channel[ch] = seg->c[0].channel[ch] * (1 - pos)
2218 + seg->c[1].channel[ch] * pos;
2219 }
2220}
2221
2222/*
2223=item hue_up_cinterp(put, pos, seg)
2224
2225Calculates the fountain color based on scaling a HSV value. The hue
2226increases as the fill parameter increases.
2227
2228=cut
2229*/
2230static void
2231hue_up_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
2232 int ch;
2233 for (ch = 0; ch < MAXCHANNELS; ++ch) {
2234 out->channel[ch] = seg->c[0].channel[ch] * (1 - pos)
2235 + seg->c[1].channel[ch] * pos;
2236 }
2237 i_hsv_to_rgbf(out);
2238}
2239
2240/*
2241=item hue_down_cinterp(put, pos, seg)
2242
2243Calculates the fountain color based on scaling a HSV value. The hue
2244decreases as the fill parameter increases.
2245
2246=cut
2247*/
2248static void
2249hue_down_cinterp(i_fcolor *out, double pos, i_fountain_seg *seg) {
2250 int ch;
2251 for (ch = 0; ch < MAXCHANNELS; ++ch) {
2252 out->channel[ch] = seg->c[0].channel[ch] * (1 - pos)
2253 + seg->c[1].channel[ch] * pos;
2254 }
2255 i_hsv_to_rgbf(out);
2256}
2257
2258/*
2259=item simple_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
2260
2261Simple grid-based super-sampling.
2262
2263=cut
2264*/
2265static int
f1ac5027 2266simple_ssample(i_fcolor *out, double x, double y, struct fount_state *state) {
6607600c 2267 i_fcolor *work = state->ssample_data;
8d14daab 2268 i_img_dim dx, dy;
f1ac5027 2269 int grid = state->parm;
6607600c
TC
2270 double base = -0.5 + 0.5 / grid;
2271 double step = 1.0 / grid;
2272 int ch, i;
2273 int samp_count = 0;
2274
2275 for (dx = 0; dx < grid; ++dx) {
2276 for (dy = 0; dy < grid; ++dy) {
2277 if (fount_getat(work+samp_count, x + base + step * dx,
f1ac5027 2278 y + base + step * dy, state)) {
6607600c
TC
2279 ++samp_count;
2280 }
2281 }
2282 }
2283 for (ch = 0; ch < MAXCHANNELS; ++ch) {
2284 out->channel[ch] = 0;
2285 for (i = 0; i < samp_count; ++i) {
2286 out->channel[ch] += work[i].channel[ch];
2287 }
2288 /* we divide by 4 rather than samp_count since if there's only one valid
2289 sample it should be mostly transparent */
2290 out->channel[ch] /= grid * grid;
2291 }
2292 return samp_count;
2293}
2294
2295/*
2296=item random_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
2297
2298Random super-sampling.
2299
2300=cut
2301*/
2302static int
f1ac5027
TC
2303random_ssample(i_fcolor *out, double x, double y,
2304 struct fount_state *state) {
6607600c
TC
2305 i_fcolor *work = state->ssample_data;
2306 int i, ch;
f1ac5027 2307 int maxsamples = state->parm;
6607600c
TC
2308 double rand_scale = 1.0 / RAND_MAX;
2309 int samp_count = 0;
2310 for (i = 0; i < maxsamples; ++i) {
2311 if (fount_getat(work+samp_count, x - 0.5 + rand() * rand_scale,
f1ac5027 2312 y - 0.5 + rand() * rand_scale, state)) {
6607600c
TC
2313 ++samp_count;
2314 }
2315 }
2316 for (ch = 0; ch < MAXCHANNELS; ++ch) {
2317 out->channel[ch] = 0;
2318 for (i = 0; i < samp_count; ++i) {
2319 out->channel[ch] += work[i].channel[ch];
2320 }
2321 /* we divide by maxsamples rather than samp_count since if there's
2322 only one valid sample it should be mostly transparent */
2323 out->channel[ch] /= maxsamples;
2324 }
2325 return samp_count;
2326}
2327
2328/*
2329=item circle_ssample(out, parm, x, y, state, ffunc, rpfunc, segs, count)
2330
2331Super-sampling around the circumference of a circle.
2332
2333I considered saving the sin()/cos() values and transforming step-size
2334around the circle, but that's inaccurate, though it may not matter
2335much.
2336
2337=cut
2338 */
2339static int
f1ac5027
TC
2340circle_ssample(i_fcolor *out, double x, double y,
2341 struct fount_state *state) {
6607600c
TC
2342 i_fcolor *work = state->ssample_data;
2343 int i, ch;
f1ac5027 2344 int maxsamples = state->parm;
6607600c
TC
2345 double angle = 2 * PI / maxsamples;
2346 double radius = 0.3; /* semi-random */
2347 int samp_count = 0;
2348 for (i = 0; i < maxsamples; ++i) {
2349 if (fount_getat(work+samp_count, x + radius * cos(angle * i),
f1ac5027 2350 y + radius * sin(angle * i), state)) {
6607600c
TC
2351 ++samp_count;
2352 }
2353 }
2354 for (ch = 0; ch < MAXCHANNELS; ++ch) {
2355 out->channel[ch] = 0;
2356 for (i = 0; i < samp_count; ++i) {
2357 out->channel[ch] += work[i].channel[ch];
2358 }
2359 /* we divide by maxsamples rather than samp_count since if there's
2360 only one valid sample it should be mostly transparent */
2361 out->channel[ch] /= maxsamples;
2362 }
2363 return samp_count;
2364}
2365
2366/*
2367=item fount_r_none(v)
2368
2369Implements no repeats. Simply clamps the fill value.
2370
2371=cut
2372*/
2373static double
2374fount_r_none(double v) {
2375 return v < 0 ? 0 : v > 1 ? 1 : v;
2376}
2377
2378/*
2379=item fount_r_sawtooth(v)
2380
2381Implements sawtooth repeats. Clamps negative values and uses fmod()
2382on others.
2383
2384=cut
2385*/
2386static double
2387fount_r_sawtooth(double v) {
2388 return v < 0 ? 0 : fmod(v, 1.0);
2389}
2390
2391/*
2392=item fount_r_triangle(v)
2393
2394Implements triangle repeats. Clamps negative values, uses fmod to get
2395a range 0 through 2 and then adjusts values > 1.
2396
2397=cut
2398*/
2399static double
2400fount_r_triangle(double v) {
2401 if (v < 0)
2402 return 0;
2403 else {
2404 v = fmod(v, 2.0);
2405 return v > 1.0 ? 2.0 - v : v;
2406 }
2407}
2408
2409/*
2410=item fount_r_saw_both(v)
2411
2412Implements sawtooth repeats in the both postive and negative directions.
2413
2414Adjusts the value to be postive and then just uses fmod().
2415
2416=cut
2417*/
2418static double
2419fount_r_saw_both(double v) {
2420 if (v < 0)
2421 v += 1+(int)(-v);
2422 return fmod(v, 1.0);
2423}
2424
2425/*
2426=item fount_r_tri_both(v)
2427
2428Implements triangle repeats in the both postive and negative directions.
2429
2430Uses fmod on the absolute value, and then adjusts values > 1.
2431
2432=cut
2433*/
2434static double
2435fount_r_tri_both(double v) {
2436 v = fmod(fabs(v), 2.0);
2437 return v > 1.0 ? 2.0 - v : v;
2438}
2439
f1ac5027
TC
2440/*
2441=item fill_fountf(fill, x, y, width, channels, data)
2442
2443The fill function for fountain fills.
2444
2445=cut
2446*/
2447static void
50c75381
TC
2448fill_fountf(i_fill_t *fill, i_img_dim x, i_img_dim y, i_img_dim width,
2449 int channels, i_fcolor *data) {
f1ac5027 2450 i_fill_fountain_t *f = (i_fill_fountain_t *)fill;
efdc2568 2451
43c5dacb
TC
2452 while (width--) {
2453 i_fcolor c;
2454 int got_one;
2455
2456 if (f->state.ssfunc)
2457 got_one = f->state.ssfunc(&c, x, y, &f->state);
2458 else
2459 got_one = fount_getat(&c, x, y, &f->state);
98747309
TC
2460
2461 if (got_one)
2462 *data++ = c;
43c5dacb
TC
2463
2464 ++x;
f1ac5027
TC
2465 }
2466}
2467
2468/*
2469=item fount_fill_destroy(fill)
2470
2471=cut
2472*/
2473static void
2474fount_fill_destroy(i_fill_t *fill) {
2475 i_fill_fountain_t *f = (i_fill_fountain_t *)fill;
2476 fount_finish_state(&f->state);
2477}
2478
6607600c
TC
2479/*
2480=back
2481
2482=head1 AUTHOR
2483
2484Arnar M. Hrafnkelsson <addi@umich.edu>
2485
2486Tony Cook <tony@develop-help.com> (i_fountain())
2487
2488=head1 SEE ALSO
2489
2490Imager(3)
2491
2492=cut
2493*/