001    package railo.runtime.img;
002    
003    public class NeuQuant {
004    
005            protected static final int netsize = 256; /* number of colours used */
006    
007            /* four primes near 500 - assume no image has a length so large */
008            /* that it is divisible by all four primes */
009            protected static final int prime1 = 499;
010            protected static final int prime2 = 491;
011            protected static final int prime3 = 487;
012            protected static final int prime4 = 503;
013    
014            protected static final int minpicturebytes = (3 * prime4);
015            /* minimum size for input image */
016    
017            /* Program Skeleton
018               ----------------
019               [select samplefac in range 1..30]
020               [read image from input file]
021               pic = (unsigned char*) malloc(3*width*height);
022               initnet(pic,3*width*height,samplefac);
023               learn();
024               unbiasnet();
025               [write output image header, using writecolourmap(f)]
026               inxbuild();
027               write output image using inxsearch(b,g,r)      */
028    
029            /* Network Definitions
030               ------------------- */
031    
032            protected static final int maxnetpos = (netsize - 1);
033            protected static final int netbiasshift = 4; /* bias for colour values */
034            protected static final int ncycles = 100; /* no. of learning cycles */
035    
036            /* defs for freq and bias */
037            protected static final int intbiasshift = 16; /* bias for fractions */
038            protected static final int intbias = (( 1) << intbiasshift);
039            protected static final int gammashift = 10; /* gamma = 1024 */
040            protected static final int gamma = (( 1) << gammashift);
041            protected static final int betashift = 10;
042            protected static final int beta = (intbias >> betashift); /* beta = 1/1024 */
043            protected static final int betagamma =
044                    (intbias << (gammashift - betashift));
045    
046            /* defs for decreasing radius factor */
047            protected static final int initrad = (netsize >> 3); /* for 256 cols, radius starts */
048            protected static final int radiusbiasshift = 6; /* at 32.0 biased by 6 bits */
049            protected static final int radiusbias = (( 1) << radiusbiasshift);
050            protected static final int initradius = (initrad * radiusbias); /* and decreases by a */
051            protected static final int radiusdec = 30; /* factor of 1/30 each cycle */
052    
053            /* defs for decreasing alpha factor */
054            protected static final int alphabiasshift = 10; /* alpha starts at 1.0 */
055            protected static final int initalpha = (1 << alphabiasshift);
056    
057            protected int alphadec; /* biased by 10 bits */
058    
059            /* radbias and alpharadbias used for radpower calculation */
060            protected static final int radbiasshift = 8;
061            protected static final int radbias = (1 << radbiasshift);
062            protected static final int alpharadbshift = (alphabiasshift + radbiasshift);
063            protected static final int alpharadbias = (( 1) << alpharadbshift);
064    
065            /* Types and Global Variables
066            -------------------------- */
067    
068            protected byte[] thepicture; /* the input image itself */
069            protected int lengthcount; /* lengthcount = H*W*3 */
070    
071            protected int samplefac; /* sampling factor 1..30 */
072    
073            //   typedef int pixel[4];                /* BGRc */
074            protected int[][] network; /* the network itself - [netsize][4] */
075    
076            protected int[] netindex = new int[256];
077            /* for network lookup - really 256 */
078    
079            protected int[] bias = new int[netsize];
080            /* bias and freq arrays for learning */
081            protected int[] freq = new int[netsize];
082            protected int[] radpower = new int[initrad];
083            /* radpower for precomputation */
084    
085            /* Initialise network in range (0,0,0) to (255,255,255) and set parameters
086               ----------------------------------------------------------------------- */
087            public NeuQuant(byte[] thepic, int len, int sample) {
088    
089                    int i;
090                    int[] p;
091    
092                    thepicture = thepic;
093                    lengthcount = len;
094                    samplefac = sample;
095    
096                    network = new int[netsize][];
097                    for (i = 0; i < netsize; i++) {
098                            network[i] = new int[4];
099                            p = network[i];
100                            p[0] = p[1] = p[2] = (i << (netbiasshift + 8)) / netsize;
101                            freq[i] = intbias / netsize; /* 1/netsize */
102                            bias[i] = 0;
103                    }
104            }
105            
106            public byte[] colorMap() {
107                    byte[] map = new byte[3 * netsize];
108                    int[] index = new int[netsize];
109                    for (int i = 0; i < netsize; i++)
110                            index[network[i][3]] = i;
111                    int k = 0;
112                    for (int i = 0; i < netsize; i++) {
113                            int j = index[i];
114                            map[k++] = (byte) (network[j][0]);
115                            map[k++] = (byte) (network[j][1]);
116                            map[k++] = (byte) (network[j][2]);
117                    }
118                    return map;
119            }
120            
121            /* Insertion sort of network and building of netindex[0..255] (to do after unbias)
122               ------------------------------------------------------------------------------- */
123            public void inxbuild() {
124    
125                    int i, j, smallpos, smallval;
126                    int[] p;
127                    int[] q;
128                    int previouscol, startpos;
129    
130                    previouscol = 0;
131                    startpos = 0;
132                    for (i = 0; i < netsize; i++) {
133                            p = network[i];
134                            smallpos = i;
135                            smallval = p[1]; /* index on g */
136                            /* find smallest in i..netsize-1 */
137                            for (j = i + 1; j < netsize; j++) {
138                                    q = network[j];
139                                    if (q[1] < smallval) { /* index on g */
140                                            smallpos = j;
141                                            smallval = q[1]; /* index on g */
142                                    }
143                            }
144                            q = network[smallpos];
145                            /* swap p (i) and q (smallpos) entries */
146                            if (i != smallpos) {
147                                    j = q[0];
148                                    q[0] = p[0];
149                                    p[0] = j;
150                                    j = q[1];
151                                    q[1] = p[1];
152                                    p[1] = j;
153                                    j = q[2];
154                                    q[2] = p[2];
155                                    p[2] = j;
156                                    j = q[3];
157                                    q[3] = p[3];
158                                    p[3] = j;
159                            }
160                            /* smallval entry is now in position i */
161                            if (smallval != previouscol) {
162                                    netindex[previouscol] = (startpos + i) >> 1;
163                                    for (j = previouscol + 1; j < smallval; j++)
164                                            netindex[j] = i;
165                                    previouscol = smallval;
166                                    startpos = i;
167                            }
168                    }
169                    netindex[previouscol] = (startpos + maxnetpos) >> 1;
170                    for (j = previouscol + 1; j < 256; j++)
171                            netindex[j] = maxnetpos; /* really 256 */
172            }
173            
174            /* Main Learning Loop
175               ------------------ */
176            public void learn() {
177    
178                    int i, j, b, g, r;
179                    int radius, rad, alpha, step, delta, samplepixels;
180                    byte[] p;
181                    int pix, lim;
182    
183                    if (lengthcount < minpicturebytes)
184                            samplefac = 1;
185                    alphadec = 30 + ((samplefac - 1) / 3);
186                    p = thepicture;
187                    pix = 0;
188                    lim = lengthcount;
189                    samplepixels = lengthcount / (3 * samplefac);
190                    delta = samplepixels / ncycles;
191                    alpha = initalpha;
192                    radius = initradius;
193    
194                    rad = radius >> radiusbiasshift;
195                    if (rad <= 1)
196                            rad = 0;
197                    for (i = 0; i < rad; i++)
198                            radpower[i] =
199                                    alpha * (((rad * rad - i * i) * radbias) / (rad * rad));
200    
201                    //fprintf(stderr,"beginning 1D learning: initial radius=%d\n", rad);
202    
203                    if (lengthcount < minpicturebytes)
204                            step = 3;
205                    else if ((lengthcount % prime1) != 0)
206                            step = 3 * prime1;
207                    else {
208                            if ((lengthcount % prime2) != 0)
209                                    step = 3 * prime2;
210                            else {
211                                    if ((lengthcount % prime3) != 0)
212                                            step = 3 * prime3;
213                                    else
214                                            step = 3 * prime4;
215                            }
216                    }
217    
218                    i = 0;
219                    while (i < samplepixels) {
220                            b = (p[pix + 0] & 0xff) << netbiasshift;
221                            g = (p[pix + 1] & 0xff) << netbiasshift;
222                            r = (p[pix + 2] & 0xff) << netbiasshift;
223                            j = contest(b, g, r);
224    
225                            altersingle(alpha, j, b, g, r);
226                            if (rad != 0)
227                                    alterneigh(rad, j, b, g, r); /* alter neighbours */
228    
229                            pix += step;
230                            if (pix >= lim)
231                                    pix -= lengthcount;
232    
233                            i++;
234                            if (delta == 0)
235                                    delta = 1;
236                            if (i % delta == 0) {
237                                    alpha -= alpha / alphadec;
238                                    radius -= radius / radiusdec;
239                                    rad = radius >> radiusbiasshift;
240                                    if (rad <= 1)
241                                            rad = 0;
242                                    for (j = 0; j < rad; j++)
243                                            radpower[j] =
244                                                    alpha * (((rad * rad - j * j) * radbias) / (rad * rad));
245                            }
246                    }
247                    //fprintf(stderr,"finished 1D learning: final alpha=%f !\n",((float)alpha)/initalpha);
248            }
249            
250            /* Search for BGR values 0..255 (after net is unbiased) and return colour index
251               ---------------------------------------------------------------------------- */
252            public int map(int b, int g, int r) {
253    
254                    int i, j, dist, a, bestd;
255                    int[] p;
256                    int best;
257    
258                    bestd = 1000; /* biggest possible dist is 256*3 */
259                    best = -1;
260                    i = netindex[g]; /* index on g */
261                    j = i - 1; /* start at netindex[g] and work outwards */
262    
263                    while ((i < netsize) || (j >= 0)) {
264                            if (i < netsize) {
265                                    p = network[i];
266                                    dist = p[1] - g; /* inx key */
267                                    if (dist >= bestd)
268                                            i = netsize; /* stop iter */
269                                    else {
270                                            i++;
271                                            if (dist < 0)
272                                                    dist = -dist;
273                                            a = p[0] - b;
274                                            if (a < 0)
275                                                    a = -a;
276                                            dist += a;
277                                            if (dist < bestd) {
278                                                    a = p[2] - r;
279                                                    if (a < 0)
280                                                            a = -a;
281                                                    dist += a;
282                                                    if (dist < bestd) {
283                                                            bestd = dist;
284                                                            best = p[3];
285                                                    }
286                                            }
287                                    }
288                            }
289                            if (j >= 0) {
290                                    p = network[j];
291                                    dist = g - p[1]; /* inx key - reverse dif */
292                                    if (dist >= bestd)
293                                            j = -1; /* stop iter */
294                                    else {
295                                            j--;
296                                            if (dist < 0)
297                                                    dist = -dist;
298                                            a = p[0] - b;
299                                            if (a < 0)
300                                                    a = -a;
301                                            dist += a;
302                                            if (dist < bestd) {
303                                                    a = p[2] - r;
304                                                    if (a < 0)
305                                                            a = -a;
306                                                    dist += a;
307                                                    if (dist < bestd) {
308                                                            bestd = dist;
309                                                            best = p[3];
310                                                    }
311                                            }
312                                    }
313                            }
314                    }
315                    return (best);
316            }
317            public byte[] process() {
318                    learn();
319                    unbiasnet();
320                    inxbuild();
321                    return colorMap();
322            }
323            
324            /* Unbias network to give byte values 0..255 and record position i to prepare for sort
325               ----------------------------------------------------------------------------------- */
326            public void unbiasnet() {
327    
328                    int i;
329    
330                    for (i = 0; i < netsize; i++) {
331                            network[i][0] >>= netbiasshift;
332                            network[i][1] >>= netbiasshift;
333                            network[i][2] >>= netbiasshift;
334                            network[i][3] = i; /* record colour no */
335                    }
336            }
337            
338            /* Move adjacent neurons by precomputed alpha*(1-((i-j)^2/[r]^2)) in radpower[|i-j|]
339               --------------------------------------------------------------------------------- */
340            protected void alterneigh(int rad, int i, int b, int g, int r) {
341    
342                    int j, k, lo, hi, a, m;
343                    int[] p;
344    
345                    lo = i - rad;
346                    if (lo < -1)
347                            lo = -1;
348                    hi = i + rad;
349                    if (hi > netsize)
350                            hi = netsize;
351    
352                    j = i + 1;
353                    k = i - 1;
354                    m = 1;
355                    while ((j < hi) || (k > lo)) {
356                            a = radpower[m++];
357                            if (j < hi) {
358                                    p = network[j++];
359                                    try {
360                                            p[0] -= (a * (p[0] - b)) / alpharadbias;
361                                            p[1] -= (a * (p[1] - g)) / alpharadbias;
362                                            p[2] -= (a * (p[2] - r)) / alpharadbias;
363                                    } catch (Exception e) {
364                                    } // prevents 1.3 miscompilation
365                            }
366                            if (k > lo) {
367                                    p = network[k--];
368                                    try {
369                                            p[0] -= (a * (p[0] - b)) / alpharadbias;
370                                            p[1] -= (a * (p[1] - g)) / alpharadbias;
371                                            p[2] -= (a * (p[2] - r)) / alpharadbias;
372                                    } catch (Exception e) {
373                                    }
374                            }
375                    }
376            }
377            
378            /* Move neuron i towards biased (b,g,r) by factor alpha
379               ---------------------------------------------------- */
380            protected void altersingle(int alpha, int i, int b, int g, int r) {
381    
382                    /* alter hit neuron */
383                    int[] n = network[i];
384                    n[0] -= (alpha * (n[0] - b)) / initalpha;
385                    n[1] -= (alpha * (n[1] - g)) / initalpha;
386                    n[2] -= (alpha * (n[2] - r)) / initalpha;
387            }
388            
389            /* Search for biased BGR values
390               ---------------------------- */
391            protected int contest(int b, int g, int r) {
392    
393                    /* finds closest neuron (min dist) and updates freq */
394                    /* finds best neuron (min dist-bias) and returns position */
395                    /* for frequently chosen neurons, freq[i] is high and bias[i] is negative */
396                    /* bias[i] = gamma*((1/netsize)-freq[i]) */
397    
398                    int i, dist, a, biasdist, betafreq;
399                    int bestpos, bestbiaspos, bestd, bestbiasd;
400                    int[] n;
401    
402                    bestd = ~(( 1) << 31);
403                    bestbiasd = bestd;
404                    bestpos = -1;
405                    bestbiaspos = bestpos;
406    
407                    for (i = 0; i < netsize; i++) {
408                            n = network[i];
409                            dist = n[0] - b;
410                            if (dist < 0)
411                                    dist = -dist;
412                            a = n[1] - g;
413                            if (a < 0)
414                                    a = -a;
415                            dist += a;
416                            a = n[2] - r;
417                            if (a < 0)
418                                    a = -a;
419                            dist += a;
420                            if (dist < bestd) {
421                                    bestd = dist;
422                                    bestpos = i;
423                            }
424                            biasdist = dist - ((bias[i]) >> (intbiasshift - netbiasshift));
425                            if (biasdist < bestbiasd) {
426                                    bestbiasd = biasdist;
427                                    bestbiaspos = i;
428                            }
429                            betafreq = (freq[i] >> betashift);
430                            freq[i] -= betafreq;
431                            bias[i] += (betafreq << gammashift);
432                    }
433                    freq[bestpos] += beta;
434                    bias[bestpos] -= betagamma;
435                    return (bestbiaspos);
436            }
437    }