added chapter06

Chapter06/Chapter06_mcode/functions/distribute_points.m

+function C_new = distribute_points(C)
+% DISTRIBUTE_POINTS   Distributes snake points equidistantly
+%   Author: vand@dtu.dk
+
+C = C([1:end,1],:); % closing the curve
+N = size(C,1); % number of points (+ 1, due to closing)
+dist = sqrt(sum(diff(C).^2,2)); % edge segment lengths
+t = [0;cumsum(dist)]; % current positions along the snake
+tq = linspace(0,t(end),N)'; % equidistant positions
+C_new(:,1) = interp1(t,C(:,1),tq); % distributed x
+C_new(:,2) = interp1(t,C(:,2),tq); % distributed y
+C_new = C_new(1:end-1,:); % opening the curve again

Chapter06/Chapter06_mcode/functions/remove_intersections.m

+function S = remove_intersections(S)
+%REMOVE_INTERSECTIONS   Reorder snake points to remove self-intersections.
+% S = remove_intersection(S)
+% Input: snake S represented by a N-by-2 matrix 
+% Output: snake S represented by a N-by-2 matrix 
+
+S1 = [S;S(1,:)];
+n1 = size(S1,1);
+n = n1-1;
+
+for i = 1:n1-3
+    for j = i+2:n1-1
+        if ( is_crossing(S1(i,:), S1(i+1,:), S1(j,:), S1(j+1,:)) )
+            f = i+1;
+            t = j;
+            if ( j-i > n/2 )
+                f = j+1;
+                t = i+n;
+            end
+            while ( f < t )
+                idF = mod(f,n);
+                if ( idF == 0 )
+                    idF = n;
+                end
+                f = f + 1;
+                
+                idT = mod(t,n);
+                if ( idT == 0 )
+                    idT = n;
+                end
+                t = t - 1;
+                tmp = S1(idF,:);
+                S1(idF,:) = S1(idT,:);
+                S1(idT,:) = tmp;
+            
+            end
+            S1(end,:) = S1(1,:);
+        end
+    end
+end
+
+S = S1(1:end-1,:);
+
+function is_cross = is_crossing(p1, p2, p3, p4)
+% detects crossings between a two line segments: p1 to p2 and p3 to p4
+
+is_cross = false;
+
+p21 = p2-p1;
+p34 = p3-p4;
+p31 = p3-p1;
+
+alpha = (p34(2)*p31(1)-p31(2)*p34(1))/(p34(2)*p21(1)-p21(2)*p34(1));
+if alpha>0 && alpha<1
+    beta = (p21(2)*p31(1)-p31(2)*p21(1))/(p21(2)*p34(1)-p34(2)*p21(1));
+    if beta>0 && beta<1
+        is_cross = true;
+    end
+end
\ No newline at end of file

Chapter06/snake_functions.py

+"""
+Helping functions for snakes.
+"""
+
+import numpy as np
+import scipy.interpolate
+
+def distribute_points(snake):
+    """ Distributes snake points equidistantly."""
+    N = snake.shape[1]
+    
+    # Compute length of line segments.
+    d = np.sqrt(np.sum((np.roll(snake, -1, axis=1) - snake)**2, axis=0)) 
+    f = scipy.interpolate.interp1d(np.hstack((0, np.cumsum(d))), 
+                                   np.hstack((snake, snake[:,0:1])))
+    return(f(sum(d) * np.arange(N) / N))
+
+def is_crossing(p1, p2, p3, p4):
+    """ Check if the line segments (p1, p2) and (p3, p4) cross."""
+    crossing = False
+    d21 = p2 - p1
+    d43 = p4 - p3
+    d31 = p3 - p1
+    det = d21[0]*d43[1] - d21[1]*d43[0] # Determinant
+    if det != 0.0 and d21[0] != 0.0 and d21[1] != 0.0:
+        a = d43[0]/d21[0] - d43[1]/d21[1]
+        b = d31[1]/d21[1] - d31[0]/d21[0]
+        if a != 0.0:
+            u = b/a
+            if d21[0] > 0:
+                t = (d43[0]*u + d31[0])/d21[0]
+            else:
+                t = (d43[1]*u + d31[1])/d21[1]
+            crossing = 0 < u < 1 and 0 < t < 1         
+    return crossing
+
+def is_counterclockwise(snake):
+    """ Check if points are ordered counterclockwise."""
+    return np.dot(snake[0, 1:] - snake[0, :-1],
+                  snake[1, 1:] + snake[1, :-1]) < 0
+
+def remove_intersections(snake):
+    """ Reorder snake points to remove self-intersections.
+        Arguments: snake represented by a 2-by-N array.
+        Returns: snake.
+    """
+    pad_snake = np.append(snake, snake[:,0].reshape(2,1), axis=1)
+    pad_n = pad_snake.shape[1]
+    n = pad_n - 1 
+    
+    for i in range(pad_n - 3):
+        for j in range(i + 2, pad_n - 1):
+            pts = pad_snake[:, [i, i + 1, j, j + 1]]
+            if is_crossing(pts[:, 0], pts[:, 1], pts[:, 2], pts[:, 3]):
+                # Reverse vertices of smallest loop
+                rb = i + 1 # Reverse begin
+                re = j     # Reverse end
+                if j - i > n // 2:
+                    # Other loop is smallest
+                    rb = j + 1
+                    re = i + n                    
+                while rb < re:
+                    ia = rb % n
+                    rb = rb + 1                    
+                    ib = re % n
+                    re = re - 1                    
+                    pad_snake[:, [ia, ib]] = pad_snake[:, [ib, ia]]                    
+                pad_snake[:,-1] = pad_snake[:,0]                
+    snake = pad_snake[:, :-1]
+    if is_counterclockwise(snake):
+        return snake
+    else:
+        return np.flip(snake, axis=1)
+