diff --git a/skimage/morphology/data/3/pngs/1.png b/skimage/morphology/data/3/pngs/1.png new file mode 100644 index 00000000..da940f11 Binary files /dev/null and b/skimage/morphology/data/3/pngs/1.png differ diff --git a/skimage/morphology/data/3/pngs/2.png b/skimage/morphology/data/3/pngs/2.png new file mode 100644 index 00000000..ada13a9e Binary files /dev/null and b/skimage/morphology/data/3/pngs/2.png differ diff --git a/skimage/morphology/data/3/pngs/3.png b/skimage/morphology/data/3/pngs/3.png new file mode 100644 index 00000000..00ee7d1f Binary files /dev/null and b/skimage/morphology/data/3/pngs/3.png differ diff --git a/skimage/morphology/data/3/pngs/5.png b/skimage/morphology/data/3/pngs/5.png new file mode 100644 index 00000000..8378036e Binary files /dev/null and b/skimage/morphology/data/3/pngs/5.png differ diff --git a/skimage/morphology/data/3/pngs/6.png b/skimage/morphology/data/3/pngs/6.png new file mode 100644 index 00000000..3c6051ae Binary files /dev/null and b/skimage/morphology/data/3/pngs/6.png differ diff --git a/skimage/morphology/data/3/pngs/7.png b/skimage/morphology/data/3/pngs/7.png new file mode 100644 index 00000000..61eb7001 Binary files /dev/null and b/skimage/morphology/data/3/pngs/7.png differ diff --git a/skimage/morphology/data/3/stack.tif b/skimage/morphology/data/3/stack.tif new file mode 100644 index 00000000..28752135 Binary files /dev/null and b/skimage/morphology/data/3/stack.tif differ diff --git a/skimage/morphology/data/3/stack_fiji.tif b/skimage/morphology/data/3/stack_fiji.tif new file mode 100644 index 00000000..93737675 Binary files /dev/null and b/skimage/morphology/data/3/stack_fiji.tif differ diff --git a/skimage/morphology/data/4/1_fiji.png b/skimage/morphology/data/4/1_fiji.png new file mode 100644 index 00000000..03a7b79a Binary files /dev/null and b/skimage/morphology/data/4/1_fiji.png differ diff --git a/skimage/morphology/data/4/pngs/1.png b/skimage/morphology/data/4/pngs/1.png new file mode 100644 index 00000000..d03660d2 Binary files /dev/null and b/skimage/morphology/data/4/pngs/1.png differ diff --git a/skimage/morphology/data/4/pngs/2.png b/skimage/morphology/data/4/pngs/2.png new file mode 100644 index 00000000..efa5c104 Binary files /dev/null and b/skimage/morphology/data/4/pngs/2.png differ diff --git a/skimage/morphology/data/4/pngs/3.png b/skimage/morphology/data/4/pngs/3.png new file mode 100644 index 00000000..28100eee Binary files /dev/null and b/skimage/morphology/data/4/pngs/3.png differ diff --git a/skimage/morphology/data/4/pngs/4.png b/skimage/morphology/data/4/pngs/4.png new file mode 100644 index 00000000..28100eee Binary files /dev/null and b/skimage/morphology/data/4/pngs/4.png differ diff --git a/skimage/morphology/data/4/pngs/5.png b/skimage/morphology/data/4/pngs/5.png new file mode 100644 index 00000000..28100eee Binary files /dev/null and b/skimage/morphology/data/4/pngs/5.png differ diff --git a/skimage/morphology/data/4/pngs/6.png b/skimage/morphology/data/4/pngs/6.png new file mode 100644 index 00000000..28100eee Binary files /dev/null and b/skimage/morphology/data/4/pngs/6.png differ diff --git a/skimage/morphology/data/4/stack.tif b/skimage/morphology/data/4/stack.tif new file mode 100644 index 00000000..a4f9f6b7 Binary files /dev/null and b/skimage/morphology/data/4/stack.tif differ diff --git a/skimage/morphology/data/4/stack_fiji.tif b/skimage/morphology/data/4/stack_fiji.tif new file mode 100644 index 00000000..c0d182fd Binary files /dev/null and b/skimage/morphology/data/4/stack_fiji.tif differ diff --git a/skimage/morphology/data/cross.png b/skimage/morphology/data/cross.png new file mode 100644 index 00000000..bd6db21c Binary files /dev/null and b/skimage/morphology/data/cross.png differ diff --git a/skimage/morphology/data/cross.txt b/skimage/morphology/data/cross.txt new file mode 100644 index 00000000..2118f7f1 --- /dev/null +++ b/skimage/morphology/data/cross.txt @@ -0,0 +1,24 @@ +0 0 0 0 0 255 255 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 255 255 255 255 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 255 255 255 255 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 255 255 255 255 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 255 255 255 255 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 255 255 255 255 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 0 0 0 +0 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 0 0 0 +0 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 0 0 +0 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 0 0 +0 0 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 0 0 +0 0 0 0 0 255 255 255 255 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 255 255 255 255 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 255 255 255 255 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 255 255 255 255 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 255 255 255 255 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 255 255 255 255 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 255 255 255 255 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 255 255 255 255 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 255 255 255 255 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 255 255 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 diff --git a/skimage/morphology/data/cross_fiji.txt b/skimage/morphology/data/cross_fiji.txt new file mode 100644 index 00000000..c4ea69a9 --- /dev/null +++ b/skimage/morphology/data/cross_fiji.txt @@ -0,0 +1,24 @@ +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 255 255 255 255 0 255 255 255 255 255 255 255 255 255 255 255 255 0 0 0 0 +0 0 0 0 0 0 0 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 diff --git a/skimage/morphology/data/loop.png b/skimage/morphology/data/loop.png new file mode 100644 index 00000000..05fdc259 Binary files /dev/null and b/skimage/morphology/data/loop.png differ diff --git a/skimage/morphology/data/loop.txt b/skimage/morphology/data/loop.txt new file mode 100644 index 00000000..277e1b47 --- /dev/null +++ b/skimage/morphology/data/loop.txt @@ -0,0 +1,24 @@ +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 255 255 255 255 255 255 255 255 255 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 255 255 255 255 255 255 255 255 255 255 255 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 255 255 255 255 255 255 255 0 0 0 255 255 0 0 0 0 0 +0 0 0 0 0 0 0 0 255 255 255 255 0 0 0 0 0 0 255 255 0 0 0 0 0 +0 0 0 0 0 0 0 255 255 255 255 0 0 0 0 0 0 0 255 255 0 0 0 0 0 +0 0 0 0 0 0 255 255 255 255 0 0 0 0 0 0 255 255 255 255 0 0 0 0 0 +0 0 0 0 0 0 255 255 255 255 0 0 0 0 0 255 255 255 255 0 0 0 0 0 0 +0 0 0 0 0 0 0 255 255 255 0 0 0 255 255 255 255 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 255 255 0 255 255 255 255 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 255 255 255 255 255 255 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 255 255 255 255 255 255 255 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 255 255 255 255 255 255 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 255 255 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 255 255 255 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 255 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 diff --git a/skimage/morphology/data/loop_fiji.txt b/skimage/morphology/data/loop_fiji.txt new file mode 100644 index 00000000..4d38bfac --- /dev/null +++ b/skimage/morphology/data/loop_fiji.txt @@ -0,0 +1,24 @@ +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 255 255 255 255 255 255 255 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 255 0 0 0 0 0 0 0 255 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 255 0 0 0 0 0 0 0 0 255 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 255 0 0 0 0 0 0 0 0 0 255 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 255 0 0 0 0 0 0 0 0 0 255 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 255 0 0 0 0 0 0 0 255 255 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 255 0 0 0 0 0 0 255 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 255 0 0 0 0 255 255 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 255 0 0 0 255 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 255 255 255 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 255 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 diff --git a/skimage/morphology/data/strip.png b/skimage/morphology/data/strip.png new file mode 100644 index 00000000..74e9b3e4 Binary files /dev/null and b/skimage/morphology/data/strip.png differ diff --git a/skimage/morphology/data/strip.txt b/skimage/morphology/data/strip.txt new file mode 100644 index 00000000..abc45a08 --- /dev/null +++ b/skimage/morphology/data/strip.txt @@ -0,0 +1,12 @@ +0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 255 255 255 0 0 0 0 +0 0 0 0 255 255 255 0 0 0 0 +0 0 0 0 255 255 255 0 0 0 0 +0 0 0 0 255 255 255 0 0 0 0 +0 0 0 0 255 255 255 0 0 0 0 +0 0 0 0 255 255 255 0 0 0 0 +0 0 0 0 255 255 255 0 0 0 0 +0 0 0 0 255 255 255 0 0 0 0 +0 0 0 0 255 255 255 0 0 0 0 +0 0 0 0 255 255 255 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 diff --git a/skimage/morphology/data/strip_fiji.txt b/skimage/morphology/data/strip_fiji.txt new file mode 100644 index 00000000..9ff3049c --- /dev/null +++ b/skimage/morphology/data/strip_fiji.txt @@ -0,0 +1,12 @@ +0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 255 0 0 0 0 0 +0 0 0 0 0 255 0 0 0 0 0 +0 0 0 0 0 255 0 0 0 0 0 +0 0 0 0 0 255 0 0 0 0 0 +0 0 0 0 0 255 0 0 0 0 0 +0 0 0 0 0 255 0 0 0 0 0 +0 0 0 0 0 255 0 0 0 0 0 +0 0 0 0 0 255 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 diff --git a/skimage/morphology/data/two-hole.png b/skimage/morphology/data/two-hole.png new file mode 100644 index 00000000..8754e2be Binary files /dev/null and b/skimage/morphology/data/two-hole.png differ diff --git a/skimage/morphology/data/two-hole.txt b/skimage/morphology/data/two-hole.txt new file mode 100644 index 00000000..680e7beb --- /dev/null +++ b/skimage/morphology/data/two-hole.txt @@ -0,0 +1,16 @@ +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 255 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 255 255 255 255 255 0 0 0 0 0 +0 0 0 0 0 255 255 255 255 255 255 255 255 255 0 255 0 +0 0 0 255 255 255 255 255 255 255 255 255 255 255 255 255 0 +0 0 255 255 255 255 0 0 255 255 255 0 0 255 255 255 0 +0 0 255 255 255 0 0 0 255 255 255 0 0 255 255 255 0 +0 0 0 255 255 0 0 0 255 255 255 0 0 255 255 255 0 +0 0 0 255 255 255 0 0 255 255 255 0 0 255 255 255 0 +0 0 0 0 255 255 0 0 255 255 255 0 0 255 255 255 0 +0 0 0 0 255 255 255 255 255 255 255 255 255 255 255 255 0 +0 0 0 0 0 255 255 255 255 255 255 255 255 255 255 255 0 +0 0 0 0 0 255 255 255 255 255 255 255 255 255 255 255 0 +0 0 0 0 0 0 0 255 255 255 255 255 255 255 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 diff --git a/skimage/morphology/data/two-hole_fiji.txt b/skimage/morphology/data/two-hole_fiji.txt new file mode 100644 index 00000000..6c16d3f7 --- /dev/null +++ b/skimage/morphology/data/two-hole_fiji.txt @@ -0,0 +1,16 @@ +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 255 255 0 255 255 0 0 0 0 0 +0 0 0 0 0 255 255 0 0 255 0 0 255 0 0 0 0 +0 0 0 0 255 0 0 0 0 255 0 0 0 255 0 0 0 +0 0 0 255 0 0 0 0 0 255 0 0 0 0 255 0 0 +0 0 0 255 0 0 0 0 0 255 0 0 0 0 255 0 0 +0 0 0 0 255 0 0 0 0 255 0 0 0 0 255 0 0 +0 0 0 0 255 0 0 0 0 255 0 0 0 0 255 0 0 +0 0 0 0 0 255 0 0 0 255 0 0 0 0 255 0 0 +0 0 0 0 0 0 255 255 0 255 0 0 0 255 0 0 0 +0 0 0 0 0 0 0 0 255 0 255 255 255 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 diff --git a/skimage/morphology/skel.py b/skimage/morphology/skel.py new file mode 100644 index 00000000..e53c7ee6 --- /dev/null +++ b/skimage/morphology/skel.py @@ -0,0 +1,678 @@ +from __future__ import division, print_function, absolute_import + +import numpy as np + + +def _prepare_image(img_in): + """Convert to a binary image, pad the it w/ zeros, and ensure it's 3D. + """ + if img_in.ndim < 2 or img_in.ndim > 3: + raise ValueError('expect 2D, got ndim = %s' % img_in.ndim) + + img = img_in.copy() + + if img.ndim == 2: + img = img.reshape((1,) + img.shape) + + # normalize to binary + img[img != 0] = 1 + + # pad w/ zeros to simplify dealing w/ neighborhood of a pixel + img_o = np.zeros(tuple(s + 2 for s in img.shape), + dtype=np.uint8) + img_o[1:-1, 1:-1, 1:-1] = img.astype(np.uint8) + + return img_o + + +def _postprocess_image(img_o): + """Clip the image (padding is an implementation detail), convert to b/w. + """ + img_oo = img_o[1:-1, 1:-1, 1:-1] + img_oo = img_oo.squeeze() + img_oo *= 255 + return img_oo + + +def get_neighborhood(img, p, r, c): + """Get the neighborhood of a pixel. + + Assume zero boundary conditions. Image is already padded, so no + out-of-bounds checking. + """ + neighborhood = np.zeros(27, dtype=np.uint8) + + neighborhood[0] = img[p-1, r-1, c-1] + neighborhood[1] = img[p-1, r, c-1] + neighborhood[2] = img[p-1, r+1, c-1] + + neighborhood[ 3] = img[p-1, r-1, c] + neighborhood[ 4] = img[p-1, r, c] + neighborhood[ 5] = img[p-1, r+1, c] + + neighborhood[ 6] = img[p-1, r-1, c+1] + neighborhood[ 7] = img[p-1, r, c+1] + neighborhood[ 8] = img[p-1, r+1, c+1] + + neighborhood[ 9] = img[p, r-1, c-1] + neighborhood[10] = img[p, r, c-1] + neighborhood[11] = img[p, r+1, c-1] + + neighborhood[12] = img[p, r-1, c] + neighborhood[13] = img[p, r, c] + neighborhood[14] = img[p, r+1, c] + + neighborhood[15] = img[p, r-1, c+1] + neighborhood[16] = img[p, r, c+1] + neighborhood[17] = img[p, r+1, c+1] + + neighborhood[18] = img[p+1, r-1, c-1] + neighborhood[19] = img[p+1, r, c-1] + neighborhood[20] = img[p+1, r+1, c-1] + + neighborhood[21] = img[p+1, r-1, c] + neighborhood[22] = img[p+1, r, c] + neighborhood[23] = img[p+1, r+1, c] + + neighborhood[24] = img[p+1, r-1, c+1] + neighborhood[25] = img[p+1, r, c+1] + neighborhood[26] = img[p+1, r+1, c+1] + + return neighborhood + + +###### look-up tables +def fill_numpoints_LUT(n=256): + p = int(np.log2(n) + 1) + return np.sum(np.arange(n)[:, None] & (1 << np.arange(p)) != 0, axis=1) + +NUMPOINTS_LUT = fill_numpoints_LUT() + + +def fill_Euler_LUT(): + LUT = np.zeros(256, dtype=int) + + LUT[1] = 1 + LUT[3] = -1 + LUT[5] = -1 + LUT[7] = 1 + LUT[9] = -3 + LUT[11] = -1 + LUT[13] = -1 + LUT[15] = 1 + LUT[17] = -1 + LUT[19] = 1 + LUT[21] = 1 + LUT[23] = -1 + LUT[25] = 3 + LUT[27] = 1 + LUT[29] = 1 + LUT[31] = -1 + LUT[33] = -3 + LUT[35] = -1 + LUT[37] = 3 + LUT[39] = 1 + LUT[41] = 1 + LUT[43] = -1 + LUT[45] = 3 + LUT[47] = 1 + LUT[49] = -1 + LUT[51] = 1 + + LUT[53] = 1 + LUT[55] = -1 + LUT[57] = 3 + LUT[59] = 1 + LUT[61] = 1 + LUT[63] = -1 + LUT[65] = -3 + LUT[67] = 3 + LUT[69] = -1 + LUT[71] = 1 + LUT[73] = 1 + LUT[75] = 3 + LUT[77] = -1 + LUT[79] = 1 + LUT[81] = -1 + LUT[83] = 1 + LUT[85] = 1 + LUT[87] = -1 + LUT[89] = 3 + LUT[91] = 1 + LUT[93] = 1 + LUT[95] = -1 + LUT[97] = 1 + LUT[99] = 3 + LUT[101] = 3 + LUT[103] = 1 + + LUT[105] = 5 + LUT[107] = 3 + LUT[109] = 3 + LUT[111] = 1 + LUT[113] = -1 + LUT[115] = 1 + LUT[117] = 1 + LUT[119] = -1 + LUT[121] = 3 + LUT[123] = 1 + LUT[125] = 1 + LUT[127] = -1 + LUT[129] = -7 + LUT[131] = -1 + LUT[133] = -1 + LUT[135] = 1 + LUT[137] = -3 + LUT[139] = -1 + LUT[141] = -1 + LUT[143] = 1 + LUT[145] = -1 + LUT[147] = 1 + LUT[149] = 1 + LUT[151] = -1 + LUT[153] = 3 + LUT[155] = 1 + + LUT[157] = 1 + LUT[159] = -1 + LUT[161] = -3 + LUT[163] = -1 + LUT[165] = 3 + LUT[167] = 1 + LUT[169] = 1 + LUT[171] = -1 + LUT[173] = 3 + LUT[175] = 1 + LUT[177] = -1 + LUT[179] = 1 + LUT[181] = 1 + LUT[183] = -1 + LUT[185] = 3 + LUT[187] = 1 + LUT[189] = 1 + LUT[191] = -1 + LUT[193] = -3 + LUT[195] = 3 + LUT[197] = -1 + LUT[199] = 1 + LUT[201] = 1 + LUT[203] = 3 + LUT[205] = -1 + LUT[207] = 1 + + LUT[209] = -1 + LUT[211] = 1 + LUT[213] = 1 + LUT[215] = -1 + LUT[217] = 3 + LUT[219] = 1 + LUT[221] = 1 + LUT[223] = -1 + LUT[225] = 1 + LUT[227] = 3 + LUT[229] = 3 + LUT[231] = 1 + LUT[233] = 5 + LUT[235] = 3 + LUT[237] = 3 + LUT[239] = 1 + LUT[241] = -1 + LUT[243] = 1 + LUT[245] = 1 + LUT[247] = -1 + LUT[249] = 3 + LUT[251] = 1 + LUT[253] = 1 + LUT[255] = -1 + return LUT + +LUT = fill_Euler_LUT() + + +### Octants (indexOctantXXX functions) +OCTANTS = tuple(range(8)) +NEB, NWB, SEB, SWB, NEU, NWU, SEU, SWU = OCTANTS + +neib_idx = np.empty((8, 7), dtype=int) +neib_idx[NEB, ...] = [2, 1, 11, 10, 5, 4, 14] +neib_idx[NWB, ...] = [0, 9, 3, 12, 1, 10, 4] +neib_idx[SEB, ...] = [8, 7, 17, 16, 5, 4, 14] +neib_idx[SWB, ...] = [6, 15, 7, 16, 3, 12, 4] +neib_idx[NEU, ...] = [20, 23, 19, 22, 11, 14, 10] +neib_idx[NWU, ...] = [18, 21, 9, 12, 19, 22, 10] +neib_idx[SEU, ...] = [26, 23, 17, 14, 25, 22, 16] +neib_idx[SWU, ...] = [24, 25, 15, 16, 21, 22, 12] + +def index_octants(octant, neighbors): + n = 1 + for j, idx in enumerate(neib_idx[octant]): + if neighbors[idx] == 1: + n |= 2**(7 - j) + return n + + +def is_surfacepoint(neighbors, points_LUT): + for octant in OCTANTS: + n = index_octants(octabt, neighbors) + if n not in (240, 165, 170) and points_LUT[n] > 2: + return False + return True + + +def is_Euler_invariant(neighbors): + """Check if a point is Euler invariant. + + Calculate Euler characteristc for each octant and sum up. + + Parameters + ---------- + neighbors : ndarray, shape (27,) + neighbors of a point + + Returns + ------- + bool + + """ + euler_char = 0 + for octant in OCTANTS: + n = index_octants(octant, neighbors) + euler_char += LUT[n] + return euler_char == 0 + + +def is_simple_point(neighbors): + """Check is a point is a Simple Point. + + This method is named 'N(v)_labeling' in [Lee94]. + Outputs the number of connected objects in a neighborhood of a point + after this point would have been removed. + + Parameters + ---------- + neighbors : ndarray, shape(27,) + neighbors of the point + + Returns + ------- + bool + Whether the point is simple or not. + + """ + # copy neighbors for labeling + # ignore center pixel (i=13) when counting (see [Lee94]) + cube = np.r_[neighbors[:13], neighbors[14:]] + + # set initial label + label = 2 + + # for all point in the neighborhood + for i in range(26): + if cube[i] == 1: + # voxel has not been labeled yet + # start recursion with any octant that contains the point i + if i in (0, 1, 3, 4, 9, 10, 12): + octree_labeling(1, label, cube) + elif i in (2, 5, 11, 13): + octree_labeling(2, label, cube) + elif i in (6, 7, 14, 15): + octree_labeling(3, label, cube) + elif i in (8, 16): + octree_labeling(4, label, cube) + elif i in (17, 18, 20, 21): + octree_labeling(5, label, cube) + elif i in (19, 22): + octree_labeling(6, label, cube) + elif i in (23, 24): + octree_labeling(7, label, cube) + elif i == 25: + octree_labeling(8, label, cube) + else: + raise ValueError("Never be here. i = %s" % i) + label += 1 + if label - 2 >= 2: + return False + return True + + +def octree_labeling(octant, label, cube): + """This is a recursive method that calculates the number of connected + components in the 3D neighborhood after the center pixel would + have been removed. + + Parameters + ---------- + octant : int + octant index + label : int + the current label of the center point + cube : ndarray, shape(26,) + local neighborhood of the point + + """ + # check if there are points in the octant with value 1 + if octant == 1: + # set points in this octant to current label + # and recursive labeling of adjacent octants + if cube[0] == 1: + cube[0] = label + if cube[1] == 1: + cube[1] = label + octree_labeling(2, label, cube) + if cube[3] == 1: + cube[3] = label + octree_labeling(3, label, cube) + if cube[4] == 1: + cube[4] = label + octree_labeling(2, label, cube) + octree_labeling(3, label, cube) + octree_labeling(4, label, cube) + if cube[9] == 1: + cube[9] = label + octree_labeling(5, label, cube) + if cube[10] == 1: + cube[10] = label + octree_labeling(2, label, cube) + octree_labeling(5, label, cube) + octree_labeling(6, label, cube) + if cube[12] == 1: + cube[12] = label + octree_labeling(3, label, cube) + octree_labeling(5, label, cube) + octree_labeling(7, label, cube) + + if octant == 2: + if cube[1] == 1: + cube[1] = label + octree_labeling(1, label, cube) + if cube[4] == 1: + cube[4] = label + octree_labeling(1, label, cube) + octree_labeling(3, label, cube) + octree_labeling(4, label, cube) + if cube[10] == 1: + cube[10] = label + octree_labeling(1, label, cube) + octree_labeling(5, label, cube) + octree_labeling(6, label, cube) + if cube[2] == 1: + cube[2] = label + if cube[5] == 1: + cube[5] = label + octree_labeling(4, label, cube) + if cube[11] == 1: + cube[11] = label + octree_labeling(6, label, cube) + if cube[13] == 1: + cube[13] = label + octree_labeling(4, label, cube) + octree_labeling(6, label, cube) + octree_labeling(8, label, cube) + + if octant ==3: + if cube[3] == 1: + cube[3] = label + octree_labeling(1, label, cube) + if cube[4] == 1: + cube[4] = label + octree_labeling(1, label, cube) + octree_labeling(2, label, cube) + octree_labeling(4, label, cube) + if cube[12] == 1: + cube[12] = label + octree_labeling(1, label, cube) + octree_labeling(5, label, cube) + octree_labeling(7, label, cube) + if cube[6] == 1: + cube[6] = label + if cube[7] == 1: + cube[7] = label + octree_labeling(4, label, cube) + if cube[14] == 1: + cube[14] = label + octree_labeling(7, label, cube) + if cube[15] == 1: + cube[15] = label + octree_labeling(4, label, cube) + octree_labeling(7, label, cube) + octree_labeling(8, label, cube) + + if octant == 4: + if cube[4] == 1: + cube[4] = label + octree_labeling(1, label, cube) + octree_labeling(2, label, cube) + octree_labeling(3, label, cube) + if cube[5] == 1: + cube[5] = label + octree_labeling(2, label, cube) + if cube[13] == 1: + cube[13] = label + octree_labeling(2, label, cube) + octree_labeling(6, label, cube) + octree_labeling(8, label, cube) + if cube[7] == 1: + cube[7] = label + octree_labeling(3, label, cube) + if cube[15] == 1: + cube[15] = label + octree_labeling(3, label, cube) + octree_labeling(7, label, cube) + octree_labeling(8, label, cube) + if cube[8] == 1: + cube[8] = label + if cube[16] == 1: + cube[16] = label + octree_labeling(8, label, cube) + + if octant == 5: + if cube[9] == 1: + cube[9] = label + octree_labeling(1, label, cube) + if cube[10] == 1: + cube[10] = label + octree_labeling(1, label, cube) + octree_labeling(2, label, cube) + octree_labeling(6, label, cube) + if cube[12] == 1: + cube[12] = label + octree_labeling(1, label, cube) + octree_labeling(3, label, cube) + octree_labeling(7, label, cube) + if cube[17] == 1: + cube[17] = label + if cube[18] == 1: + cube[18] = label + octree_labeling(6, label, cube) + if cube[20] == 1: + cube[20] = label + octree_labeling(7, label, cube) + if cube[21] == 1: + cube[21] = label + octree_labeling(6, label, cube) + octree_labeling(7, label, cube) + octree_labeling(8, label, cube) + + if octant == 6: + if cube[10] == 1: + cube[10] = label + octree_labeling(1, label, cube) + octree_labeling(2, label, cube) + octree_labeling(5, label, cube) + if cube[11] == 1: + cube[11] = label + octree_labeling(2, label, cube) + if cube[13] == 1: + cube[13] = label + octree_labeling(2, label, cube) + octree_labeling(4, label, cube) + octree_labeling(8, label, cube) + if cube[18] == 1: + cube[18] = label + octree_labeling(5, label, cube) + if cube[21] == 1: + cube[21] = label + octree_labeling(5, label, cube) + octree_labeling(7, label, cube) + octree_labeling(8, label, cube) + if cube[19] == 1: + cube[19] = label + if cube[22] == 1: + cube[22] = label + octree_labeling(8, label, cube) + + if octant == 7: + if cube[12] == 1: + cube[12] = label + octree_labeling(1, label, cube) + octree_labeling(3, label, cube) + octree_labeling(5, label, cube) + if cube[14] == 1: + cube[14] = label + octree_labeling(3, label, cube) + if cube[15] == 1: + cube[15] = label + octree_labeling(3, label, cube) + octree_labeling(4, label, cube) + octree_labeling(8, label, cube) + if cube[20] == 1: + cube[20] = label + octree_labeling(5, label, cube) + if cube[21] == 1: + cube[21] = label + octree_labeling(5, label, cube) + octree_labeling(6, label, cube) + octree_labeling(8, label, cube) + if cube[23] == 1: + cube[23] = label + if cube[24] == 1: + cube[24] = label + octree_labeling(8, label, cube) + + if octant == 8: + if cube[13] == 1: + cube[13] = label + octree_labeling(2, label, cube) + octree_labeling(4, label, cube) + octree_labeling(6, label, cube) + if cube[15] == 1: + cube[15] = label + octree_labeling(3, label, cube) + octree_labeling(4, label, cube) + octree_labeling(7, label, cube) + if cube[16] == 1: + cube[16] = label + octree_labeling(4, label, cube) + if cube[21] == 1: + cube[21] = label + octree_labeling(5, label, cube) + octree_labeling(6, label, cube) + octree_labeling(7, label, cube) + if cube[22] == 1: + cube[22] = label + octree_labeling(6, label, cube) + if cube[24] == 1: + cube[24] = label + octree_labeling(7, label, cube) + if cube[25] == 1: + cube[25] = label + + +def _loop_through(img, curr_border): + """Inner loop of compute_thin_image. + + return simple_border_points as a list to be rechecked sequentially. + """ + # loop through the image + # NB: each loop is from 1 to size-1: img is padded from all sides + simple_border_points = [] + + ### XXX: 2D images + ### if the original is 2D, img.shape[0] == 3, the algorithm removes too much + ### because all points are considered 'boundary' in the 3rd direction. + ### Hence just bail out + if img.shape[0] == 3 and curr_border in (5, 6): + print("skipping curr_border = ", curr_border) + return [] + + for p in range(1, img.shape[0] - 1): + for r in range(1, img.shape[1] - 1): + for c in range(1, img.shape[2] - 1): + + # check if pixel is foreground + if img[p, r, c] != 1: + continue + + is_border_pt = (curr_border == 1 and img[p, r, c-1] <= 0 or #N + curr_border == 2 and img[p, r, c+1] <= 0 or #S + curr_border == 3 and img[p, r+1, c] <= 0 or #E + curr_border == 4 and img[p, r-1, c] <= 0 or #W + curr_border == 5 and img[p+1, r, c] <= 0 or #U + curr_border == 6 and img[p-1, r, c] <= 0) #B + if not is_border_pt: + # current point is not deletable + continue + + neighborhood = get_neighborhood(img, p, r, c) + + # check if (p, r, c) is an endpoint. An endpoint has exactly + # one neighbor in the 26-neighborhood. + # The center pixel is counted, thus r.h.s. is 2 + if neighborhood.sum() == 2: + continue + + # check if point is Euler invariant (condition 1 in [Lee94]) + # if it is not, it's not deletable + if not is_Euler_invariant(neighborhood): + continue + + # check if point is simple (i.e., deletion does not + # change connectivity in the 3x3x3 neighborhood) + # this are conditions 2 and 3 in [Lee94] + if not is_simple_point(neighborhood): + continue + + # ok, add (p, r, c) to the list of simple border points + simple_border_points.append((p, r, c)) + return simple_border_points + + +def compute_thin_image(img_in): + + ### prepare + img = _prepare_image(img_in) + + ### compute + unchanged_borders = 0 + + # loop through the image several times until there is no change for all + # the six border types + while unchanged_borders < 6: + unchanged_borders = 0 + for curr_border in (4, 3, 2, 1, 5, 6): + + simple_border_points = _loop_through(img, curr_border) + print(curr_border, " : ", simple_border_points, '\n') + + # sequential re-checking to preserve connectivity when deleting + # in a parallel way + no_change = True + for pt in simple_border_points: + p, r, c = pt + neighb = get_neighborhood(img, p, r, c) + if is_simple_point(neighb): + img[p, r, c] = 0 + no_change = False + else: + print(" *** ", pt, is_simple_point(neighb)) + + if no_change: + unchanged_borders += 1 + simple_border_points = [] + + img = _postprocess_image(img) + return img + + +if __name__ == "__main__": + pass diff --git a/skimage/morphology/test_simple_2d.py b/skimage/morphology/test_simple_2d.py new file mode 100644 index 00000000..c4a5cd2b --- /dev/null +++ b/skimage/morphology/test_simple_2d.py @@ -0,0 +1,84 @@ +from __future__ import division, print_function, absolute_import + +import numpy as np +from numpy.testing import assert_equal + +from skimage import io + +from skel import compute_thin_image + +import matplotlib.pyplot as plt +import matplotlib.ticker as ticker + +# nose test generators: + +# 2D images +def test_simple_2d_images(): + for fname in ("strip", "loop", "cross", "two-hole"): + yield check_skel, fname + +# trivial 3D images +def test_simple_3d(): + for fname in ['3/stack', '4/stack']: + yield check_skel_3d, fname + + +def check_skel(fname, viz=False): + # compute the thin image and compare the result to that of ImageJ + img = np.loadtxt('data/' + fname + '.txt', dtype=np.uint8) + + if viz: + ax = _viz(img, **dict(marker='s', color='b', s=99, alpha=0.2)) + + # compute + img1_2d = compute_thin_image(img) + + if viz: + ax = _viz(img1_2d, ax, **dict(marker='o', color='r', + s=80, alpha=0.7, label='us')) + + # compare to FIJI + img_f = np.loadtxt('data/' + fname + '_fiji.txt', dtype=np.uint8) + + if not viz: + # actually compare images + assert_equal(img1_2d, img_f) + else: + ax = _viz(img_f, ax, **dict(marker='o', color='g', s=45, label='fiji')) + + ax.legend() + ax.grid(True) + + def yformatter(val, pos): + return int(img.shape[1] - val + 1) + def xformatter(val, pos): + return int(val + 1) + ax.xaxis.set_major_formatter(ticker.FuncFormatter(xformatter)) + ax.yaxis.set_major_formatter(ticker.FuncFormatter(yformatter)) + + plt.show() + + +def _viz(img, ax=None, **kwds): + if ax is None: + import matplotlib.pyplot as plt + fix, ax = plt.subplots() + + x, y = np.nonzero(img) + ax.scatter(y, img.shape[1] - x, **kwds) + return ax + + +def check_skel_3d(fname): + img = io.imread('data/' + fname + '.tif') + img_f = io.imread('data/' + fname + '_fiji.tif') + + img_s = compute_thin_image(img) + assert_equal(img_s, img_f) + + +if __name__ == "__main__": + import sys + if len(sys.argv) < 2: + sys.exit("Expect an image name from the data/ directory.") + check_skel(sys.argv[1], True) diff --git a/skimage/morphology/trivial.py b/skimage/morphology/trivial.py new file mode 100644 index 00000000..46656d93 --- /dev/null +++ b/skimage/morphology/trivial.py @@ -0,0 +1,73 @@ +import numpy as np +import matplotlib.pyplot as plt + +from skel import prepare_image, compute_thin_image + + +def get_rhombus(n=64, L=22, width=7): + img = np.zeros((n, n), dtype=np.int8) + + x = np.arange(L, dtype=int) + y = L - x + + for w in range(width): + img[x + n//2, y + n//2 + w] = 1 + img[-x + n//2, y + n//2 + w] = 1 + img[x + n//2, -y + n//2 + w] = 1 + img[-x + n//2, -y + n//2 + w] = 1 + + return img + +def get_strip(): + img = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=np.uint8) + return img + + +def get_square(n=64, L=16, width=4): + img = np.zeros((n, n), dtype=np.int8) + + x = np.arange(-L//2 + 1, L//2) + for w in range(width): + img[n//2 + x, n//2 + L//2 + w] = 1 + img[n//2 + x, n//2 - L//2 - w] = 1 + + img[n//2 + L//2 - w, n//2 + x] = 1 + img[n//2 + x, n//2 - L//2 - w] = 1 + + return img + + +def get_loop(): + img = np.loadtxt('Untitled.txt', dtype=np.uint8) + return img + + +if __name__ == "__main__": + +## img = get_rhombus() +## img = get_strip() + img = get_loop() + + x, y = np.nonzero(img) + plt.scatter(x, y, marker='s', color='b', s=40, alpha=0.3) + + # skeletonize + img1 = prepare_image(img) + img1 = compute_thin_image(img1) + + img1_2d = img1[1, 1:, 1:] + x, y = np.nonzero(img1_2d) + + plt.scatter(x, y, marker='o', color='r') + plt.show()