How to convert the following trigonometric function into a python function for boundary curvature calculation?

Viewed 72

I am trying to understand a few slides from this source

Specifically, this example at slide 59:
enter image description here
The part I do not understand is how to go from the chain-code to the curvature.

I believe the formula is given in slide 56:
enter image description here

But if I try to implement this in python I get different results. For example:

import matplotlib.pyplot as plt
# Dataset
x = [0, 1, 2, 2, 3, 4, 5, 6, 6, 7, 8]
y = [0, 0, 0, 1, 1, 2, 2, 1, 0, 0, 0]
# Show data
plt.scatter(x, y)
plt.plot(x, y)
plt.axis('equal')
plt.show()

enter image description here

import math
i = 4  # Taking the 5th point, at index 4, with supposed curvature of 1 from the slide
k = 1
a = math.atan((y[i+k]-y[i])/(x[i+k]-x[i]))
b = math.atan((y[i]-y[i-k])/(x[i]-x[i-k]))
result = (a - b) % (2 * math.pi)  # = 0.7853981633974483

So clearly I a missing something, but what?

1 Answers

The "curvature" in the first image is the difference between two subsequent "chain-codes" modulo 8. So for example for chain codes 0 0 2 0 1 0 7 6 0 0 the 4th entry in curvature is 1-0 = 1 while the sixth is 7-0 = 7 = -1 (mod 8). In Python you can calculate it like this:

>>> def mod8(x):
...     m = x % 8
...     return m if m < 4 else m - 8
... 
>>> cc = [0, 0, 2, 0, 1, 0, 7, 6, 0, 0]
>>> [mod8(a - b) for (a, b) in zip(cc[1:], cc[:-1])]
[0, 2, -2, 1, -1, -1, -1, 2, 0]

If you compare this with the formula that uses atan, what the formula is missing is the conversion of the angles from radians to the units where 1 is 45 degrees (pi/4). Your result 0.7853981633974483 is correct according to the formula, but if you expected to get 1.0 you would have to divide the result by math.pi/4.

Related