Description of select3d

0001 function [pout, vout, viout, facevout, faceiout]  = select3d(obj)
0002 %SELECT3D(H) Determines the selected point in 3-D data space.
0003 %  P = SELECT3D determines the point, P, in data space corresponding
0004 %  to the current selection position. P is a point on the first
0005 %  patch or surface face intersected along the selection ray. If no
0006 %  face is encountered along the selection ray, P returns empty.
0007 %
0008 %  P = SELECT3D(H) constrains selection to graphics handle H and,
0009 %  if applicable, any of its children. H can be a figure, axes,
0010 %  patch, or surface object.
0011 %
0012 %  [P V] = SELECT3D(...), V is the closest face or line vertex
0013 %  selected based on the figure's current object.
0014 %
0015 %  [P V VI] = SELECT3D(...), VI is the index into the object's
0016 %  x,y,zdata properties corresponding to V, the closest face vertex
0017 %  selected.
0018 %
0019 %  [P V VI FACEV] = SELECT3D(...), FACE is an array of vertices
0020 %  corresponding to the face polygon containing P and V.
0021 %
0022 %  [P V VI FACEV FACEI] = SELECT3D(...), FACEI is the row index into
0023 %  the object's face array corresponding to FACE. For patch
0024 %  objects, the face array can be obtained by doing
0025 %  get(mypatch,'faces'). For surface objects, the face array
0026 %  can be obtained from the output of SURF2PATCH (see
0027 %  SURF2PATCH for more information).
0028 %
0029 %  RESTRICTIONS:
0030 %  SELECT3D supports surface, patch, or line object primitives. For surface
0031 %  and patches, the algorithm assumes non-self-intersecting planar faces.
0032 %  For line objects, the algorithm always returns P as empty, and V will
0033 %  be the closest vertex relative to the selection point.
0034 %
0035 %  Example:
0036 %
0037 %  h = surf(peaks);
0038 %  zoom(10);
0039 %  disp('Click anywhere on the surface, then hit return')
0040 %  pause
0041 %  [p v vi face facei] = select3d;
0042 %  marker1 = line('xdata',p(1),'ydata',p(2),'zdata',p(3),'marker','o',...
0043 %                 'erasemode','xor','markerfacecolor','k');
0044 %  marker2 = line('xdata',v(1),'ydata',v(2),'zdata',v(3),'marker','o',...
0045 %                 'erasemode','xor','markerfacecolor','k');
0046 %  marker2 = line('erasemode','xor','xdata',face(1,:),'ydata',face(2,:),...
0047 %                 'zdata',face(3,:),'linewidth',10);
0048 %  disp(sprintf('\nYou clicked at\nX: %.2f\nY: %.2f\nZ: %.2f',p(1),p(2),p(3)'))
0049 %  disp(sprintf('\nThe nearest vertex is\nX: %.2f\nY: %.2f\nZ: %.2f',v(1),v(2),v(3)'))
0050 %
0051 %  Version 1.2 2-15-02
0052 %  Copyright Joe Conti 2002
0053 %  Send comments to jconti@mathworks.com
0054 %
0055 %  See also GINPUT, GCO.
0056 
0057 % Output variables
0058 pout = [];
0059 vout = [];
0060 viout = [];
0061 facevout = [];
0062 faceiout = [];
0063 
0064 % other variables
0065 ERRMSG = 'Input argument must be a valid graphics handle';
0066 isline = logical(0);
0067 isperspective = logical(0);
0068 
0069 % Parse input arguments
0070 if nargin<1
0071    obj = gco;
0072 end
0073 
0074 if isempty(obj) | ~ishandle(obj) | length(obj)~=1
0075     error(ERRMSG);
0076 end
0077 
0078 % if obj is a figure
0079 if strcmp(get(obj,'type'),'figure')
0080     fig = obj;
0081     ax = get(fig,'currentobject');
0082     currobj = get(fig,'currentobject');
0083     
0084     % bail out if not a child of the axes
0085     if ~strcmp(get(get(currobj,'parent'),'type'),'axes')
0086         return;
0087     end
0088     
0089 % if obj is an axes
0090 elseif strcmp(get(obj,'type'),'axes')
0091     ax = obj;
0092     fig = get(ax,'parent');
0093     currobj = get(fig,'currentobject');
0094     currax = get(currobj,'parent');
0095     
0096     % Bail out if current object is under an unspecified axes
0097     if ~isequal(ax,currax)
0098         return;
0099     end
0100 
0101 % if obj is child of axes
0102 elseif strcmp(get(get(obj,'parent'),'type'),'axes')
0103     currobj = obj;
0104     ax = get(obj,'parent');
0105     fig = get(ax,'parent');
0106     
0107 % Bail out
0108 else 
0109     return
0110 end
0111 
0112 axchild = currobj;
0113 obj_type = get(axchild,'type');
0114 is_perspective = strcmp(get(ax,'projection'),'perspective');
0115 
0116 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0117 %% Get projection transformation %%
0118 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0119 
0120 % syntax not supported in old versions of MATLAB
0121 [a b] = view(ax); 
0122 xform = viewmtx(a,b);
0123 if is_perspective
0124     warning(sprintf('%s does not support perspective axes projection.',mfilename));
0125     d = norm(camtarget(ax)-campos(ax))
0126     P = [1 0 0 0; 
0127          0 1 0 0;
0128          0 0 1 0;
0129          0 0 -1/d 1];
0130      xform = P*xform;
0131 end
0132 
0133 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0134 %% Get vertex, face, and current point data %%
0135 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0136 cp = get(ax,'currentpoint')';
0137 
0138 % If surface object
0139 if strcmp(obj_type,'surface')
0140     % Get surface face and vertices
0141     fv = surf2patch(axchild);
0142     vert = fv.vertices;
0143     faces = fv.faces;
0144     
0145 % If patch object
0146 elseif strcmp(obj_type,'patch')
0147     vert = get(axchild,'vertices');
0148     faces = get(axchild,'faces');
0149     
0150 % If line object
0151 elseif strcmp(obj_type,'line')
0152     xdata = get(axchild,'xdata');
0153     ydata = get(axchild,'ydata');
0154     zdata = get(axchild,'zdata');
0155     vert = [xdata', ydata',zdata'];
0156     faces = []; 
0157     isline = logical(1);
0158 
0159 % Ignore all other objects
0160 else     
0161    return;
0162 end
0163     
0164 % Add z if empty
0165 if size(vert,2)==2
0166    vert(:,3) = zeros(size(vert(:,2)));
0167    if isline
0168        zdata = vert(:,3);
0169    end
0170 end
0171 
0172 % NaN and Inf check
0173 nan_inf_test1 = isnan(faces) | isinf(faces);
0174 nan_inf_test2 = isnan(vert) | isinf(vert);
0175 if any(nan_inf_test1(:)) | any(nan_inf_test2(:))
0176     warning(sprintf('%s does not support NaNs or Infs in face/vertex data.',mfilename));
0177 end
0178 
0179 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0180 %% Normalize for data aspect ratio %%
0181 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0182 dar = get(ax,'DataAspectRatio');
0183 
0184 ncp(1,:) = cp(1,:)./dar(1);
0185 ncp(2,:) = cp(2,:)./dar(2);
0186 ncp(3,:) = cp(3,:)./dar(3);
0187 ncp(4,:) = ones(size(ncp(3,:)));  
0188 
0189 nvert(:,1) = vert(:,1)./dar(1);
0190 nvert(:,2) = vert(:,2)./dar(2);
0191 nvert(:,3) = vert(:,3)./dar(3);
0192 nvert(:,4) = ones(size(nvert(:,3)));  
0193 
0194 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0195 %% Transform data to view space %%
0196 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0197 xvert = xform*nvert';  
0198 xcp = xform*ncp; 
0199 
0200 if is_perspective % normalize 4th dimension
0201     xcp(1,:) = xcp(1,:)./xcp(4,:);
0202     xcp(2,:) = xcp(2,:)./xcp(4,:);
0203     xcp(3,:) = xcp(3,:)./xcp(4,:);
0204     xcp(4,:) = xcp(4,:)./xcp(4,:);
0205 
0206     xvert(1,:) = xvert(1,:)./xvert(4,:);
0207     xvert(2,:) = xvert(2,:)./xvert(4,:);
0208     xvert(3,:) = xvert(3,:)./xvert(4,:);
0209     xvert(4,:) = xvert(4,:)./xvert(4,:);
0210 end
0211 
0212 % Ignore 3rd & 4th dimensions for crossing test
0213 xvert(4,:) = [];  
0214 xvert(3,:) = [];
0215 xcp(4,:) = []; 
0216 xcp(3,:) = [];
0217 
0218 % For debugging
0219 % if 0
0220 %     ax1 = getappdata(ax,'testselect3d');
0221 %     if isempty(ax1) | ~ishandle(ax1)
0222 %         fig = figure;
0223 %         ax1 = axes;
0224 %         axis(ax1,'equal');
0225 %         setappdata(ax,'testselect3d',ax1);
0226 %     end
0227 %     cla(ax1);
0228 %     patch('parent',ax1,'faces',faces,'vertices',xvert','facecolor','none','edgecolor','k');
0229 %     line('parent',ax1,'xdata',xcp(1,2),'ydata',xcp(2,2),'zdata',0,'marker','o','markerfacecolor','r','erasemode','xor');
0230 % end
0231 
0232 % Translate vertices so that the selection point is at the origin.
0233 xvert(1,:) = xvert(1,:) - xcp(1,2);
0234 xvert(2,:) = xvert(2,:) - xcp(2,2);
0235 
0236 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0237 %% simple algorithm (almost naive algorithm!) for line objects %%
0238 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0239 if isline
0240 
0241     % Ignoring line width and marker attributes, find closest
0242     % vertex in 2-D view space.
0243     d = xvert(1,:).*xvert(1,:) + xvert(2,:).*xvert(2,:);
0244     [val i] = min(d);
0245     i = i(1); % enforce only one output
0246     
0247     % Assign output
0248     vout = [ xdata(i) ydata(i) zdata(i)];
0249     viout = i;
0250     
0251     return % Bail out early
0252 end
0253    
0254 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0255 %% Perform 2-D crossing test (Jordan Curve Theorem) %%
0256 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0257 
0258 % Find all vertices that have y components less than zero
0259 vert_with_negative_y = zeros(size(faces));
0260 face_y_vert = xvert(2,faces);
0261 ind_vert_with_negative_y = find(face_y_vert<0); 
0262 vert_with_negative_y(ind_vert_with_negative_y) = logical(1);
0263 
0264 % Find all the line segments that span the x axis
0265 is_line_segment_spanning_x = abs(diff([vert_with_negative_y, vert_with_negative_y(:,1)],1,2));
0266 
0267 % Find all the faces that have line segments that span the x axis
0268 ind_is_face_spanning_x = find(any(is_line_segment_spanning_x,2));
0269 
0270 % Ignore data that doesn't span the x axis
0271 candidate_faces = faces(ind_is_face_spanning_x,:);
0272 vert_with_negative_y = vert_with_negative_y(ind_is_face_spanning_x,:);
0273 is_line_segment_spanning_x = is_line_segment_spanning_x(ind_is_face_spanning_x,:);
0274 
0275 % Create line segment arrays
0276 pt1 = candidate_faces;
0277 pt2 = [candidate_faces(:,2:end), candidate_faces(:,1)];
0278 
0279 % Point 1
0280 x1 = reshape(xvert(1,pt1),size(pt1));
0281 y1 = reshape(xvert(2,pt1),size(pt1));
0282 
0283 % Point 2
0284 x2 = reshape(xvert(1,pt2),size(pt2));
0285 y2 = reshape(xvert(2,pt2),size(pt2));
0286 
0287 % Cross product of vector to origin with line segment
0288 cross_product_test = -x1.*(y2-y1) > -y1.*(x2-x1);
0289 
0290 % Find all line segments that cross the positive x axis
0291 crossing_test = (cross_product_test==vert_with_negative_y) & is_line_segment_spanning_x;
0292 
0293 % If the number of line segments is odd, then we intersected the polygon
0294 s = sum(crossing_test,2);
0295 s = mod(s,2);
0296 ind_intersection_test = find(s~=0);
0297 
0298 % Bail out early if no faces were hit
0299 if isempty(ind_intersection_test)
0300    return;    
0301 end
0302 
0303 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0304 %% Plane/ray intersection test %%
0305 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0306 % Perform plane/ray intersection with the faces that passed
0307 % the polygon intersection tests. Grab the only the first
0308 % three vertices since that is all we need to define a plane).
0309 % assuming planar polygons.
0310 candidate_faces = candidate_faces(ind_intersection_test,1:3);
0311 candidate_faces = reshape(candidate_faces',1,prod(size(candidate_faces)));
0312 vert = vert';
0313 candidate_facev = vert(:,candidate_faces);
0314 candidate_facev = reshape(candidate_facev,3,3,length(ind_intersection_test));
0315 
0316 % Get three contiguous vertices along polygon
0317 v1 = squeeze(candidate_facev(:,1,:));
0318 v2 = squeeze(candidate_facev(:,2,:));
0319 v3 = squeeze(candidate_facev(:,3,:));
0320 
0321 % Get normal to face plane
0322 vec1 = [v2-v1];
0323 vec2 = [v3-v2];
0324 crs = cross(vec1,vec2);
0325 mag = sqrt(sum(crs.*crs));
0326 nplane(1,:) = crs(1,:)./mag;
0327 nplane(2,:) = crs(2,:)./mag;
0328 nplane(3,:) = crs(3,:)./mag;
0329 
0330 % Compute intersection between plane and ray
0331 cp1 = cp(:,1);  
0332 cp2 = cp(:,2);  
0333 d = cp2-cp1;
0334 dp = dot(-nplane,v1);
0335 
0336 %A = dot(nplane,d);
0337 A(1,:) = nplane(1,:).*d(1);
0338 A(2,:) = nplane(2,:).*d(2);
0339 A(3,:) = nplane(3,:).*d(3);
0340 A = sum(A,1);
0341 
0342 %B = dot(nplane,pt1)
0343 B(1,:) = nplane(1,:).*cp1(1);
0344 B(2,:) = nplane(2,:).*cp1(2);
0345 B(3,:) = nplane(3,:).*cp1(3);
0346 B = sum(B,1);
0347 
0348 % Distance to intersection point
0349 t = (-dp-B)./A;
0350 
0351 % Find "best" distance (smallest)
0352 [tbest ind_best] = min(t);
0353 
0354 % Determine intersection point
0355 pout = cp1 + tbest .* d;
0356 
0357 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0358 %% Assign additional output variables %%
0359 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0360 if nargout>1
0361     
0362     % Get face index and vertices
0363     faceiout = ind_is_face_spanning_x(ind_intersection_test(ind_best));
0364     facevout = vert(:,faces(faceiout,:));
0365     
0366     % Determine index of closest face vertex intersected
0367     facexv = xvert(:,faces(faceiout,:));
0368     dist = sqrt(facexv(1,:).*facexv(1,:) +  facexv(2,:).*facexv(2,:));
0369     min_dist = min(dist);
0370     min_index = find(dist==min_dist);
0371     
0372     % Get closest vertex index and vertex
0373     viout = faces(faceiout,min_index);
0374     vout = vert(:,viout);
0375 end
select3d

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
