/* 
Lightness computation - Andrew Blake, Dept. CS, Edinburgh University.
adjusted by Gavin Brelstaff 16/11/84.

Altered to allow for a SQUARE GRID OF SKEWED TRIANGLES rather than
equilaterals.

  WITH NON-FLOATING POINT VARIATION OF THE OVERRELAXTION PARAMETER,
  implemented by Gavin Brelstaff 21/11/84.

Extended to pictures greater than 128x128 by Gavin Brelstaff 14/12/84.

WITH LOGARITHMIC INPUT CONVERSION
WITH OUTPUT OF MEAN ZERO            both by Gavin Brelstaff   14/1/85
 
With threshold percentage option 4/2/85 Gavin Brelstaff.

picture format:
  first byte is image size (square image) if less than 128, else use -x option,
  then image grey data in left to right, top to bottom format.

to run, 

  see vision:light.doc

*/
 
 
#include <stdio.h>
#include <math.h>
#define progname "light"

#define MAXSIZE 128 /* max. image size */
#define MAXAREA MAXSIZE*MAXSIZE

short nbits = 0;                /* picture scaling: nbits +8 bits */
long x[MAXAREA], x0[MAXAREA] ;  /* pics: relaxing, initial */
char conflag[MAXAREA];          /* flags indicating connectivity of array
                                   and results of thresholding (top 2 bits) */
char diag[MAXAREA];             /* number of surrounding equilateral triangles*/

short asize ;    /* pic dimensions */
long area;
short header=1;/* set if pic file contains a header char */
     
float th_percentage=6.0;         /* input parameter the threshold percentage*/
int th=100;                      /* threshold value for (gradient*gradient) */
#define MAXITER 2000
int iter=MAXITER;               /* maximum no of iterations */
int prec=16;                    /* constant specifying precision of integer
                                   computation */

FILE *fp, *fq = NULL;            /* input,output picture channels */
 
short JORpara  = 8;               /* The relaxtion parameter w as in YOUNG p72*/
                                  /* w = 1 + JORpara/16. For over-relaxtion   */
                                  /* w usually lies in the range 1.0 <w <2.0  */
                                  /* N.B. Blakes lambda = 1 - w .             */

char tflag=0;   /* set if threshold specified */
char mflag=1;   /* set if mean zero option for output */
char bflag=0;   /* set if beeps for attention */
char lflag=1;   /* set if log conversion used */
char mess_flag=1;   /* set if messages reported to screen */

typedef long (*PFL)();
PFL fnp;                         /* global pointer to function that multiplies*/
                                 /* values (x) by the JORpara/16 using shifts */
                                 /* to avoid slow floating point calculations */

/* functions that multiply by sixteenths using shifting op.s */
long fn1(x) long x;{return(x>>4);}
long fn2(x) long x;{return(x>>3);}
long fn3(x) long x;{return((x>>3)+(x>>4));}
long fn4(x) long x;{return(x>>2);}
long fn5(x) long x;{return((x>>2)+(x>>4));}
long fn6(x) long x;{return((x>>2)+(x>>3));}
long fn7(x) long x;{return((x>>2)+(x>>3)+(x>>4));}
long fn8(x) long x;{return(x>>1);}
long fn9(x) long x;{return((x>>1)+(x>>4));}
long fn10(x) long x;{return((x>>1)+(x>>3));}   /* e.g. returns 10x/16  */
long fn11(x) long x;{return((x>>1)+(x>>3)+(x>>4));}
long fn12(x) long x;{return((x>>1)+(x>>2));}
long fn13(x) long x;{return((x>>1)+(x>>2)+(x>>4));}
long fn14(x) long x;{return((x>>1)+(x>>2)+(x>>3));}
long fn15(x) long x;{return((x>>1)+(x>>2)+(x>>3)+(x>>4));}
long fn16(x) long x;{return(x);}
 
setfnpointer(l)   /* sets function pointer to approipriate shifting function */
short l;
{
   switch(l) {
     case(1): fnp=fn1;break;
     case(2): fnp=fn2;break;
     case(3): fnp=fn3;break;
     case(4): fnp=fn4;break;
     case(5): fnp=fn5;break;
     case(6): fnp=fn6;break;
     case(7): fnp=fn7;break;
     case(8): fnp=fn8;break;
     case(9): fnp=fn9;break;
     case(10): fnp=fn10;break;
     case(11): fnp=fn11;break;
     case(12): fnp=fn12;break;
     case(13): fnp=fn13;break;
     case(14): fnp=fn14;break;
     case(15): fnp=fn15;break;
     case(16): fnp=fn16;break;
  }
}
 
 
main(argc,argv)
char **argv;
{

  argv++;
  while ((argc>1)&&((*argv)[0] == '-'))
  {

    switch((*argv)[1]) {
    case 'j':
      sscanf(*argv+2,"%hd",&JORpara);
      if(JORpara<1 || JORpara>16) error("JOR parameter out of range 1 .. 16");
      break;
    case 'x':
      sscanf(*argv+2,"%hd",&asize);
      if(asize<0||asize>MAXSIZE) error("picture to big");
      header=0;
      break;
    case 'b':  /* beep flag */
      bflag=1;
      break;
    case 's':  /* supress messages to screen flag */
      mess_flag=0;
      break;
    case 'l':  /* turn off logarithmic input conversion option */
      lflag=0;
      break;
    case 'm':  /* turn off mean-zero output option */
      mflag=0;
      break;
    case 't':
      tflag=1;
      sscanf(*argv+2,"%f",&th_percentage);
      if(th_percentage<0.0||100.0<th_percentage)
        error("bad threshold percentage");
      break;
    case 'i':
      sscanf(*argv+2,"%d",&iter);
      if(iter<0||MAXITER<iter) error("bad number of iterations");
      break;
    case 'p':
      sscanf(*argv+2,"%d",&prec);
      if(prec<0||16<prec) error("bad precision specified");
      break;
    default:
      error("unknown flag");
      break;
    }
    --argc;argv++;
  }
 
  if (argc<3)
    error("missing arguments");
  else
  {
    if((fp = fopen(*argv++,"r")) == NULL)
      error("can't access input file");
    --argc;

    if((fq = fopen(*argv,"w")) == NULL)
      error("can't access output file");
    fclose(fq);

  }

  setfnpointer(JORpara);
  readpic();
  fclose(fp);
  log_scale(lflag,tflag);

  setupcon();      /* and thresholds */
  scale(prec);
  iterate();
  scale(0);
  if(!mflag)  makepos(x);
  else make_meanzero(x);

  if((fq = fopen(*argv++,"w")) == NULL)
    error("can't reaccess output file");
  writepic(x);
  fclose(fq);
}
 

makepos(y)  /* ensure pic is everywhere positive and in range */
long *y;
{
long *p, max, min, *pend;

  min=100000; max = -100000;
  for (p=y,pend=y+area; p<pend; p++)
  {
    if (*p > max) max = (*p);
    if (*p < min) min = (*p);
  }

  for (p=y,pend=y+area; p<pend; p++)
    *p = (*p) -min;

  max=max-min;
  while (max>255) {
    scale(nbits-1);
    max = ((max+1)>>1);
  }
}

#define TMAX 256
#define TOP_INPUT_INTENSITY 256

log_scale(lflag,tflag)
char lflag,tflag;
{
long log_table[TMAX];  /* temporary log conversion look up table */
double CONST;
register long  *p;
register long i;

    if(!lflag){
      nbits=0;
      if(tflag)
         th = ((int) th_percentage);  /* set threshold option */
      return(0);  
    }

    nbits = 6;     /* 11 + 3 - 8 */    /* global variable */
    if(tflag)  th=((int) ( (TOP_INPUT_INTENSITY<<nbits)*
                            log( (double) (1.0 + th_percentage/100.0))/
                            log( (double) TOP_INPUT_INTENSITY) ) );
    th=th*th;

    log_table[0]=0;
    CONST = 2048.0/log((double) 2.0) ;  /* 2^11/loge2 */
    for(i=1,p=log_table+1;i<TMAX;i++) 
      *p++ = ( (long) (log((double) i)*CONST));

    for (p=x,i=0;i<area;i++ )
      *p++ = log_table[*p];

/*
    for(i=0,p=log_table;i<TMAX;i++) 
      if(mess_flag) fprintf(stderr,"%d %ld\n",i,*p++);
*/
   
}

error(s)
char *s;
{
  if(mess_flag) fprintf(stderr,"%s: %s\n",progname,s);
  exit(1);
}
 
readpic()       /*read from fp, set up x as integer arrays*/
{
 
register long *p;
long i, *pend;
 
 
  if (header)    asize=getc(fp)&0377;
  if (/*(asize<8)||*/(asize>MAXSIZE))
    error("array size out of range");
  area = asize*asize;
 
  for (p=x,pend=x+area; p<pend;)
    *p++ = 0377& getc(fp);
}
 
 
writepic(y)      /* send y out on fq. low bytes only */
long *y;
{
 
long *p;
long *pend,i;
 
  if(header)putc(asize,fq);
 
  for (p=y, pend=y+area; p<pend;)
    putc(*p++,fq);
}


make_meanzero(y)  
long *y;
{
long *p , *pend,total,mean;

  for (p=y,pend=y+area,total=0; p<pend; p++)
    total += *p;

  mean= total/area;

  for (p=y,pend=y+area; p<pend; p++) {
    *p = (*p) -mean;
    if(*p>127) *p =  127;
    else if(*p<-128) *p = -128;
  }
}
 

scale(n)  /* scale x,x0 to n+8 bits */
{
 
register long s, *ptop;
register long *p,*q;
 
 
  s=n-nbits;
  nbits=n;
 
  if (s>0)
  {
    for (p=x,q=x0,ptop=x+area; p<ptop;)
    {
      *p++ = (*p<<s);
      *q++ = (*q<<s);
    }
  }
  else if (s<0)
  {
    s= -s;
    for (p=x,q=x0,ptop=x+area; p<ptop;)
    {
      *p++ = (*p>>s);
      *q++ = (*q>>s);
    }
  }
}


/* constants for setupcon and adjust*/
#define dirNNE 01
#define dirNEE 02
#define dirES  04
#define dirSSW 010
#define dirSWW 020
#define dirWN  040

#define offN  (-asize)
#define offNE (1-asize)
#define offE   1
#define offS   asize
#define offSW (asize-1)
#define offW  (-1)

#define flagU 0100
#define flagD 0200
 
 
/*macros for setupcon and adjust */

#define setru2(dir,off) {if (*conp&(dir)) {\
  if (*conp&flagU) *q = *q+ *p - *(p+(off)) ;\
  (*d)++; } }

#define setrd2(dir,off) {if (*conp&(dir)) {\
  if (*conp&flagD) *q = *q+ *p - *(p+(off)) ;\
  (*d)++; } }

#define setru3(dir,off,offU) {if (*conp&(dir)) {\
  if (*(conp+(offU))&flagU) *q = *q+ *p - *(p+(off)) ;\
  (*d)++; } }

#define setrd3(dir,off,offD) {if (*conp&(dir)) {\
  if (*(conp+(offD))&flagD) *q = *q+ *p - *(p+(off)) ;\
  (*d)++; } }

#define setru4(dir,off1,off2,offU) {if (*conp&(dir)) {\
  if (*(conp+(offU))&flagU) *q = *q+ ((*p<<1)) - *(p+(off1)) - *(p+(off2));\
  (*d)+=2; } }

#define setrd4(dir,off1,off2,offD) {if (*conp&(dir)) {\
  if (*(conp+(offD))&flagD) *q = *q+ ((*p<<1)) - *(p+(off1)) - *(p+(off2));\
  (*d)+=2; } }

#define gu (gval(*p, *(p-asize), *(p-asize+1))  >th)

#define gd (gval(*p, *(p+asize), *(p+asize-1))  >th)

#define adj2(dir,off)  {if (*conp&dir) inc+=(*(p+(off)));}

#define adj3(dir,off1,off2)  {if (*conp&dir) inc+=(*(p+(off1)) + *(p+(off2)));}



long gval(a,b,c)
long a,b,c;
{
register long t;
long g;

  t=a-b; g=t*t;
  t=b-c; g += t*t;
/*  t=c-a; g += t*t;     SKEW VERSION */

  return(g);
}


setupcon()
{
  long *p,*q,*ptop;
  register char *conp,*d;
  long i;
  

/* initialise conflag. all 6 adjacent equilateral exist. */

for (i=0,conp=conflag;i++<area; )
    *conp++=077;


/* set up edge conditions - remove non-existent triangles over the edge */

  for (i=0, conp=conflag; i++<asize; conp+=asize) {
    *conp &= ~(dirSSW|dirSWW|dirWN);
    *(conp+asize-1) &= ~(dirNNE|dirNEE|dirES);
  }

  for (i=0, conp=conflag; i++<asize; conp++) {
    *conp &= ~(dirWN|dirNNE|dirNEE);
    *(conp+area-asize) &= ~(dirES|dirSSW|dirSWW);
  }


/* set up upper/lower triangle gradient flags
   - set flag if gradient in triangle exceeds threshold */

for (i=0,p=x,conp=conflag;i++<area;p++,conp++) {
    if ((*conp&dirNNE)&&(gu))
      *conp |= flagU;
    if ((*conp&dirSSW)&&(gd))
      *conp |= flagD;
}


/* setup r=div(E) in x0  and diag array */

  for (i=0,p=x,q=x0,d=diag,conp=conflag;i++<area;p++,q++,d++,conp++) {
    *q=0;*d=0;
    setru2(dirNNE,offN);
    setrd3(dirNEE,offE ,offNE);
    setru4(dirES ,offE ,offS ,offS );
    setrd2(dirSSW,offS);
    setru3(dirSWW,offW ,offSW);
    setrd4(dirWN ,offW ,offN ,offN );
  }

}  /* end of setupcon */





/* perform one iteration of relaxation */

long adjust()
{
  long *p,*q, inc, v;
  char *conp,*d;
  long change;
  long i;

  change=0;
  for(i=0,p=x,q=x0,d=diag,conp=conflag;i++<area;p++,conp++) {
    inc=0;
    adj2(dirNNE,offN);
    adj2(dirNEE,offE );
    adj3(dirES ,offE ,offS );
    adj2(dirSSW,offS);
    adj2(dirSWW,offW );
    adj3(dirWN ,offW ,offN );

    v= (*p);
    inc = ((inc+ (*q++))/(*d++));

/* over relaxation */
    inc = inc + (*fnp)(inc);    /* N.B. pointer to function used */
    *p  = inc - (*fnp)(v);
    inc = (*p) - v;
    if (inc<0)
      inc = -inc;
    if (change<inc)
      change=inc;
  }/* end of for */

  return(change>>4);
}


#define BEEP   7

iterate()
{
  short i,j;
  long a;

  if(mess_flag) fprintf(stderr,"threshold = %3.2f, iteration limit = %d, precision = %d size = %d\n"
                                               ,th_percentage,iter,prec,asize);
  for (i=0, a=1; (i<iter)&&(a>0); i+=10)
  {
    if (!(i%100))
      if(mess_flag) fprintf(stderr,"\n");
    for (j=0; j++ <10;)
      a=adjust();
    if(mess_flag) fprintf(stderr,"%ld ", a);if(bflag) putchar(BEEP);
  }
  if(mess_flag) fprintf(stderr,"total iterations: %d\n",iter);
}