From a85a18dd3af2384240ef849df592ce02eec904f8 Mon Sep 17 00:00:00 2001 From: nick black Date: Sat, 13 Mar 2021 14:46:21 -0500 Subject: [PATCH] kohonen neural net experiment --- src/lib/kohonen.c | 344 +++++++++++++++++++++++++++++++++++++++++++++ src/lib/neuquant.h | 39 +++++ src/lib/sixel.c | 153 +++++--------------- 3 files changed, 417 insertions(+), 119 deletions(-) create mode 100644 src/lib/kohonen.c create mode 100644 src/lib/neuquant.h diff --git a/src/lib/kohonen.c b/src/lib/kohonen.c new file mode 100644 index 000000000..27387d8c4 --- /dev/null +++ b/src/lib/kohonen.c @@ -0,0 +1,344 @@ +#include +#include +#include "neuquant.h" + +// kohonen neural net adapted from dekker's "kohonen neural networks for +// optimal colour quantization" (2009) + +#define netsize 256 +#define maxnetpos (netsize-1) +#define netbiasshift 4 /* bias for color values */ +#define ncycles 100 /* no. of learning cycles */ + +/* defs for freq and bias */ +#define intbiasshift 16 /* bias for fractions */ +#define intbias (((int) 1)<>betashift) /* beta = 1/1024 */ +#define betagamma (intbias<<(gammashift-betashift)) + +/* defs for decreasing radius factor */ +#define initrad (netsize>>3) /* for 256 cols, radius starts */ +#define radiusbiasshift 6 /* at 32.0 biased by 6 bits */ +#define radiusbias (((int) 1)<network[i]; + p[0] = p[1] = p[2] = (i << (netbiasshift+8))/netsize; + p[3] = 0; + freq[i] = intbias/netsize; /* 1/netsize */ + bias[i] = 0; + } + ret->data = data; + ret->linesize = linesize; + ret->samplefac = sample; + ret->leny = leny; + ret->lenx = lenx; + } + return ret; +} + + +/* Unbias network to give byte values 0..255 and record position i to prepare for sort + ----------------------------------------------------------------------------------- */ + +void unbiasnet(kohonenctx* kctx){ + int i,j,temp; + + for (i=0; inetwork[i][j] + (1 << (netbiasshift - 1))) >> netbiasshift; + if (temp > 255) temp = 255; + kctx->network[i][j] = temp; + } + kctx->network[i][3] = i; /* record color no */ + } +} + + +void netcolor(const kohonenctx* kctx, int color, unsigned char rgb[static 3]){ + rgb[0] = kctx->network[color][0]; + rgb[1] = kctx->network[color][1]; + rgb[2] = kctx->network[color][2]; +} + +void inxbuild(kohonenctx* kctx){ + int i,j,smallpos,smallval; + int *p,*q; + int previouscol,startpos; + + previouscol = 0; + startpos = 0; + for (i=0; inetwork[i]; + smallpos = i; + smallval = p[1]; /* index on g */ + /* find smallest in i..netsize-1 */ + for (j=i+1; jnetwork[j]; + if (q[1] < smallval) { /* index on g */ + smallpos = j; + smallval = q[1]; /* index on g */ + } + } + q = kctx->network[smallpos]; + /* swap p (i) and q (smallpos) entries */ + if (i != smallpos) { + j = q[0]; q[0] = p[0]; p[0] = j; + j = q[1]; q[1] = p[1]; p[1] = j; + j = q[2]; q[2] = p[2]; p[2] = j; + j = q[3]; q[3] = p[3]; p[3] = j; + } + /* smallval entry is now in position i */ + if (smallval != previouscol) { + netindex[previouscol] = (startpos+i)>>1; + for (j=previouscol+1; j>1; + for (j=previouscol+1; j<256; j++) netindex[j] = maxnetpos; /* really 256 */ +} + + +int inxsearch(kohonenctx* kctx, int r, int g, int b){ + int i,j,dist,a,bestd; + int *p; + int best; + + bestd = 1000; /* biggest possible dist is 256*3 */ + best = -1; + i = netindex[g]; /* index on g */ + j = i-1; /* start at netindex[g] and work outwards */ + + while ((i=0)) { + if (inetwork[i]; + dist = p[1] - g; /* inx key */ + if (dist >= bestd) i = netsize; /* stop iter */ + else { + i++; + if (dist<0) dist = -dist; + a = p[0] - r; if (a<0) a = -a; + dist += a; + if (dist=0) { + p = kctx->network[j]; + dist = g - p[1]; /* inx key - reverse dif */ + if (dist >= bestd) j = -1; /* stop iter */ + else { + j--; + if (dist<0) dist = -dist; + a = p[0] - r; if (a<0) a = -a; + dist += a; + if (distnetwork[i]; + dist = n[0] - r; if (dist<0) dist = -dist; + a = n[1] - g; if (a<0) a = -a; + dist += a; + a = n[2] - b; if (a<0) a = -a; + dist += a; + if (dist>(intbiasshift-netbiasshift)); + if (biasdist> betashift); + *f++ -= betafreq; + *p++ += (betafreq<network[i]; /* alter hit neuron */ + *n -= (alpha*(*n - r)) / initalpha; + n++; + *n -= (alpha*(*n - g)) / initalpha; + n++; + *n -= (alpha*(*n - b)) / initalpha; +} + + +/* Move adjacent neurons by precomputed alpha*(1-((i-j)^2/[r]^2)) in radpower[|i-j|] + --------------------------------------------------------------------------------- */ + +static void +alterneigh(kohonenctx* kctx, int rad, int i, int r, int g, int b){ + int j, k, lo, hi, a; + int *p, *q; + + lo = i - rad; if (lo<-1) lo=-1; + hi = i + rad; if (hi>netsize) hi=netsize; + + j = i + 1; + k = i - 1; + q = radpower; + while ((j < hi) || (k > lo)) { + a = (*(++q)); + if (jnetwork[j]; + *p -= (a*(*p - r)) / alpharadbias; + p++; + *p -= (a*(*p - g)) / alpharadbias; + p++; + *p -= (a*(*p - b)) / alpharadbias; + j++; + } + if (k>lo) { + p = kctx->network[k]; + *p -= (a*(*p - r)) / alpharadbias; + p++; + *p -= (a*(*p - g)) / alpharadbias; + p++; + *p -= (a*(*p - b)) / alpharadbias; + k--; + } + } +} + + +/* Main Learning Loop + ------------------ */ + +static inline +const unsigned char* get_pixel(const kohonenctx* kctx, int pixel){ + int line = pixel / kctx->lenx; + int offx = pixel % kctx->lenx; + return kctx->data + kctx->linesize * line + offx * 4; +} + +void learn(kohonenctx* kctx){ + int radius,rad,alpha,step,delta,samplepixels; + int i,j; + + alphadec = 30 + ((kctx->samplefac - 1) / 3); + int pixel = 0; + int lim = kctx->leny * kctx->lenx; + samplepixels = lim / (3 * kctx->samplefac); + delta = samplepixels / ncycles; + alpha = initalpha; + radius = initradius; + + rad = radius >> radiusbiasshift; + if (rad <= 1) rad = 0; + for (i=0; i= lim){ + pixel -= lim; + } + + i++; + if(i % delta == 0){ + alpha -= alpha / alphadec; + radius -= radius / radiusdec; + rad = radius >> radiusbiasshift; + if(rad <= 1){ + rad = 0; + } + for(j=0; j 256)) - -// first pass: extract up to 256 sixelspace colors over arbitrarily many sixels -// sixelspace is 0..100 corresponding to 0..255, lame =[ -typedef struct colortable { - int colors; - int sixelcount; - unsigned char table[CENTSIZE * MAXCOLORS]; // components + dtable index -} colortable; // second pass: construct data for extracted colors over the sixels typedef struct sixeltable { - colortable* ctab; unsigned char* data; // |colors|x|sixelcount|-byte arrays + kohonenctx* kctx; + int sixelcount; } sixeltable; -static inline int -ctable_to_dtable(const unsigned char* ctable){ - if(MAXCOLORS > 256){ - return ctable[RGBSIZE] * 256 + ctable[RGBSIZE + 1]; - }else{ - return ctable[RGBSIZE]; - } -} - -static inline void -dtable_to_ctable(int dtable, unsigned char* ctable){ - if(MAXCOLORS > 256){ - ctable[RGBSIZE] = dtable / 256; - ctable[RGBSIZE + 1] = dtable % 256; - }else{ - ctable[RGBSIZE] = dtable; - } -} - -// returns the index at which the provided color can be found *in the -// dtable*, possibly inserting it into the ctable. returns -1 if the -// color is not in the table and the table is full. -static int -find_color(colortable* ctab, unsigned char comps[static RGBSIZE]){ - int i; - if(ctab->colors){ - int l, r; - l = 0; - r = ctab->colors - 1; - do{ - i = l + (r - l) / 2; -//fprintf(stderr, "%02x%02x%02x L %d R %d m %d\n", comps[0], comps[1], comps[2], l, r, i); - int cmp = memcmp(ctab->table + i * CENTSIZE, comps, RGBSIZE); - if(cmp == 0){ - return ctable_to_dtable(ctab->table + i * CENTSIZE); - } - if(cmp < 0){ - l = i + 1; - }else{ // key is smaller - r = i - 1; - } -//fprintf(stderr, "BCMP: %d L %d R %d m: %d\n", cmp, l, r, i); - }while(l <= r); - if(r < 0){ - i = 0; - }else if(l == ctab->colors){ - i = ctab->colors; - }else{ - i = l; - } - if(ctab->colors == MAXCOLORS){ - return -1; - } - if(i < ctab->colors){ - memmove(ctab->table + (i + 1) * CENTSIZE, ctab->table + i * CENTSIZE, - (ctab->colors - i) * CENTSIZE); - } - }else{ - i = 0; - } - memcpy(ctab->table + i * CENTSIZE, comps, RGBSIZE); - dtable_to_ctable(ctab->colors, ctab->table + i * CENTSIZE); - ++ctab->colors; - return ctab->colors - 1; - //return ctable_to_dtable(ctab->table + i * CENTSIZE); -} - // rather inelegant preprocess of the entire image. colors are converted to the // 100x100x100 sixel colorspace, and built into a table. static int @@ -104,19 +20,17 @@ extract_ctable_inner(const uint32_t* data, int linesize, int begy, int begx, for(int visy = begy ; visy < (begy + leny) ; visy += 6){ for(int visx = begx ; visx < (begx + lenx) ; visx += 1){ for(int sy = visy ; sy < (begy + leny) && sy < visy + 6 ; ++sy){ - const uint32_t* rgb = (const uint32_t*)(data + (linesize / CENTSIZE * sy) + visx); + const uint32_t* rgb = (const uint32_t*)(data + (linesize / 4 * sy) + visx); if(rgba_trans_p(ncpixel_a(*rgb))){ continue; } - unsigned char comps[RGBSIZE]; - break_sixel_comps(comps, *rgb); - int c = find_color(stab->ctab, comps); + int c = inxsearch(stab->kctx, ncpixel_r(*rgb), ncpixel_g(*rgb), ncpixel_b(*rgb)); if(c < 0){ //fprintf(stderr, "FAILED FINDING COLOR AUGH\n"); return -1; } - stab->data[c * stab->ctab->sixelcount + pos] |= (1u << (sy - visy)); -//fprintf(stderr, "color %d pos %d: 0x%x\n", c, pos, stab->data[c * stab->ctab->sixelcount + pos]); + stab->data[c * stab->sixelcount + pos] |= (1u << (sy - visy)); +//fprintf(stderr, "color %d pos %d: 0x%x\n", c, pos, stab->data[c * stab->sixelcount + pos]); } ++pos; } @@ -126,8 +40,7 @@ extract_ctable_inner(const uint32_t* data, int linesize, int begy, int begx, static inline void initialize_stable(sixeltable* stab){ - stab->ctab->colors = 0; - memset(stab->data, 0, stab->ctab->sixelcount * MAXCOLORS); + memset(stab->data, 0, stab->sixelcount * MAXCOLORS); } // Use as many of the original colors as we can, but not more than will fit @@ -182,42 +95,43 @@ write_sixel_data(FILE* fp, int lenx, sixeltable* stab){ // unspecified pixels, which we do not want. //fprintf(fp, "\"1;1;%d;%d", lenx, leny); - for(int i = 0 ; i < stab->ctab->colors ; ++i){ - const unsigned char* rgb = stab->ctab->table + i * CENTSIZE; - fprintf(fp, "#%d;2;%u;%u;%u", i, rgb[0], rgb[1], rgb[2]); + for(int i = 0 ; i < MAXCOLORS ; ++i){ + unsigned char rgb[3]; + netcolor(stab->kctx, i, rgb); +fprintf(stderr, "RGB[%03d]: %3d %3d %3d\n", i, rgb[0], rgb[1], rgb[2]); + fprintf(fp, "#%d;2;%u;%u;%u", i, rgb[0] * 100 / 255, rgb[1] * 100 / 255, rgb[2] * 100 / 255); } int p = 0; - while(p < stab->ctab->sixelcount){ - for(int i = 0 ; i < stab->ctab->colors ; ++i){ + while(p < stab->sixelcount){ + for(int i = 0 ; i < MAXCOLORS ; ++i){ int printed = 0; int seenrle = 0; // number of repetitions unsigned char crle = 0; // character being repeated - int idx = ctable_to_dtable(stab->ctab->table + i * CENTSIZE); - for(int m = p ; m < stab->ctab->sixelcount && m < p + lenx ; ++m){ -//fprintf(stderr, "%d ", idx * stab->ctab->sixelcount + m); -//fputc(stab->data[idx * stab->ctab->sixelcount + m] + 63, stderr); + for(int m = p ; m < stab->sixelcount && m < p + lenx ; ++m){ +//fprintf(stderr, "%d ", i * stab->sixelcount + m); +//fputc(stab->data[i * stab->sixelcount + m] + 63, stderr); if(seenrle){ - if(stab->data[idx * stab->ctab->sixelcount + m] == crle){ + if(stab->data[i * stab->sixelcount + m] == crle){ ++seenrle; }else{ write_rle(&printed, i, fp, seenrle, crle); seenrle = 1; - crle = stab->data[idx * stab->ctab->sixelcount + m]; + crle = stab->data[i * stab->sixelcount + m]; } }else{ seenrle = 1; - crle = stab->data[idx * stab->ctab->sixelcount + m]; + crle = stab->data[i * stab->sixelcount + m]; } } if(crle){ write_rle(&printed, i, fp, seenrle, crle); } - if(i + 1 < stab->ctab->colors){ + if(i + 1 < MAXCOLORS){ if(printed){ fputc('$', fp); } }else{ - if(p + lenx < stab->ctab->sixelcount){ + if(p + lenx < stab->sixelcount){ fputc('-', fp); } } @@ -264,26 +178,27 @@ int sixel_blit_inner(ncplane* nc, int placey, int placex, int lenx, int sixel_blit(ncplane* nc, int placey, int placex, int linesize, const void* data, int begy, int begx, int leny, int lenx, unsigned cellpixx){ - colortable* ctab = malloc(sizeof(*ctab)); - if(ctab == NULL){ - return -1; - } - ctab->sixelcount = (lenx - begx) * ((leny - begy + 5) / 6); + int sixelcount = (lenx - begx) * ((leny - begy + 5) / 6); sixeltable stable = { - .ctab = ctab, - .data = malloc(MAXCOLORS * ctab->sixelcount), + .data = malloc(MAXCOLORS * sixelcount), + .sixelcount = sixelcount, }; if(stable.data == NULL){ - free(ctab); return -1; } + if((stable.kctx = initnet(data, leny, linesize, lenx, 1)) == NULL){ + free(stable.data); + } + learn(stable.kctx); + unbiasnet(stable.kctx); + inxbuild(stable.kctx); if(extract_color_table(data, linesize, begy, begx, leny, lenx, &stable)){ - free(ctab); free(stable.data); + freenet(stable.kctx); return -1; } int r = sixel_blit_inner(nc, placey, placex, lenx, &stable, cellpixx); free(stable.data); - free(ctab); + freenet(stable.kctx); return r; }