/*
Description: This program is a plugin for the Blender sequence editor.
It applies an unsharp mask to one strip.
Author: Andrew McCargar (andrew@outsideworld.org)
Website: http://www.outsideworld.org/projects/blender_plugins.html
code used from the Gimp unsharp filter:
 *  http://www.steppe.com/~winston/gimp/unsharp.html
 *
 * Copyright (C) 1999 Winston Chang
 *                    <wchang3@students.wisc.edu>
 *                    <winston@steppe.com>
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.

Last Modified: Thus Dec 23 23:41:35 EST 2009
*/

#include "math.h"
#include "plugin.h"

char name[24]= "Unsharp mask";

/* structure for buttons,
 * butcode name default min max 0
 */

VarStruct varstr[]= {
   LABEL, 	"Unsharp mask", 0.0, 	0.0, 	0.0, 	"",
   NUMSLI|FLO,	"Radius ",		5.0,	0.0,	20.0, "Size of the edges to be enhanced", 
   NUMSLI|FLO,	"Amount ",		0.5,	0.0,	1.0, "magnatude of edge contrast", 
   NUMSLI|INT,	"Threshold ",		0.0,	0.0,	255.0, "minimum brightness change that will be sharpened", 
   TOG|INT,	"Y filter",		1.0,	0.0,	1.0, "filter along y as well as x axis (slower but better)", 
};


/* The cast struct is for input in the main doit function
 Varstr and Cast must have the same variables in the same order */

typedef struct Cast {
 int dummy; /* because of the 'label' button */
 float radius;
 float amount;
 int threshold;
 int yfilter;
} Cast;

/* cfra: the current frame */

float cfra;

void plugin_seq_doit(Cast *, float, float, int, int,
 ImBuf *, ImBuf *, ImBuf *, ImBuf *);

int plugin_seq_getversion(void) {
 return B_PLUGIN_VERSION;
}

void plugin_but_changed(int but) {
}

void plugin_init() {
}

void plugin_getinfo(PluginInfo *info) {
 info->name= name;
 info->nvars= sizeof(varstr)/sizeof(VarStruct);
 info->cfra= &cfra;

 info->varstr= varstr;

 info->init= plugin_init;
 info->seq_doit= (SeqDoit) plugin_seq_doit;
 info->callback= plugin_but_changed;
}

/* local function prototypes */
static inline void blur_line
  (double* ctable,
   double* cmatrix,
   int cmatrix_length,
   char* cur_col,
   char* dest_col,
   int y,
   long bytes);
static double* gen_lookup_table(double* cmatrix, int cmatrix_length);
static int gen_convolve_matrix(double radius, double** cmatrix_p);
static inline int round2int(double i);
void merge_bufs(char *source, char *dest, float amount, int threshold, int height, int width);
/* end of local function prototypes */

void plugin_seq_doit(Cast *cast, float facf0, float facf1, int width, int height,
 ImBuf *ibuf1, ImBuf *ibuf2, ImBuf *outbuf, ImBuf *use) {
 char *in1= (char *)ibuf1->rect;
 char *out= (char *)outbuf->rect;
 int length = (width * 4);
 char *cur_row = (char *) malloc (length);
 char *dest_row = (char *) malloc (length);
 char *cur_col = (char *) malloc (height);
 char *dest_col = (char *) malloc (height);
 int x, y; /* counter variables */

 double *cmatrix = NULL;
 int cmatrix_length;
 double *ctable;

 /* generate convolution matrix and make sure it's smaller than each dimension */
 cmatrix_length = gen_convolve_matrix(cast->radius, &cmatrix);
 /* generate lookup table */
 ctable = gen_lookup_table(cmatrix, cmatrix_length);

 /* blur the rows */
 for (y=0;y<height;y++) {
   for (x=0;x<length;x++) {
     cur_row[x] = in1[y*length+x];
   }
   blur_line(ctable, cmatrix, cmatrix_length, cur_row, dest_row, width, 4);
   for (x=0;x<length;x++) {
     out[y*length+x] = dest_row[x];
   }
 }
 /* then blur the collums (optional)*/
 if (cast->yfilter) {
   for (x=0;x<length;x++) {
     for (y=0;y<height;y++) {
       cur_col[y] = out[y*length+x];
     }
     blur_line(ctable, cmatrix, cmatrix_length, cur_col, dest_col, height, 1);
     for (y=0;y<height;y++) {
       out[y*length+x] = dest_col[y];
     }
   }
 } 

 /* Finally merge the two buffers*/
 merge_bufs((char *)ibuf1->rect, (char *)outbuf->rect, (cast->amount * facf0), cast->threshold, height, width);

 /* free the memory we took */
 free(cur_row);
 free(dest_row);
 free(cur_col);
 free(dest_col);
 free(cmatrix);
 free(ctable);
} 

static inline void blur_line
  (double* ctable,
   double* cmatrix,
   int cmatrix_length,
   char* cur_col,
   char* dest_col,
   int y,
   long bytes)
{
  double scale;
  double sum;
  int i=0, j=0;
  int row;
  int cmatrix_middle = cmatrix_length/2;

  double *cmatrix_p;
  char  *cur_col_p;
  char  *cur_col_p1;
  char  *dest_col_p;
  double *ctable_p;


  /* this first block is the same as the non-optimized version --
   * it is only used for very small pictures, so speed isn't a
   * big concern.
   */
  if (cmatrix_length > y) {

    for (row = 0; row < y ; row++) {
      scale=0;
      /* find the scale factor */
      for (j = 0; j < y ; j++) {
        /* if the index is in bounds, add it to the scale counter */
        if ((j + cmatrix_length/2 - row >= 0) &&
            (j + cmatrix_length/2 - row < cmatrix_length))
          scale += cmatrix[j + cmatrix_length/2 - row];
      }
      for (i = 0; i<bytes; i++) {
        sum = 0;
        for (j = 0; j < y; j++) {
          if ( (j >= row - cmatrix_length/2) && (j <= row + cmatrix_length/2) )
            sum += cur_col[j*bytes + i] * cmatrix[j];
        }
        dest_col[row*bytes + i] = (char)round2int(sum / scale);
      }
    }
  }

  else {
    /* for the edge condition, we only use available info and scale to one */
    for (row = 0; row < cmatrix_middle; row++) {
      /* find scale factor */
      scale=0;
      for (j = cmatrix_middle - row; j<cmatrix_length; j++)
        scale += cmatrix[j];
      for (i = 0; i<bytes; i++) {
        sum = 0;
        for (j = cmatrix_middle - row; j<cmatrix_length; j++) {
          sum += cur_col[(row + j-cmatrix_middle)*bytes + i] * cmatrix[j];
        }
        dest_col[row*bytes + i] = (char)round2int(sum / scale);
      }
    }
    /* go through each pixel in each col */
    dest_col_p = dest_col + row*bytes;
    for ( ; row < y-cmatrix_middle; row++) {
      cur_col_p = (row - cmatrix_middle) * bytes + cur_col;
      for (i = 0; i<bytes; i++) {
        sum = 0;
        cmatrix_p = cmatrix;
        cur_col_p1 = cur_col_p;
        ctable_p = ctable;
        for (j = cmatrix_length; j>0; j--) {
          sum += *(ctable_p + *cur_col_p1);
          cur_col_p1 += bytes;
          ctable_p += 256;
        }
        cur_col_p++;
        *(dest_col_p++) = round2int(sum);
      }
    }

    /* for the edge condition , we only use available info, and scale to one */
    for ( ; row < y; row++) {
      /* find scale factor */
      scale=0;
      for (j = 0; j< y-row + cmatrix_middle; j++)
        scale += cmatrix[j];
      for (i = 0; i<bytes; i++) {
        sum = 0;
        for (j = 0; j<y-row + cmatrix_middle; j++) {
          sum += cur_col[(row + j-cmatrix_middle)*bytes + i] * cmatrix[j];
        }
        dest_col[row*bytes + i] = (char)round2int(sum / scale);
      }
    }
  }
}

static int gen_convolve_matrix(double radius, double** cmatrix_p) {
  int matrix_length;
  int matrix_midpoint;
  double* cmatrix;
  int i,j;
  double std_dev;
  double sum;

  /* we want to generate a matrix that goes out a certain radius
   * from the center, so we have to go out ceil(rad-0.5) pixels,
   * inlcuding the center pixel.  Of course, that's only in one direction,
   * so we have to go the same amount in the other direction, but not count
   * the center pixel again.  So we double the previous result and subtract
   * one.
   * The radius parameter that is passed to this function is used as
   * the standard deviation, and the radius of effect is the
   * standard deviation * 2.  It's a little confusing.
   */
  radius = fabs(radius) + 1.0;

  std_dev = radius;
  radius = std_dev * 2;

  /* go out 'radius' in each direction */
  matrix_length = 2 * ceil(radius-0.5) + 1;
  if (matrix_length <= 0) matrix_length = 1;
  matrix_midpoint = matrix_length/2 + 1;
  *cmatrix_p = (double*)(malloc(matrix_length * sizeof(double)));
  cmatrix = *cmatrix_p;

  /*  Now we fill the matrix by doing a numeric integration approximation
   * from -2*std_dev to 2*std_dev, sampling 50 points per pixel.
   * We do the bottom half, mirror it to the top half, then compute the
   * center point, which has 51 sample points.
   *  Otherwise asymmetric quantization errors will occur.
   *  The formula to integrate is e^-(x^2/2s^2).
   */

  /* first we do the top (right) half of matrix */
  for (i=matrix_length/2 + 1; i<matrix_length; i++) {
    double base_x = i - floor(matrix_length/2) - 0.5;
    sum = 0;
    for (j=1; j<=50; j++) {
      if ( base_x+0.02*j <= radius )
        sum += exp( -(base_x+0.02*j)*(base_x+0.02*j) /
                         (2*std_dev*std_dev) );
    }
    cmatrix[i] = sum/50;
  }

  /* mirror the thing to the bottom half */
  for (i=0; i<=matrix_length/2; i++) {
    cmatrix[i] = cmatrix[matrix_length-1-i];
  }

  /* find center val -- calculate an odd number of quanta to make it symmetric,
   * even if the center point is weighted slightly higher than others. */
  sum = 0;
  for (j=0; j<=50; j++) {
    sum += exp( -(0.5+0.02*j)*(0.5+0.02*j) /
                     (2*std_dev*std_dev) );
  }
  cmatrix[matrix_length/2] = sum/51;


  /* normalize the distribution by scaling the total sum to one */
  sum=0;
  for (i=0; i<matrix_length; i++) sum += cmatrix[i];
  for (i=0; i<matrix_length; i++) cmatrix[i] = cmatrix[i] / sum;

  return matrix_length;
}

static double* gen_lookup_table(double* cmatrix, int cmatrix_length) {
  int i, j;
  double* lookup_table = malloc(cmatrix_length * 256 * sizeof(double));

  double* lookup_table_p = lookup_table;
  double* cmatrix_p      = cmatrix;
  for (i=0; i<cmatrix_length; i++) {
    for (j=0; j<256; j++) {
      *(lookup_table_p++) = *cmatrix_p * (double)j;
    }
    cmatrix_p++;
  }

  return lookup_table;
}

static inline int round2int(double i) {
  return (int)(i + 0.5);
}

void merge_bufs(char *source, char *dest, float amount, int threshold, int height, int width) {
  /* merge the source and destination (which currently contains
     the blurred version) images */
  int value, diff;
  int x, y, z; /*counter variables*/

  for (y = 0; y < height; y++) {
    value = 0;
    /* combine the two rows*/
    for (x = 0; x < width; x++) {
      for (z = 0; z < 3; z++) {
        diff = (*source - *dest);
        /* do tresholding */
        if ( abs(2*diff) < threshold )
          diff = 0;

        value = *source++ + amount * diff;

        if (value < 0) *dest++ = 0;
        else if (value > 255) *dest++ = 255;
        else  *dest++ = value; 
      }
      *dest++ = *source++; /*copy over the alpha channel */
    }
  }
}