--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/image/src/deriche.cc	Fri Dec 08 06:43:30 2006 +0000
@@ -0,0 +1,308 @@
+ /* $Id$ */
+ #include <octave/oct.h>
+ /****************************************************************************
+  * (C)opyright Christian Kotz 2006
+  * This code has no warranty whatsover. Do what you like with this code 
+  *  as long as you leave this copyright in place.
+  ****************************************************************************
+  * author: Christian Kotz 
+  * date:     $Date$
+  * version: $Revision$
+  *
+  * (email: christian dot kotz at gmx dot net)
+  *
+  * History: 
+  * $Log$
+  * Revision 1.1  2006/12/08 06:43:30  cocus
+  * fast c implementation to replace m file deriche.m
+  *
+  */
+/*  
+              "-*- texinfo -*-\n\
+ @deftypefn{Loadable Function} {@var{b}} = deriche(@var{a}, @var{n}, @var{m})\n\
+ \n\
+ @cindex deriche edge detector\n\
+ \n\
+ Return edge detector image of @var{a} image according to an algorithm by Rachid Deriche. \n\
+ Matrix @var{a} is a real matrix, and @var{n} a non-negative real kernel scaling parameter (default 1.0).\
+ Specify @var{m} of 0 for a gradient magnitude result (default), @var{m} of 1 for a vector\
+ gradient result.\n @var{n} and @var{m} are optional arguments.\
+ Processing time is independent on var{n}. see Klette, Zameroni: Handbuch der\
+ Operatoren fuer die Bildverarbeitung, vieweg 2. ed. 1995 pp. 224--229. for\
+ details.\n\
+ Original paper: Deriche, R.: Fast algorithms for low-level vision: IEEE Trans PAMI-12 (1990) pp 78--87\n\
+"
+*/
+ 
+ static void dericheAbs(const double  *p, double *q, unsigned w, 	unsigned h, unsigned linLen, double alpha);
+ static void dericheVec(const double  *p, double *q, unsigned w,  unsigned h, unsigned linLen, double alpha);
+ 
+ DEFUN_DLD(deriche, args, ,
+            "-*- texinfo -*-\n\
+@deftypefn{Loadable Function} {@var{b}} = deriche(@var{a}, @var{n}, @var{m})\n\
+ \n\
+ @cindex deriche edge detector\n\
+ Return edge detector image of @var{a} image according to an algorithm by Rachid Deriche. \n\
+ Matrix @var{a} is a real matrix, and @var{n} a non-negative real kernel scaling parameter (default 1.0).\
+ Specify @var{m} = 0 for a gradient magnitude result (default), @var{m} = 1 for a vector\
+ gradient result.\n @var{n} and @var{m} are optional arguments.\
+ \n\n\
+ Processing time is independent on @var{n}.\n\
+ see for details: Klette, Zameroni: Handbuch der Operatoren fuer die Bildverarbeitung, vieweg 2. ed. 1995 pp. 224--229.\n\
+ Original paper: Deriche, R.: Fast algorithms for low-level vision: IEEE Trans PAMI-12 (1990) pp 78--87.\
+ \n\n\
+ Example:\
+  @example\n\
+  a = double(imread('myimg.png'));\n\
+  b = deriche(a, 1.0, 1);\n\
+  imshow(b(:,:,1));\n\
+  imshow(b(:,:,2));\n\
+  @end example\n\
+ @end deftypefn\
+ ")
+
+ {  
+     enum Method { absgrad, vecgrad, polargrad };
+     const int nargin = args.length();
+     
+     if (nargin < 2 || nargin > 3){
+        error("call to deriche needs 1 or 2 arguments supplied.");
+     }       
+     
+     const double alpha = (nargin <  2) ? 1.0: args(1).double_value();  
+     Method method = absgrad;
+     if (args.length() >  2){
+	int m = (int)(args(2).double_value());
+        switch(m){
+	case 0: break;
+	case 1: method = vecgrad; break;
+        case 2: method = polargrad;		
+            error("not yet implemented. Use builtin 'card2pol' after method 2 (cartesian vector grad).");		
+	   break;
+        default: error("unknown method parameter.");		
+	}	     
+     }
+  
+     Matrix p(args(0).matrix_value());
+     const int h = p.rows();
+     const int w = p.columns();
+     switch (method){
+     case absgrad:{
+        Matrix b(h, w);	
+        dericheAbs(p.fortran_vec(), b.fortran_vec(), h, w, h, alpha);
+        return octave_value(b);     
+     }
+     case vecgrad:{
+        ArrayN<double> b(dim_vector(h,w,2));
+        dericheVec(p.fortran_vec(), b.fortran_vec(), h, w, h, alpha);
+        return octave_value(b);
+     }
+     default:
+	error("method not yet implemented.");
+        return octave_value_list();
+     }     
+ }
+ 
+ // q has to be dense gapless, for w and liLen may differ
+ static void dericheAbs(const double  *p, double *q, unsigned w, unsigned h, unsigned linLen, double alpha){
+  double a(1.0-exp(-alpha));
+  a = - (a*a);
+  double b1(-2.0 * exp(-alpha));
+  double b2(exp(-2.0*alpha));
+  double a0(-a/(1.0-a*b1-b2));
+  double a1(a0*(alpha-1)*exp(-alpha));
+  double a2(a1-a0*b1);
+  double a3(-a0*b2);
+  double *tmp = 0;
+  unsigned const sz = h*w;	 
+  try {
+    tmp = new double[2*h*w + 2*w];
+	  if (!tmp) error("alloc error");
+    memset(tmp, 0, 2*h*w+2*linLen * sizeof(double));
+    double* B1 = tmp;
+    double* B2 = B1 + h *w;
+    double* Z3 = B2 + h * w;
+    double* Z2 = Z3 + w;
+  
+    const double  *ze; // int8
+    double  *za; // int8
+    double *Ba1;
+    double *Ba2;
+
+  // Berechnung von H
+  int y;
+  for(y=2; y < h; y++){  // (i)
+    ze = p + linLen*y;
+    Ba1 = B1 + w*y;  
+    for(int x=0;x < w; x++)
+      Ba1[x] = ze[x] - b1* *(Ba1 + x - w) - b2 * *(Ba1 + x -w -w);   
+  };
+ 
+  for(y = h-3 ; y >= 0 ; y--){       // (ii)
+     ze = p + (y+1) * linLen;
+     Ba1 = B1 + w*y;
+     Ba2 = B2 + w*y;
+    int x;
+    for(x=0; x < w; x++){
+      Ba2[x] = ze[x] - b1 * Ba2[x+w] - b2 * Ba2[x+w+w];
+      Ba1[x] = a * (Ba1[x] - Ba2[x]);
+    };
+  };
+  
+  for(y=0;y<h;y++){ // (iii, iv)
+     Ba1 = B1 + w*y; // Ba1 ist Z1 im Buch
+     int x;
+     for(x=2;x<w;x++)
+       Z2[x] = a0 * Ba1[x] + a1 * *(Ba1 + x - 1) - b1 * *(Z2 + x -1) - b2 * *(Z2 + x-2);
+     for(x = w-3; x >= 0 ; x--)
+       Z3[x] = a2 * Ba1[x+1] + a3 * Ba1[x+2] - b1 * Z3[x+1] - b2 * Z3[x+2];
+     for(x=0;x<w;x++){
+         q[y*w+x] = Z2[x] + Z3[x];
+     };
+   }
+  
+    // Berechnung von V
+    memset (Z2, 0, w*sizeof(double));
+    memset (Z3, 0, w*sizeof(double));
+  
+  for(y=0; y < h; y++){  // (v, vi)
+      ze = p + linLen*y;
+      Ba1 = B1 + w*y;
+      int x;
+      for(x=2;x < w; x++)
+        Z2[x] = *(ze+x-1) - b1 * *(Z2+x-1) - b2 * *(Z2+x-2);
+      for(x=w-3; x >=0 ; x--)
+        Z3[x] = ze[x+1] - b1 * Z3[x+1] - b2 * Z3[x+2];
+      for(x=0; x < w; x++)
+        Ba1[x] = a * (Z2[x] - Z3[x]);
+    };
+    for(y = 2 ; y < h ; y++){       // (vii)
+       Ba2 = B2 + w*y;
+       Ba1 = B1 + w*y;
+       int x;
+       for(x=0; x < w; x++)
+          Ba2[x] = (a0 + a1) * Ba1[x] - b1 * *(Ba2+x-w) - b2 * *(Ba2+x-w-w);
+    };
+    for(y = h - 3 ; y >= 0 ; y--){  // (viii)
+       Ba1 = B1 + y * w;
+       Ba2 = B2 + y * w;
+       memcpy(Z2, Ba2, w * sizeof(double)); // save contents of row in Z2
+       int x;
+       for(x= 0; x < w; x++){
+          Ba2[x] = a2 * Ba1[x+w] + a3 * Ba1[x+w+w]
+            - b1 * Ba2[x+w] - b2 * Ba2[x+w+w];
+       };
+       for(x= 0; x < w; x++){//  memset (B1, 0, h*w*sizeof(double));
+          double z1 = Ba2[x] + Z2[x];
+          double z2 = q[y*w+x];
+          q[y*w+x] = sqrt(z1 * z1 + z2 * z2);
+       };
+    }  
+  }catch(...){
+	delete [] tmp;
+	throw;
+  }
+  delete[] tmp;  
+  }
+
+  // q has to be dense gapless, for w and liLen may differ
+  static void dericheVec(const double  *p, double *q, unsigned w, unsigned h, unsigned linLen, double alpha){
+  double a(1.0-exp(-alpha));
+  a = - (a*a);
+  double b1(-2.0 * exp(-alpha));
+  double b2(exp(-2.0*alpha));
+  double a0(-a/(1.0-a*b1-b2));
+  double a1(a0*(alpha-1)*exp(-alpha));
+  double a2(a1-a0*b1);
+  double a3(-a0*b2);
+  double *tmp = 0;
+  double *r=q+h*w;
+  unsigned const sz = h*w;	 
+  try {
+    tmp = new double[2*h*w + 2*w];
+	  if (!tmp) error("alloc error");
+    memset(tmp, 0, 2*h*w+2*linLen * sizeof(double));
+    double* B1 = tmp;
+    double* B2 = B1 + h *w;
+    double* Z3 = B2 + h * w;
+    double* Z2 = Z3 + w;
+  
+    const double  *ze; // int8
+    double  *za; // int8
+    double *Ba1;
+    double *Ba2;
+
+  // Berechnung von H
+  int y;
+  for(y=2; y < h; y++){  // (i)
+    ze = p + linLen*y;
+    Ba1 = B1 + w*y;  
+    for(int x=0;x < w; x++)
+      Ba1[x] = ze[x] - b1* *(Ba1 + x - w) - b2 * *(Ba1 + x -w -w);   
+  };
+ 
+  for(y = h-3 ; y >= 0 ; y--){       // (ii)
+     ze = p + (y+1) * linLen;
+     Ba1 = B1 + w*y;
+     Ba2 = B2 + w*y;
+    int x;
+    for(x=0; x < w; x++){
+      Ba2[x] = ze[x] - b1 * Ba2[x+w] - b2 * Ba2[x+w+w];
+      Ba1[x] = a * (Ba1[x] - Ba2[x]);
+    };
+  };
+  
+  for(y=0;y<h;y++){ // (iii, iv)
+     Ba1 = B1 + w*y; // Ba1 ist Z1 im Buch
+     int x;
+     for(x=2;x<w;x++)
+       Z2[x] = a0 * Ba1[x] + a1 * *(Ba1 + x - 1) - b1 * *(Z2 + x -1) - b2 * *(Z2 + x-2);
+     for(x = w-3; x >= 0 ; x--)
+       Z3[x] = a2 * Ba1[x+1] + a3 * Ba1[x+2] - b1 * Z3[x+1] - b2 * Z3[x+2];
+     for(x=0;x<w;x++){
+       q[y*w+x] =  Z2[x] + Z3[x];
+     };
+   }
+  
+    // Berechnung von V
+    memset (Z2, 0, w*sizeof(double));
+    memset (Z3, 0, w*sizeof(double));
+  
+  for(y=0; y < h; y++){  // (v, vi)
+      ze = p + linLen*y;
+      Ba1 = B1 + w*y;
+      int x;
+      for(x=2;x < w; x++)
+        Z2[x] = *(ze+x-1) - b1 * *(Z2+x-1) - b2 * *(Z2+x-2);
+      for(x=w-3; x >=0 ; x--)
+        Z3[x] = ze[x+1] - b1 * Z3[x+1] - b2 * Z3[x+2];
+      for(x=0; x < w; x++)
+        Ba1[x] = a * (Z2[x] - Z3[x]);
+    };
+    for(y = 2 ; y < h ; y++){       // (vii)
+       Ba2 = B2 + w*y;
+       Ba1 = B1 + w*y;
+       int x;
+       for(x=0; x < w; x++)
+          Ba2[x] = (a0 + a1) * Ba1[x] - b1 * *(Ba2+x-w) - b2 * *(Ba2+x-w-w);
+    };
+    for(y = h - 3 ; y >= 0 ; y--){  // (viii)
+       Ba1 = B1 + y * w;
+       Ba2 = B2 + y * w;
+       memcpy(Z2, Ba2, w * sizeof(double)); // save contents of row in Z2
+       int x;
+       for(x= 0; x < w; x++){
+          Ba2[x] = a2 * Ba1[x+w] + a3 * Ba1[x+w+w]
+            - b1 * Ba2[x+w] - b2 * Ba2[x+w+w];
+       };
+       for(x= 0; x < w; x++){//  memset (B1, 0, h*w*sizeof(double));
+          r[y*w+x]  = Ba2[x] + Z2[x];
+       };
+    }  
+  }catch(...){
+	delete [] tmp;
+	throw;
+  }
+  delete[] tmp;  
+  }
+