% General testing program % Calls various projection methods, forms histograms of bad cases clear all; warning off all; Window = 1; % size of rectangle containing all data points n = 2; % number of data points histCols = 20; % number of columns in histograms (maximum recordable precision) histRows = 6; % number of rows in histograms (number of projection methods) Emin = 10.^(-histCols); % min allowed error hist1=zeros(histRows,histCols); hist2=zeros(histRows,histCols); % files for histograms N=0; N0=0; N1=0; N2=0; N3=0; Nmiss=0; % various counters of randomly generated conics Times = zeros(histRows,1); % for timing E1 = zeros(histRows,n); E1 = zeros(histRows,n); XYbest = zeros(n,2); XYprojAll = zeros(n,2,histRows); while(1) % main loop over random conics and points A = randn(6,1); % random conic parameters % the following line makes the conic nearly-parabolic (ill-conditioned) % this line can be commented out when generating typical conics %e=1e-13; A(3) = A(2)^2/A(1) + e*randn; % make it a parabola (if desired) A = A/norm(A); % normalization (unnecessary) % discard the conic if it is of the wrong type: if (IsCrossingWindow(A,Window) == 0) % discard the conic if it doed not Nmiss = Nmiss + 1; % cross the data window continue; end; [ParG,code] = AtoG(A); % ParG is vector of geometric parameters of the conic if (code > 3) % discard degenerate conics N0 = N0 + 1; continue; end; N = N + 1; if (code==1), N1 = N1 + 1; end; % count ellipses if (code==2), N2 = N2 + 1; end; % count hyperbolas if (code==3), N3 = N3 + 1; end; % count parabolas XY = 2*Window*(rand(n,2)-0.5*ones(n,2)); % random points Dmin = realmax*ones(n,1); tic XYproj = ProjectPointsOntoConicByWEPqz(XY,A); % 1st projection method Times(1) = Times(1) + toc; E2(1,:) = ProjectionQuality(XY,XYproj,A); D = DistanceToConicApprx(XY,XYproj,A); for i=1:n if Dmin(i) > D(i) Dmin(i) = D(i); XYbest(i,:) = XYproj(i,:); end end XYprojAll(:,:,1) = XYproj; tic XYproj = ProjectPointsOntoConicByWEPinv(XY,A); % 2nd projection method Times(2) = Times(2) + toc; E2(2,:) = ProjectionQuality(XY,XYproj,A); D = DistanceToConicApprx(XY,XYproj,A); for i=1:n if Dmin(i) > D(i) Dmin(i) = D(i); XYbest(i,:) = XYproj(i,:); end end XYprojAll(:,:,2) = XYproj; tic XYproj = ProjectPointsOntoConicByWEPeq3Schur(XY,A); % 3rd projection method Times(3) = Times(3) + toc; E2(3,:) = ProjectionQuality(XY,XYproj,A); D = DistanceToConicApprx(XY,XYproj,A); for i=1:n if Dmin(i) > D(i) Dmin(i) = D(i); XYbest(i,:) = XYproj(i,:); end end XYprojAll(:,:,3) = XYproj; tic XYproj = ProjectPointsOntoConicByEberlyModified(XY,A); % 4th projection method Times(4) = Times(4) + toc; E2(4,:) = ProjectionQuality(XY,XYproj,A); D = DistanceToConicApprx(XY,XYproj,A); for i=1:n if Dmin(i) > D(i) Dmin(i) = D(i); XYbest(i,:) = XYproj(i,:); end end XYprojAll(:,:,4) = XYproj; tic XYproj = ProjectPointsOntoConicByEberlyOrig(XY,A); % 5th projection method Times(5) = Times(5) + toc; E2(5,:) = ProjectionQuality(XY,XYproj,A); D = DistanceToConicApprx(XY,XYproj,A); for i=1:n if Dmin(i) > D(i) Dmin(i) = D(i); XYbest(i,:) = XYproj(i,:); end end XYprojAll(:,:,5) = XYproj; tic XYproj = ProjectPointsOntoConicByQuarticEquation(XY,A); % 6th projection method Times(6) = Times(6) + toc; E2(6,:) = ProjectionQuality(XY,XYproj,A); D = DistanceToConicApprx(XY,XYproj,A); for i=1:n if Dmin(i) > D(i) Dmin(i) = D(i); XYbest(i,:) = XYproj(i,:); end end XYprojAll(:,:,6) = XYproj; % Record the results of all projection methods % first, for each projected point compute the distance to the best projection for i=1:n for j=1:histRows E1(j,i) = norm(XYprojAll(i,:,j)-XYbest(i,:)); E1(j,i) = max(E1(j,i),Emin); E2(j,i) = max(E2(j,i),Emin); end end % next for each projected point record the errors computed in two ways: % E1 stores the distance to the best projected point % E2 stores the quality of the projection (E in the article) % update the histograms for i=1:n for j=1:histRows k = floor(-log10(E1(j,i))); if (k<1), k=1; end if (k<=histCols), hist1(j,k) = hist1(j,k) + 1;; end k = floor(-log10(E2(j,i))); if (k<1), k=1; end if (k<=histCols), hist2(j,k) = hist2(j,k) + 1;; end end end % print out the results if mod(N*n,10000)==0 fprintf(1,'\n Histogram of distances to the best projection point:\n\n'); for j=1:histRows for k=1:histCols fprintf(1,' %5d',hist1(j,k)); end fprintf(1,'\n'); end fprintf(1,'\n Histogram of qualities of projections:\n\n'); for j=1:histRows for k=1:histCols fprintf(1,' %5d',hist2(j,k)); end fprintf(1,'\n'); end fprintf(1,'\n Timing:') for j=1:histRows fprintf(1,' %8.5f',Times(j)/N); end fprintf(1,'\n\n conics used: %d (Ell: %d Hyp: %d Par: %d)\n',N,N1,N2,N3) fprintf(1,'\n conics removed: %d degenerate and %d missing the window\n',N0,Nmiss) end if (N*n>=100000000) return; end end