new file mode 100644
--- /dev/null
+++ b/inst/bwdist.m
@@ -0,0 +1,50 @@
+## Copyright (C) 2006  Søren Hauberg
+## 
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+## 
+## This program is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details. 
+## 
+## You should have received a copy of the GNU General Public License
+## along with this file.  If not, write to the Free Software Foundation,
+## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} @var{d} = bwdist(@var{bw}, @var{method})
+## Computes the distance transform of the binary image @var{bw}.
+## The result @var{d} is a matrix of the same size as @var{bw}, where
+## each value is the shortest distance to a non-zero pixel in @var{bw}.
+## 
+## @var{method} changes the used distance function. Currently
+## the Euclidian distance is the only supported distance function.
+## @end deftypefn
+
+function D = bwdist(bw, method = "euclidian")
+  ## Check input
+  if (nargin == 0)
+    print_usage();
+  endif
+  
+  if (!ismatrix(bw) || ndims(bw) != 2)
+    error("bwdist: input must be a 2-dimensional matrix");
+  endif
+  
+  if (!ischar(method))
+    error("bwdist: method name must be a string");
+  endif
+
+  ## Do the work
+  bw = (bw != 0);
+  switch (lower(method(1)))
+    case "e" 
+      ## Euclidian distance transform
+      D = __bwdist(bw);
+    otherwise
+      error("bwdist: unsupported method '%s'", method);
+  endswitch
+endfunction

--- a/src/Makefile
+++ b/src/Makefile
@@ -13,7 +13,7 @@
 endif
 
 all: cordflt2.oct bwlabel.oct bwfill.oct rotate_scale.oct \
-	houghtf.oct graycomatrix.oct deriche.oct \
+	houghtf.oct graycomatrix.oct deriche.oct __bwdist.oct \
 	$(JPEG) $(PNG) $(IMAGEMAGICK)
 
 jpgread.oct: jpgread.cc

new file mode 100644
--- /dev/null
+++ b/src/__bwdist.cc
@@ -0,0 +1,148 @@
+/*
+Copyright (C) 2006 Pedro Felzenszwalb
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program; if not, write to the Free Software
+Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+*/
+
+/* 
+ * Remark by Søren Hauberg, december 28th 2006.
+ *
+ * This code was mainly written by Pedro Felzenszwalb and published
+ * on http://people.cs.uchicago.edu/~pff/dt/
+ * Pedro kindly released his code under the GPL and I ported it to
+ * Octave.
+ */
+
+#include <octave/oct.h>
+
+#define INF 1E20
+
+template <class T>
+inline T square(const T &x) { return x*x; };
+
+/* dt of 1d function using squared distance */
+static float *dt(float *f, int n) {
+  float *d = new float[n];
+  int *v = new int[n];
+  float *z = new float[n+1];
+  int k = 0;
+  v[0] = 0;
+  z[0] = -INF;
+  z[1] = +INF;
+  for (int q = 1; q <= n-1; q++) {
+    float s  = ((f[q]+square(q))-(f[v[k]]+square(v[k])))/(2*q-2*v[k]);
+    while (s <= z[k]) {
+      k--;
+      s  = ((f[q]+square(q))-(f[v[k]]+square(v[k])))/(2*q-2*v[k]);
+    }
+    k++;
+    v[k] = q;
+    z[k] = s;
+    z[k+1] = +INF;
+  }
+
+  k = 0;
+  for (int q = 0; q <= n-1; q++) {
+    while (z[k+1] < q)
+      k++;
+    d[q] = square(q-v[k]) + f[v[k]];
+  }
+
+  delete [] v;
+  delete [] z;
+  return d;
+}
+
+/* dt of 2d function using squared distance */
+static void dt(NDArray &im) {
+  const int width = im.dim2();
+  const int height = im.dim1();
+  float *f = new float[std::max(width,height)];
+
+  // transform along columns
+  for (int x = 0; x < width; x++) {
+    for (int y = 0; y < height; y++) {
+      f[y] = im(x, y);
+    }
+    float *d = dt(f, height);
+    for (int y = 0; y < height; y++) {
+      im(x, y) = d[y];
+    }
+    delete [] d;
+  }
+
+  // transform along rows
+  for (int y = 0; y < height; y++) {
+    for (int x = 0; x < width; x++) {
+      f[x] = im(x, y);
+    }
+    float *d = dt(f, width);
+    for (int x = 0; x < width; x++) {
+      im(x, y) = d[x];
+    }
+    delete [] d;
+  }
+
+  delete f;
+}
+
+
+/* dt of binary image using squared distance */
+void dt(boolNDArray &im, NDArray &out) {
+  const int width = im.dim2();
+  const int height = im.dim1();
+
+  for (int y = 0; y < height; y++) {
+    for (int x = 0; x < width; x++) {
+      if (im(x, y))
+        out(x, y) = 0;
+      else
+        out(x, y) = INF;
+    }
+  }
+
+  dt(out);
+}
+
+DEFUN_DLD(__bwdist, args, nargout, "\
+-*- texinfo -*-\n\
+@deftypefn {Function File} __bwdist (@var{bw})\n\
+Computes the Euclidian Distance Transform for a binary image @var{bw}.\n\
+You should not call this function directly, instead call 'bwdist'.\n\
+@seealso{bwdist}\n\
+@end deftypefn\n\
+")
+{
+  octave_value_list retval;
+  
+  boolNDArray bw = args(0).bool_array_value();
+  if (error_state) {
+    error("__bwdist: input must be a boolean matrix");
+    return retval;
+  }
+  
+  const dim_vector dims = bw.dims();
+  if (dims.length() != 2) {
+    error("__bwdist: currently only binary images are supported");
+    return retval;
+  }
+  
+  NDArray *out = new NDArray(dims);
+  
+  dt(bw, *out);
+  
+  retval.append(*out);
+  return retval;
+}