Mercurial > forge

changeset 11218:5e36bb3c4688 octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
deriche: removing function due to non-standard license
author carandraug
date Thu, 08 Nov 2012 15:54:10 +0000
parents 6b967074572d
children 2427efc9c2fc
files main/image/NEWS main/image/devel/README main/image/devel/deriche.cc main/image/src/Makefile main/image/src/deriche.cc
diffstat 5 files changed, 298 insertions(+), 289 deletions(-) [+]
line wrap: on
line diff
--- a/main/image/NEWS	Thu Nov 08 11:40:58 2012 +0000
+++ b/main/image/NEWS	Thu Nov 08 15:54:10 2012 +0000
@@ -17,6 +17,8 @@
 
       imginfo
 
+ ** The function `deriche' has been removed.
+
  ** The complete set of functions to work with Analyze 7.5 files has been
     implemented. See `analyze75info', `analyze75read' and `analyze75write'.
 
--- a/main/image/devel/README	Thu Nov 08 11:40:58 2012 +0000
+++ b/main/image/devel/README	Thu Nov 08 15:54:10 2012 +0000
@@ -4,3 +4,10 @@
 
 They are helpul, though, to know how functions have been developped
 and for debugging.
+
+The function deriche has been removed from the package due to its non-standard
+license. Its original author and copyright owner does not reply to e-mails
+christian.kotz@gmx.net or facebook messages http://www.facebook.com/KoraxKorakos
+
+The deriche function would ideally be implemented again as an edge option,
+possible as part of the canny option since they are very similar.
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/image/devel/deriche.cc	Thu Nov 08 15:54:10 2012 +0000
@@ -0,0 +1,288 @@
+//Copyright (C) 2006 Christian Kotz <christian.kotz@gmx.net>
+//
+// This code has no warranty whatsover. Do what you like with this code 
+// as long as you leave this copyright in place.
+
+#include <octave/oct.h>
+
+
+static void dericheAbs(const double  *p, double *q, int w, int h, int linLen, double alpha);
+static void dericheVec(const double  *p, double *q, int w, int h, int 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\
+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\
+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 < 1 || nargin > 2){
+    error("call to deriche needs 1 or 2 arguments supplied.");
+    return octave_value_list ();
+  }
+  
+  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).");
+      return octave_value_list ();
+    default:
+      error("unknown method parameter.");
+      return octave_value_list ();
+    }
+  }
+
+  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:{
+      NDArray 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, int w, int h, int 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;
+  //const int sz = h*w; // unused 
+  try {
+    tmp = new double[2*h*w + 2*w];
+    if (!tmp) {
+      error("alloc error");
+      return;
+    }
+    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 // unused
+    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, int w, int h, int 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;
+  //const int sz = h*w;  // unused
+  try {
+    tmp = new double[2*h*w + 2*w];
+    if (!tmp) {
+      error("alloc error");
+      return;
+    }
+    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 // unused
+    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;
+}
--- a/main/image/src/Makefile	Thu Nov 08 11:40:58 2012 +0000
+++ b/main/image/src/Makefile	Thu Nov 08 15:54:10 2012 +0000
@@ -2,7 +2,7 @@
 
 all: __spatial_filtering__.oct __bilateral__.oct __custom_gaussian_smoothing__.oct \
 	__boundary__.oct bwlabel.oct bwfill.oct rotate_scale.oct hough_line.oct \
-	graycomatrix.oct deriche.oct __bwdist__.oct nonmax_supress.oct bwlabeln.oct
+	graycomatrix.oct __bwdist__.oct nonmax_supress.oct bwlabeln.oct
 
 %.oct: %.cc
 	$(MKOCTFILE) $<
--- a/main/image/src/deriche.cc	Thu Nov 08 11:40:58 2012 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,288 +0,0 @@
-//Copyright (C) 2006 Christian Kotz <christian.kotz@gmx.net>
-//
-// This code has no warranty whatsover. Do what you like with this code 
-// as long as you leave this copyright in place.
-
-#include <octave/oct.h>
-
-
-static void dericheAbs(const double  *p, double *q, int w, int h, int linLen, double alpha);
-static void dericheVec(const double  *p, double *q, int w, int h, int 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\
-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\
-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 < 1 || nargin > 2){
-    error("call to deriche needs 1 or 2 arguments supplied.");
-    return octave_value_list ();
-  }
-  
-  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).");
-      return octave_value_list ();
-    default:
-      error("unknown method parameter.");
-      return octave_value_list ();
-    }
-  }
-
-  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:{
-      NDArray 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, int w, int h, int 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;
-  //const int sz = h*w; // unused 
-  try {
-    tmp = new double[2*h*w + 2*w];
-    if (!tmp) {
-      error("alloc error");
-      return;
-    }
-    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 // unused
-    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, int w, int h, int 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;
-  //const int sz = h*w;  // unused
-  try {
-    tmp = new double[2*h*w + 2*w];
-    if (!tmp) {
-      error("alloc error");
-      return;
-    }
-    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 // unused
-    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;
-}