First commit - moving from LuaForge to SourceForge

1 files changed, 503 insertions, 0 deletions

diff --git a/src/process/im_arithmetic_bin.cpp b/src/process/im_arithmetic_bin.cpp
new file mode 100644
index 0000000..74fe010
--- /dev/null
+++ b/src/process/im_arithmetic_bin.cpp
@@ -0,0 +1,503 @@
+/** \file
+ * \brief Binary Arithmetic Operations
+ *
+ * See Copyright Notice in im_lib.h
+ * $Id: im_arithmetic_bin.cpp,v 1.1 2008/10/17 06:16:33 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 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
+        DoBinaryOp((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;
+    }
+  }
+}
+
+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;
+  }
+}
+
+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
+        DoBinaryConstOp((imbyte*)src_image1->data[i], (imushort)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)
+        DoBinaryConstOp((imushort*)src_image1->data[i], (imushort)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)
+        DoBinaryConstOp((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];
+    imProcessUnArithmeticOp(image, acum_image, IM_UN_INC);
+  }
+
+  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
+    imProcessUnArithmeticOp(aux_image, dst_image, IM_UN_INC);
+  }
+
+  // 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++;
+  }
+}