/** \file
 * \brief Binary Arithmetic Operations
 *
 * See Copyright Notice in im_lib.h
 * $Id: im_arithmetic_bin.cpp,v 1.3 2010/01/08 03:49:05 scuri Exp $
 */


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

#include "im_process_pon.h"
#include "im_math_op.h"

#include <stdlib.h>
#include <memory.h>


template <class T1, class T2, class T3> 
static void DoBinaryOp(T1 *map1, T2 *map2, T3 *map, int count, int op)
{
  int i;

  switch(op)
  {
  case IM_BIN_ADD:
    for (i = 0; i < count; i++)
      map[i] = add_op((T3)map1[i], (T3)map2[i]);
    break;
  case IM_BIN_SUB:
    for (i = 0; i < count; i++)
      map[i] = sub_op((T3)map1[i], (T3)map2[i]);
    break;
  case IM_BIN_MUL:
    for (i = 0; i < count; i++)
      map[i] = mul_op((T3)map1[i], (T3)map2[i]);
    break;
  case IM_BIN_DIV:
    for (i = 0; i < count; i++)
      map[i] = div_op((T3)map1[i], (T3)map2[i]);
    break;
  case IM_BIN_DIFF:
    for (i = 0; i < count; i++)
      map[i] = diff_op((T3)map1[i], (T3)map2[i]);
    break;
  case IM_BIN_MIN:
    for (i = 0; i < count; i++)
      map[i] = min_op((T3)map1[i], (T3)map2[i]);
    break;
  case IM_BIN_MAX:
    for (i = 0; i < count; i++)
      map[i] = max_op((T3)map1[i], (T3)map2[i]);
    break;
  case IM_BIN_POW:
    for (i = 0; i < count; i++)
      map[i] = pow_op((T3)map1[i], (T3)map2[i]);
    break;
  }
}

static void DoBinaryOpByte(imbyte *map1, imbyte *map2, imbyte *map, int count, int op)
{
  int i;

  switch(op)
  {
  case IM_BIN_ADD:
    for (i = 0; i < count; i++)
      map[i] = (imbyte)crop_byte(add_op((int)map1[i], (int)map2[i]));
    break;
  case IM_BIN_SUB:
    for (i = 0; i < count; i++)
      map[i] = (imbyte)crop_byte(sub_op((int)map1[i], (int)map2[i]));
    break;
  case IM_BIN_MUL:
    for (i = 0; i < count; i++)
      map[i] = (imbyte)crop_byte(mul_op((int)map1[i], (int)map2[i]));
    break;
  case IM_BIN_DIV:
    for (i = 0; i < count; i++)
      map[i] = (imbyte)crop_byte(div_op((int)map1[i], (int)map2[i]));
    break;
  case IM_BIN_DIFF:
    for (i = 0; i < count; i++)
      map[i] = (imbyte)crop_byte(diff_op((int)map1[i], (int)map2[i]));
    break;
  case IM_BIN_MIN:
    for (i = 0; i < count; i++)
      map[i] = (imbyte)crop_byte(min_op((int)map1[i], (int)map2[i]));
    break;
  case IM_BIN_MAX:
    for (i = 0; i < count; i++)
      map[i] = (imbyte)crop_byte(max_op((int)map1[i], (int)map2[i]));
    break;
  case IM_BIN_POW:
    for (i = 0; i < count; i++)
      map[i] = (imbyte)crop_byte(pow_op((int)map1[i], (int)map2[i]));
    break;
  }
}

static void DoBinaryOpCpxReal(imcfloat *map1, float *map2, imcfloat *map, int count, int op)
{
  int i;

  switch(op)
  {
  case IM_BIN_ADD:
    for (i = 0; i < count; i++)
      map[i] = add_op(map1[i], map2[i]);
    break;
  case IM_BIN_SUB:
    for (i = 0; i < count; i++)
      map[i] = sub_op(map1[i], map2[i]);
    break;
  case IM_BIN_MUL:
    for (i = 0; i < count; i++)
      map[i] = mul_op(map1[i], map2[i]);
    break;
  case IM_BIN_DIV:
    for (i = 0; i < count; i++)
      map[i] = div_op(map1[i], (imcfloat)map2[i]);
    break;
  case IM_BIN_DIFF:
    for (i = 0; i < count; i++)
      map[i] = diff_op(map1[i], map2[i]);
    break;
  case IM_BIN_MIN:
    for (i = 0; i < count; i++)
      map[i] = min_op(map1[i], map2[i]);
    break;
  case IM_BIN_MAX:
    for (i = 0; i < count; i++)
      map[i] = max_op(map1[i], map2[i]);
    break;
  case IM_BIN_POW:
    for (i = 0; i < count; i++)
      map[i] = pow_op(map1[i], map2[i]);
    break;
  }
}

void imProcessArithmeticOp(const imImage* src_image1, const imImage* src_image2, imImage* dst_image, int op)
{
  int count = src_image1->count;

  for (int i = 0; i < src_image1->depth; i++)
  {
    switch(src_image1->data_type)
    {
    case IM_BYTE:
      if (dst_image->data_type == IM_FLOAT)
        DoBinaryOp((imbyte*)src_image1->data[i], (imbyte*)src_image2->data[i], (float*)dst_image->data[i], count, op);
      else if (dst_image->data_type == IM_USHORT)
        DoBinaryOp((imbyte*)src_image1->data[i], (imbyte*)src_image2->data[i], (imushort*)dst_image->data[i], count, op);
      else if (dst_image->data_type == IM_INT)
        DoBinaryOp((imbyte*)src_image1->data[i], (imbyte*)src_image2->data[i], (int*)dst_image->data[i], count, op);
      else
        DoBinaryOpByte((imbyte*)src_image1->data[i], (imbyte*)src_image2->data[i], (imbyte*)dst_image->data[i], count, op);
      break;
    case IM_USHORT:
      if (dst_image->data_type == IM_FLOAT)
        DoBinaryOp((imushort*)src_image1->data[i], (imushort*)src_image2->data[i], (float*)dst_image->data[i], count, op);
      else if (dst_image->data_type == IM_INT)
        DoBinaryOp((imushort*)src_image1->data[i], (imushort*)src_image2->data[i], (int*)dst_image->data[i], count, op);
      else
        DoBinaryOp((imushort*)src_image1->data[i], (imushort*)src_image2->data[i], (imushort*)dst_image->data[i], count, op);
      break;
    case IM_INT:
      if (dst_image->data_type == IM_FLOAT)
        DoBinaryOp((int*)src_image1->data[i], (int*)src_image2->data[i], (float*)dst_image->data[i], count, op);
      else
        DoBinaryOp((int*)src_image1->data[i], (int*)src_image2->data[i], (int*)dst_image->data[i], count, op);
      break;
    case IM_FLOAT:
      DoBinaryOp((float*)src_image1->data[i], (float*)src_image2->data[i], (float*)dst_image->data[i], count, op);
      break;
    case IM_CFLOAT:
      if (src_image2->data_type == IM_FLOAT)
        DoBinaryOpCpxReal((imcfloat*)src_image1->data[i], (float*)src_image2->data[i], (imcfloat*)dst_image->data[i], count, op);
      else
        DoBinaryOp((imcfloat*)src_image1->data[i], (imcfloat*)src_image2->data[i], (imcfloat*)dst_image->data[i], count, op);
      break;
    }
  }
}

template <class T>
static inline T blend_op(const T& v1, const T& v2, const float& alpha)
{
  return (T)(alpha*v1 + (1.0f - alpha)*v2);
}

template <class T> 
static void DoBlendConst(T *map1, T *map2, T *map, int count, float alpha)
{
  for (int i = 0; i < count; i++)
    map[i] = blend_op(map1[i], map2[i], alpha);
}

void imProcessBlendConst(const imImage* src_image1, const imImage* src_image2, imImage* dst_image, float alpha)
{
  int count = src_image1->count;

  for (int i = 0; i < src_image1->depth; i++)
  {
    switch(src_image1->data_type)
    {
    case IM_BYTE:
      DoBlendConst((imbyte*)src_image1->data[i], (imbyte*)src_image2->data[i], (imbyte*)dst_image->data[i], count, alpha);
      break;
    case IM_USHORT:
      DoBlendConst((imushort*)src_image1->data[i], (imushort*)src_image2->data[i], (imushort*)dst_image->data[i], count, alpha);
      break;
    case IM_INT:
      DoBlendConst((int*)src_image1->data[i], (int*)src_image2->data[i], (int*)dst_image->data[i], count, alpha);
      break;
    case IM_FLOAT:
      DoBlendConst((float*)src_image1->data[i], (float*)src_image2->data[i], (float*)dst_image->data[i], count, alpha);
      break;
    case IM_CFLOAT:
      DoBlendConst((imcfloat*)src_image1->data[i], (imcfloat*)src_image2->data[i], (imcfloat*)dst_image->data[i], count, alpha);
      break;
    }
  }
}

template <class T, class TA> 
static void DoBlend(T *map1, T *map2, TA *alpha, T *map, int count, TA max)
{
  for (int i = 0; i < count; i++)
    map[i] = blend_op(map1[i], map2[i], ((float)alpha[i])/max);
}

void imProcessBlend(const imImage* src_image1, const imImage* src_image2, const imImage* alpha, imImage* dst_image)
{
  int count = src_image1->count;

  for (int i = 0; i < src_image1->depth; i++)
  {
    switch(src_image1->data_type)
    {
    case IM_BYTE:
      DoBlend((imbyte*)src_image1->data[i], (imbyte*)src_image2->data[i], (imbyte*)alpha->data[0], (imbyte*)dst_image->data[i], count, (imbyte)255);
      break;
    case IM_USHORT:
      DoBlend((imushort*)src_image1->data[i], (imushort*)src_image2->data[i], (imushort*)alpha->data[0], (imushort*)dst_image->data[i], count, (imushort)65535);
      break;
    case IM_INT:
      DoBlend((int*)src_image1->data[i], (int*)src_image2->data[i], (int*)alpha->data[0], (int*)dst_image->data[i], count, (int)2147483647);
      break;
    case IM_FLOAT:
      DoBlend((float*)src_image1->data[i], (float*)src_image2->data[i], (float*)alpha->data[0], (float*)dst_image->data[i], count, 1.0f);
      break;
    case IM_CFLOAT:
      DoBlend((imcfloat*)src_image1->data[i], (imcfloat*)src_image2->data[i], (float*)alpha->data[0], (imcfloat*)dst_image->data[i], count, 1.0f);
      break;
    }
  }
}

#define COMPOSE_OVER(_SRC, _SRC_ALPHA, _DST, _TMP_MULTI, _TMP_ALPHA) (T)(((_SRC_ALPHA)*(_SRC) + (_TMP_MULTI)*(_DST)) / (_TMP_ALPHA))
#define ALPHA_BLEND(_src,_dst,_alpha) (T)(((_src) * (_alpha) + (_dst) * (max - (_alpha))) / max)

template <class T, class TA>
static inline T compose_op(const T& v1, const T& v2, const T& alpha1, const T& alpha2, const TA& max)
{
  if (alpha1 != max)   /* some transparency */                                                                     
  {                                                                                                                    
    if (alpha1 != 0) /* source not full transparent */                                                             
    {                                                                                                                  
      if (alpha2 == 0) /* destiny full transparent */                                                            
      {                                                                                                                
        return v1;                                                                                           
      }                                                                                                                
      else if (alpha2 == max) /* destiny opaque */                                                               
      {                                                                                                                
        return ALPHA_BLEND(v1, v2, alpha1);                                                   
      }                                                                                                                
      else /* (0<alpha2<max && 0<alpha1<max) destiny and source are semi-transparent */                      
      {                                                                                                                
        /* Closed Compositing SRC over DST  (see smith95a.pdf)        */                                               
        /* Colors NOT Premultiplied by Alpha                          */                                               
        /* DST = SRC * SRC_ALPHA + DST * DST_ALPHA * (1 - SRC_ALPHA)  */                                               
        /* DST_ALPHA = SRC_ALPHA + DST_ALPHA * (1 - SRC_ALPHA)        */                                               
        /* DST /= DST_ALPHA */                                                                                         
        TA _tmp_multi = alpha2 * (max - alpha1);                                                            
        TA _tmp_src_alpha = alpha1*max;                                                                           
        TA _tmp_alpha = _tmp_src_alpha + _tmp_multi;                                                                  
        return COMPOSE_OVER(v1, _tmp_src_alpha, v2, _tmp_multi, _tmp_alpha);                     
      }                                                                                                                
    }                                                                                                                  
    else  /* (alpha1 == 0) source full transparent */                                                              
    {                                                                                                                  
      return v2;                                                                                           
    }                                                                                                                  
  }                                                                                                                    
  else  /* (alpha1 == max) source has no alpha = opaque */                                                         
  {                                                                                                                    
    return v1;                                                                                               
  }                                                                                                                    
}

template <class T, class TA>
static inline T compose_alpha_op(const T& alpha1, const T& alpha2, const TA& max)
{
  if (alpha1 != max)   /* some transparency */                                                                     
  {                                                                                                                    
    if (alpha1 != 0) /* source not full transparent */                                                             
    {                                                                                                                  
      if (alpha2 == 0) /* destiny full transparent */                                                            
      {                                                                                                                
        return alpha1;                                                                                     
      }                                                                                                                
      else if (alpha2 == max) /* destiny opaque */                                                               
      {                                                                                                                
        /* alpha2 is not changed */                                                                              
        return alpha2;
      }                                                                                                                
      else /* (0<alpha2<max && 0<alpha1<max) destiny and source are semi-transparent */                      
      {                                                                                                                
        /* Closed Compositing SRC over DST  (see smith95a.pdf)        */                                               
        /* Colors NOT Premultiplied by Alpha                          */                                               
        /* DST = SRC * SRC_ALPHA + DST * DST_ALPHA * (1 - SRC_ALPHA)  */                                               
        /* DST_ALPHA = SRC_ALPHA + DST_ALPHA * (1 - SRC_ALPHA)        */                                               
        /* DST /= DST_ALPHA */                                                                                         
        TA _tmp_multi = alpha2 * (max - alpha1);                                                            
        TA _tmp_src_alpha = alpha1*max;                                                                           
        TA _tmp_alpha = _tmp_src_alpha + _tmp_multi;                                                                  
        return (T)(_tmp_alpha / max);
      }                                                                                                                
    }                                                                                                                  
    else  /* (alpha1 == 0) source full transparent */                                                              
    {                                                                                                                  
      /* alpha2 is not changed */                                                                                
      return alpha2;
    }                                                                                                                  
  }                                                                                                                    
  else  /* (alpha1 == max) source has no alpha = opaque */                                                         
  {                                                                                                                    
    return (unsigned char)max;   /* set destiny as opaque */                                                   
  }                                                                                                                    
}

template <class T, class TA> 
static void DoCompose(T *map1, T *map2, T *alpha1, T *alpha2, T *map, int count, TA max)
{
  for (int i = 0; i < count; i++)
    map[i] = compose_op(map1[i], map2[i], alpha1[i], alpha2[i], max);
}

template <class T, class TA> 
static void DoComposeAlpha(T *alpha1, T *alpha2, T *dst_alpha, int count, TA max)
{
  for (int i = 0; i < count; i++)
    dst_alpha[i] = compose_alpha_op(alpha1[i], alpha2[i], max);
}

void imProcessCompose(const imImage* src_image1, const imImage* src_image2, imImage* dst_image)
{
  int count = src_image1->count, 
      src_alpha = src_image1->depth;

  if (!src_image1->has_alpha || !src_image2->has_alpha || !dst_image->has_alpha)
    return;

  for (int i = 0; i < src_image1->depth; i++)
  {
    switch(src_image1->data_type)
    {
    case IM_BYTE:
      DoCompose((imbyte*)src_image1->data[i], (imbyte*)src_image2->data[i], (imbyte*)src_image1->data[src_alpha], (imbyte*)src_image2->data[src_alpha], (imbyte*)dst_image->data[i], count, (int)255);
      break;
    case IM_USHORT:
      DoCompose((imushort*)src_image1->data[i], (imushort*)src_image2->data[i], (imushort*)src_image1->data[src_alpha], (imushort*)src_image2->data[src_alpha], (imushort*)dst_image->data[i], count, (int)65535);
      break;
    case IM_INT:
      DoCompose((int*)src_image1->data[i], (int*)src_image2->data[i], (int*)src_image1->data[src_alpha], (int*)src_image2->data[src_alpha], (int*)dst_image->data[i], count, (int)2147483647);
      break;
    case IM_FLOAT:
      DoCompose((float*)src_image1->data[i], (float*)src_image2->data[i], (float*)src_image1->data[src_alpha], (float*)src_image2->data[src_alpha], (float*)dst_image->data[i], count, 1.0f);
      break;
    }
  }

  /* one more for the alpha channel */
  switch(src_image1->data_type)
  {
  case IM_BYTE:
    DoComposeAlpha((imbyte*)src_image1->data[src_alpha], (imbyte*)src_image2->data[src_alpha], (imbyte*)dst_image->data[src_alpha], count, (int)255);
    break;
  case IM_USHORT:
    DoComposeAlpha((imushort*)src_image1->data[src_alpha], (imushort*)src_image2->data[src_alpha], (imushort*)dst_image->data[src_alpha], count, (int)65535);
    break;
  case IM_INT:
    DoComposeAlpha((int*)src_image1->data[src_alpha], (int*)src_image2->data[src_alpha], (int*)dst_image->data[src_alpha], count, (int)2147483647);
    break;
  case IM_FLOAT:
    DoComposeAlpha((float*)src_image1->data[src_alpha], (float*)src_image2->data[src_alpha], (float*)dst_image->data[src_alpha], count, 1.0f);
    break;
  }
}

static void DoBinaryConstOpCpxReal(imcfloat *map1, float value, imcfloat *map, int count, int op)
{
  int i;

  switch(op)
  {
  case IM_BIN_ADD:
    for (i = 0; i < count; i++)
      map[i] = add_op(map1[i], value);
    break;
  case IM_BIN_SUB:
    for (i = 0; i < count; i++)
      map[i] = sub_op(map1[i], value);
    break;
  case IM_BIN_MUL:
    for (i = 0; i < count; i++)
      map[i] = mul_op(map1[i], value);
    break;
  case IM_BIN_DIV:
    for (i = 0; i < count; i++)
      map[i] = div_op(map1[i], (imcfloat)value);
    break;
  case IM_BIN_DIFF:
    for (i = 0; i < count; i++)
      map[i] = diff_op(map1[i], value);
    break;
  case IM_BIN_MIN:
    for (i = 0; i < count; i++)
      map[i] = min_op(map1[i], value);
    break;
  case IM_BIN_MAX:
    for (i = 0; i < count; i++)
      map[i] = max_op(map1[i], value);
    break;
  case IM_BIN_POW:
    for (i = 0; i < count; i++)
      map[i] = pow_op(map1[i], value);
    break;
  }
}

template <class T1, class T2, class T3> 
static void DoBinaryConstOp(T1 *map1, T2 value, T3 *map, int count, int op)
{
  int i;

  switch(op)
  {
  case IM_BIN_ADD:
    for (i = 0; i < count; i++)
      map[i] = (T3)add_op((T2)map1[i], value);
    break;
  case IM_BIN_SUB:
    for (i = 0; i < count; i++)
      map[i] = (T3)sub_op((T2)map1[i], value);
    break;
  case IM_BIN_MUL:
    for (i = 0; i < count; i++)
      map[i] = (T3)mul_op((T2)map1[i], value);
    break;
  case IM_BIN_DIV:
    for (i = 0; i < count; i++)
      map[i] = (T3)div_op((T2)map1[i], value);
    break;
  case IM_BIN_DIFF:
    for (i = 0; i < count; i++)
      map[i] = (T3)diff_op((T2)map1[i], value);
    break;
  case IM_BIN_MIN:
    for (i = 0; i < count; i++)
      map[i] = (T3)min_op((T2)map1[i], value);
    break;
  case IM_BIN_MAX:
    for (i = 0; i < count; i++)
      map[i] = (T3)max_op((T2)map1[i], value);
    break;
  case IM_BIN_POW:
    for (i = 0; i < count; i++)
      map[i] = (T3)pow_op((T2)map1[i], value);
    break;
  }
}

template <class T1> 
static void DoBinaryConstOpByte(T1 *map1, int value, imbyte *map, int count, int op)
{
  int i;

  switch(op)
  {
  case IM_BIN_ADD:
    for (i = 0; i < count; i++)
      map[i] = (imbyte)crop_byte(add_op((int)map1[i], value));
    break;
  case IM_BIN_SUB:
    for (i = 0; i < count; i++)
      map[i] = (imbyte)crop_byte(sub_op((int)map1[i], value));
    break;
  case IM_BIN_MUL:
    for (i = 0; i < count; i++)
      map[i] = (imbyte)crop_byte(mul_op((int)map1[i], value));
    break;
  case IM_BIN_DIV:
    for (i = 0; i < count; i++)
      map[i] = (imbyte)crop_byte(div_op((int)map1[i], value));
    break;
  case IM_BIN_DIFF:
    for (i = 0; i < count; i++)
      map[i] = (imbyte)crop_byte(diff_op((int)map1[i], value));
    break;
  case IM_BIN_MIN:
    for (i = 0; i < count; i++)
      map[i] = (imbyte)crop_byte(min_op((int)map1[i], value));
    break;
  case IM_BIN_MAX:
    for (i = 0; i < count; i++)
      map[i] = (imbyte)crop_byte(max_op((int)map1[i], value));
    break;
  case IM_BIN_POW:
    for (i = 0; i < count; i++)
      map[i] = (imbyte)crop_byte(pow_op((int)map1[i], value));
    break;
  }
}

void imProcessArithmeticConstOp(const imImage* src_image1, float value, imImage* dst_image, int op)
{
  int count = src_image1->count;

  for (int i = 0; i < src_image1->depth; i++)
  {
    switch(src_image1->data_type)
    {
    case IM_BYTE:
      if (dst_image->data_type == IM_FLOAT)
        DoBinaryConstOp((imbyte*)src_image1->data[i], (float)value, (float*)dst_image->data[i], count, op);
      else if (dst_image->data_type == IM_USHORT)
        DoBinaryConstOp((imbyte*)src_image1->data[i], (imushort)value, (imushort*)dst_image->data[i], count, op);
      else if (dst_image->data_type == IM_INT)
        DoBinaryConstOp((imbyte*)src_image1->data[i], (int)value, (int*)dst_image->data[i], count, op);
      else
        DoBinaryConstOpByte((imbyte*)src_image1->data[i], (int)value, (imbyte*)dst_image->data[i], count, op);
      break;
    case IM_USHORT:
      if (dst_image->data_type == IM_FLOAT)
        DoBinaryConstOp((imushort*)src_image1->data[i], (float)value, (float*)dst_image->data[i], count, op);
      else if (dst_image->data_type == IM_INT)
        DoBinaryConstOp((imushort*)src_image1->data[i], (int)value, (int*)dst_image->data[i], count, op);
      else if (dst_image->data_type == IM_BYTE)
        DoBinaryConstOpByte((imushort*)src_image1->data[i], (int)value, (imbyte*)dst_image->data[i], count, op);
      else
        DoBinaryConstOp((imushort*)src_image1->data[i], (imushort)value, (imushort*)dst_image->data[i], count, op);
      break;
    case IM_INT:
      if (dst_image->data_type == IM_FLOAT)
        DoBinaryConstOp((int*)src_image1->data[i], (float)value, (float*)dst_image->data[i], count, op);
      else if (dst_image->data_type == IM_USHORT)
        DoBinaryConstOp((int*)src_image1->data[i], (int)value, (imushort*)dst_image->data[i], count, op);
      else if (dst_image->data_type == IM_BYTE)
        DoBinaryConstOpByte((int*)src_image1->data[i], (int)value, (imbyte*)dst_image->data[i], count, op);
      else
        DoBinaryConstOp((int*)src_image1->data[i], (int)value, (int*)dst_image->data[i], count, op);
      break;
    case IM_FLOAT:
      DoBinaryConstOp((float*)src_image1->data[i], (float)value, (float*)dst_image->data[i], count, op);
      break;
    case IM_CFLOAT:
      DoBinaryConstOpCpxReal((imcfloat*)src_image1->data[i], (float)value, (imcfloat*)dst_image->data[i], count, op);
      break;
    }
  }
}

void imProcessMultipleMean(const imImage** src_image_list, int src_image_count, imImage* dst_image)
{
  const imImage* image1 = src_image_list[0];

  int data_type = image1->data_type;
  if (image1->data_type == IM_BYTE)
    data_type = IM_USHORT;

  imImage *acum_image = imImageCreate(image1->width, image1->height, image1->color_space, data_type);
  if (!acum_image)
    return;

  for(int i = 0; i < src_image_count; i++)
  {
    const imImage *image = src_image_list[i];
    imProcessArithmeticOp(image, acum_image, acum_image, IM_BIN_ADD);  /* acum_image += image */

  }

  imProcessArithmeticConstOp(acum_image, float(src_image_count), dst_image, IM_BIN_DIV);

  imImageDestroy(acum_image);
}

void imProcessMultipleStdDev(const imImage** src_image_list, int src_image_count, const imImage *mean_image, imImage* dst_image)
{
  imImage* aux_image = imImageClone(dst_image);
  if (!aux_image)
    return;

  // sdtdev = sqrt( sum(sqr(x - m)) / N)

  // a = sum(sqr(x - m))
  for(int i = 0; i < src_image_count; i++)
  {
    // aux_image = image - mean_image
    imProcessArithmeticOp(src_image_list[i], mean_image, aux_image, IM_BIN_SUB);

    // aux_image = aux_image * aux_image
    imProcessUnArithmeticOp(aux_image, aux_image, IM_UN_SQR);

    // dst_image += aux_image
    imProcessArithmeticOp(aux_image, dst_image, dst_image, IM_BIN_ADD);
  }

  // dst_image = dst_image / src_image_count;
  imProcessArithmeticConstOp(dst_image, float(src_image_count), dst_image, IM_BIN_DIV);

  // dst_image = sqrt(dst_image);
  imProcessUnArithmeticOp(dst_image, dst_image, IM_UN_SQRT);

  imImageDestroy(aux_image);
}

template <class DT> 
static float AutoCovCalc(int width, int height, DT *src_map, DT *mean_map, int x, int y, float count)
{
  float value = 0;
  int ni = height - y;
  int nj = width - x;
  int offset, offset1;
  int next = width*y + x;

  for (int i = 0; i < ni; i++)
  {
    for (int j = 0; j < nj; j++)
    {
      offset = width*i + j;
      offset1 = offset + next;
      value += float(src_map[offset] - mean_map[offset]) * float(src_map[offset1] - mean_map[offset1]);
    }
  }

  return (value/count);
}

template <class DT> 
static int AutoCov(int width, int height, DT *src_map, DT *mean_map, float *dst_map, int counter)
{
  int count = width*height;

  for (int y = 0; y < height; y++)
  {
    for (int x = 0; x < width; x++)
    {
      *dst_map = AutoCovCalc(width, height, src_map, mean_map, x, y, (float)count);
      dst_map++;
    }

    if (!imCounterInc(counter))
      return 0;
  }

  return 1;
}

int imProcessAutoCovariance(const imImage* image, const imImage* mean_image, imImage* dst_image)
{
  int ret = 0;

  int counter = imCounterBegin("Auto Convariance");
  imCounterTotal(counter, image->depth*image->height, "Processing...");

  for (int i = 0; i < image->depth; i++)
  {
    switch(image->data_type)
    {
    case IM_BYTE:
      ret = AutoCov(image->width, image->height, (imbyte*)image->data[i], (imbyte*)mean_image->data[i], (float*)dst_image->data[i], counter);
      break;
    case IM_USHORT:
      ret = AutoCov(image->width, image->height, (imushort*)image->data[i], (imushort*)mean_image->data[i], (float*)dst_image->data[i], counter);
      break;
    case IM_INT:
      ret = AutoCov(image->width, image->height, (int*)image->data[i], (int*)mean_image->data[i], (float*)dst_image->data[i], counter);
      break;
    case IM_FLOAT:
      ret = AutoCov(image->width, image->height, (float*)image->data[i], (float*)mean_image->data[i], (float*)dst_image->data[i], counter);
      break;
    }

    if (!ret)
      break;
  }

  imCounterEnd(counter);

  return ret;
}

void imProcessMultiplyConj(const imImage* image1, const imImage* image2, imImage* NewImage)
{
  int total_count = image1->count*image1->depth;

  imcfloat* map = (imcfloat*)NewImage->data[0];
  imcfloat* map1 = (imcfloat*)image1->data[0];
  imcfloat* map2 = (imcfloat*)image2->data[0];
  imcfloat tmp; // this will allow an in-place operation

  for (int i = 0; i < total_count; i++)
  {
    tmp.real = map1->real * map2->real + map1->imag * map2->imag; 
    tmp.imag = map1->real * map2->imag - map1->imag * map2->real;
    *map = tmp;

    map++;
    map1++;
    map2++;
  }
}