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