1 /*
2 * FIG : Facility for Interactive Generation of figures
3 * Copyright (c) 1985-1988 by Supoj Sutanthavibul
4 * Parts Copyright (c) 1989-2015 by Brian V. Smith
5 * Parts Copyright (c) 1991 by Paul King
6 * Parts Copyright (c) 2016-2020 by Thomas Loimer
7 *
8 * Copyright (c) 1994 Anthony Dekker
9 *
10 * [NEUQUANT Neural-Net quantization algorithm by Anthony Dekker, 1994.
11 * See "Kohonen neural networks for optimal colour quantization"
12 * in "Network: Computation in Neural Systems" Vol. 5 (1994) pp 351-367.
13 * for a discussion of the algorithm.]
14
15 * Any party obtaining a copy of these files is granted, free of charge, a
16 * full and unrestricted irrevocable, world-wide, paid up, royalty-free,
17 * nonexclusive right and license to deal in this software and documentation
18 * files (the "Software"), including without limitation the rights to use,
19 * copy, modify, merge, publish, distribute, sublicense and/or sell copies of
20 * the Software, and to permit persons who receive copies from any such
21 * party to do so, with the only requirement being that the above copyright
22 * and this permission notice remain intact.
23 *
24 */
25
26 /*
27 * Neural-Net quantization algorithm based on work of Anthony Dekker
28 */
29
30 #include "f_neuclrtab.h"
31
32 #include <stdio.h>
33 #include <stdlib.h>
34 #include <string.h>
35 #include <limits.h> /* INT_MAX */
36
37
38 static void initnet(void);
39 static void learn(void);
40 static void unbiasnet(void);
41 static void cpyclrtab(void);
42 static void inxbuild(void);
43 static int contest(register int b, register int g, register int r);
44 static void alterneigh(int rad, int i, register int b, register int g, register int r);
45 static void altersingle(register int alpha, register int i, register int b, register int g, register int r);
46 static int inxsearch(register int b, register int g, register int r);
47
48 #define MAXNETSIZE 256
49
50 /* our color table (global) */
51 BYTE clrtab[256][3];
52
53 static int netsize;
54
55 #ifndef DEFSMPFAC
56 # ifdef SPEED
57 # define DEFSMPFAC (240/SPEED+3)
58 # else
59 # define DEFSMPFAC 30 /* only sample every 30th pixel */
60 # endif /* SPEED */
61 #endif /* DEFSMPFAC */
62
63 int samplefac = DEFSMPFAC; /* sampling factor */
64
65 /* Samples array starts off holding spacing between adjacent
66 * samples, and ends up holding actual BGR sample values.
67 */
68 static BYTE *thesamples;
69 static int nsamples;
70 static BYTE *cursamp;
71 static long skipcount;
72
73 #define MAXSKIP (1<<24-1)
74
75 #define nskip(sp) ((long)(sp)[0]<<16|(long)(sp)[1]<<8|(long)(sp)[2])
76
77 #define setskip(sp,n) ((sp)[0]=(n)>>16,(sp)[1]=((n)>>8)&255,(sp)[2]=(n)&255)
78
79 /* npixels is the total number of pixels defined */
80
81
82 void neu_dith_colrs (register BYTE *bs, register COLR (*cs), int n);
83
84 int
neu_init(long int npixels)85 neu_init(long int npixels) /* initialize our sample array */
86
87 {
88 if (npixels < MIN_NEU_SAMPLES)
89 samplefac = 1;
90 else
91 samplefac = DEFSMPFAC;
92 return neu_init2(npixels);
93 }
94
95 int
neu_init2(long int npixels)96 neu_init2(long int npixels) /* initialize our sample array */
97
98 {
99 register int nsleft;
100 register long sv;
101 double rval, cumprob;
102 long npleft;
103
104 nsamples = npixels/samplefac;
105 if (nsamples < MIN_NEU_SAMPLES)
106 return((int)-MIN_NEU_SAMPLES/nsamples-1); /* THIS IS DIFFERENT FROM THE ORIGINAL -1 */
107 thesamples = (BYTE *)malloc(nsamples*3);
108 if (thesamples == NULL)
109 return(-1);
110 cursamp = thesamples;
111 npleft = npixels;
112 nsleft = nsamples;
113 while (nsleft) {
114 rval = drand48(); /* random distance to next sample */
115 sv = 0;
116 cumprob = 0.;
117 while ((cumprob += (1.-cumprob)*nsleft/(npleft-sv)) < rval)
118 sv++;
119 if (nsleft == nsamples)
120 skipcount = sv;
121 else {
122 setskip(cursamp, sv);
123 cursamp += 3;
124 }
125 npleft -= sv+1;
126 nsleft--;
127 }
128 setskip(cursamp, npleft); /* tag on end to skip the rest */
129 cursamp = thesamples;
130 return(0);
131 }
132
133
134 void
neu_pixel(register BYTE * col)135 neu_pixel(register BYTE *col) /* add pixel to our samples */
136 {
137 if (!skipcount--) {
138 skipcount = nskip(cursamp);
139 cursamp[0] = col[N_BLU];
140 cursamp[1] = col[N_GRN];
141 cursamp[2] = col[N_RED];
142 cursamp += 3;
143 }
144 }
145
146
147 void
neu_colrs(register COLR (* cs),register int n)148 neu_colrs(register COLR (*cs), register int n) /* add a scanline to our samples */
149 {
150 while (n > skipcount) {
151 cs += skipcount;
152 n -= skipcount+1;
153 skipcount = nskip(cursamp);
154 cursamp[0] = cs[0][N_BLU];
155 cursamp[1] = cs[0][N_GRN];
156 cursamp[2] = cs[0][N_RED];
157 cs++;
158 cursamp += 3;
159 }
160 skipcount -= n;
161 }
162
163
164 int
neu_clrtab(int ncolors)165 neu_clrtab(int ncolors) /* make new color table using ncolors */
166 {
167 netsize = ncolors;
168 if (netsize > MAXNETSIZE) netsize = MAXNETSIZE;
169 initnet();
170 learn();
171 unbiasnet();
172 cpyclrtab();
173 inxbuild();
174 /* we're done with our samples */
175 free((char *)thesamples);
176 /* reset dithering function */
177 neu_dith_colrs((BYTE *)NULL, (COLR *)NULL, 0);
178 /* return new color table size */
179 return(netsize);
180 }
181
182
183 int
neu_map_pixel(register BYTE * col)184 neu_map_pixel(register BYTE *col) /* get pixel for color */
185 {
186 return(inxsearch(col[N_BLU],col[N_GRN],col[N_RED]));
187 }
188
189
neu_map_colrs(register BYTE * bs,register COLR (* cs),register int n)190 void neu_map_colrs(register BYTE *bs, register COLR (*cs), register int n) /* convert a scanline to color index values */
191 {
192 while (n-- > 0) {
193 *bs++ = inxsearch(cs[0][N_BLU],cs[0][N_GRN],cs[0][N_RED]);
194 cs++;
195 }
196 }
197
198
neu_dith_colrs(register BYTE * bs,register COLR (* cs),int n)199 void neu_dith_colrs(register BYTE *bs, register COLR (*cs), int n) /* convert scanline to dithered index values */
200 {
201 static short (*cerr)[3] = NULL;
202 static int N = 0;
203 int err[3], errp[3];
204 register int x, i;
205
206 if (n != N) { /* get error propogation array */
207 if (N) {
208 free((char *)cerr);
209 cerr = NULL;
210 }
211 if (n)
212 cerr = (short (*)[3])malloc(3*n*sizeof(short));
213 if (cerr == NULL) {
214 N = 0;
215 neu_map_colrs(bs, cs, n);
216 return;
217 }
218 N = n;
219 memset(cerr, 0, 3*N*sizeof(short));
220 }
221 err[0] = err[1] = err[2] = 0;
222 for (x = 0; x < n; x++) {
223 for (i = 0; i < 3; i++) { /* dither value */
224 errp[i] = err[i];
225 err[i] += cerr[x][i];
226 #ifdef MAXERR
227 if (err[i] > MAXERR) err[i] = MAXERR;
228 else if (err[i] < -MAXERR) err[i] = -MAXERR;
229 #endif /* MAXERR */
230 err[i] += cs[x][i];
231 if (err[i] < 0) err[i] = 0;
232 else if (err[i] > 255) err[i] = 255;
233 }
234 bs[x] = inxsearch(err[N_BLU],err[N_GRN],err[N_RED]);
235 for (i = 0; i < 3; i++) { /* propagate error */
236 err[i] -= clrtab[bs[x]][i];
237 err[i] /= 3;
238 cerr[x][i] = err[i] + errp[i];
239 }
240 }
241 }
242
243 /* The following was adapted and modified slightly from the original */
244
245 #define bool int
246 #define false 0
247 #define true 1
248
249 /* network defs */
250 #define maxnetpos (netsize-1)
251 #define netbiasshift 4 /* bias for colour values */
252 #define ncycles 100 /* no. of learning cycles */
253
254 /* defs for freq and bias */
255 #define intbiasshift 16 /* bias for fractions */
256 #define intbias (((int) 1)<<intbiasshift)
257 #define gammashift 10 /* gamma = 1024 */
258 #define gamma (((int) 1)<<gammashift)
259 #define betashift 10
260 #define beta (intbias>>betashift) /* beta = 1/1024 */
261 #define betagamma (intbias<<(gammashift-betashift))
262
263 /* defs for decreasing radius factor */
264 #define initrad (netsize>>3) /* for 256 cols, radius starts */
265 #define radiusbiasshift 6 /* at 32.0 biased by 6 bits */
266 #define radiusbias (((int) 1)<<radiusbiasshift)
267 #define initradius (initrad*radiusbias) /* and decreases by a */
268 #define radiusdec 30 /* factor of 1/30 each cycle */
269
270 /* defs for decreasing alpha factor */
271 #define alphabiasshift 10 /* alpha starts at 1.0 */
272 #define initalpha (((int) 1)<<alphabiasshift)
273 int alphadec; /* biased by 10 bits */
274
275 /* radbias and alpharadbias used for radpower calculation */
276 #define radbiasshift 8
277 #define radbias (((int) 1)<<radbiasshift)
278 #define alpharadbshift (alphabiasshift+radbiasshift)
279 #define alpharadbias (((int) 1)<<alpharadbshift)
280
281 /* four primes near 500 - assume no image has a length so large */
282 /* that it is divisible by all four primes */
283 #define prime1 499
284 #define prime2 491
285 #define prime3 487
286 #define prime4 503
287
288 /* hardwire these */
289
290 #define thepicture thesamples
291 #define lengthcount (nsamples*3)
292 #define samplefac 1 /* this MUST be 1 */
293
294 typedef int pixel[4]; /* BGRc */
295 pixel network[MAXNETSIZE];
296
297 int netindex[256]; /* for network lookup - really 256 */
298
299 int bias [MAXNETSIZE]; /* bias and freq arrays for learning */
300 int freq [MAXNETSIZE];
301 int radpower[MAXNETSIZE>>3]; /* radpower for precomputation */
302
303
304 /* initialise network in range (0,0,0) to (255,255,255) */
305
306 static void
initnet(void)307 initnet(void)
308 {
309 register int i;
310 register int *p;
311
312 for (i=0; i<netsize; i++) {
313 p = network[i];
314 p[0] = p[1] = p[2] = (i << (netbiasshift+8))/netsize;
315 freq[i] = intbias/netsize; /* 1/netsize */
316 bias[i] = 0;
317 }
318 }
319
320
321 /* do after unbias - insertion sort of network and build netindex[0..255] */
322
323 static void
inxbuild(void)324 inxbuild(void)
325 {
326 register int i,j,smallpos,smallval;
327 register int *p,*q;
328 int previouscol,startpos;
329
330 previouscol = 0;
331 startpos = 0;
332 for (i=0; i<netsize; i++) {
333 p = network[i];
334 smallpos = i;
335 smallval = p[1]; /* index on g */
336 /* find smallest in i..netsize-1 */
337 for (j=i+1; j<netsize; j++) {
338 q = network[j];
339 if (q[1] < smallval) { /* index on g */
340 smallpos = j;
341 smallval = q[1]; /* index on g */
342 }
343 }
344 q = network[smallpos];
345 /* swap p (i) and q (smallpos) entries */
346 if (i != smallpos) {
347 j = q[0]; q[0] = p[0]; p[0] = j;
348 j = q[1]; q[1] = p[1]; p[1] = j;
349 j = q[2]; q[2] = p[2]; p[2] = j;
350 j = q[3]; q[3] = p[3]; p[3] = j;
351 }
352 /* smallval entry is now in position i */
353 if (smallval != previouscol) {
354 netindex[previouscol] = (startpos+i)>>1;
355 for (j=previouscol+1; j<smallval; j++) netindex[j] = i;
356 previouscol = smallval;
357 startpos = i;
358 }
359 }
360 netindex[previouscol] = (startpos+maxnetpos)>>1;
361 for (j=previouscol+1; j<256; j++) netindex[j] = maxnetpos; /* really 256 */
362 }
363
364 static int
inxsearch(register int b,register int g,register int r)365 inxsearch(register int b, register int g, register int r) /* accepts real BGR values after net is unbiased */
366
367 {
368 register int i,j,dist,a,bestd;
369 register int *p;
370 int best;
371
372 bestd = 1000; /* biggest possible dist is 256*3 */
373 best = -1;
374 i = netindex[g]; /* index on g */
375 j = i-1; /* start at netindex[g] and work outwards */
376
377 while ((i<netsize) || (j>=0)) {
378 if (i<netsize) {
379 p = network[i];
380 dist = p[1] - g; /* inx key */
381 if (dist >= bestd) i = netsize; /* stop iter */
382 else {
383 i++;
384 if (dist<0) dist = -dist;
385 a = p[0] - b; if (a<0) a = -a;
386 dist += a;
387 if (dist<bestd) {
388 a = p[2] - r; if (a<0) a = -a;
389 dist += a;
390 if (dist<bestd) {bestd=dist; best=p[3];}
391 }
392 }
393 }
394 if (j>=0) {
395 p = network[j];
396 dist = g - p[1]; /* inx key - reverse dif */
397 if (dist >= bestd) j = -1; /* stop iter */
398 else {
399 j--;
400 if (dist<0) dist = -dist;
401 a = p[0] - b; if (a<0) a = -a;
402 dist += a;
403 if (dist<bestd) {
404 a = p[2] - r; if (a<0) a = -a;
405 dist += a;
406 if (dist<bestd) {bestd=dist; best=p[3];}
407 }
408 }
409 }
410 }
411 return(best);
412 }
413
414
415 /* finds closest neuron (min dist) and updates freq */
416 /* finds best neuron (min dist-bias) and returns position */
417 /* for frequently chosen neurons, freq[i] is high and bias[i] is negative */
418 /* bias[i] = gamma*((1/netsize)-freq[i]) */
419
420 static int
contest(register int b,register int g,register int r)421 contest(register int b, register int g, register int r) /* accepts biased BGR values */
422
423 {
424 register int i,dist,a,biasdist,betafreq;
425 int bestpos,bestbiaspos,bestd,bestbiasd;
426 register int *p,*f, *n;
427
428 bestd = INT_MAX;
429 bestbiasd = bestd;
430 bestpos = -1;
431 bestbiaspos = bestpos;
432 p = bias;
433 f = freq;
434
435 for (i=0; i<netsize; i++) {
436 n = network[i];
437 dist = n[0] - b; if (dist<0) dist = -dist;
438 a = n[1] - g; if (a<0) a = -a;
439 dist += a;
440 a = n[2] - r; if (a<0) a = -a;
441 dist += a;
442 if (dist<bestd) {bestd=dist; bestpos=i;}
443 biasdist = dist - ((*p)>>(intbiasshift-netbiasshift));
444 if (biasdist<bestbiasd) {bestbiasd=biasdist; bestbiaspos=i;}
445 betafreq = (*f >> betashift);
446 *f++ -= betafreq;
447 *p++ += (betafreq<<gammashift);
448 }
449 freq[bestpos] += beta;
450 bias[bestpos] -= betagamma;
451 return(bestbiaspos);
452 }
453
454
455 /* move neuron i towards (b,g,r) by factor alpha */
456
457 static void
altersingle(register int alpha,register int i,register int b,register int g,register int r)458 altersingle(register int alpha, register int i, register int b, register int g, register int r) /* accepts biased BGR values */
459
460 {
461 register int *n;
462
463 n = network[i]; /* alter hit neuron */
464 *n -= (alpha*(*n - b)) / initalpha;
465 n++;
466 *n -= (alpha*(*n - g)) / initalpha;
467 n++;
468 *n -= (alpha*(*n - r)) / initalpha;
469 }
470
471
472 /* move neurons adjacent to i towards (b,g,r) by factor */
473 /* alpha*(1-((i-j)^2/[r]^2)) precomputed as radpower[|i-j|]*/
474
475 static void
alterneigh(int rad,int i,register int b,register int g,register int r)476 alterneigh(int rad, int i, register int b, register int g, register int r) /* accents biased BGR values */
477
478
479 {
480 register int j,k,lo,hi,a;
481 register int *p, *q;
482
483 lo = i-rad; if (lo<-1) lo = -1;
484 hi = i+rad; if (hi>netsize) hi=netsize;
485
486 j = i+1;
487 k = i-1;
488 q = radpower;
489 while ((j<hi) || (k>lo)) {
490 a = (*(++q));
491 if (j<hi) {
492 p = network[j];
493 *p -= (a*(*p - b)) / alpharadbias;
494 p++;
495 *p -= (a*(*p - g)) / alpharadbias;
496 p++;
497 *p -= (a*(*p - r)) / alpharadbias;
498 j++;
499 }
500 if (k>lo) {
501 p = network[k];
502 *p -= (a*(*p - b)) / alpharadbias;
503 p++;
504 *p -= (a*(*p - g)) / alpharadbias;
505 p++;
506 *p -= (a*(*p - r)) / alpharadbias;
507 k--;
508 }
509 }
510 }
511
512
513 static void
learn(void)514 learn(void)
515 {
516 register int i,j,b,g,r;
517 int radius,rad,alpha,step,delta,samplepixels;
518 register unsigned char *p;
519 unsigned char *lim;
520
521 alphadec = 30 + ((samplefac-1)/3);
522 p = thepicture;
523 lim = thepicture + lengthcount;
524 samplepixels = lengthcount/(3*samplefac);
525 delta = samplepixels/ncycles;
526 alpha = initalpha;
527 radius = initradius;
528
529 rad = radius >> radiusbiasshift;
530 if (rad <= 1) rad = 0;
531 for (i=0; i<rad; i++)
532 radpower[i] = alpha*(((rad*rad - i*i)*radbias)/(rad*rad));
533
534 if ((lengthcount%prime1) != 0) step = 3*prime1;
535 else {
536 if ((lengthcount%prime2) !=0) step = 3*prime2;
537 else {
538 if ((lengthcount%prime3) !=0) step = 3*prime3;
539 else step = 3*prime4;
540 }
541 }
542
543 i = 0;
544 while (i < samplepixels) {
545 b = p[0] << netbiasshift;
546 g = p[1] << netbiasshift;
547 r = p[2] << netbiasshift;
548 j = contest(b,g,r);
549
550 altersingle(alpha,j,b,g,r);
551 if (rad) alterneigh(rad,j,b,g,r); /* alter neighbours */
552
553 p += step;
554 if (p >= lim) p -= lengthcount;
555
556 i++;
557 if (i%delta == 0) {
558 alpha -= alpha / alphadec;
559 radius -= radius / radiusdec;
560 rad = radius >> radiusbiasshift;
561 if (rad <= 1) rad = 0;
562 for (j=0; j<rad; j++)
563 radpower[j] = alpha*(((rad*rad - j*j)*radbias)/(rad*rad));
564 }
565 }
566 }
567
568 /* unbias network to give 0..255 entries */
569 /* which can then be used for colour map */
570 /* and record position i to prepare for sort */
571
572 static void
unbiasnet(void)573 unbiasnet(void)
574 {
575 int i,j;
576
577 for (i=0; i<netsize; i++) {
578 for (j=0; j<3; j++)
579 network[i][j] >>= netbiasshift;
580 network[i][3] = i; /* record colour no */
581 }
582 }
583
584 /* Don't do this until the network has been unbiased */
585
586 static void
cpyclrtab(void)587 cpyclrtab(void)
588 {
589 register int i,j,k;
590
591 for (j=0; j<netsize; j++) {
592 k = network[j][3];
593 for (i = 0; i < 3; i++)
594 clrtab[k][i] = network[j][2-i];
595 }
596 }
597