~ubuntu-branches/ubuntu/wily/pngnq/wily : revision 3

1

/* NeuQuant Neural-Net Quantization Algorithm

2

* ------------------------------------------

3

*

4

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

* See also http://members.ozemail.com.au/~dekker/NEUQUANT.HTML

11

*

12

* Any party obtaining a copy of these files from the author, directly or

13

* indirectly, is granted, free of charge, a full and unrestricted irrevocable,

14

* world-wide, paid up, royalty-free, nonexclusive right and license to deal

15

* in this software and documentation files (the "Software"), including without

16

* limitation the rights to use, copy, modify, merge, publish, distribute, sublicense,

17

* and/or sell copies of the Software, and to permit persons who receive

18

* copies from any such party to do so, with the only requirement being

19

* that this copyright notice remain intact.

20

*

21

*

22

* Modified to process 32bit RGBA images.

23

* Stuart Coyle July 2004

24

*/

25

26

27

#include "neuquant32.h"

28

29

30

/* Network Definitions

31

------------------- */

32

33

#define maxnetpos (netsize-1)

34

#define netbiasshift 4 /* bias for colour values */

35

#define ncycles 100 /* no. of learning cycles */

36

37

/* defs for freq and bias */

38

#define intbiasshift 16 /* bias for fractions */

39

#define intbias (((int) 1)<<intbiasshift)

40

#define gammashift 10 /* gamma = 1024 */

41

#define gamma (((int) 1)<<gammashift)

42

#define betashift 10

43

#define beta (intbias>>betashift) /* beta = 1/1024 */

44

#define betagamma (intbias<<(gammashift-betashift))

45

46

/* defs for decreasing radius factor */

47

#define initrad (netsize>>3) /* for 256 cols, radius starts */

48

#define radiusbiasshift 6 /* at 32.0 biased by 6 bits */

49

#define radiusbias (((int) 1)<<radiusbiasshift)

50

#define initradius (initrad*radiusbias) /* and decreases by a */

51

#define radiusdec 30 /* factor of 1/30 each cycle */

52

53

/* defs for decreasing alpha factor */

54

#define alphabiasshift 10 /* alpha starts at 1.0 */

55

#define initalpha (((int) 1)<<alphabiasshift)

56

int alphadec; /* biased by 10 bits */

57

58

/* radbias and alpharadbias used for radpower calculation */

59

#define radbiasshift 8

60

#define radbias (((int) 1)<<radbiasshift)

61

#define alpharadbshift (alphabiasshift+radbiasshift)

62

#define alpharadbias (((int) 1)<<alpharadbshift)

63

64

65

/* Types and Global Variables

66

-------------------------- */

67

68

static unsigned char *thepicture; /* the input image itself */

69

static int lengthcount; /* lengthcount = H*W*3 */

70

71

static int samplefac; /* sampling factor 1..30 */

72

73

74

typedef int nq_pixel[5]; /* ABGRc */

75

static nq_pixel network[netsize]; /* the network itself */

76

77

static int netindex[256]; /* for network lookup - really 256 */

78

79

static int bias [netsize]; /* bias and freq arrays for learning */

80

static int freq [netsize];

81

static int radpower[initrad]; /* radpower for precomputation */

82

83

84

/* Initialise network in range (0,0,0,0) to (255,255,255,255) and set parameters

85

----------------------------------------------------------------------- */

86

87

void initnet(thepic,len,sample)

88

unsigned char *thepic;

89

int len;

90

int sample;

91

{

92

register int i;

93

register int *p;

94

95

thepicture = thepic;

96

lengthcount = len;

97

samplefac = sample;

98

99

for (i=0; i<netsize; i++) {

100

p = network[i];

101

p[0] = p[1] = p[2] = p[3] = (i << (netbiasshift+8))/netsize;

102

freq[i] = intbias/netsize; /* 1/netsize */

103

bias[i] = 0;

104

}

105

}

106

107

108

/* Unbias network to give byte values 0..255 and record position i to prepare for sort

109

----------------------------------------------------------------------------------- */

110

111

void unbiasnet()

112

{

113

int i,j,temp;

114

115

for (i=0; i<netsize; i++) {

116

for (j=0; j<4; j++) {

117

/* OLD CODE: network[i][j] >>= netbiasshift; */

118

/* Fix based on bug report by Juergen Weigert jw@suse.de */

119

temp = (network[i][j] + (1 << (netbiasshift - 1))) >> netbiasshift;

120

if (temp > 255) temp = 255;

121

network[i][j] = temp;

122

}

123

network[i][4] = i; /* record colour no */

124

}

125

}

126

127

128

/* Output colour map

129

----------------- */

130

131

void writecolourmap(f)

132

FILE *f;

133

{

134

int i,j;

135

136

for (i=3; i>=0; i--)

137

for (j=0; j<netsize; j++)

138

putc(network[j][i], f);

139

}

140

141

142

/* Output colormap to unsigned char ptr in RGBA format */

143

void getcolormap(map)

144

unsigned char *map;

145

{

146

int i,j;

147

for(j=0; j<netsize; j++){

148

for (i=3; i>=0; i--){

149

*map = network[j][i];

150

map++;

151

}

152

}

153

}

154

155

156

/* Insertion sort of network and building of netindex[0..255] (to do after unbias)

157

------------------------------------------------------------------------------- */

158

159

void inxbuild()

160

{

161

register int i,j,smallpos,smallval;

162

register int *p,*q;

163

int previouscol,startpos;

164

165

previouscol = 0;

166

startpos = 0;

167

for (i=0; i<netsize; i++) {

168

p = network[i];

169

smallpos = i;

170

smallval = p[2]; /* index on g */

171

/* find smallest in i..netsize-1 */

172

for (j=i+1; j<netsize; j++) {

173

q = network[j];

174

if (q[2] < smallval) { /* index on g */

175

smallpos = j;

176

smallval = q[2]; /* index on g */

177

}

178

}

179

q = network[smallpos];

180

/* swap p (i) and q (smallpos) entries */

181

if (i != smallpos) {

182

j = q[0]; q[0] = p[0]; p[0] = j;

183

j = q[1]; q[1] = p[1]; p[1] = j;

184

j = q[2]; q[2] = p[2]; p[2] = j;

185

j = q[3]; q[3] = p[3]; p[3] = j;

186

j = q[4]; q[4] = p[4]; p[4] = j;

187

}

188

/* smallval entry is now in position i */

189

if (smallval != previouscol) {

190

netindex[previouscol] = (startpos+i)>>1;

191

for (j=previouscol+1; j<smallval; j++) netindex[j] = i;

192

previouscol = smallval;

193

startpos = i;

194

}

195

}

196

netindex[previouscol] = (startpos+maxnetpos)>>1;

197

for (j=previouscol+1; j<256; j++) netindex[j] = maxnetpos; /* really 256 */

198

}

199

200

201

/* Search for ABGR values 0..255 (after net is unbiased) and return colour index

202

---------------------------------------------------------------------------- */

203

int inxsearch(al,b,g,r)

204

register int al,b,g,r;

205

{

206

register int i,j,dist,a,bestd;

207

register int *p;

208

int best;

209

210

bestd = 1000; /* biggest possible dist is 256*3 */

211

best = -1;

212

i = netindex[g]; /* index on g */

213

j = i-1; /* start at netindex[g] and work outwards */

214

215

while ((i<netsize) || (j>=0)) {

216

if (i<netsize) {

217

p = network[i];

218

dist = p[2] - g; /* inx key */

219

if (dist >= bestd) i = netsize; /* stop iter */

220

else {

221

i++;

222

if (dist<0) dist = -dist;

223

a = p[1] - b; if (a<0) a = -a;

224

dist += a;

225

if (dist<bestd) {

226

a = p[3] - r; if (a<0) a = -a;

227

dist += a;

228

}

229

if(dist<bestd) {

230

a = p[0] - al; if (a<0) a = -a;

231

dist += a;

232

}

233

if (dist<bestd) {bestd=dist; best=p[4];}

234

}

235

}

236

237

if (j>=0) {

238

p = network[j];

239

dist = g - p[2]; /* inx key - reverse dif */

240

if (dist >= bestd) j = -1; /* stop iter */

241

else {

242

j--;

243

if (dist<0) dist = -dist;

244

a = p[1] - b; if (a<0) a = -a;

245

dist += a;

246

if (dist<bestd) {

247

a = p[3] - r; if (a<0) a = -a;

248

dist += a;

249

}

250

if(dist<bestd) {

251

a = p[0] - al; if (a<0) a = -a;

252

dist += a;

253

}

254

if (dist<bestd) {bestd=dist; best=p[4];}

255

}

256

}

257

}

258

259

return(best);

260

}

261

262

263

264

265

/* Search for biased ABGR values

266

---------------------------- */

267

268

int contest(al,b,g,r)

269

register int al,b,g,r;

270

{

271

/* finds closest neuron (min dist) and updates freq */

272

/* finds best neuron (min dist-bias) and returns position */

273

/* for frequently chosen neurons, freq[i] is high and bias[i] is negative */

274

/* bias[i] = gamma*((1/netsize)-freq[i]) */

275

276

register int i,dist,a,biasdist,betafreq;

277

int bestpos,bestbiaspos,bestd,bestbiasd;

278

register int *p,*f, *n;

279

280

bestd = ~(((int) 1)<<31);

281

bestbiasd = bestd;

282

bestpos = -1;

283

bestbiaspos = bestpos;

284

p = bias;

285

f = freq;

286

287

for (i=0; i<netsize; i++) {

288

n = network[i];

289

dist = n[0] - al; if (dist<0) dist = -dist;

290

a = n[1] - b; if (a<0) a = -a;

291

dist += a;

292

a = n[2] - g; if (a<0) a = -a;

293

dist += a;

294

a = n[3] - r; if (a<0) a = -a;

295

dist += a;

296

if (dist<bestd) {bestd=dist; bestpos=i;}

297

biasdist = dist - ((*p)>>(intbiasshift-netbiasshift));

298

if (biasdist<bestbiasd) {bestbiasd=biasdist; bestbiaspos=i;}

299

betafreq = (*f >> betashift);

300

*f++ -= betafreq;

301

*p++ += (betafreq<<gammashift);

302

}

303

freq[bestpos] += beta;

304

bias[bestpos] -= betagamma;

305

return(bestbiaspos);

306

}

307

308

309

/* Move neuron i towards biased (a,b,g,r) by factor alpha

310

---------------------------------------------------- */

311

312

void altersingle(alpha,i,al,b,g,r)

313

register int alpha,i,b,g,r;

314

{

315

register int *n;

316

317

n = network[i]; /* alter hit neuron */

318

*n -= (alpha*(*n - al)) / initalpha;

319

n++;

320

*n -= (alpha*(*n - b)) / initalpha;

321

n++;

322

*n -= (alpha*(*n - g)) / initalpha;

323

n++;

324

*n -= (alpha*(*n - r)) / initalpha;

325

}

326

327

328

/* Move adjacent neurons by precomputed alpha*(1-((i-j)^2/[r]^2)) in radpower[|i-j|]

329

--------------------------------------------------------------------------------- */

330

331

void alterneigh(rad,i,al,b,g,r)

332

int rad,i;

333

register int al,b,g,r;

334

{

335

register int j,k,lo,hi,a;

336

register int *p, *q;

337

338

lo = i-rad; if (lo<-1) lo=-1;

339

hi = i+rad; if (hi>netsize) hi=netsize;

340

341

j = i+1;

342

k = i-1;

343

q = radpower;

344

while ((j<hi) || (k>lo)) {

345

a = (*(++q));

346

if (j<hi) {

347

p = network[j];

348

*p -= (a*(*p - al)) / alpharadbias;

349

p++;

350

*p -= (a*(*p - b)) / alpharadbias;

351

p++;

352

*p -= (a*(*p - g)) / alpharadbias;

353

p++;

354

*p -= (a*(*p - r)) / alpharadbias;

355

j++;

356

}

357

if (k>lo) {

358

p = network[k];

359

*p -= (a*(*p - al)) / alpharadbias;

360

p++;

361

*p -= (a*(*p - b)) / alpharadbias;

362

p++;

363

*p -= (a*(*p - g)) / alpharadbias;

364

p++;

365

*p -= (a*(*p - r)) / alpharadbias;

366

k--;

367

}

368

}

369

}

370

371

372

/* Main Learning Loop

373

------------------ */

374

375

void learn(int verbose) /* Stu: N.B. added parameter so that main() could control verbosity. */

376

{

377

register int i,j,al,b,g,r;

378

int radius,rad,alpha,step,delta,samplepixels;

379

register unsigned char *p;

380

unsigned char *lim;

381

382

alphadec = 30 + ((samplefac-1)/3);

383

p = thepicture;

384

lim = thepicture + lengthcount;

385

samplepixels = lengthcount/(3*samplefac);

386

delta = samplepixels/ncycles; /* here's a problem with small images: samplepixels < ncycles => delta = 0 */

387

if(delta==0) delta = 1; /* kludge to fix */

388

alpha = initalpha;

389

radius = initradius;

390

391

rad = radius >> radiusbiasshift;

392

if (rad <= 1) rad = 0;

393

for (i=0; i<rad; i++)

394

radpower[i] = alpha*(((rad*rad - i*i)*radbias)/(rad*rad));

395

396

if(verbose) fprintf(stderr,"beginning 1D learning: initial radius=%d\n", rad);

397

398

if ((lengthcount%prime1) != 0) step = 4*prime1;

399

else {

400

if ((lengthcount%prime2) !=0) step = 4*prime2;

401

else {

402

if ((lengthcount%prime3) !=0) step = 4*prime3;

403

else step = 4*prime4;

404

}

405

}

406

407

i = 0;

408

while (i < samplepixels) {

409

al = p[3] << netbiasshift;

410

b = p[2] << netbiasshift;

411

g = p[1] << netbiasshift;

412

r = p[0] << netbiasshift;

413

j = contest(al,b,g,r);

414

415

altersingle(alpha,j,al,b,g,r);

416

if (rad) alterneigh(rad,j,al,b,g,r); /* alter neighbours */

417

418

p += step;

419

if (p >= lim) p -= lengthcount;

420

421

i++;

422

if (i%delta == 0) { /* FPE here if delta=0*/

423

alpha -= alpha / alphadec;

424

radius -= radius / radiusdec;

425

rad = radius >> radiusbiasshift;

426

if (rad <= 1) rad = 0;

427

for (j=0; j<rad; j++)

428

radpower[j] = alpha*(((rad*rad - j*j)*radbias)/(rad*rad));

429

}

430

}

431

if(verbose) fprintf(stderr,"finished 1D learning: final alpha=%f !\n",((float)alpha)/initalpha);

432

}