1 /* 2 * $RCSfile: NeuQuantOpImage.java,v $ 3 * 4 * Copyright (c) 2005 Sun Microsystems, Inc. All rights reserved. 5 * 6 * Use is subject to license terms. 7 * 8 * $Revision: 1.2 $ 9 * $Date: 2005/05/10 01:03:23 $ 10 * $State: Exp $ 11 */ 12 package com.lightcrafts.media.jai.opimage; 13 import java.awt.Rectangle; 14 import java.awt.image.RenderedImage; 15 import java.util.Map; 16 import com.lightcrafts.mediax.jai.ImageLayout; 17 import com.lightcrafts.mediax.jai.LookupTableJAI; 18 import com.lightcrafts.mediax.jai.PlanarImage; 19 import com.lightcrafts.mediax.jai.iterator.RandomIter; 20 import com.lightcrafts.mediax.jai.iterator.RandomIterFactory; 21 import com.lightcrafts.mediax.jai.ROI; 22 23 /** 24 * An <code>OpImage</code> implementing the "ColorQuantizer" operation as 25 * described in <code>com.lightcrafts.mediax.jai.operator.ExtremaDescriptor</code> 26 * based on the median-cut algorithm. 27 * 28 * This is based on a java-version of Anthony Dekker's implementation of 29 * NeuQuant Neural-Net Quantization Algorithm 30 * 31 * NEUQUANT Neural-Net quantization algorithm by Anthony Dekker, 1994. 32 * See "Kohonen neural networks for optimal colour quantization" 33 * in "Network: Computation in Neural Systems" Vol. 5 (1994) pp 351-367. 34 * for a discussion of the algorithm. 35 * 36 * Any party obtaining a copy of these files from the author, directly or 37 * indirectly, is granted, free of charge, a full and unrestricted irrevocable, 38 * world-wide, paid up, royalty-free, nonexclusive right and license to deal 39 * in this software and documentation files (the "Software"), including without 40 * limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, 41 * and/or sell copies of the Software, and to permit persons who receive 42 * copies from any such party to do so, with the only requirement being 43 * that this copyright notice remain intact. 44 * 45 * @see com.lightcrafts.mediax.jai.operator.ExtremaDescriptor 46 * @see ExtremaCRIF 47 */ 48 public class NeuQuantOpImage extends ColorQuantizerOpImage { 49 /** four primes near 500 - assume no image has a length so large 50 * that it is divisible by all four primes 51 */ 52 protected static final int prime1 = 499; 53 protected static final int prime2 = 491; 54 protected static final int prime3 = 487; 55 protected static final int prime4 = 503; 56 57 /* minimum size for input image */ 58 protected static final int minpicturebytes = (3 * prime4); 59 60 /** The size of the histogram. */ 61 private int ncycles; 62 63 /* Program Skeleton 64 ---------------- 65 [select samplefac in range 1..30] 66 [read image from input file] 67 pic = (unsigned char*) malloc(3*width*height); 68 initnet(pic,3*width*height,samplefac); 69 learn(); 70 unbiasnet(); 71 [write output image header, using writecolourmap(f)] 72 inxbuild(); 73 write output image using inxsearch(b,g,r) */ 74 75 /* Network Definitions 76 ------------------- */ 77 78 private final int maxnetpos = maxColorNum - 1; 79 private final int netbiasshift = 4; /* bias for colour values */ 80 81 /* defs for freq and bias */ 82 private final int intbiasshift = 16; /* bias for fractions */ 83 private final int intbias = 1 << intbiasshift; 84 private final int gammashift = 10; /* gamma = 1024 */ 85 private final int gamma = 1 << gammashift; 86 private final int betashift = 10; 87 private final int beta = intbias >> betashift; /* beta = 1/1024 */ 88 private final int betagamma = intbias << (gammashift - betashift); 89 90 /* defs for decreasing radius factor */ 91 private final int initrad = maxColorNum >> 3; 92 private final int radiusbiasshift = 6; /* at 32.0 biased by 6 bits */ 93 private final int radiusbias = 1 << radiusbiasshift; 94 private final int initradius = initrad * radiusbias; /* and decreases by a */ 95 private final int radiusdec = 30; /* factor of 1/30 each cycle */ 96 97 /* defs for decreasing alpha factor */ 98 private final int alphabiasshift = 10; /* alpha starts at 1.0 */ 99 private final int initalpha = 1 << alphabiasshift; 100 101 private int alphadec; /* biased by 10 bits */ 102 103 /* radbias and alpharadbias used for radpower calculation */ 104 private final int radbiasshift = 8; 105 private final int radbias = 1 << radbiasshift; 106 private final int alpharadbshift = alphabiasshift + radbiasshift; 107 private final int alpharadbias = 1 << alpharadbshift; 108 109 // typedef int pixel[4]; /* BGRc */ 110 private int[][] network; /* the network itself - [maxColorNum][4] */ 111 112 private int[] netindex = new int[256]; /* for network lookup - really 256 */ 113 114 private int[] bias = new int[maxColorNum]; /* bias and freq arrays for learning */ 115 private int[] freq = new int[maxColorNum]; 116 private int[] radpower = new int[initrad]; /* radpower for precomputation */ 117 118 /** 119 * Constructs an <code>NeuQuantOpImage</code>. 120 * 121 * @param source The source image. 122 */ NeuQuantOpImage(RenderedImage source, Map config, ImageLayout layout, int maxColorNum, int upperBound, ROI roi, int xPeriod, int yPeriod)123 public NeuQuantOpImage(RenderedImage source, 124 Map config, 125 ImageLayout layout, 126 int maxColorNum, 127 int upperBound, 128 ROI roi, 129 int xPeriod, 130 int yPeriod) { 131 super(source, config, layout, maxColorNum, roi, xPeriod, yPeriod); 132 133 colorMap = null; 134 this.ncycles = upperBound; 135 } 136 train()137 protected synchronized void train() { 138 139 // intialize the network 140 network = new int[maxColorNum][]; 141 for (int i = 0; i < maxColorNum; i++) { 142 network[i] = new int[4]; 143 int[] p = network[i]; 144 p[0] = p[1] = p[2] = (i << (netbiasshift + 8)) / maxColorNum; 145 freq[i] = intbias / maxColorNum; /* 1/maxColorNum */ 146 bias[i] = 0; 147 } 148 149 PlanarImage source = getSourceImage(0); 150 Rectangle rect = source.getBounds(); 151 152 if (roi != null) 153 rect = roi.getBounds(); 154 155 RandomIter iterator = RandomIterFactory.create(source, rect); 156 157 int samplefac = xPeriod * yPeriod; 158 int startX = rect.x / xPeriod; 159 int startY = rect.y / yPeriod; 160 int offsetX = rect.x % xPeriod; 161 int offsetY = rect.y % yPeriod; 162 int pixelsPerLine = (rect.width - 1) / xPeriod + 1; 163 int numSamples = 164 pixelsPerLine * ((rect.height - 1) / yPeriod + 1); 165 166 if (numSamples < minpicturebytes) 167 samplefac = 1; 168 169 alphadec = 30 + ((samplefac - 1) / 3); 170 int pix = 0; 171 172 int delta = numSamples / ncycles; 173 int alpha = initalpha; 174 int radius = initradius; 175 176 int rad = radius >> radiusbiasshift; 177 if (rad <= 1) 178 rad = 0; 179 for (int i = 0; i < rad; i++) 180 radpower[i] = alpha * (((rad * rad - i * i) * radbias) / (rad * rad)); 181 182 int step; 183 if (numSamples < minpicturebytes) 184 step = 3; 185 else if ((numSamples % prime1) != 0) 186 step = 3 * prime1; 187 else { 188 if ((numSamples % prime2) != 0) 189 step = 3 * prime2; 190 else { 191 if ((numSamples % prime3) != 0) 192 step = 3 * prime3; 193 else 194 step = 3 * prime4; 195 } 196 } 197 198 int[] pixel = new int[3]; 199 200 for (int i = 0; i < numSamples;) { 201 int y = (pix / pixelsPerLine + startY) * yPeriod + offsetY; 202 int x = (pix % pixelsPerLine + startX) * xPeriod + offsetX; 203 204 try { 205 iterator.getPixel(x, y, pixel); 206 } catch (Exception e) { 207 continue; 208 } 209 210 int b = pixel[2] << netbiasshift; 211 int g = pixel[1] << netbiasshift; 212 int r = pixel[0] << netbiasshift; 213 214 int j = contest(b , g, r); 215 216 altersingle(alpha, j, b , g, r); 217 if (rad != 0) 218 alterneigh(rad, j, b , g, r); /* alter neighbours */ 219 220 pix += step; 221 if (pix >= numSamples) 222 pix -= numSamples; 223 224 i++; 225 if (i % delta == 0) { 226 alpha -= alpha / alphadec; 227 radius -= radius / radiusdec; 228 rad = radius >> radiusbiasshift; 229 if (rad <= 1) 230 rad = 0; 231 for (j = 0; j < rad; j++) 232 radpower[j] = alpha * (((rad * rad - j * j) * radbias) / (rad * rad)); 233 } 234 } 235 236 unbiasnet(); 237 inxbuild(); 238 createLUT(); 239 setProperty("LUT", colorMap); 240 setProperty("JAI.LookupTable", colorMap); 241 } 242 createLUT()243 private void createLUT() { 244 colorMap = new LookupTableJAI(new byte[3][maxColorNum]); 245 byte[][] map = colorMap.getByteData(); 246 int[] index = new int[maxColorNum]; 247 for (int i = 0; i < maxColorNum; i++) 248 index[network[i][3]] = i; 249 for (int i = 0; i < maxColorNum; i++) { 250 int j = index[i]; 251 map[2][i] = (byte) (network[j][0]); 252 map[1][i] = (byte) (network[j][1]); 253 map[0][i] = (byte) (network[j][2]); 254 } 255 } 256 257 /** Insertion sort of network and building of netindex[0..255] 258 * (to do after unbias) 259 */ inxbuild()260 private void inxbuild() { 261 int previouscol = 0; 262 int startpos = 0; 263 for (int i = 0; i < maxColorNum; i++) { 264 int[] p = network[i]; 265 int smallpos = i; 266 int smallval = p[1]; /* index on g */ 267 /* find smallest in i..maxColorNum-1 */ 268 int j; 269 for (j = i + 1; j < maxColorNum; j++) { 270 int[] q = network[j]; 271 if (q[1] < smallval) { /* index on g */ 272 smallpos = j; 273 smallval = q[1]; /* index on g */ 274 } 275 } 276 int[] q = network[smallpos]; 277 /* swap p (i) and q (smallpos) entries */ 278 if (i != smallpos) { 279 j = q[0]; q[0] = p[0]; p[0] = j; 280 j = q[1]; q[1] = p[1]; p[1] = j; 281 j = q[2]; q[2] = p[2]; p[2] = j; 282 j = q[3]; q[3] = p[3]; p[3] = j; 283 } 284 /* smallval entry is now in position i */ 285 if (smallval != previouscol) { 286 netindex[previouscol] = (startpos + i) >> 1; 287 for (j = previouscol + 1; j < smallval; j++) 288 netindex[j] = i; 289 previouscol = smallval; 290 startpos = i; 291 } 292 } 293 netindex[previouscol] = (startpos + maxnetpos) >> 1; 294 for (int j = previouscol + 1; j < 256; j++) 295 netindex[j] = maxnetpos; /* really 256 */ 296 } 297 298 /** Search for BGR values 0..255 (after net is unbiased) and 299 * return colour index 300 */ findNearestEntry(int r, int g, int b)301 protected byte findNearestEntry(int r, int g, int b) { 302 int bestd = 1000; /* biggest possible dist is 256*3 */ 303 int best = -1; 304 int i = netindex[g]; /* index on g */ 305 int j = i - 1; /* start at netindex[g] and work outwards */ 306 307 while (i < maxColorNum || j >= 0) { 308 if (i < maxColorNum) { 309 int[] p = network[i]; 310 int dist = p[1] - g; /* inx key */ 311 if (dist >= bestd) 312 i = maxColorNum; /* stop iter */ 313 else { 314 i++; 315 if (dist < 0) 316 dist = -dist; 317 int a = p[0] - b; 318 if (a < 0) 319 a = -a; 320 dist += a; 321 if (dist < bestd) { 322 a = p[2] - r; 323 if (a < 0) 324 a = -a; 325 dist += a; 326 if (dist < bestd) { 327 bestd = dist; 328 best = p[3]; 329 } 330 } 331 } 332 } 333 334 if (j >= 0) { 335 int[] p = network[j]; 336 int dist = g - p[1]; /* inx key - reverse dif */ 337 if (dist >= bestd) 338 j = -1; /* stop iter */ 339 else { 340 j--; 341 if (dist < 0) 342 dist = -dist; 343 int a = p[0] - b; 344 if (a < 0) 345 a = -a; 346 dist += a; 347 if (dist < bestd) { 348 a = p[2] - r; 349 if (a < 0) 350 a = -a; 351 dist += a; 352 if (dist < bestd) { 353 bestd = dist; 354 best = p[3]; 355 } 356 } 357 } 358 } 359 } 360 return (byte)best; 361 } 362 363 /** Unbias network to give byte values 0..255 and record 364 * position i to prepare for sort. 365 */ unbiasnet()366 private void unbiasnet() { 367 for (int i = 0; i < maxColorNum; i++) { 368 network[i][0] >>= netbiasshift; 369 network[i][1] >>= netbiasshift; 370 network[i][2] >>= netbiasshift; 371 network[i][3] = i; /* record colour no */ 372 } 373 } 374 375 /** Move adjacent neurons by precomputed 376 * alpha*(1-((i-j)^2/[r]^2)) in radpower[|i-j|] 377 */ 378 alterneigh(int rad, int i, int b, int g, int r)379 private void alterneigh(int rad, int i, int b, int g, int r) { 380 int lo = i - rad; 381 if (lo < -1) 382 lo = -1; 383 int hi = i + rad; 384 if (hi > maxColorNum) 385 hi = maxColorNum; 386 387 int j = i + 1; 388 int k = i - 1; 389 int m = 1; 390 while ((j < hi) || (k > lo)) { 391 int a = radpower[m++]; 392 if (j < hi) { 393 int[] p = network[j++]; 394 // try { 395 p[0] -= (a * (p[0] - b)) / alpharadbias; 396 p[1] -= (a * (p[1] - g)) / alpharadbias; 397 p[2] -= (a * (p[2] - r)) / alpharadbias; 398 // } catch (Exception e) {} // prevents 1.3 miscompilation 399 } 400 if (k > lo) { 401 int[] p = network[k--]; 402 // try { 403 p[0] -= (a * (p[0] - b)) / alpharadbias; 404 p[1] -= (a * (p[1] - g)) / alpharadbias; 405 p[2] -= (a * (p[2] - r)) / alpharadbias; 406 // } catch (Exception e) {} 407 } 408 } 409 } 410 411 412 /** Move neuron i towards biased (b,g,r) by factor alpha. */ altersingle(int alpha, int i, int b, int g, int r)413 private void altersingle(int alpha, int i, int b, int g, int r) { 414 /* alter hit neuron */ 415 int[] n = network[i]; 416 n[0] -= (alpha * (n[0] - b)) / initalpha; 417 n[1] -= (alpha * (n[1] - g)) / initalpha; 418 n[2] -= (alpha * (n[2] - r)) / initalpha; 419 } 420 421 422 /** Search for biased BGR values. */ contest(int b, int g, int r)423 private int contest(int b, int g, int r) { 424 /* finds closest neuron (min dist) and updates freq */ 425 /* finds best neuron (min dist-bias) and returns position */ 426 /* for frequently chosen neurons, freq[i] is high and bias[i] is negative */ 427 /* bias[i] = gamma*((1/maxColorNum)-freq[i]) */ 428 int bestd = ~(((int) 1) << 31); 429 int bestbiasd = bestd; 430 int bestpos = -1; 431 int bestbiaspos = bestpos; 432 433 for (int i = 0; i < maxColorNum; i++) { 434 int[] n = network[i]; 435 int dist = n[0] - b; 436 if (dist < 0) 437 dist = -dist; 438 int a = n[1] - g; 439 if (a < 0) 440 a = -a; 441 dist += a; 442 a = n[2] - r; 443 if (a < 0) 444 a = -a; 445 dist += a; 446 if (dist < bestd) { 447 bestd = dist; 448 bestpos = i; 449 } 450 int biasdist = dist - ((bias[i]) >> (intbiasshift - netbiasshift)); 451 if (biasdist < bestbiasd) { 452 bestbiasd = biasdist; 453 bestbiaspos = i; 454 } 455 int betafreq = (freq[i] >> betashift); 456 freq[i] -= betafreq; 457 bias[i] += (betafreq << gammashift); 458 } 459 freq[bestpos] += beta; 460 bias[bestpos] -= betagamma; 461 return (bestbiaspos); 462 } 463 } 464