Linia Krzywa 2D

Funkcja LineCurvature2D oblicza krzywiznę linii 2D. Najpierw dopasowuje wielokąty do punktów, a następnie oblicza krzywiznę analityczną na podstawie wielokątów.

Zmienne wejściowe:

1. `Vertices`:

Macierz M x 2
Opis: Lista punktów tworzących linię, gdzie każda z M wierszy zawiera współrzędne x i y danego punktu.
Przykład: Vertices = [x1, y1; x2, y2; ...; xM, yM]

2. `Lines` (opcjonalnie):

Macierz N x 2
Opis: Lista odcinków linii zdefiniowana przez indeksy punktów z Vertices. Każdy wiersz zawiera dwa indeksy, które definiują odcinek linii.
Przykład: Lines = [index1, index2; index2, index3; ...; indexN-1, indexN]
Domyślnie: Jeśli Lines nie jest podane, przyjmuje się, że Lines = [1 2; 2 3; ...; M-1 M], co oznacza, że punkty są połączone w kolejności, w jakiej są podane.

Zmienna wyjściowa:

1. `k`:

Wektor M x 1
Opis: Wartości krzywizny dla każdego z M punktów linii. Każdy element wektora k odpowiada krzywiźnie w odpowiednim punkcie Vertices.
Przykład: k = [k1; k2; ...; kM], gdzie ki to krzywizna w punkcie i.

Przykład użycia:

Jeśli mamy punkty linii:

Vertices = [0, 0; 1, 1; 2, 0; 3, -1]; i chcemy obliczyć krzywiznę, możemy wywołać funkcję: k = LineCurvature2D(Vertices);

Jeśli chcemy określić własne odcinki linii:

Lines = [1, 2; 2, 3; 3, 4];
k = LineCurvature2D(Vertices, Lines);

Przykłady:

Koło

1. Generowanie losowych kątów:

r = sort(rand(15,1)) * 2 * pi;

Generujemy 15 losowych liczb w zakresie od 0 do 1 (rand(15,1)).
Sortujemy te liczby rosnąco (sort).
Przekształcamy je na kąty w zakresie od 0 do 2π poprzez mnożenie przez 2π.

2. Tworzenie punktów okręgu:

Vertices = [sin(r) cos(r)] * 10;

Używamy funkcji trygonometrycznych sin i cos, aby przekształcić kąty na współrzędne punktów na okręgu.
Skalujemy współrzędne przez 10, aby uzyskać punkty na okręgu o promieniu 10.
Vertices to macierz 15x2, gdzie każda para [sin(r_i) cos(r_i)] * 10 reprezentuje współrzędne punktu.

3. Tworzenie linii:

Lines = [(1:size(Vertices,1))' (2:size(Vertices,1)+1']; Lines(end,2) = 1;

Tworzymy macierz Lines, która definiuje odcinki linii przez indeksy punktów w Vertices.
Każdy wiersz Lines łączy punkt i z punktem i+1.
Ostatnia para jest specjalnie ustawiona (Lines(end,2) = 1), aby zamknąć okrąg, łącząc ostatni punkt z pierwszym.

4. Obliczenie krzywizny:

k = LineCurvature2D(Vertices, Lines);

Wywołujemy funkcję LineCurvature2D z punktami Vertices i liniami Lines, aby obliczyć krzywiznę.
k to wektor wartości krzywizny dla każdego punktu.

figure, hold on;

Tworzy nowe okno wykresu.
Włącza tryb hold on, co pozwala na nałożenie wielu wykresów na jednej figurze bez ich nadpisywania.

5. Obliczenie normalnych do linii:

N = LineNormals2D(Vertices,Lines);

Wywołuje funkcję LineNormals2D, która oblicza normalne do linii w punktach Vertices na podstawie Lines.
N to macierz zawierająca wektory normalne dla każdego punktu.

6. Skalowanie wartości krzywizny:

k = k * 100;

Skaluje wartości krzywizny k przez mnożenie przez 100, aby lepiej widoczne były na wykresie.
N to macierz zawierająca wektory normalne dla każdego punktu.

7. Wykreślanie normalnych do linii:

plot([Vertices(:,1) Vertices(:,1)+k.*N(:,1)]',[Vertices(:,2) Vertices(:,2)+k.*N(:,2)]','g');

Wykreśla wektory normalne w każdym punkcie Vertices.
Vertices(:,1) to współrzędne x punktów początkowych.
Vertices(:,1) + k.*N(:,1) to współrzędne x końców wektorów normalnych.
Vertices(:,2) to współrzędne y punktów początkowych.
Vertices(:,2) + k.*N(:,2) to współrzędne y końców wektorów normalnych.
Kolor wykresu: zielony ('g').

8. Wykreślanie linii:

plot([Vertices(Lines(:,1),1) Vertices(Lines(:,2),1)]',[Vertices(Lines(:,1),2) Vertices(Lines(:,2),2)]','b');

Wykreśla linie łączące punkty z Vertices zgodnie z Lines.
Vertices(Lines(:,1),1) i Vertices(Lines(:,2),1) to współrzędne x punktów początkowych i końcowych odcinków linii.
Vertices(Lines(:,1),2) i Vertices(Lines(:,2),2) to współrzędne y punktów początkowych i końcowych odcinków linii.
Kolor wykresu: niebieski ('b').

9. Wykreślanie idealnego okręgu:

plot(sin(0:0.01:2*pi)*10, cos(0:0.01:2*pi)*10, 'r.');

Wykreśla idealny okrąg o promieniu 10.
sin(0:0.01:2*pi)*10 i cos(0:0.01:2*pi)*10 to współrzędne x i y punktów na okręgu.
Kolor wykresu: czerwony ('r').

10. Ustawienie równej skali osi:

axis equal;

Ustawia równe jednostki na osiach x i y, co sprawia, że okrąg wygląda proporcjonalnie.

Przykład, Ręka:

1. Wczytanie danych z pliku:

load('testdata');

Wczytuje dane z pliku testdata.mat, które zawierają zmienne Vertices (współrzędne punktów) i Lines (lista odcinków linii).

2. Obliczenie krzywizny linii:

k = LineCurvature2D(Vertices, Lines);

Wywołuje funkcję LineCurvature2D, która oblicza wartości krzywizny dla punktów zdefiniowanych w Vertices i Lines.
k to wektor zawierający wartości krzywizny dla każdego punktu.

3. Utworzenie nowej figury i włączenie trybu hold on:

figure, hold on;

Tworzy nowe okno wykresu.
Włącza tryb hold on, co pozwala na nałożenie wielu wykresów na jednej figurze bez ich nadpisywania.

4. Obliczenie normalnych do linii:

N = LineNormals2D(Vertices, Lines);

Wywołuje funkcję LineNormals2D, która oblicza wektory normalne dla punktów z Vertices na podstawie Lines.
N to macierz zawierająca wektory normalne dla każdego punktu.

5. Skalowanie wartości krzywizny:

k = k * 100;

Skaluje wartości krzywizny k przez mnożenie przez 100, aby lepiej widoczne były na wykresie.

6. Wykreślanie normalnych do linii:

plot([Vertices(:,1) Vertices(:,1) + k .* N(:,1)]', [Vertices(:,2) Vertices(:,2) + k .* N(:,2)]', 'g');

Wykreśla wektory normalne w każdym punkcie Vertices.
Vertices(:,1) to współrzędne x punktów początkowych.
Vertices(:,1) + k .* N(:,1) to współrzędne x końców wektorów normalnych.
Vertices(:,2) to współrzędne y punktów początkowych.
Vertices(:,2) + k .* N(:,2) to współrzędne y końców wektorów normalnych.
Kolor wykresu: zielony ('g').

7. Wykreślanie linii:

plot([Vertices(Lines(:,1),1) Vertices(Lines(:,2),1)]', [Vertices(Lines(:,1),2) Vertices(Lines(:,2),2)]', 'b'); `

Wykreśla linie łączące punkty z Vertices zgodnie z Lines.
Vertices(Lines(:,1),1) i Vertices(Lines(:,2),1) to współrzędne x punktów początkowych i końcowych odcinków linii.
Vertices(Lines(:,1),2) i Vertices(Lines(:,2),2) to współrzędne y punktów początkowych i końcowych odcinków linii.
Kolor wykresu: niebieski ('b').

8. Wykreślanie punktów:

plot(Vertices(:,1), Vertices(:,2), 'r.');

Wykreśla punkty Vertices.
Kolor wykresu: czerwony ('r').

9. Ustawienie równej skali osi:

axis equal;

Ustawia równe jednostki na osiach x i y, co zapewnia proporcjonalne wyświetlanie.

10. Lista zmiennych i ich opis:

Vertices:

Typ: Macierz M x 2
Opis: Współrzędne punktów (wierzchołków) linii.

Lines:

Typ: Macierz N x 2
Opis: Lista odcinków linii, zdefiniowana przez indeksy punktów w Vertices.

k:

Typ: Wektor M x 1
Opis: Wartości krzywizny dla punktów Vertices, skalowane przez 100.

N:

Typ: Macierz M x 2
Opis: Wektory normalne dla punktów z Vertices, obliczone przez LineNormals2D.

Funkcja

Obliczenie krzywizny linii na podstawie zadanych wierzchołków (Vertices) i odcinków linii (Lines). Obejmuje on kroki takie jak uzupełnianie brakujących sąsiadów, przeliczanie wektorów normalnych oraz obliczanie krzywizny na podstawie dopasowania wielomianu do wierzchołków.

Kod:

% Function is written by D.Kroon University of Twente (August 2011)

% If no line-indices, assume a x(1) connected with x(2), x(3) with x(4) ...
if(nargin<2)
    Lines=[(1:(size(Vertices,1)-1))' (2:size(Vertices,1))'];
end

% Get left and right neighbor of each points
Na=zeros(size(Vertices,1),1); Nb=zeros(size(Vertices,1),1);
Na(Lines(:,1))=Lines(:,2); Nb(Lines(:,2))=Lines(:,1);

% Check for end of line points, without a left or right neighbor
checkNa=Na==0; checkNb=Nb==0;
Naa=Na; Nbb=Nb;
Naa(checkNa)=find(checkNa); Nbb(checkNb)=find(checkNb);

% If no left neighbor use two right neighbors, and the same for right... 
Na(checkNa)=Nbb(Nbb(checkNa)); Nb(checkNb)=Naa(Naa(checkNb));

% Correct for sampeling differences
Ta=-sqrt(sum((Vertices-Vertices(Na,:)).^2,2));
Tb=sqrt(sum((Vertices-Vertices(Nb,:)).^2,2)); 

% If no left neighbor use two right neighbors, and the same for right... 
Ta(checkNa)=-Ta(checkNa); Tb(checkNb)=-Tb(checkNb);

% Fit a polygons to the vertices 
% x=a(3)*t^2 + a(2)*t + a(1) 
% y=b(3)*t^2 + b(2)*t + b(1) 
% we know the x,y of every vertice and set t=0 for the vertices, and
% t=Ta for left vertices, and t=Tb for right vertices,  
x = [Vertices(Na,1) Vertices(:,1) Vertices(Nb,1)];
y = [Vertices(Na,2) Vertices(:,2) Vertices(Nb,2)];
M = [ones(size(Tb)) -Ta Ta.^2 ones(size(Tb)) zeros(size(Tb)) zeros(size(Tb)) ones(size(Tb)) -Tb Tb.^2];
invM=inverse3(M);
a(:,1)=invM(:,1,1).*x(:,1)+invM(:,2,1).*x(:,2)+invM(:,3,1).*x(:,3);
a(:,2)=invM(:,1,2).*x(:,1)+invM(:,2,2).*x(:,2)+invM(:,3,2).*x(:,3);
a(:,3)=invM(:,1,3).*x(:,1)+invM(:,2,3).*x(:,2)+invM(:,3,3).*x(:,3);
b(:,1)=invM(:,1,1).*y(:,1)+invM(:,2,1).*y(:,2)+invM(:,3,1).*y(:,3);
b(:,2)=invM(:,1,2).*y(:,1)+invM(:,2,2).*y(:,2)+invM(:,3,2).*y(:,3);
b(:,3)=invM(:,1,3).*y(:,1)+invM(:,2,3).*y(:,2)+invM(:,3,3).*y(:,3);

% Calculate the curvature from the fitted polygon
k = 2*(a(:,2).*b(:,3)-a(:,3).*b(:,2)) ./ ((a(:,2).^2+b(:,2).^2).^(3/2));

function  Minv  = inverse3(M)
% This function does inv(M) , but then for an array of 3x3 matrices
adjM(:,1,1)=  M(:,5).*M(:,9)-M(:,8).*M(:,6);
adjM(:,1,2)=  -(M(:,4).*M(:,9)-M(:,7).*M(:,6));
adjM(:,1,3)=  M(:,4).*M(:,8)-M(:,7).*M(:,5);
adjM(:,2,1)=  -(M(:,2).*M(:,9)-M(:,8).*M(:,3));
adjM(:,2,2)=  M(:,1).*M(:,9)-M(:,7).*M(:,3);
adjM(:,2,3)=  -(M(:,1).*M(:,8)-M(:,7).*M(:,2));
adjM(:,3,1)=  M(:,2).*M(:,6)-M(:,5).*M(:,3);
adjM(:,3,2)=  -(M(:,1).*M(:,6)-M(:,4).*M(:,3));
adjM(:,3,3)=  M(:,1).*M(:,5)-M(:,4).*M(:,2);
detM=M(:,1).*M(:,5).*M(:,9)-M(:,1).*M(:,8).*M(:,6)-M(:,4).*M(:,2).*M(:,9)+M(:,4).*M(:,8).*M(:,3)+M(:,7).*M(:,2).*M(:,6)-M(:,7).*M(:,5).*M(:,3);
Minv=bsxfun(@rdivide,adjM,detM);

Kroki działania:

1. Sprawdzenie i przypisanie domyślnych indeksów linii:

if(nargin<2) Lines=[(1:(size(Vertices,1)-1))' (2:size(Vertices,1))']; end

Jeśli nie podano Lines, domyślnie zakłada, że punkty są połączone sekwencyjnie.

2. Znalezienie lewego i prawego sąsiada dla każdego punktu:

Na=zeros(size(Vertices,1),1); Nb=zeros(size(Vertices,1),1); Na(Lines(:,1))=Lines(:,2); Nb(Lines(:,2))=Lines(:,1);

Na i Nb przechowują indeksy lewych i prawych sąsiadów dla każdego punktu.

3. Sprawdzenie punktów bez sąsiadów i uzupełnienie braków:

checkNa=Na==0; checkNb=Nb==0;
Naa=Na; Nbb=Nb;
Naa(checkNa)=find(checkNa); Nbb(checkNb)=find(checkNb);

Na(checkNa)=Nbb(Nbb(checkNa)); Nb(checkNb)=Naa(Naa(checkNb));

Jeśli punkt nie ma lewego lub prawego sąsiada, używa dwóch sąsiadów z przeciwnej strony.

4. Korekta różnic próbkowania:

Ta=-sqrt(sum((Vertices-Vertices(Na,:)).^2,2));
Tb=sqrt(sum((Vertices-Vertices(Nb,:)).^2,2)); 

Ta(checkNa)=-Ta(checkNa); Tb(checkNb)=-Tb(checkNb);

Oblicza odległości do sąsiadów, uwzględniając brakujące wartości.

5. Dopasowanie wielomianu do wierzchołków:

x = [Vertices(Na,1) Vertices(:,1) Vertices(Nb,1)]; y = [Vertices(Na,2) Vertices(:,2) Vertices(Nb,2)]; M = [ones(size(Tb)) -Ta Ta.^2 ones(size(Tb)) zeros(size(Tb)) zeros(size(Tb)) ones(size(Tb)) -Tb Tb.^2]; invM=inverse3(M); a(:,1)=invM(:,1,1).*x(:,1)+invM(:,2,1).*x(:,2)+invM(:,3,1).*x(:,3); a(:,2)=invM(:,1,2).*x(:,1)+invM(:,2,2).*x(:,2)+invM(:,3,2).*x(:,3); a(:,3)=invM(:,1,3).*x(:,1)+invM(:,2,3).*x(:,2)+invM(:,3,3).*x(:,3); b(:,1)=invM(:,1,1).*y(:,1)+invM(:,2,1).*y(:,2)+invM(:,3,1).*y(:,3); b(:,2)=invM(:,1,2).*y(:,1)+invM(:,2,2).*y(:,2)+invM(:,3,2).*y(:,3); b(:,3)=invM(:,1,3).*y(:,1)+invM(:,2,3).*y(:,2)+invM(:,3,3).*y(:,3);

Dopasowuje wielomiany do współrzędnych punktów, biorąc pod uwagę sąsiadów.

6. Obliczenie krzywizny:

k = 2*(a(:,2).*b(:,3)-a(:,3).*b(:,2)) ./ ((a(:,2).^2+b(:,2).^2).^(3/2));

7. Odwracanie macierzy 3x3:

function  Minv  = inverse3(M)
adjM(:,1,1)=  M(:,5).*M(:,9)-M(:,8).*M(:,6);
adjM(:,1,2)=  -(M(:,4).*M(:,9)-M(:,7).*M(:,6));
adjM(:,1,3)=  M(:,4).*M(:,8)-M(:,7).*M(:,5);
adjM(:,2,1)=  -(M(:,2).*M(:,9)-M(:,8).*M(:,3));
adjM(:,2,2)=  M(:,1).*M(:,9)-M(:,7).*M(:,3);
adjM(:,2,3)=  -(M(:,1).*M(:,8)-M(:,7).*M(:,2));
adjM(:,3,1)=  M(:,2).*M(:,6)-M(:,5).*M(:,3);
adjM(:,3,2)=  -(M(:,1).*M(:,6)-M(:,4).*M(:,3));
adjM(:,3,3)=  M(:,1).*M(:,5)-M(:,4).*M(:,2);
detM=M(:,1).*M(:,5).*M(:,9)-M(:,1).*M(:,8).*M(:,6)-M(:,4).*M(:,2).*M(:,9)+M(:,4).*M(:,8).*M(:,3)+M(:,7).*M(:,2).*M(:,6)-M(:,7).*M(:,5).*M(:,3);
Minv=bsxfun(@rdivide,adjM,detM);

Lista zmiennych i ich opis: