source: trunk/Lucide/SOURCE/gui/neuquant.cpp @ 172

Last change on this file since 172 was 172, checked in by Eugene Romanenko, 15 years ago

neuquant optimization

File size: 13.0 KB
Line 
1/* NeuQuant Neural-Net Quantization Algorithm
2 * ------------------------------------------
3 *
4 * Copyright (c) 1994 Anthony Dekker
5 *
6 * NEUQUANT Neural-Net quantization algorithm by Anthony Dekker, 1994.
7 * See "Kohonen neural networks for optimal colour quantization"
8 * in "Network: Computation in Neural Systems" Vol. 5 (1994) pp 351-367.
9 * for a discussion of the algorithm.
10 *
11 * Any party obtaining a copy of these files from the author, directly or
12 * indirectly, is granted, free of charge, a full and unrestricted irrevocable,
13 * world-wide, paid up, royalty-free, nonexclusive right and license to deal
14 * in this software and documentation files (the "Software"), including without
15 * limitation the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 * and/or sell copies of the Software, and to permit persons who receive
17 * copies from any such party to do so, with the only requirement being
18 * that this copyright notice remain intact.
19 */
20
21#include <os2.h>
22
23#include <lupixbuf.xh>
24
25#include "neuquant.h"
26#include "globals.h"
27
28
29// Four primes near 500 - assume no image has a length so large
30// that it is divisible by all four primes
31// ==========================================================
32
33#define prime1      499
34#define prime2      491
35#define prime3      487
36#define prime4      503
37
38// ----------------------------------------------------------------
39
40NeuQuantizer::NeuQuantizer( LuPixbuf *pb, int PaletteSize )
41{
42    netsize = PaletteSize;
43    maxnetpos = netsize - 1;
44    initrad = netsize < 8 ? 1 : (netsize >> 3);
45    initradius = (initrad * radiusbias);
46
47    network = NULL;
48
49    network = (pixel *)malloc(netsize * sizeof(pixel));
50    bias = (int *)malloc(netsize * sizeof(int));
51    freq = (int *)malloc(netsize * sizeof(int));
52    radpower = (int *)malloc(initrad * sizeof(int));
53
54    pixbuf = pb;
55    bpp        = pixbuf->getBpp(ev);
56    img_width  = pixbuf->getWidth(ev);
57    img_height = pixbuf->getHeight(ev);
58    // DIB line length in bytes (should be equal to 3 x W)
59    img_line   = img_width * bpp;
60    img_line_len = pixbuf->getRowSize(ev);
61    pixbits    = (unsigned char *)pixbuf->getDataPtr(ev);
62}
63
64NeuQuantizer::~NeuQuantizer()
65{
66    if(network) free(network);
67    if(bias) free(bias);
68    if(freq) free(freq);
69    if(radpower) free(radpower);
70}
71
72///////////////////////////////////////////////////////////////////////////
73// Initialise network in range (0,0,0) to (255,255,255) and set parameters
74// -----------------------------------------------------------------------
75
76void NeuQuantizer::initnet()
77{
78    int i, *p;
79
80    for (i = 0; i < netsize; i++)
81    {
82        p = network[i];
83        p[FI_RGBA_BLUE] = p[FI_RGBA_GREEN] = p[FI_RGBA_RED] = (i << (netbiasshift+8))/netsize;
84        freq[i] = intbias/netsize;  /* 1/netsize */
85        bias[i] = 0;
86    }
87}
88
89///////////////////////////////////////////////////////////////////////////////////////
90// Unbias network to give byte values 0..255 and record position i to prepare for sort
91// ------------------------------------------------------------------------------------
92
93void NeuQuantizer::unbiasnet()
94{
95    int i, j, temp;
96
97    for (i = 0; i < netsize; i++)
98    {
99        for (j = 0; j < 3; j++)
100        {
101            // OLD CODE: network[i][j] >>= netbiasshift;
102            // Fix based on bug report by Juergen Weigert jw@suse.de
103            temp = (network[i][j] + (1 << (netbiasshift - 1))) >> netbiasshift;
104            if (temp > 255) temp = 255;
105            network[i][j] = temp;
106        }
107        network[i][3] = i;          // record colour no
108    }
109}
110
111//////////////////////////////////////////////////////////////////////////////////
112// Insertion sort of network and building of netindex[0..255] (to do after unbias)
113// -------------------------------------------------------------------------------
114
115void NeuQuantizer::inxbuild()
116{
117    int i,j,smallpos,smallval;
118    int *p,*q;
119    int previouscol,startpos;
120
121    previouscol = 0;
122    startpos = 0;
123    for (i = 0; i < netsize; i++)
124    {
125        p = network[i];
126        smallpos = i;
127        smallval = p[FI_RGBA_GREEN];            // index on g
128        // find smallest in i..netsize-1
129        for (j = i+1; j < netsize; j++)
130        {
131            q = network[j];
132            if (q[FI_RGBA_GREEN] < smallval) {  // index on g
133                smallpos = j;
134                smallval = q[FI_RGBA_GREEN];    // index on g
135            }
136        }
137        q = network[smallpos];
138        // swap p (i) and q (smallpos) entries
139        if (i != smallpos)
140        {
141            j = q[FI_RGBA_BLUE];  q[FI_RGBA_BLUE]  = p[FI_RGBA_BLUE];  p[FI_RGBA_BLUE]  = j;
142            j = q[FI_RGBA_GREEN]; q[FI_RGBA_GREEN] = p[FI_RGBA_GREEN]; p[FI_RGBA_GREEN] = j;
143            j = q[FI_RGBA_RED];   q[FI_RGBA_RED]   = p[FI_RGBA_RED];   p[FI_RGBA_RED]   = j;
144            j = q[3];   q[3] = p[3];   p[3] = j;
145        }
146        // smallval entry is now in position i
147        if (smallval != previouscol)
148        {
149            netindex[previouscol] = (startpos+i)>>1;
150            for (j = previouscol+1; j < smallval; j++) {
151                netindex[j] = i;
152            }
153            previouscol = smallval;
154            startpos = i;
155        }
156    }
157    netindex[previouscol] = (startpos+maxnetpos)>>1;
158    for (j = previouscol+1; j < 256; j++)
159        netindex[j] = maxnetpos; // really 256
160}
161
162///////////////////////////////////////////////////////////////////////////////
163// Search for BGR values 0..255 (after net is unbiased) and return colour index
164// ----------------------------------------------------------------------------
165
166int NeuQuantizer::inxsearch(int b, int g, int r)
167{
168    int i, j, dist, a, bestd;
169    int *p;
170    int best;
171
172    bestd = 1000;       // biggest possible dist is 256*3
173    best = -1;
174    i = netindex[g];    // index on g
175    j = i-1;            // start at netindex[g] and work outwards
176
177    while ((i < netsize) || (j >= 0))
178    {
179        if (i < netsize)
180        {
181            p = network[i];
182            dist = p[FI_RGBA_GREEN] - g;                // inx key
183            if (dist >= bestd)
184                i = netsize;    // stop iter
185            else
186            {
187                i++;
188                if (dist < 0)
189                    dist = -dist;
190                a = p[FI_RGBA_BLUE] - b;
191                if (a < 0)
192                    a = -a;
193                dist += a;
194                if (dist < bestd)
195                {
196                    a = p[FI_RGBA_RED] - r;
197                    if (a<0)
198                        a = -a;
199                    dist += a;
200                    if (dist < bestd) {
201                        bestd = dist;
202                        best = p[3];
203                    }
204                }
205            }
206        }
207        if (j >= 0)
208        {
209            p = network[j];
210            dist = g - p[FI_RGBA_GREEN];            // inx key - reverse dif
211            if (dist >= bestd)
212                j = -1; // stop iter
213            else
214            {
215                j--;
216                if (dist < 0)
217                    dist = -dist;
218                a = p[FI_RGBA_BLUE] - b;
219                if (a<0)
220                    a = -a;
221                dist += a;
222                if (dist < bestd)
223                {
224                    a = p[FI_RGBA_RED] - r;
225                    if (a<0)
226                        a = -a;
227                    dist += a;
228                    if (dist < bestd) {
229                        bestd = dist;
230                        best = p[3];
231                    }
232                }
233            }
234        }
235    }
236    return best;
237}
238
239///////////////////////////////
240// Search for biased BGR values
241// ----------------------------
242
243int NeuQuantizer::contest(int b, int g, int r)
244{
245    // finds closest neuron (min dist) and updates freq
246    // finds best neuron (min dist-bias) and returns position
247    // for frequently chosen neurons, freq[i] is high and bias[i] is negative
248    // bias[i] = gamma*((1/netsize)-freq[i])
249
250    int i,dist,a,biasdist,betafreq;
251    int bestpos,bestbiaspos,bestd,bestbiasd;
252    int *p,*f, *n;
253
254    bestd = ~(((int) 1)<<31);
255    bestbiasd = bestd;
256    bestpos = -1;
257    bestbiaspos = bestpos;
258    p = bias;
259    f = freq;
260
261    for (i = 0; i < netsize; i++)
262    {
263        n = network[i];
264        dist = n[FI_RGBA_BLUE] - b;
265        if (dist < 0)
266            dist = -dist;
267        a = n[FI_RGBA_GREEN] - g;
268        if (a < 0)
269            a = -a;
270        dist += a;
271        a = n[FI_RGBA_RED] - r;
272        if (a < 0)
273            a = -a;
274        dist += a;
275        if (dist < bestd) {
276            bestd = dist;
277            bestpos = i;
278        }
279        biasdist = dist - ((*p)>>(intbiasshift-netbiasshift));
280        if (biasdist < bestbiasd) {
281            bestbiasd = biasdist;
282            bestbiaspos = i;
283        }
284        betafreq = (*f >> betashift);
285        *f++ -= betafreq;
286        *p++ += (betafreq << gammashift);
287    }
288    freq[bestpos] += beta;
289    bias[bestpos] -= betagamma;
290    return bestbiaspos;
291}
292
293///////////////////////////////////////////////////////
294// Move neuron i towards biased (b,g,r) by factor alpha
295// ----------------------------------------------------
296
297void NeuQuantizer::altersingle(int alpha, int i, int b, int g, int r)
298{
299    int *n;
300
301    n = network[i];             // alter hit neuron
302    n[FI_RGBA_BLUE]  -= (alpha * (n[FI_RGBA_BLUE]  - b)) / initalpha;
303    n[FI_RGBA_GREEN] -= (alpha * (n[FI_RGBA_GREEN] - g)) / initalpha;
304    n[FI_RGBA_RED]   -= (alpha * (n[FI_RGBA_RED]   - r)) / initalpha;
305}
306
307////////////////////////////////////////////////////////////////////////////////////
308// Move adjacent neurons by precomputed alpha*(1-((i-j)^2/[r]^2)) in radpower[|i-j|]
309// ---------------------------------------------------------------------------------
310
311void NeuQuantizer::alterneigh(int rad, int i, int b, int g, int r)
312{
313    int j, k, lo, hi, a;
314    int *p, *q;
315
316    lo = i - rad;   if (lo < -1) lo = -1;
317    hi = i + rad;   if (hi > netsize) hi = netsize;
318
319    j = i+1;
320    k = i-1;
321    q = radpower;
322    while ((j < hi) || (k > lo))
323    {
324        a = (*(++q));
325        if (j < hi)
326        {
327            p = network[j];
328            p[FI_RGBA_BLUE]  -= (a * (p[FI_RGBA_BLUE] - b)) / alpharadbias;
329            p[FI_RGBA_GREEN] -= (a * (p[FI_RGBA_GREEN] - g)) / alpharadbias;
330            p[FI_RGBA_RED]   -= (a * (p[FI_RGBA_RED] - r)) / alpharadbias;
331            j++;
332        }
333        if (k > lo)
334        {
335            p = network[k];
336            p[FI_RGBA_BLUE]  -= (a * (p[FI_RGBA_BLUE] - b)) / alpharadbias;
337            p[FI_RGBA_GREEN] -= (a * (p[FI_RGBA_GREEN] - g)) / alpharadbias;
338            p[FI_RGBA_RED]   -= (a * (p[FI_RGBA_RED] - r)) / alpharadbias;
339            k--;
340        }
341    }
342}
343
344/////////////////////
345// Main Learning Loop
346// ------------------
347
348void NeuQuantizer::learn(int sampling_factor)
349{
350    int i, j, b, g, r;
351    int radius, rad, alpha, step, delta, samplepixels;
352    int alphadec; // biased by 10 bits
353    long pos, lengthcount;
354
355    // image size as viewed by the scan algorithm
356    lengthcount = img_width * img_height * bpp;
357
358    // number of samples used for the learning phase
359    samplepixels = lengthcount / (3 * sampling_factor);
360
361    // decrease learning rate after delta pixel presentations
362    delta = samplepixels / ncycles;
363    if(delta == 0) {
364        // avoid a 'divide by zero' error with very small images
365        delta = 1;
366    }
367
368    // initialize learning parameters
369    alphadec = 30 + ((sampling_factor - 1) / 3);
370    alpha = initalpha;
371    radius = initradius;
372
373    rad = radius >> radiusbiasshift;
374    if (rad <= 1) rad = 0;
375    for (i = 0; i < rad; i++)
376        radpower[i] = alpha*(((rad*rad - i*i)*radbias)/(rad*rad));
377
378    // initialize pseudo-random scan
379    if ((lengthcount % prime1) != 0)
380        step = bpp*prime1;
381    else {
382        if ((lengthcount % prime2) != 0)
383            step = bpp*prime2;
384        else {
385            if ((lengthcount % prime3) != 0)
386                step = bpp*prime3;
387            else
388                step = bpp*prime4;
389        }
390    }
391
392    i = 0;      // iteration
393    pos = 0;    // pixel position
394
395    while (i < samplepixels)
396    {
397        // get next learning sample
398        getSample(pos, &b, &g, &r);
399
400        // find winning neuron
401        j = contest(b, g, r);
402
403        // alter winner
404        altersingle(alpha, j, b, g, r);
405
406        // alter neighbours
407        if (rad) alterneigh(rad, j, b, g, r);
408
409        // next sample
410        pos += step;
411        while (pos >= lengthcount) pos -= lengthcount;
412
413        i++;
414        if (i % delta == 0)
415        {
416            // decrease learning rate and also the neighborhood
417            alpha -= alpha / alphadec;
418            radius -= radius / radiusdec;
419            rad = radius >> radiusbiasshift;
420            if (rad <= 1) rad = 0;
421            for (j = 0; j < rad; j++)
422                radpower[j] = alpha * (((rad*rad - j*j) * radbias) / (rad*rad));
423        }
424    }
425}
426
427
Note: See TracBrowser for help on using the repository browser.