#include <stddef.h>
#include <imapio.h>
#include <libimg.h>

/*-----------------------------------------------------------------
     Hurst operator

     1998.11.04
     filename: hurst.lc
     author  : W.J.Bokhove
----------------------------------------------------------------*/

/*
M_DOCDEF	hurst

NAME   hurst

SYNOPSIS
    #include <libimg.h>

       void hurst(src, dst, lines)
       sep uchar src[];   /* source image *//*
       sep uchar dst[];   /* destination image *//*
       uint lines;        /* number of lines to work on *//*

MEMORY CONSUMPTION
   This function consumes about 274 lines of IMAP memory

DESCRIPTION
   Performs the Hurst operator. Calculates in each area the difference
   between the brightest and darkest pixel. Then it fits a line through
   the log(distance from central pixel to area) and log(range). The slope
   of this line is used to represent the roughness of the area.

       h g h
     f e d e f
   h e c b c e h
   g d b a b d g
   h e c b c e h
     f e d e f
       h g h

M_DOCEND
*/

void hurst(sep uchar src[], sep uchar dest[], uint lines)
{
  //     h g h
  //   f e d e f
  // h e c b c e h
  // g d b a b d g
  // h e c b c e h
  //   f e d e f
  //     h g h

  #define swapsep(a,b,t)  mif((a)>(b)) {t=b;b=a;a=t;}

sep uchar ln_lut =:(0,32,51,64,74,82,90,96,101,106,110,114,118,121,125,128,130,133,135,138,140,142,144,146,148,150,152,153
,155,156,158,159,161,162,164,165,166,167,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,183,184,185,186
,187,188,188,189,190,191,191,192,193,193,194,195,195,196,197,197,198,199,199,200,200,201,202,202,203,203,204,204,205
,205,206,206,207,207,208,208,209,209,210,210,211,211,212,212,213,213,214,214,215,215,215,216,216,217,217,217,218,218
,219,219,219,220,220,221,221,221,222,222,222,223,223,224,224,224,225,225,225,226,226,226,227,227,227,228,228,228,229
,229,229,230,230,230,230,231,231,231,232,232,232,233,233,233,233,234,234,234,235,235,235,235,236,236,236,237,237,237
,237,238,238,238,238,239,239,239,239,240,240,240,240,241,241,241,241,242,242,242,242,243,243,243,243,243,244,244,244
,244,245,245,245,245,246,246,246,246,246,247,247,247,247,247,248,248,248,248,249,249,249,249,249,250,250,250,250,250
,251,251,251,251,251,252,252,252,252,252,252,253,253,253,253,253,254,254,254,254,254,255,255,255,255:);

  int i, j;
  sep uchar wrk[256];
  sep uchar max1, min1, max2, min2, max3, min3, tmp;
  sep uchar mintmp1, maxtmp1, mintmp2, maxtmp2;
  sep uint sumdif;
  sep ulong sumdifdst;
  
  for (i=1;i<256;i++) {
   wrk[i] = ln_lut:[i:];
  }

  dest[0] = 0;
  dest[1] = 0;
  dest[2] = 0;
  for (i=3;i<lines-3;i++) {
    /* 7 groups */
    max1 = src[i-1];          min1 = src[i+1];
    swapsep(min1, max1, tmp);
    max2 = src[i-2];          min2 = src[i+2];
    swapsep(min2, max2, tmp);
    max3 = src[i-3];          min3 = src[i+3];
    swapsep(min3, max3, tmp);

  // group b/1
    maxtmp1 = :<src[i];                 mintmp1 = :>src[i];
    swapsep(mintmp1, maxtmp1, tmp);
    maxtmp1 = max(maxtmp1, max1);       mintmp1 = min(mintmp1, min1);
    tmp = (maxtmp1 - mintmp1);
    sumdif = wrk[tmp];
  // group c/2
    maxtmp1 = max(:>max1, :<max1);      mintmp1 = min(:>min1, :<min1);
    tmp = (maxtmp1 - mintmp1);
    sumdif += wrk[tmp];
    sumdifdst = (wrk[tmp]<<4);
  // group d/3
    maxtmp1 = (src[i]:>2);              mintmp1 = (src[i]:<2);
    swapsep(mintmp1, maxtmp1, tmp);
    maxtmp1 = max(maxtmp1, max2);       mintmp1 = min(mintmp1, min2);
    tmp = (maxtmp1 - mintmp1);
    sumdif += wrk[tmp];
    sumdifdst += (wrk[tmp]<<5);
  // group e/4
    maxtmp1 = max(:>max2, :<max2);      mintmp1 = min(:>min2, :<min2);
    maxtmp2 = max((max1:>2), (max1:<2));
    mintmp2 = min((min1:>2), (min1:<2));
    maxtmp1 = max(maxtmp1, maxtmp2);    mintmp1 = min(mintmp1, mintmp2);
    tmp = (maxtmp1 - mintmp1);
    sumdif += wrk[tmp];
    sumdifdst += (wrk[tmp]*37);
  // group f/5
    maxtmp1 = max((max2:>2), (max2:<2));
    mintmp1 = min((min2:>2), (min2:<2));
    tmp = (maxtmp1 - mintmp1);
    sumdif += wrk[tmp];
    sumdifdst += (wrk[tmp]*48);
  // group g/6
    maxtmp1 = (src[i]:>3);              mintmp1 = (src[i]:<3);
    swapsep (mintmp1, maxtmp1, tmp);
    maxtmp1 = max(maxtmp1, max3);       mintmp1 = min(mintmp1, min3);
    tmp = (maxtmp1 - mintmp1);
    sumdif += wrk[tmp];
    sumdifdst += (wrk[tmp]*51);
  // group h/7
    maxtmp1 = max(:>max3, :<max3);      mintmp1 = min(:>min3, :<min3);
    maxtmp2 = max((max1:>3), (max1:<3));
    mintmp2 = min((min1:>3), (min1:<3));
    maxtmp1 = max(maxtmp1, maxtmp2);    mintmp1 = min(mintmp1, mintmp2);
    tmp = (maxtmp1 - mintmp1);
    sumdif += wrk[tmp];
    sumdifdst += (wrk[tmp]*53);

    dest[i] = ((7*sumdifdst) - (236*sumdif)) / 16313;
  }
  dest[lines-3] = 0;
  dest[lines-2] = 0;
  dest[lines-1] = 0;
}