0001 function [Indx ,ddmin] = vb_find_nearest_point(Vold,Vnew,Rmax,Nstep,Mode,Disp)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031 if ~exist('Nstep','var')|isempty(Nstep), Nstep = 100; end;
0032 if ~exist('Disp','var'), Disp = 0; end;
0033
0034
0035 Npoint = size(Vnew,1);
0036 Indx = zeros(Npoint,1);
0037 ddmin = zeros(Npoint,1);
0038 inew = zeros(Npoint,1);
0039
0040 Nold = size(Vold,1);
0041 dd = zeros(Nold,1);
0042 iold = zeros(Nold,1);
0043 Vold2 = zeros(Nold,3);
0044
0045
0046 Zmax = max(Vnew(:,3));
0047 Zmin = min(Vnew(:,3));
0048 Zstep = (Zmax-Zmin)/Nstep;
0049
0050 if Zstep == 0
0051 Zstep = 0.001;
0052 Zlist = [Zmin , Zmin + Zstep/2];
0053 Nstep = 1;
0054 else
0055 Zlist = Zmin:Zstep:Zmax;
0056 Zlist(Nstep+1) = Zmax + Zstep/2;
0057 end
0058
0059 if ~exist('Rmax','var')|isempty(Rmax), Rmax = Zstep; end;
0060
0061 if Disp == 1,
0062 prg = 0;
0063 prg_all = Nstep;
0064 h = waitbar(0,'Nearest point search');
0065 vb_disp_nonl(sprintf('%3d %% processed',ceil(100*(prg/prg_all))));
0066 end
0067
0068
0069
0070 for n=1:Nstep,
0071
0072 inew = find( Vnew(:,3) >= Zlist(n) ...
0073 & Vnew(:,3) < Zlist(n+1) );
0074 iold = find( Vold(:,3) >= (Zlist(n) - Rmax) ...
0075 & Vold(:,3) < (Zlist(n+1) + Rmax) );
0076
0077 Nnew = length(inew);
0078 Nold2 = length(iold);
0079 Vold2 = Vold(iold,:);
0080
0081 switch Disp
0082 case 1,
0083 waitbar(n/Nstep);
0084 for ii=1:15; vb_disp_nonl(sprintf('\b')); end
0085 vb_disp_nonl(sprintf('%3d %% processed',ceil(100*(prg/prg_all))));
0086 prg = prg+1;
0087 case 2,
0088 fprintf('Searching for %d vs %d points in %d-th slice \n',Nnew,Nold2,n)
0089 end
0090
0091 if Nold2 == 0, ddmin(inew) = Rmax; continue; end;
0092
0093 for i=1:Nnew,
0094 k = inew(i);
0095 Vnew2 = Vnew( k,:);
0096
0097
0098 dd = (Vold2(:,1) - Vnew2(1)).^2 ...
0099 + (Vold2(:,2) - Vnew2(2)).^2 ...
0100 + (Vold2(:,3) - Vnew2(3)).^2;
0101
0102
0103 [md,ix] = min(dd);
0104 Indx(k) = iold(ix);
0105 ddmin(k) = md;
0106 end;
0107
0108 end
0109
0110 ddmin = sqrt(ddmin);
0111
0112 if Disp == 1,
0113 close(h);
0114 drawnow;
0115 vb_disp_nonl(sprintf('\n'));
0116 end
0117
0118 if exist('Mode','var') & Mode == 0; return; end;
0119
0120
0121 ixmax = find( ddmin >= Rmax | Indx == 0 );
0122 Nmax = length(ixmax);
0123
0124 if Nmax==0, return; end;
0125
0126 switch Disp
0127 case 1,
0128 fprintf('Searching large distance %d vs %d points\n',Nmax,Nold)
0129 h = waitbar(0,'Distant point search');
0130 prg = 0;
0131 prg_all = Nmax;
0132 vb_disp_nonl(sprintf('%3d %% processed',ceil(100*(prg/prg_all))));
0133 case 2,
0134 fprintf('Searching large distance %d vs %d points\n',Nmax,Nold)
0135 end
0136
0137
0138 for i=1:Nmax,
0139 k = ixmax(i);
0140 Vnew2 = Vnew( k,:);
0141
0142
0143 dd = (Vold(:,1) - Vnew2(1)).^2 ...
0144 + (Vold(:,2) - Vnew2(2)).^2 ...
0145 + (Vold(:,3) - Vnew2(3)).^2;
0146
0147
0148 [md,ix] = min(dd);
0149 Indx(k) = ix;
0150 ddmin(k) = md;
0151
0152 if Disp == 1,
0153 waitbar(i/Nmax);
0154 for ii=1:15; vb_disp_nonl(sprintf('\b')); end
0155 vb_disp_nonl(sprintf('%3d %% processed',ceil(100*(prg/prg_all))));
0156 prg = prg+1;
0157 end
0158 end
0159
0160 ddmin(ixmax) = sqrt(ddmin(ixmax));
0161
0162 if Disp == 1,
0163 close(h);
0164 drawnow;
0165 vb_disp_nonl(sprintf('\n'));
0166 end
0167
0168 clear inew dd iold Vold2
0169
0170 return;