/** \file
 * \brief Rank Convolution Operations
 *
 * See Copyright Notice in im_lib.h
 * $Id: Exp $
 */


#include <im.h>
#include <im_util.h>
#include <im_counter.h>

#include "im_process.h"

#include <stdlib.h>
#include <stdio.h>
#include <memory.h>
#include <string.h>
#include <math.h>


template <class T, class DT> 
static int DoConvolveRankFunc(T *map, DT* new_map, int width, int height, int kw, int kh, T (*func)(T* value, int ks), int counter)
{
  T* value = new T[kw*kh];
  int offset, new_offset, i, j, x, y, v;
  int kh2, kw2, ks;

  kh2 = kh/2;
  kw2 = kw/2;
  ks = kh*kw;

  for(j = 0; j < height; j++)
  {
    new_offset = j * width;

    for(i = 0; i < width; i++)
    {
      v = 0;
    
      for(y = -kh2; y <= kh2; y++)
      {
        if (j + y < 0)             // pass the bottom border
          offset = -(y + j + 1) * width;
        else if (j + y >= height)  // pass the top border
          offset = (2*height - 1 - (j + y)) * width;
        else
          offset = (j + y) * width;

        for(x = -kw2; x <= kw2; x++)
        {
          if (i + x < 0)            // pass the left border
            value[v] = map[offset - (i + x + 1)];
          else if (i + x >= width)  // pass the right border
            value[v] = map[offset + 2*width - 1 - (i + x)];
          else
            value[v] = map[offset + (i + x)];

          v++;
        }
      }
      
      new_map[new_offset + i] = (DT)func(value, ks);
    }    

    if (!imCounterInc(counter))
    {
      delete[] value;
      return 0;
    }
  }

  delete[] value;
  return 1;
}

static int compare_imReal(const void *elem1, const void *elem2) 
{
  float* v1 = (float*)elem1;
  float* v2 = (float*)elem2;

  if (*v1 < *v2)
    return -1;

  if (*v1 > *v2)
    return 1;

  return 0;
}

static int compare_imInt(const void *elem1, const void *elem2) 
{
  int* v1 = (int*)elem1;
  int* v2 = (int*)elem2;

  if (*v1 < *v2)
    return -1;

  if (*v1 > *v2)
    return 1;

  return 0;
}

static int compare_imUShort(const void *elem1, const void *elem2) 
{
  imushort* v1 = (imushort*)elem1;
  imushort* v2 = (imushort*)elem2;

  if (*v1 < *v2)
    return -1;

  if (*v1 > *v2)
    return 1;

  return 0;
}

static int compare_imByte(const void *elem1, const void *elem2) 
{
  imbyte* v1 = (imbyte*)elem1;
  imbyte* v2 = (imbyte*)elem2;

  if (*v1 < *v2)
    return -1;

  if (*v1 > *v2)
    return 1;

  return 0;
}

static imbyte median_op_byte(imbyte* value, int ks)
{
  qsort(value, ks, sizeof(imbyte), compare_imByte);
  return value[ks/2];
}

static imushort median_op_ushort(imushort* value, int ks)
{
  qsort(value, ks, sizeof(imushort), compare_imUShort);
  return value[ks/2];
}

static int median_op_int(int* value, int ks)
{
  qsort(value, ks, sizeof(int), compare_imInt);
  return value[ks/2];
}

static float median_op_real(float* value, int ks)
{
  qsort(value, ks, sizeof(float), compare_imReal);
  return value[ks/2];
}

int imProcessMedianConvolve(const imImage* src_image, imImage* dst_image, int ks)
{
  int i, ret = 0;
  int counter;

  counter = imCounterBegin("Median Filter");
  imCounterTotal(counter, src_image->depth*src_image->height, "Filtering...");

  for (i = 0; i < src_image->depth; i++)
  {
    switch(src_image->data_type)
    {
    case IM_BYTE:
      ret = DoConvolveRankFunc((imbyte*)src_image->data[i], (imbyte*)dst_image->data[i], 
                               src_image->width, src_image->height, ks, ks, median_op_byte, counter);
      break;                                                                                
    case IM_USHORT:                                                                           
      ret = DoConvolveRankFunc((imushort*)src_image->data[i], (imushort*)dst_image->data[i], 
                               src_image->width, src_image->height, ks, ks, median_op_ushort, counter);
      break;                                                                                
    case IM_INT:                                                                           
      ret = DoConvolveRankFunc((int*)src_image->data[i], (int*)dst_image->data[i], 
                               src_image->width, src_image->height, ks, ks, median_op_int, counter);
      break;                                                                                
    case IM_FLOAT:                                                                           
      ret = DoConvolveRankFunc((float*)src_image->data[i], (float*)dst_image->data[i], 
                               src_image->width, src_image->height, ks, ks, median_op_real, counter);
      break;                                                                                
    }
    
    if (!ret) 
      break;
  }

  imCounterEnd(counter);

  return ret;
}

template <class T> 
static inline void local_minmax(T* value, int ks, T &min, T &max)
{
  max=*value;
  min=*value;
  *value++;

  for (int i = 1; i < ks; i++)
  {
    if (*value < min)
      min = *value;

    if (*value > max)
      max = *value;

    value++;
  }
}

static imbyte range_op_byte(imbyte* value, int ks)
{
  imbyte min, max;
  local_minmax(value, ks, min, max);
  return max-min;
}

static imushort range_op_ushort(imushort* value, int ks)
{
  imushort min, max;
  local_minmax(value, ks, min, max);
  return max-min;
}

static int range_op_int(int* value, int ks)
{
  int min, max;
  local_minmax(value, ks, min, max);
  return max-min;
}

static float range_op_real(float* value, int ks)
{
  float min, max;
  local_minmax(value, ks, min, max);
  return max-min;
}

int imProcessRangeConvolve(const imImage* src_image, imImage* dst_image, int ks)
{
  int i, ret = 0;
  int counter;

  counter = imCounterBegin("Range Filter");
  imCounterTotal(counter, src_image->depth*src_image->height, "Filtering...");

  for (i = 0; i < src_image->depth; i++)
  {
    switch(src_image->data_type)
    {
    case IM_BYTE:
      ret = DoConvolveRankFunc((imbyte*)src_image->data[i], (imbyte*)dst_image->data[i], 
                               src_image->width, src_image->height, ks, ks, range_op_byte, counter);
      break;                                                                                
    case IM_USHORT:                                                                           
      ret = DoConvolveRankFunc((imushort*)src_image->data[i], (imushort*)dst_image->data[i], 
                               src_image->width, src_image->height, ks, ks, range_op_ushort, counter);
      break;                                                                                
    case IM_INT:                                                                           
      ret = DoConvolveRankFunc((int*)src_image->data[i], (int*)dst_image->data[i], 
                               src_image->width, src_image->height, ks, ks, range_op_int, counter);
      break;                                                                                
    case IM_FLOAT:                                                                           
      ret = DoConvolveRankFunc((float*)src_image->data[i], (float*)dst_image->data[i], 
                               src_image->width, src_image->height, ks, ks, range_op_real, counter);
      break;                                                                                
    }

    if (!ret) 
      break;
  }

  imCounterEnd(counter);

  return ret;
}

/*
Local variable threshold by the method of Bernsen.

Description:	
    If the difference between the largest and the smallest
		pixel value within the 'dx'*'dy' window is greater than
		or equal to 'cmin' (local contrast threshold), the average
		of the two values is used as threshold.

		Pixels in homogenous areas (difference below 'cmin') 
    are assumed to be below the threshold.

Reference:	
    Bernsen, J: "Dynamic thresholding of grey-level images"
		Proc. of the 8th ICPR, Paris, Oct 1986, 1251-1255.

Author:         Oivind Due Trier

Copyright 1990, Blab, UiO
Image processing lab, Department of Informatics
University of Oslo
*/

static int thresAux = 0;

static imbyte contrast_thres_op_byte(imbyte* value, int ks)
{
  int c, t;
  imbyte v = value[ks/2], min, max;

  local_minmax(value, ks, min, max);

  c = max-min;

  if (c < thresAux) 
    return 0;
  else
  { 
    t = ((int)max + (int)min) / 2;

    if (v >= t)
      return 1;
    else
      return 0;
  }
}

static imushort contrast_thres_op_ushort(imushort* value, int ks)
{
  int c, t;
  imushort v = value[ks/2], min, max;

  local_minmax(value, ks, min, max);

  c = max-min;

  if (c < thresAux) 
    return 0;
  else
  { 
    t = ((int)max + (int)min) / 2;

    if (v >= t)
      return 1;
    else
      return 0;
  }
}

static int contrast_thres_op_int(int* value, int ks)
{
  int c, t;
  int v = value[ks/2], min, max;

  local_minmax(value, ks, min, max);

  c = max-min;

  if (c < thresAux) 
    return 0;
  else
  { 
    t = ((int)max + (int)min) / 2;

    if (v >= t)
      return 1;
    else
      return 0;
  }
}

int imProcessRangeContrastThreshold(const imImage* src_image, imImage* dst_image, int ks, int min_range)
{
  int ret = 0;
  int counter = imCounterBegin("Range Contrast Threshold");
  imCounterTotal(counter, src_image->depth*src_image->height, "Filtering...");

  thresAux = min_range;

  switch(src_image->data_type)
  {
  case IM_BYTE:
    ret = DoConvolveRankFunc((imbyte*)src_image->data[0], (imbyte*)dst_image->data[0], 
                             src_image->width, src_image->height, ks, ks, contrast_thres_op_byte, counter);
    break;                                                                                
  case IM_USHORT:                                                                           
    ret = DoConvolveRankFunc((imushort*)src_image->data[0], (imbyte*)dst_image->data[0], 
                             src_image->width, src_image->height, ks, ks, contrast_thres_op_ushort, counter);
    break;                                                                                
  case IM_INT:                                                                           
    ret = DoConvolveRankFunc((int*)src_image->data[0], (imbyte*)dst_image->data[0], 
                             src_image->width, src_image->height, ks, ks, contrast_thres_op_int, counter);
    break;                                                                                
  }

  imCounterEnd(counter);

  return ret;
}

static imbyte max_thres_op_byte(imbyte* value, int ks)
{
  imbyte v = value[ks/2], min, max;

  if (v < thresAux) 
    return 0;

  local_minmax(value, ks, min, max);

  if (v < max)
    return 0;

  return 1;
}

static imushort max_thres_op_ushort(imushort* value, int ks)
{
  imushort v = value[ks/2], min, max;

  if (v < thresAux) 
    return 0;

  local_minmax(value, ks, min, max);

  if (v < max)
    return 0;

  return 1;
}

static int max_thres_op_int(int* value, int ks)
{
  int v = value[ks/2], min, max;

  if (v < thresAux) 
    return 0;

  local_minmax(value, ks, min, max);

  if (v < max)
    return 0;

  return 1;
}

int imProcessLocalMaxThreshold(const imImage* src_image, imImage* dst_image, int ks, int min_thres)
{
  int ret = 0;
  int counter = imCounterBegin("Local Max Threshold");
  imCounterTotal(counter, src_image->depth*src_image->height, "Filtering...");

  thresAux = min_thres;

  switch(src_image->data_type)
  {
  case IM_BYTE:
    ret = DoConvolveRankFunc((imbyte*)src_image->data[0], (imbyte*)dst_image->data[0], 
                             src_image->width, src_image->height, ks, ks, max_thres_op_byte, counter);
    break;                                                                                
  case IM_USHORT:                                                                           
    ret = DoConvolveRankFunc((imushort*)src_image->data[0], (imbyte*)dst_image->data[0], 
                             src_image->width, src_image->height, ks, ks, max_thres_op_ushort, counter);
    break;                                                                                
  case IM_INT:                                                                           
    ret = DoConvolveRankFunc((int*)src_image->data[0], (imbyte*)dst_image->data[0], 
                             src_image->width, src_image->height, ks, ks, max_thres_op_int, counter);
    break;                                                                                
  }

  imCounterEnd(counter);

  return ret;
}

static imbyte rank_closest_op_byte(imbyte* value, int ks)
{
  imbyte v = value[ks/2];
  imbyte min, max;

  local_minmax(value, ks, min, max);

  if (v - min < max - v) 
    return min;
  else
    return max;
}

static imushort rank_closest_op_ushort(imushort* value, int ks)
{
  imushort v = value[ks/2];
  imushort min, max;

  local_minmax(value, ks, min, max);

  if (v - min < max - v) 
    return min;
  else
    return max;
}

static int rank_closest_op_int(int* value, int ks)
{
  int v = value[ks/2];
  int min, max;

  local_minmax(value, ks, min, max);

  if (v - min < max - v) 
    return min;
  else
    return max;
}

static float rank_closest_op_real(float* value, int ks)
{
  float v = value[ks/2];
  float min, max;

  local_minmax(value, ks, min, max);

  if (v - min < max - v) 
    return min;
  else
    return max;
}


int imProcessRankClosestConvolve(const imImage* src_image, imImage* dst_image, int ks)
{
  int i, ret = 0;
  int counter;

  counter = imCounterBegin("Rank Closest");
  imCounterTotal(counter, src_image->depth*src_image->height, "Filtering...");

  for (i = 0; i < src_image->depth; i++)
  {
    switch(src_image->data_type)
    {
    case IM_BYTE:
      ret = DoConvolveRankFunc((imbyte*)src_image->data[i], (imbyte*)dst_image->data[i], 
                               src_image->width, src_image->height, ks, ks, rank_closest_op_byte, counter);
      break;                                                                                
    case IM_USHORT:                                                                           
      ret = DoConvolveRankFunc((imushort*)src_image->data[i], (imushort*)dst_image->data[i], 
                               src_image->width, src_image->height, ks, ks, rank_closest_op_ushort, counter);
      break;                                                                                
    case IM_INT:                                                                           
      ret = DoConvolveRankFunc((int*)src_image->data[i], (int*)dst_image->data[i], 
                               src_image->width, src_image->height, ks, ks, rank_closest_op_int, counter);
      break;                                                                                
    case IM_FLOAT:                                                                           
      ret = DoConvolveRankFunc((float*)src_image->data[i], (float*)dst_image->data[i], 
                               src_image->width, src_image->height, ks, ks, rank_closest_op_real, counter);
      break;                                                                                
    }
    
    if (!ret) 
      break;
  }

  imCounterEnd(counter);

  return ret;
}

static imbyte rank_max_op_byte(imbyte* value, int ks)
{
  imbyte min, max;
  local_minmax(value, ks, min, max);
  return max;
}

static imushort rank_max_op_ushort(imushort* value, int ks)
{
  imushort min, max;
  local_minmax(value, ks, min, max);
  return max;
}

static int rank_max_op_int(int* value, int ks)
{
  int min, max;
  local_minmax(value, ks, min, max);
  return max;
}

static float rank_max_op_real(float* value, int ks)
{
  float min, max;
  local_minmax(value, ks, min, max);
  return max;
}

int imProcessRankMaxConvolve(const imImage* src_image, imImage* dst_image, int ks)
{
  int i, ret = 0;
  int counter;

  counter = imCounterBegin("Rank Max");
  imCounterTotal(counter, src_image->depth*src_image->height, "Filtering...");

  for (i = 0; i < src_image->depth; i++)
  {
    switch(src_image->data_type)
    {
    case IM_BYTE:
      ret = DoConvolveRankFunc((imbyte*)src_image->data[i], (imbyte*)dst_image->data[i], 
                               src_image->width, src_image->height, ks, ks, rank_max_op_byte, counter);
      break;                                                                                
    case IM_USHORT:                                                                           
      ret = DoConvolveRankFunc((imushort*)src_image->data[i], (imushort*)dst_image->data[i], 
                               src_image->width, src_image->height, ks, ks, rank_max_op_ushort, counter);
      break;                                                                                
    case IM_INT:                                                                           
      ret = DoConvolveRankFunc((int*)src_image->data[i], (int*)dst_image->data[i], 
                               src_image->width, src_image->height, ks, ks, rank_max_op_int, counter);
      break;                                                                                
    case IM_FLOAT:                                                                           
      ret = DoConvolveRankFunc((float*)src_image->data[i], (float*)dst_image->data[i], 
                               src_image->width, src_image->height, ks, ks, rank_max_op_real, counter);
      break;                                                                                
    }
    
    if (!ret) 
      break;
  }

  imCounterEnd(counter);

  return ret;
}

static imbyte rank_min_op_byte(imbyte* value, int ks)
{
  imbyte min, max;
  local_minmax(value, ks, min, max);
  return min;
}

static imushort rank_min_op_ushort(imushort* value, int ks)
{
  imushort min, max;
  local_minmax(value, ks, min, max);
  return min;
}

static int rank_min_op_int(int* value, int ks)
{
  int min, max;
  local_minmax(value, ks, min, max);
  return min;
}

static float rank_min_op_real(float* value, int ks)
{
  float min, max;
  local_minmax(value, ks, min, max);
  return min;
}

int imProcessRankMinConvolve(const imImage* src_image, imImage* dst_image, int ks)
{
  int i, ret = 0;
  int counter;

  counter = imCounterBegin("Rank Min");
  imCounterTotal(counter, src_image->depth*src_image->height, "Filtering...");

  for (i = 0; i < src_image->depth; i++)
  {
    switch(src_image->data_type)
    {
    case IM_BYTE:
      ret = DoConvolveRankFunc((imbyte*)src_image->data[i], (imbyte*)dst_image->data[i], 
                               src_image->width, src_image->height, ks, ks, rank_min_op_byte, counter);
      break;                                                                                
    case IM_USHORT:                                                                           
      ret = DoConvolveRankFunc((imushort*)src_image->data[i], (imushort*)dst_image->data[i], 
                               src_image->width, src_image->height, ks, ks, rank_min_op_ushort, counter);
      break;                                                                                
    case IM_INT:                                                                           
      ret = DoConvolveRankFunc((int*)src_image->data[i], (int*)dst_image->data[i], 
                               src_image->width, src_image->height, ks, ks, rank_min_op_int, counter);
      break;                                                                                
    case IM_FLOAT:                                                                           
      ret = DoConvolveRankFunc((float*)src_image->data[i], (float*)dst_image->data[i], 
                               src_image->width, src_image->height, ks, ks, rank_min_op_real, counter);
      break;                                                                                
    }
    
    if (!ret) 
      break;
  }

  imCounterEnd(counter);

  return ret;
}
