mirror of
https://github.com/wassname/scikit-image.git
synced 2026-07-20 12:40:31 +08:00
REF: Convert find_contours to Cython.
This commit is contained in:
Stefan van der Walt
@@ -1,262 +0,0 @@
|
||||
#include <Python.h>
|
||||
#include "numpy/arrayobject.h"
|
||||
|
||||
static char _find_contours_doc[] =
|
||||
"This module defines C helper functions for find_contours";
|
||||
|
||||
|
||||
static char iterate_and_store_doc[] =
|
||||
"iterate_and_store(array, level, vertex_connect_high)\n\
|
||||
\n\
|
||||
Iterate across the given array in a marching-squares fashion, looking for\n\
|
||||
segments that cross 'level'. If such a segment is found, its coordinates are\n\
|
||||
added to a growing list of segments, which is returned by the function.\n\
|
||||
if vertex_connect_high is nonzero, high-values pixels are considered to be\n\
|
||||
face+vertex connected into objects; otherwise low-valued pixels are.";
|
||||
|
||||
|
||||
// Nasty macros to inline the inner loop of interpolating the position of
|
||||
// the contour and add an appropriate tuple to the output list.
|
||||
// These macros define blocks of code that can only be used within the inner
|
||||
// loop of 'iterate_and_store' because they use variables therefrom.
|
||||
|
||||
#define GET_FRACTION(from_value, to_value) \
|
||||
((level - from_value) / (to_value - from_value))
|
||||
|
||||
#define TOP { \
|
||||
output0 = coords[0]; \
|
||||
output1 = coords[1] + GET_FRACTION(*ul_ptr, *ur_ptr); \
|
||||
}
|
||||
|
||||
#define BOTTOM { \
|
||||
output0 = coords[0] + 1; \
|
||||
output1 = coords[1] + GET_FRACTION(*ll_ptr, *lr_ptr); \
|
||||
}
|
||||
|
||||
#define LEFT { \
|
||||
output0 = coords[0] + GET_FRACTION(*ul_ptr, *ll_ptr); \
|
||||
output1 = coords[1]; \
|
||||
}
|
||||
|
||||
#define RIGHT { \
|
||||
output0 = coords[0] + GET_FRACTION(*ur_ptr, *lr_ptr); \
|
||||
output1 = coords[1] + 1; \
|
||||
}
|
||||
|
||||
#define ADD_TUPLE { \
|
||||
PyObject* tuple = Py_BuildValue("(dd)", output0, output1); \
|
||||
if (!tuple) { \
|
||||
Py_DECREF(double_array); \
|
||||
Py_DECREF(arc_list); \
|
||||
return NULL; \
|
||||
} \
|
||||
char res = PyList_Append(arc_list, tuple); \
|
||||
Py_DECREF(tuple); \
|
||||
if (res < 0) { \
|
||||
Py_DECREF(double_array); \
|
||||
Py_DECREF(arc_list); \
|
||||
return NULL; \
|
||||
} \
|
||||
}
|
||||
|
||||
#define ADD_SEGMENT(START, END) { \
|
||||
double output0, output1; \
|
||||
START \
|
||||
ADD_TUPLE \
|
||||
END \
|
||||
ADD_TUPLE \
|
||||
}
|
||||
|
||||
static PyObject*
|
||||
iterate_and_store(PyObject *self, PyObject *args)
|
||||
{
|
||||
PyObject* array;
|
||||
double level;
|
||||
int vertex_connect_high;
|
||||
if (!PyArg_ParseTuple(args, "Odi:iterate_and_store", &array, &level,
|
||||
&vertex_connect_high)) {
|
||||
return NULL;
|
||||
}
|
||||
|
||||
PyObject* double_array = PyArray_FromAny(array,
|
||||
PyArray_DescrFromType(NPY_DOUBLE), 2, 2, NPY_CONTIGUOUS | NPY_ALIGNED,
|
||||
NULL);
|
||||
if (!double_array) {
|
||||
return NULL;
|
||||
}
|
||||
|
||||
npy_intp *dims = PyArray_DIMS(double_array);
|
||||
if (dims[0] < 2 || dims[1] < 2) {
|
||||
Py_DECREF(double_array);
|
||||
PyErr_SetString(PyExc_ValueError, "Input array must be at least 2x2.");
|
||||
return NULL;
|
||||
}
|
||||
|
||||
// The plan is to iterate a 2x2 square across the input array. This means
|
||||
// that the upper-left corner of the square needs to iterate across a
|
||||
// sub-array that's one-less-large in each direction (so that the square
|
||||
// never steps out of bounds). The square is represented by four pointers:
|
||||
// ul, ur, ll, and lr (for 'upper left', etc.). We also maintain the current
|
||||
// 2D coordinates for the position of the upper-left pointer. Note that we
|
||||
// ensured that the array is of type 'double' and is C-contiguous (last
|
||||
// index varies the fastest).
|
||||
|
||||
// Current coords start at 0,0.
|
||||
npy_intp coords[2] = {0,0};
|
||||
// Precompute the size of the array minus 2 in each direction, so we'll know
|
||||
// when to update the coordinates and double-increment the square pointers
|
||||
// so that the upper-left pointer never visits the last column.
|
||||
npy_intp dims_m2[2];
|
||||
dims_m2[0] = dims[0] - 2;
|
||||
dims_m2[1] = dims[1] - 2;
|
||||
// Calculate the number of iterations we'll need
|
||||
npy_intp num_square_steps = (dims[0] - 1) * (dims[1] - 1);
|
||||
// and set up the square pointers.
|
||||
double* ul_ptr = PyArray_DATA(double_array);
|
||||
double* ur_ptr = ul_ptr + 1;
|
||||
double* ll_ptr = ul_ptr + dims[1];
|
||||
double* lr_ptr = ll_ptr + 1;
|
||||
|
||||
// make a list to hold the returned coordinates
|
||||
PyObject* arc_list = PyList_New(0);
|
||||
if (!arc_list) {
|
||||
Py_DECREF(double_array);
|
||||
return NULL;
|
||||
}
|
||||
while(num_square_steps--) {
|
||||
// There are sixteen different possible square types, diagramed below.
|
||||
// A + indicates that the vertex is above the contour value, and a -
|
||||
// indicates that the vertex is below or equal to the contour value.
|
||||
// The vertices of each square are:
|
||||
// ul ur
|
||||
// ll lr
|
||||
// and can be treated as a binary value with the bits in that order. Thus
|
||||
// each square case can be numbered:
|
||||
// 0-- 1+- 2-+ 3++ 4-- 5+- 6-+ 7++
|
||||
// -- -- -- -- +- +- +- +-
|
||||
//
|
||||
// 8-- 9+- 10-+ 11++ 12-- 13+- 14-+ 15++
|
||||
// -+ -+ -+ -+ ++ ++ ++ ++
|
||||
//
|
||||
// The position of the line segment that cuts through (or doesn't, in case
|
||||
// 0 and 15) each square is clear, except in cases 6 and 9. In this case,
|
||||
// where the segments are placed is determined by vertex_connect_high.
|
||||
// If vertex_connect_high is false, then lines like \\ are drawn
|
||||
// through square 6, and lines like // are drawn through square 9.
|
||||
// Otherwise, the situation is reversed.
|
||||
// Finally, recall that we draw the lines so that (moving from tail to
|
||||
// head) the lower-valued pixels are on the left of the line. So, for
|
||||
// example, case 1 entails a line slanting from the middle of the top of
|
||||
// the square to the middle of the left side of the square.
|
||||
|
||||
unsigned char square_case = 0;
|
||||
if ((*ul_ptr) > level) square_case += 1;
|
||||
if ((*ur_ptr) > level) square_case += 2;
|
||||
if ((*ll_ptr) > level) square_case += 4;
|
||||
if ((*lr_ptr) > level) square_case += 8;
|
||||
|
||||
switch(square_case)
|
||||
{
|
||||
case 0: // no line
|
||||
break;
|
||||
case 1: // top to left
|
||||
ADD_SEGMENT(TOP, LEFT);
|
||||
break;
|
||||
case 2: // right to top
|
||||
ADD_SEGMENT(RIGHT, TOP);
|
||||
break;
|
||||
case 3: // right to left
|
||||
ADD_SEGMENT(RIGHT, LEFT);
|
||||
break;
|
||||
case 4: // left to bottom
|
||||
ADD_SEGMENT(LEFT, BOTTOM);
|
||||
break;
|
||||
case 5: // top to bottom
|
||||
ADD_SEGMENT(TOP, BOTTOM);
|
||||
break;
|
||||
case 6:
|
||||
if (vertex_connect_high)
|
||||
{
|
||||
// left to top
|
||||
ADD_SEGMENT(LEFT, TOP);
|
||||
// right to bottom
|
||||
ADD_SEGMENT(RIGHT, BOTTOM);
|
||||
}
|
||||
else
|
||||
{
|
||||
// right to top
|
||||
ADD_SEGMENT(RIGHT, TOP);
|
||||
// left to bottom
|
||||
ADD_SEGMENT(LEFT, BOTTOM);
|
||||
}
|
||||
break;
|
||||
case 7: // right to bottom
|
||||
ADD_SEGMENT(RIGHT, BOTTOM);
|
||||
break;
|
||||
case 8: // bottom to right
|
||||
ADD_SEGMENT(BOTTOM, RIGHT);
|
||||
break;
|
||||
case 9:
|
||||
if (vertex_connect_high)
|
||||
{
|
||||
// top to right
|
||||
ADD_SEGMENT(TOP, RIGHT);
|
||||
// bottom to left
|
||||
ADD_SEGMENT(BOTTOM, LEFT);
|
||||
}
|
||||
else
|
||||
{
|
||||
// top to left
|
||||
ADD_SEGMENT(TOP, LEFT);
|
||||
// bottom to right
|
||||
ADD_SEGMENT(BOTTOM, RIGHT);
|
||||
}
|
||||
break;
|
||||
case 10: // bottom to top
|
||||
ADD_SEGMENT(BOTTOM, TOP);
|
||||
break;
|
||||
case 11: // bottom to left
|
||||
ADD_SEGMENT(BOTTOM, LEFT);
|
||||
break;
|
||||
case 12: // left to right
|
||||
ADD_SEGMENT(LEFT, RIGHT);
|
||||
break;
|
||||
case 13: // top to right
|
||||
ADD_SEGMENT(TOP, RIGHT);
|
||||
break;
|
||||
case 14: // left to top
|
||||
ADD_SEGMENT(LEFT, TOP);
|
||||
break;
|
||||
case 15: // no line
|
||||
break;
|
||||
} // switch square_case
|
||||
|
||||
if (coords[1] < dims_m2[1]) {
|
||||
coords[1]++;
|
||||
} else {
|
||||
coords[1] = 0;
|
||||
coords[0]++;
|
||||
// Double-increment pointers to advance them to the next row, since
|
||||
// we're skipping the last column, as far as ul_ptr is concerned.
|
||||
ul_ptr++; ur_ptr++; ll_ptr++; lr_ptr++;
|
||||
}
|
||||
ul_ptr++; ur_ptr++; ll_ptr++; lr_ptr++;
|
||||
} // iteration
|
||||
|
||||
// get rid of the double array reference that we own
|
||||
Py_DECREF(double_array);
|
||||
return arc_list;
|
||||
}
|
||||
|
||||
|
||||
static PyMethodDef _find_contours_methods[] = {
|
||||
{"iterate_and_store", iterate_and_store, METH_VARARGS, iterate_and_store_doc},
|
||||
{NULL, NULL, 0, NULL}
|
||||
};
|
||||
|
||||
PyMODINIT_FUNC
|
||||
init_find_contours(void)
|
||||
{
|
||||
Py_InitModule3("_find_contours", _find_contours_methods, _find_contours_doc);
|
||||
import_array();
|
||||
}
|
||||
@@ -0,0 +1,185 @@
|
||||
# -*- python -*-
|
||||
# cython: cdivision=True
|
||||
|
||||
import numpy as np
|
||||
cimport numpy as np
|
||||
|
||||
np.import_array()
|
||||
|
||||
cdef double _get_fraction(double from_value, double to_value, double level):
|
||||
if (to_value == from_value):
|
||||
return 0
|
||||
return ((level - from_value) / (to_value - from_value))
|
||||
|
||||
|
||||
def iterate_and_store(np.ndarray[double, ndim=2, mode='c'] array,
|
||||
double level, int vertex_connect_high):
|
||||
"""Iterate across the given array in a marching-squares fashion,
|
||||
looking for segments that cross 'level'. If such a segment is
|
||||
found, its coordinates are added to a growing list of segments,
|
||||
which is returned by the function. if vertex_connect_high is
|
||||
nonzero, high-values pixels are considered to be face+vertex
|
||||
connected into objects; otherwise low-valued pixels are.
|
||||
|
||||
"""
|
||||
if array.shape[0] < 2 or array.shape[1] < 2:
|
||||
raise ValueError("Input array must be at least 2x2.")
|
||||
|
||||
cdef list arc_list = []
|
||||
cdef int n
|
||||
|
||||
# The plan is to iterate a 2x2 square across the input array. This means
|
||||
# that the upper-left corner of the square needs to iterate across a
|
||||
# sub-array that's one-less-large in each direction (so that the square
|
||||
# never steps out of bounds). The square is represented by four pointers:
|
||||
# ul, ur, ll, and lr (for 'upper left', etc.). We also maintain the current
|
||||
# 2D coordinates for the position of the upper-left pointer. Note that we
|
||||
# ensured that the array is of type 'double' and is C-contiguous (last
|
||||
# index varies the fastest).
|
||||
|
||||
# Current coords start at 0,0.
|
||||
cdef int[2] coords
|
||||
coords[0] = 0
|
||||
coords[1] = 0
|
||||
|
||||
# Calculate the number of iterations we'll need
|
||||
cdef int num_square_steps = (array.shape[0] - 1) * (array.shape[1] - 1)
|
||||
|
||||
cdef unsigned char square_case = 0
|
||||
cdef tuple top, bottom, left, right
|
||||
cdef double ul, ur, ll, lr
|
||||
cdef int r0, r1, c0, c1
|
||||
|
||||
for n in range(num_square_steps):
|
||||
# There are sixteen different possible square types, diagramed below.
|
||||
# A + indicates that the vertex is above the contour value, and a -
|
||||
# indicates that the vertex is below or equal to the contour value.
|
||||
# The vertices of each square are:
|
||||
# ul ur
|
||||
# ll lr
|
||||
# and can be treated as a binary value with the bits in that order. Thus
|
||||
# each square case can be numbered:
|
||||
# 0-- 1+- 2-+ 3++ 4-- 5+- 6-+ 7++
|
||||
# -- -- -- -- +- +- +- +-
|
||||
#
|
||||
# 8-- 9+- 10-+ 11++ 12-- 13+- 14-+ 15++
|
||||
# -+ -+ -+ -+ ++ ++ ++ ++
|
||||
#
|
||||
# The position of the line segment that cuts through (or
|
||||
# doesn't, in case 0 and 15) each square is clear, except in
|
||||
# cases 6 and 9. In this case, where the segments are placed
|
||||
# is determined by vertex_connect_high. If
|
||||
# vertex_connect_high is false, then lines like \\ are drawn
|
||||
# through square 6, and lines like # are drawn through square
|
||||
# 9. Otherwise, the situation is reversed.
|
||||
# Finally, recall that we draw the lines so that (moving from tail to
|
||||
# head) the lower-valued pixels are on the left of the line. So, for
|
||||
# example, case 1 entails a line slanting from the middle of the top of
|
||||
# the square to the middle of the left side of the square.
|
||||
|
||||
r0, c0 = coords[0], coords[1]
|
||||
r1, c1 = r0 + 1, c0 + 1
|
||||
|
||||
ul = array[r0, c0]
|
||||
ur = array[r0, c0 + 1]
|
||||
ll = array[r0 + 1, c0]
|
||||
lr = array[r0 + 1, c0 + 1]
|
||||
|
||||
square_case = 0
|
||||
if (ul > level): square_case += 1
|
||||
if (ur > level): square_case += 2
|
||||
if (ll > level): square_case += 4
|
||||
if (lr > level): square_case += 8
|
||||
|
||||
top = coords[0], coords[1] + _get_fraction(ul, ur, level)
|
||||
bottom = coords[0] + 1, coords[1] + _get_fraction(ll, lr, level)
|
||||
left = coords[0] + _get_fraction(ul, ll, level), coords[1]
|
||||
right = coords[0] + _get_fraction(ur, lr, level), coords[1] + 1
|
||||
|
||||
if (square_case == 0):
|
||||
# no line
|
||||
pass
|
||||
elif (square_case == 1):
|
||||
# top to left
|
||||
arc_list.append(top)
|
||||
arc_list.append(left)
|
||||
elif (square_case == 2):
|
||||
# right to top
|
||||
arc_list.append(right)
|
||||
arc_list.append(top)
|
||||
elif (square_case == 3):
|
||||
# right to left
|
||||
arc_list.append(right)
|
||||
arc_list.append(left)
|
||||
elif (square_case == 4):
|
||||
# left to bottom
|
||||
arc_list.append(left)
|
||||
arc_list.append(bottom)
|
||||
elif (square_case == 5):
|
||||
# top to bottom
|
||||
arc_list.append(top)
|
||||
arc_list.append(bottom)
|
||||
elif (square_case == 6):
|
||||
if vertex_connect_high:
|
||||
arc_list.append(left)
|
||||
arc_list.append(top)
|
||||
|
||||
arc_list.append(right)
|
||||
arc_list.append(bottom)
|
||||
else:
|
||||
arc_list.append(right)
|
||||
arc_list.append(top)
|
||||
arc_list.append(left)
|
||||
arc_list.append(bottom)
|
||||
elif (square_case == 7):
|
||||
# right to bottom
|
||||
arc_list.append(right)
|
||||
arc_list.append(bottom)
|
||||
elif (square_case == 8):
|
||||
# bottom to right
|
||||
arc_list.append(bottom)
|
||||
arc_list.append(right)
|
||||
elif (square_case == 9):
|
||||
if vertex_connect_high:
|
||||
arc_list.append(top)
|
||||
arc_list.append(right)
|
||||
|
||||
arc_list.append(bottom)
|
||||
arc_list.append(left)
|
||||
else:
|
||||
arc_list.append(top)
|
||||
arc_list.append(left)
|
||||
|
||||
arc_list.append(bottom)
|
||||
arc_list.append(right)
|
||||
elif (square_case == 10):
|
||||
# bottom to top
|
||||
arc_list.append(bottom)
|
||||
arc_list.append(top)
|
||||
elif (square_case == 11):
|
||||
# bottom to left
|
||||
arc_list.append(bottom)
|
||||
arc_list.append(left)
|
||||
elif (square_case == 12):
|
||||
# lef to right
|
||||
arc_list.append(left)
|
||||
arc_list.append(right)
|
||||
elif (square_case == 13):
|
||||
# top to right
|
||||
arc_list.append(top)
|
||||
arc_list.append(right)
|
||||
elif (square_case == 14):
|
||||
# left to top
|
||||
arc_list.append(left)
|
||||
arc_list.append(top)
|
||||
elif (square_case == 15):
|
||||
# no line
|
||||
pass
|
||||
|
||||
if coords[1] < array.shape[1] - 2:
|
||||
coords[1] += 1
|
||||
else:
|
||||
coords[0] += 1
|
||||
coords[1] = 0
|
||||
|
||||
return arc_list
|
||||
@@ -84,7 +84,7 @@ def find_contours(array, level, fully_connected='low', positive_orientation='low
|
||||
around structures that are a single array element wide. Instead choose
|
||||
a middle value, as above.'''
|
||||
|
||||
array = np.asarray(array)
|
||||
array = np.asarray(array, dtype=np.double)
|
||||
if array.ndim != 2:
|
||||
raise RuntimeError('Only 2D arrays are supported.')
|
||||
level = float(level)
|
||||
|
||||
@@ -1,11 +1,17 @@
|
||||
#!/usr/bin/env python
|
||||
|
||||
from skimage._build import cython
|
||||
|
||||
import os
|
||||
base_path = os.path.abspath(os.path.dirname(__file__))
|
||||
|
||||
def configuration(parent_package='', top_path=None):
|
||||
from numpy.distutils.misc_util import Configuration, get_numpy_include_dirs
|
||||
|
||||
config = Configuration('measure', parent_package, top_path)
|
||||
config.add_data_dir('tests')
|
||||
|
||||
cython(['_find_contours.pyx'], working_path=base_path)
|
||||
|
||||
config.add_extension('_find_contours', sources=['_find_contours.c'],
|
||||
include_dirs=[get_numpy_include_dirs()])
|
||||
|
||||
Reference in New Issue
Block a user