vb_distance_min3d

 Average minimum distance from set of points {Vnew} to set of points {Vold}

  [ddsum, ddmin] = vb_distance_min3d(Vold, Vnew, Rmax, Zstep)
 --- Input
 Vnew      : New vertex point
 Vold      : Old vertex point
 Rmax      : Max radius for point search
 Zstep     : Step size to divide Z-axis
 --- Output
 ddmin(n) = Min distance from Vnew(n) to point set {Vold}
 ddsum    = sum_n ddmin(n)/Nnew

  M. Sato  2006-2-3

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function    [ddsum, ddmin] = vb_distance_min3d(Vold, Vnew, Rmax, Zstep)
0002 % Average minimum distance from set of points {Vnew} to set of points {Vold}
0003 %
0004 %  [ddsum, ddmin] = vb_distance_min3d(Vold, Vnew, Rmax, Zstep)
0005 % --- Input
0006 % Vnew      : New vertex point
0007 % Vold      : Old vertex point
0008 % Rmax      : Max radius for point search
0009 % Zstep     : Step size to divide Z-axis
0010 % --- Output
0011 % ddmin(n) = Min distance from Vnew(n) to point set {Vold}
0012 % ddsum    = sum_n ddmin(n)/Nnew
0013 %
0014 %  M. Sato  2006-2-3
0015 %
0016 % Copyright (C) 2011, ATR All Rights Reserved.
0017 % License : New BSD License(see VBMEG_LICENSE.txt)
0018 
0019 if ~exist('Zstep','var') | isempty(Zstep), Zstep = 0.002; end; % [m]
0020 
0021 % New model point number
0022 Npoint  = size(Vnew,1);
0023 ddmin    = zeros(Npoint,1);    % min distance
0024 inew    = zeros(Npoint,1);
0025 
0026 Nold    = size(Vold,1);
0027 dd        = zeros(Nold,1);
0028 iold    = zeros(Nold,1);
0029 Vold2    = zeros(Nold,3);
0030 
0031 % Slice definition
0032 Zmax  = max(Vnew(:,3));
0033 Zmin  = min(Vnew(:,3));
0034 Nstep = ceil((Zmax-Zmin)/Zstep);
0035 
0036 if Nstep == 0,
0037     Nstep = 1;
0038     Zlist = [Zmax , Zmax + Zstep/2];
0039 else
0040     Zstep = (Zmax-Zmin)/Nstep;
0041     Zlist = Zmin:Zstep:Zmax;
0042     Zlist(Nstep+1) = Zmax + Zstep/2;
0043 end
0044 
0045 if ~exist('Rmax','var') | isempty(Rmax),  Rmax  = 0.003; end; % [m]
0046 
0047 % Find nearest point from slice region
0048 % Loop for slices in z-direction
0049 for n=1:Nstep,
0050     % Find vertex points in this slice
0051     inew  = find( Vnew(:,3) >= Zlist(n) ...
0052                 & Vnew(:,3) <  Zlist(n+1) );
0053     iold  = find( Vold(:,3) >= (Zlist(n)   - Rmax) ...
0054                 & Vold(:,3) <  (Zlist(n+1) + Rmax) );
0055 
0056     Nnew  = length(inew);
0057     Nold2 = length(iold);
0058     Vold2 = Vold(iold,:);
0059     
0060 %    fprintf('Searching for %d vs %d points in %d-th slice \n',Nnew,Nold2,n)
0061 
0062     if Nold2 == 0, ddmin(inew) = Rmax; continue; end;
0063     
0064     for i=1:Nnew,
0065         k      = inew(i);            % New vertex index
0066         Vnew2 = Vnew( k,:);
0067         
0068         % Distance between new and old vertex
0069         dd    = (Vold2(:,1) - Vnew2(1)).^2 ...
0070             + (Vold2(:,2) - Vnew2(2)).^2 ...
0071             + (Vold2(:,3) - Vnew2(3)).^2;
0072             
0073         % Find nearest old vertex
0074         ddmin(k)  = min(dd); 
0075     end;
0076 
0077 end
0078 
0079 % Find min distance is larger than Rmax
0080 ixmax = find( ddmin >= Rmax );
0081 
0082 Nmax  = length(ixmax);
0083 
0084 %fprintf('Searching large distance %d vs %d points\n',Nmax,Nold)
0085 
0086 % Find nearest point from all vertex
0087 for i=1:Nmax,
0088     k      = ixmax(i);            % New vertex index
0089     Vnew2 = Vnew( k,:);
0090     
0091     % Distance between new and old vertex
0092     dd    = (Vold(:,1) - Vnew2(1)).^2 ...
0093         + (Vold(:,2) - Vnew2(2)).^2 ...
0094         + (Vold(:,3) - Vnew2(3)).^2;
0095         
0096     % Find nearest old vertex
0097     ddmin(k)  = min(dd); 
0098 end;
0099 
0100 ddsum = sum(ddmin)/Npoint;