3
$\begingroup$

I would like to know the algorithm to calculate contour line length. Suppose we have numerical data set of an function $f(x,y)$.

How could I calculate the length of the line from $(x_1,y_1)$ to $(x_2,y_2)$ on the contour $f = f_0$? I am working on this with Fortran, but the example with any other language will be OK.

Improve this question
$\endgroup$
5
  • $\begingroup$ Could you please confirm that you do not already know the contour line (i.e. a set of points along the contour). @nicoguaro's answer assumes you already know the contour. Also, there are many descriptions online for how to compute contour lines and how to compute arc length. What have you tried and what exactly have you gotten stuck on? $\endgroup$
    – Doug Lipinski
    Apr 16, 2015 at 22:40
  • $\begingroup$ Thank you guys. As @Doug Lipinski inferred, the analytic form of the contour line is not given in my case. The all we have is just the data set of ($x_i, y_i, f_{ij}$). As a first approximation, what I did was starting from $(x_0, y_0)$, search the grid point (let's say $(x_0, y_0)$) surrounding $x_0$ with the closest $f_0$. Then calculate $\Delta r= \sqrt{(x_0-x_1)^2+(y_0-y_1)^2}$, and collect all $\Delta r$. I know this algorithm is horrible, But I have no idea what should I do. $\endgroup$
    – user1048419
    Apr 17, 2015 at 18:41
  • 1
    $\begingroup$ There are two steps you'll need to do: 1) find the contour 2) compute the arc length. As nicoguaro mentioned in a comment, marching squares is a good starting point for computing the contour. Higher order approximations of where the contour crosses your grid can increase accuracy. nicoguaro's answer provides a good overview of integrating the arc length along the curve. $\endgroup$
    – Doug Lipinski
    Apr 18, 2015 at 0:30
  • $\begingroup$ I have never heard of the marching squares. Thank you. I will try that! $\endgroup$
    – user1048419
    Apr 18, 2015 at 18:16
  • $\begingroup$ I expanded the answer to cover a brief description of the steps for Marching squares. $\endgroup$
    – nicoguaro ♦
    Apr 21, 2015 at 20:32

2 Answers 2

Reset to default
3
$\begingroup$

As Doug Lipinski mentioned you need to:

1) find the contour
2) compute the arc length.

Find the contour

To find the contour you can use Marching Squares, where you rewrite

$$f(x,y) = f_0 \Leftrightarrow f(x,y) - f_0 = 0 \enspace .$$

You assign a primitive to each individual cell (each 4 cells can be considered to be the corners of a bilinear interpolation, for example). Thus, you consider the sign of the values at vertices; intersections occur on edges with sign change. There are 16 options for this, although just 4 topologically different

  • No crossing: 4 edges with same sign
  • Singlet: 1 edge with different sign
  • Double adjacent
  • Double opposite

In the double opposite case might be some ambiguities and the sign of the saddle point should be computed.

As a last step you determine the position of intersection (using interpolation along grid edges).

Compute the arc length

You want to perform the numeric integration to obtain the arc length, i.e.,

$$s = \int_{a}^{b} \sqrt { [x'(t)]^2 + [y'(t)]^2 }\, dt. $$

where $t$ is a parameter for $x(t)$ and $y(t)$.

Since your input is the numerical data for $x$ and $y$, you can directly proceed with the numerical integration. The simplest case is to connect your points with lines and add all those segments

$$s \approx \sum_{k=1}^{N} \sqrt{ (x_{k} - x_{k-1})^2 + (y_{k} - y_{k-1})^2}$$

An example in Python

import numpy as np

def arc_length(x, y):
    npts = len(x)
    arc = np.sqrt((x[1] - x[0])**2 + (y[1] - y[0])**2)
    for k in range(1, npts):
        arc = arc + np.sqrt((x[k] - x[k-1])**2 + (y[k] - y[k-1])**2)

    return arc


# Parabolic segment
npts = 1000
a = 3.
h = 5.
x = np.linspace(-a, a, npts)
y = h*(1 - x**2/a**2)
analytic = np.sqrt(a**2 + 4*h**2) + a**2/(2*h)*np.arcsinh(2*h/a)
numeric = arc_length(x, y)
print analytic, numeric

# Semicircle
npts = 1000
x = np.linspace(-1, 1, npts)
y = np.sqrt(1 - x**2)
analytic = np.pi
numeric = arc_length(x, y)
print analytic, numeric

With output

12.1673133338 12.1881924551
3.14159265359 3.20484351644

You can use higher order derivatives (interpolations or B-Splines, as suggested by Model_Math) for more accurate results, but the concept is the same.

Improve this answer
$\endgroup$
2
  • $\begingroup$ The question says that the OP has numerical data for a function $f(x,y)$, not that he/she has numerical data for the contour line (perhaps he/she could confirm that?). Computing the contour line is much more challenging than just computing arc length of a known curve. $\endgroup$
    – Doug Lipinski
    Apr 16, 2015 at 22:38
  • $\begingroup$ You are right, my bad. But you can use as an intermediate step marching squares to obtain the data for your contour. Let's wait for the response and depending on that I will add the other step. $\endgroup$
    – nicoguaro ♦
    Apr 16, 2015 at 22:53
2
$\begingroup$

Even if I knew the contour in terms of points where $f_0$ was constant to some useful approximation, I would probably still use B-Splines to adaptively refine the point sequence to converge to an arc length. (If that's gobblety gook, the reference I'm about to give will cover this among more advanced spline details.)

Consider reading (may require purchase) Contour reconstruction using recursive smoothing splines - Algorithms and experimental validation from Robotics and Autonomous Systems Volume 57, Issues 6–7, 30 June 2009, Pages 617–628

That article will either answer your question or point you in more elementary directions if they better suit you.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.