Интерполяция через точки в вашем файле данных не может быть легко выполнена. Я бы рекомендовал design B-сплайновые кривые для интерполяции точек между плоскостью 1 и 2 и плоскостью 6 и 7. Пространства между плоскостью 2 и плоскостью 6, плоскостью 7 и плоскостью 10 кажутся довольно линейными:
clear
filename = 'datafile.txt';
A = importdata(filename);
vertices = A.data(2:end);
vertices = reshape(vertices, 12, [])';
vx = vertices(:, 1:3:10);
vy = vertices(:,2:3:11);
vz = vertices(:,3:3:12);
figure; patch(vx',vy',vz',1); axis equal;
Нет простого способа выполнить такую интерполяцию, потому что вы хотите обеспечить как минимум непрерывность С1 вдоль кривой, чтобы избежать какого-либо острого края. Кривая B-сплайна здесь может быть полезна, но если вы не знакомы с ней, вам будет сложно ее запрограммировать. К счастью, я работал над проектом, требующим подгонки поверхности и кривой, и у меня есть код:
function [x,y,z] = bspline(u, ctrlp, k, knots)
U = bspbasis(u, numel(ctrlp(:,1)), k, knots);
x=U*ctrlp(:,1);
y=U*ctrlp(:,2);
z=U*ctrlp(:,3);
end
function U = bspbasis(u, nctrlp, K, knots)
nu = numel(u);
umax = max(u);
index = 1:nctrlp;
% preallocating variables
U = zeros(nu,nctrlp);
N = zeros(nctrlp+1,K);
% Calculate the denominators for basis functions (k>2). -may be useful when
% the size of point data is substantial, so the calculation is not repeated.
d1 = zeros(nctrlp,K);
d2 = d1;
for m=2:K
d1(:,m) = knots(index+m-1) - knots(index);
d2(:,m) = knots(index+m) - knots(index+1);
end
knots = knots(:);
knots1 = knots(1:nctrlp+1);
knots2 = knots(2:nctrlp+2);
knotSize = size(knots1);
knotc1 = knots(nctrlp+1);
knotc2 = knots(nctrlp+2);
for ui = 1:nu
% k = 1
u_ = u(ui)*ones(knotSize);
NA = u_>=knots1 & u_<knots2;
if u(ui) == umax && knotc1 == umax && knotc2 == umax
NA(1:nctrlp+1) = 0;
NA(end-1,1) = 1;
end
N(:,1) = NA;
% k > 2
for k = 2:K
p1 = (u(ui) - knots(index)) ./ d1(:,k) .* N(index,k-1);
p1(isnan(p1)) = 0;
p2 = (knots(index+k) - u(ui)) ./ d2(:,k) .* N(index+1,k-1);
p2(isnan(p2)) = 0;
N(index,k) = p1 + p2;
end
U(ui,:)=N(1:end-1,k);
end
end
Если вы хотите понять приведенный выше код, вы можете прочитать B-spline Wikipedia . На Youtube также есть множество учебных пособий и интерактивных инструментов, таких как этот .
Код ниже соответствует B-сплайну 3-го порядка для каждого из четырех наборов угловых вершин.
u = linspace(0,1,500);
k = 3;
for i = 1:4
ishift = (i-1) * 3 + 1;
p = vertices(:,ishift:ishift+2);
ctrlp = [p(1,:); [0 0 0]; p(2:3,:); p(5:6,:); zeros(2,3); p(7,:); p(8,:)];
ctrlp(2,:) = 5*ctrlp(3,:) - 4*ctrlp(4,:);
ctrlp(7,:) = 2*ctrlp(6,:) - ctrlp(5,:);
ctrlp(8,:) = 4*ctrlp(9,:) - 3*ctrlp(10,:);
knots = getknots(ctrlp, k);
[x_,y_,z_] = bspline(u, ctrlp, k, knots);
x(:,i) = [p(:,1); x_];
y(:,i) = [p(:,2); y_];
z(:,i) = [p(:,3); z_];
[~,I] = sort(x(:,i));
x(:,i) = x(I,i);
y(:,i) = y(I,i);
z(:,i) = z(I,i);
end
c = repmat(1:numel(x)/4,4,1)';
xx=[x;flipud(x(:,[2,3,4,1]))];
yy=[y;flipud(y(:,[2,3,4,1]))];
zz=[z;flipud(z(:,[2,3,4,1]))];
cc=[c;flipud(c)];
figure; patch(xx,yy,zz,cc);
В результате получается гладкая поверхность:
Я хотел объяснить вам, как работает код, но после того, как я боролся с ним в течение часа, я сдался ... Вместо этого я предоставлю краткое изложение некоторых ключевых моментов ниже.
В сводке используются следующие обозначения: C1, C2, ..., C10 - контрольные точки кривой B-сплайна, а V1, V2, ..., V10 - вершины, используемые для вычисления B- сплайн-кривая. На рисунке ниже показаны контрольные точки и вершины, использованные для расчета первой кривой B-сплайнов.
- Количество значений в
u определяет количество точек на конечной кривой.
- B-сплайновая кривая 3-го порядка достаточна.
- Кривая должна проходить через V1, V2, ..., V10.
- Чтобы пройти через V1 и V10, кривая B-сплайна должна быть фиксированного типа. Следовательно,
C1 = V1 и C10 = V10.
- Чтобы пройти через V2, C2 должен быть на линии, проходящей через V2, а V3 и C3 должны равняться V2. Следовательно,
C2 - V2 = d*(V2-V3) и C3 = V2.
- Чтобы пройти через V3, V4 и V5, должно быть выполнено следующее условие:
C3 = V2; C4 = V3; C5 = V5; C6 = V6.
- Чтобы пройти через V6, C7 должен быть на линии, проходящей через V5 и V6, вместе с C6, равным V6. Следовательно,
C7-V6 = d*(V6-V5).
- Чтобы пройти через V7, C8 должен быть на линии, проходящей через V7 и V8, а C9 должен равняться V7. Следовательно,
C8-V7 = d*(V7-V8); C9 = V7;.
- Поскольку
C10 = V10 и C9 = V7, кривая проходит через V8 и V9.
- Наконец, узлы могут быть либо однородными, либо оцениваться по длине хорды между контрольными точками.
Обновление: getknots
Функция getknots:
function knots = getknots(ctrlp, k)
d = sqrt(sum(diff(ctrlp).^2, 2));
ds = cumsum(d)./sum(d);
knots = [zeros(k,1); ds(k-1:end); ones(k,1)];
end
Все в одном коде для метода bspline:
clear
filename = 'datafile.txt';
A = importdata(filename);
vertices = A.data(2:end);
vertices = reshape(vertices, 12, [])';
u = linspace(0,1,500);
k = 3;
for i = 1:4
ishift = (i-1) * 3 + 1;
p = vertices(:,ishift:ishift+2);
ctrlp = [p(1,:); [0 0 0]; p(2:3,:); p(5:6,:); zeros(2,3); p(7,:); p(8,:)];
ctrlp(2,:) = 5*ctrlp(3,:) - 4*ctrlp(4,:);
ctrlp(7,:) = 2*ctrlp(6,:) - ctrlp(5,:);
ctrlp(8,:) = 4*ctrlp(9,:) - 3*ctrlp(10,:);
knots = getknots(ctrlp, k);
[x_,y_,z_] = bspline(u, ctrlp, k, knots);
x(:,i) = [p(:,1); x_];
y(:,i) = [p(:,2); y_];
z(:,i) = [p(:,3); z_];
[~,I] = sort(x(:,i));
x(:,i) = x(I,i);
y(:,i) = y(I,i);
z(:,i) = z(I,i);
end
c = repmat(1:numel(x)/4,4,1)';
xx=[x;flipud(x(:,[2,3,4,1]))];
yy=[y;flipud(y(:,[2,3,4,1]))];
zz=[z;flipud(z(:,[2,3,4,1]))];
cc=[c;flipud(c)];
figure; patch(xx,yy,zz,cc);
function [x,y,z] = bspline(u, ctrlp, k, knots)
U = bspbasis(u, numel(ctrlp(:,1)), k, knots);
x=U*ctrlp(:,1);
y=U*ctrlp(:,2);
z=U*ctrlp(:,3);
end
function U = bspbasis(u, nctrlp, K, knots)
nu = numel(u);
umax = max(u);
index = 1:nctrlp;
% preallocating variables
U = zeros(nu,nctrlp);
N = zeros(nctrlp+1,K);
% Calculate the denominators for basis functions (k>2). -may be useful when
% the size of point data is substantial, so the calculation is not repeated.
d1 = zeros(nctrlp,K);
d2 = d1;
for m=2:K
d1(:,m) = knots(index+m-1) - knots(index);
d2(:,m) = knots(index+m) - knots(index+1);
end
knots = knots(:);
knots1 = knots(1:nctrlp+1);
knots2 = knots(2:nctrlp+2);
knotSize = size(knots1);
knotc1 = knots(nctrlp+1);
knotc2 = knots(nctrlp+2);
for ui = 1:nu
% k = 1
u_ = u(ui)*ones(knotSize);
NA = u_>=knots1 & u_<knots2;
if u(ui) == umax && knotc1 == umax && knotc2 == umax
NA(1:nctrlp+1) = 0;
NA(end-1,1) = 1;
end
N(:,1) = NA;
% k > 2
for k = 2:K
p1 = (u(ui) - knots(index)) ./ d1(:,k) .* N(index,k-1);
p1(isnan(p1)) = 0;
p2 = (knots(index+k) - u(ui)) ./ d2(:,k) .* N(index+1,k-1);
p2(isnan(p2)) = 0;
N(index,k) = p1 + p2;
end
U(ui,:)=N(1:end-1,k);
end
end
function knots = getknots(ctrlp, k)
d = sqrt(sum(diff(ctrlp).^2, 2));
ds = cumsum(d)./sum(d);
knots = [zeros(k,1); ds(k-1:end); ones(k,1)];
end