📄 bounding_linear_box.m
字号:
elseif (elem>=lim_med_parallel)
test_value = 1;
test_angle = 1;
for k=1:length(nvec{nvec{i}(index)})
vec2 = CIproj(nvec{nvec{i}(index)}(k),:);
vec3 = CIproj(nvec{i}(index),:);
if ~compvec(vec2,zeros(1,length(vec2)))
if (elem < (vec2*vec3')/sqrt((vec2*vec2')*(vec3*vec3')) - GLOBAL_APPROX_PARAM.poly_epsilon)
test_value = 0;
end
factor= (vec2*vec0')/sqrt((vec2*vec2')*(vec0*vec0'));
if (lim_max_parallel - factor) > GLOBAL_APPROX_PARAM.poly_epsilon
test_angle = 0;
end
end
end
test_length = (dist(i)> edge_factor*dist(nvec{i}(index)));
if (test_value & test_angle & test_length)
j = index;
j1 = nvec{i}(j);
dist(j1) = 0;
[new_CI,new_d,new_proj] = ...
create_new_hyper(vector_elem{i}{j},CIm(i,:),CIm(j1,:),dist(i),dist(j1),CEm);
CIm(i,:) = new_CI;
dIm(i) = new_d;
CIproj(i,:) = new_proj;
CIproj(j1,:)= zeros(1,size(CIproj,2));
CIm(j1,:) = zeros(1,size(CIm,2));
dIm(j1,:) = 0;
for m=1:length(nvec)
for k=1:length(nvec{m})
if nvec{m}(k) == j1
nvec{m}(k) = i;
end
end
end
nvec{i} = [nvec{i} nvec{j1}];
vector = [];
for n=1:length(nvec{i})
if nvec{i}(n) ~= i;
vector = [vector nvec{i}(n)];
end
end
nvec{i} = unique(vector);
nvec{j1} = [];
[vector_elem{i},dist(i)]= compute_new_vertices(nvec{i},i,CIm,dIm,CEm,dEm);
%save teste1 vector_elem
vector_elem{j1} = {};
med = adjust_med(dist);
end
elseif (elem>=lim_extra_parallel)
test_value = 1;
test_angle = 1;
for k=1:length(nvec{nvec{i}(index)})
vec2 = CIproj(nvec{nvec{i}(index)}(k),:);
vec3 = CIproj(nvec{i}(index),:);
if ~compvec(vec2,zeros(1,length(vec2)))
if (elem < (vec2*vec3')/sqrt((vec2*vec2')*(vec3*vec3')) - GLOBAL_APPROX_PARAM.poly_epsilon)
test_value = 0;
end
factor= (vec2*vec0')/sqrt((vec2*vec2')*(vec0*vec0'));
if (lim_unbound - factor) > GLOBAL_APPROX_PARAM.poly_epsilon
test_angle = 0;
end
end
end
test_length = (dist(i)> edge_factor*dist(nvec{i}(index)));
teste_med_length_j1 = (dist(nvec{i}(index))<= med_length*med);
if (test_value & test_angle & teste_med_length_j1 & test_length)
j = index;
j1 = nvec{i}(j);
dist(j1) = 0;
[new_CI,new_d,new_proj] = ...
create_new_hyper(vector_elem{i}{j},CIm(i,:),CIm(j1,:),dist(i),dist(j1),CEm);
CIm(i,:) = new_CI;
dIm(i) = new_d;
CIproj(i,:) = new_proj;
CIproj(j1,:)= zeros(1,size(CIproj,2));
CIm(j1,:) = zeros(1,size(CIm,2));
dIm(j1,:) = 0;
for m=1:length(nvec)
for k=1:length(nvec{m})
if nvec{m}(k) == j1
nvec{m}(k) = i;
end
end
end
nvec{i} = [nvec{i} nvec{j1}];
vector = [];
for n=1:length(nvec{i})
if nvec{i}(n) ~= i;
vector = [vector nvec{i}(n)];
end
end
nvec{i} = unique(vector);
nvec{j1} = [];
[vector_elem{i},dist(i)]= compute_new_vertices(nvec{i},i,CIm,dIm,CEm,dEm);
vector_elem{j1} = {};
med = adjust_med(dist);
end
end
end
end
i = mod(i, length(nvec)+1) + 1;
if (i==length(nvec)+1)
if (CIproj == CIold)
END_OK = 1;
end
i=1;
end
end
% Throw away the lines of CIm that are zeroed.
CInew=[];
dInew=[];
for i = 1:length(CIproj)
if ~compvec(CIproj(i,:),zeros(1,size(CIproj,2)))
CInew(size(CInew,1)+1,:)= CIm(i,:);
dInew(length(dInew)+1)= dIm(i);
end
end
% Make sure the new mapping is contained in the invariant face.
Cface = get_linearcon_param(face_poly,'CE');
dface = get_linearcon_param(face_poly,'dE');
new_mapping = linearcon(Cface,dface,CInew,dInew');
box = and(new_mapping,face_poly);
%fprintf(1,'printing box...\n')
%plot(box);
%fprintf(1,' done. Press any key\n')
%pause;
return
% -----------------------------------------------------------------------------
function retcode = compvec(v1,v2)
if (length(v1) == length(v2))
if (v1 == v2)
retcode = 1;
else
retcode =0;
end
else
retcode =0;
end
return
% -----------------------------------------------------------------------------
function [new_CI,new_d,new_proj] = create_new_hyper(vector_elem,CI_A,CI_B,distA,distB,CE)
global GLOBAL_APPROX_PARAM
new_CI = (CI_A*distA + CI_B*distB)/(distA+distB);
if ~compvec(new_CI,zeros(1,length(new_CI)))
new_CI = new_CI/sqrt(new_CI*new_CI');
else
new_CI = CI_A;
end
for i=1:size(vector_elem,2)
d(i) = new_CI*vector_elem(:,i);
end
new_d = max(d);
vec=[];
base=[];
new_proj = new_CI;
for i=1:size(CE,1)
vec = CE(i,:);
if ~ compvec(vec,zeros(1,length(vec)))
base(i,:) = vec/sqrt(vec*vec');
new_proj = new_proj - (new_proj*base(i,:)')*base(i,:);
if abs(new_proj*base(i,:)')>GLOBAL_APPROX_PARAM.poly_epsilon
fprintf(1,'\n warning! normal vector could not be found.\n')
end
end
end
return
% -----------------------------------------------------------------------------
function [new_vertices,dist] = compute_new_vertices(neighbors,index,CI,dI,CE,dE)
new_vertices = {};
CEp =[CE;CI(index,:)];
dEp = [dE;dI(index)];
CIp = [];
dIp = [];
for o=1:length(neighbors)
CIp = [CIp ; CI(neighbors(o),:)];
dIp = [dIp ; dI(neighbors(o))];
end
new_poly = linearcon(CEp,dEp,CIp,dIp);
new_polytope = clean_up(new_poly);
new_polyhedron = polyhedron(new_polytope);
av = vertices(new_polytope);
Vm = get_param(new_polyhedron,'VI');
for j=1:length(Vm)
for k=1:length(Vm{j})
aux = Vm{j}(k);
if length(new_vertices)<j
new_vertices{j} = aux;
else
new_vertices{j} = [new_vertices{j} aux];
end
end
end
dist = 0;
for j=1:length(av)
for k=j+1:length(av)
aux_dist = sqrt((av(j) - av(k))'*(av(j) - av(k)));
if (dist <= aux_dist)
dist =aux_dist;
end
end
end
return
% -----------------------------------------------------------------------------
function med = adjust_med(dist)
x=[];
for i=1:length(dist)
if (dist(i)~=0)
x(length(x)+1)=dist(i);
end
end
med = mean(x);
return⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -