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