bounding_linear_box.m

一个matlab的将军模型

function box = bounding_linear_box(loc,mapping_k_1,face_poly)

global GLOBAL_PIHA
global GLOBAL_APPROX_PARAM

% Get resolution to be used for quantization
min_angle 				= 	GLOBAL_APPROX_PARAM.min_angle;
med_angle 				= 	GLOBAL_APPROX_PARAM.med_angle;
extra_angle				= 	GLOBAL_APPROX_PARAM.extra_angle;
max_angle 				= 	GLOBAL_APPROX_PARAM.max_angle;
unbound_angle	 		= 	GLOBAL_APPROX_PARAM.unbound_angle;
edge_factor 			= 	GLOBAL_APPROX_PARAM.edge_factor;
med_length				=  GLOBAL_APPROX_PARAM.edge_med_length ;
lim_min_parallel		=	cos(min_angle*pi/180);
lim_med_parallel		=	cos(med_angle*pi/180);
lim_extra_parallel	=	cos(extra_angle*pi/180);
lim_max_parallel		=	cos(max_angle*pi/180);
lim_unbound				=	cos(unbound_angle*pi/180);

% Extract relevant information from the mapping

Vm		= vertices;
[CEm,dEm,CIm,dIm]	= linearcon_data(mapping_k_1);
PH 	= polyhedron(mapping_k_1);
Vm		= get_param(PH,'VI');

% If dimension is less than or equal to 2, return the mapping as the result
if dim(Vm{1})<=2
   Cface 		= get_linearcon_param(face_poly,'CE');
   dface 		= get_linearcon_param(face_poly,'dE');
   CInew			= get_linearcon_param(mapping_k_1,'CI');
   dInew			= get_linearcon_param(mapping_k_1,'dI');
	new_mapping = linearcon(Cface,dface,CInew,dInew);
	box = and(new_mapping,face_poly);
   return
end

% Compute the projection of the CIm vectors  in the direction of CEm.
vec=[];
base=[];
CIproj = CIm;
for i=1:size(CEm,1)
  vec = CEm(i,:); 
  if ~ compvec(vec,zeros(1,length(vec)))
  	base(i,:) = vec/sqrt(vec*vec');
   	for j=1:size(CIm,1)
         CIproj(j,:) = CIproj(j,:) - (CIproj(j,:)*base(i,:)')*base(i,:);
         if abs(CIproj(j,:)*base(i,:)')> GLOBAL_APPROX_PARAM.poly_epsilon
            fprintf(1,'\n warning! normal vector could not be found.\n')
         end
      end
  end
end




%compute the length of each edge. The 'length' is always the biggest distance between
%two vertices composing the hyperplane. In the 3-dimension case, is the length of the edge.
dist=zeros(1,length(Vm));
for i=1:length(Vm)
   for j=1:length(Vm{i})
      for k = j+1:length(Vm{i})
         aux_dist = sqrt((Vm{i}(k) - Vm{i}(j))'*(Vm{i}(k) - Vm{i}(j)));
         if (dist(i)<=aux_dist)
            dist(i)=aux_dist;
         end
      end
   end
end
%Compute the average value of vector dist.
med =mean(dist);

% Compute the neighbors for each edge of the polyhedron
% This block scans the set of vertices for each edge, comparing each vertex with the other ones.
% If two edges have a number of common vertices greater than num_points, they are contiguous
%
neighbors = {};
vector_elem = {};
nvec = {};
num_points = dim(Vm{1}) - 2;
for i=1:length(Vm)
   
   % Finding common vertices between hyperplanes (candidates to neighbors)
   for j=1:length(Vm{i})
      for k = 1:length(Vm)
         if k~= i
            for m = 1: length(Vm{k})
%               if abs(Vm{i}(j)', - Vm{k}(m))< poly_param('epsilon')
                if ismember(Vm{i}(j)',Vm{k}(m)','rows')
						if length(neighbors)<i
                     neighbors{i} = k;
                 	else
                    	neighbors{i} = [neighbors{i} k];
                 	end
            	end
            end
         end
      end
   end
   
   % Looking at the number os vertices to decide which hyperplanes are neighbors
   for j=1:length(neighbors{i})
      str_neighbors = findstr(neighbors{i},neighbors{i}(j));
      if length(str_neighbors) >=num_points
         if length(nvec)<i
           	nvec{i}  = neighbors{i}(j);
        	else
           	nvec{i} = [nvec{i} neighbors{i}(j)];
        	end
		end
   end
   nvec{i} = unique(nvec{i});
   
   % Computing the common vertices between neighbors
   if i==25
      i
   end
   for j=1:length(nvec{i})
      for k=1:length(Vm{i})
         set_of_vertices = Vm{nvec{i}(j)}(1:length(Vm{nvec{i}(j)}));
         Vm{i}(k);
%         flag_true = 0;
%         for m=1:size(set_of_vertices,2)
%            if abs(Vm{i}(k) - set_of_vertices(:,m)) < poly_param('epsilon')
%               flag_true = 1;
%               break
%            end
%         end
         if ismember(Vm{i}(k)',set_of_vertices','rows')
         %if flag_true
            if length(vector_elem)<i
               vector_elem{i}{j} = [Vm{i}(k)];
            elseif length(vector_elem{i})<j
               vector_elem{i}{j} = [Vm{i}(k)];
            else
               aux_var1 = vector_elem{i}{j};
               vector_elem{i}{j} = [aux_var1 Vm{i}(k)];
               aux_var1 =[];
            end
         end
      end
   end
   
end
%save ws1021991 nvec vector_elem Vm PH CIm CEm dEm dIm
%Search CIproj,dIm looking for contigous edges "almost" parallel. There are the following criteria to Eliminate edges:
% If angle(v1,v2)<=min_angle
%		Find neighbors of v2 (N1,...Nn)
%		If (angle(v1,Ni)< unbound_angle with i=1:n) and 
%			(angle(v1,v2)< angle(v2,Ni)  with i=1:n) and
%			||v1|| > ||v2||
%				v2=0,replace v2 by v1 for all occurrences in nvec
%		end
%	elseIf angle(v1,v2)<=med_angle
%		Find neighbors of v2 (N1,...Nn)
%		If (angle(v1,Ni)< max_angle with i=1:n) and 
%			(angle(v1,v2)< angle(v2,Ni)  with i=1:n) and
%			||v1|| > edge_factor*||v2||
%				v2=0,replace v2 by v1 for all occurrences in nvec
%		end
%  end

i = 1;
END_OK = 0;
vec=[];
while ~END_OK
   if i==1
      CIold = CIproj;
   end
   vec0 = CIproj(i,:);
   iprod1 =[];
   elem  = 0;
   % Comparing neighbors in order to find candidates to elimination
  	if ~ compvec(vec0,zeros(1,length(vec0)))
      for j=1:length(nvec{i})
         vec1 = CIproj(nvec{i}(j),:);
         if ~compvec(vec1,zeros(1,length(vec1)))
            aux1 = vec1*vec0';
            aux2 = vec1*vec1';
            aux3 = vec0*vec0';
            iprod1(j) = aux1/sqrt(aux2*aux3);
         else
            iprod1(j)=0;
         end
   	end
      if ~isempty(iprod1)
      	[elem,index] = max(iprod1);
      	if (elem>=lim_min_parallel)
         	test_value 	= 1;
            test_angle 	= 1;
            
            % Testing angle between neighbors of the vertex to be eliminated and
            % the vector that will replace it.
        		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)> dist(nvec{i}(index)));
            if (test_value & test_angle & test_length) 
               j 				= index;
            	j1 			= nvec{i}(j);
               [new_CI,new_d,new_proj]	= ...
                  create_new_hyper(vector_elem{i}{j},CIm(i,:),CIm(j1,:),dist(i),dist(j1),CEm);
               dist(j1) 	= 0;
					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