Index: Damkjer/PointProcessing/SpatialAnalyzer/nbrspanalyze.m =================================================================== --- Damkjer/PointProcessing/SpatialAnalyzer/nbrspanalyze.m (revision 13) +++ Damkjer/PointProcessing/SpatialAnalyzer/nbrspanalyze.m (revision 14) @@ -75,9 +75,11 @@ % - maketimebar = false; - - if (~isempty(varargin) && strcmp('timebar',varargin{1})) - maketimebar = true; - end + userParams = parseInputs(varargin{:}); % Parse Inputs + +% maketimebar = false; +% +% if (~isempty(varargin) && strcmp('timebar',varargin{1})) +% maketimebar = true; +% end % Sizing metrics @@ -88,9 +90,22 @@ norms=zeros(dimensions, elements); feats=zeros(dimensions, elements); - biases=zeros(1,elements); + dists2means=zeros(1,elements); ints=zeros(1,elements); + if (~isempty(userParams.test)) + if (length(neighborhoods) ~= length(userParams.test)) + error('Damkjer:nbrspanalyze:size', ... + ['When test sets are provided, ' ... + 'there must be one for each point set.']); + end + + neighborhoods=cellfun(@(x,y) [x;y], userParams.test, neighborhoods, ... + 'UniformOutput', false); +% neighborhoods=cellfun(@(x,y) [x(1);y], userParams.test, neighborhoods, ... +% 'UniformOutput', false); + end + % This process may take a while. Display time bar while processing. - if (maketimebar) + if (userParams.maketimebar) msg='Analyzing Neighborhoods...'; tstart=tic; @@ -101,5 +116,5 @@ step=1000; - if (maketimebar) + if (userParams.maketimebar) tic; end @@ -117,6 +132,15 @@ end - [covs,bias,inty]=fastcov(cells); - + [covs,dist2mean,inty]=fastcov(cells); + +% points(:,1:10) +% NN{1:10} +% cells{1:10} +% covs{1:10} +% dist2mean{1:10} +% inty{1:10} +% +% error('stop'); + % Compute the neighborhood covariance eigenvalues. [V,D]=par_eig(covs); @@ -142,10 +166,10 @@ feats(:,elem+nbr-1)=sqrt(max(sort(D{nbr},'descend'),0)); - biases(elem+nbr-1)=bias{nbr}; + dists2means(elem+nbr-1)=dist2mean{nbr}; ints(elem+nbr-1)=inty{nbr}; end % Update the time bar. - if (maketimebar && toc > 1) + if (userParams.maketimebar && toc > 1) tic; h = timebar(elem, elements, msg, tstart, h); @@ -154,5 +178,5 @@ % Close the time bar, if still open. - if (maketimebar && all(ishghandle(h, 'figure'))) + if (userParams.maketimebar && all(ishghandle(h, 'figure'))) close(h); end @@ -160,6 +184,7 @@ % Compute Normalized Features energy=sum(feats,1); - nfeats=bsxfun(@times,feats',1./energy')'; + nfeats=bsxfun(@times,feats.',1./energy.').'; nfeats(isnan(nfeats))=1; + nfeats(nfeats==0)=1; % Compute Omnivariance @@ -170,5 +195,6 @@ % Compute Eigen Entropy entropy=-sum(nfeats.*log(nfeats))./log(dimensions); - entropy(entropy==1)=0; +% entropy(entropy==1)=0; + entropy(isnan(entropy))=0; % Compute Eigen Structure @@ -200,14 +226,16 @@ % alpha=[dims;iso]; alpha=dims; + alpha(isnan(alpha))=1; + alpha(alpha==0)=1; de=-sum(alpha.*log(alpha))./log(dimensions); - if (~isreal(de)) - disp('de...'); - disp(de); - disp('feats...'); - disp(feats); - disp('eigs...'); - disp(feats.*feats); - end +% if (~isreal(de)) +% disp('de...'); +% disp(de); +% disp('feats...'); +% disp(feats); +% disp('eigs...'); +% disp(feats.*feats); +% end % Populate feature classes @@ -222,7 +250,87 @@ classes.omnivariance=omnivar; classes.labeling=label; - classes.biases=biases; + classes.dists2means=dists2means; classes.intensity=ints; classes.de=de; end + +%****************************************************************************** +% parseInputs +%****************************************************************************** +function [userParams] = parseInputs(varargin) + % Support function to parse inputs into userParams structure + % + % Parameters: + % varargin - Variable-length input argument list. Option strings may be + % provided as abbreviations as long as they resolve to a + % unique option.Defined key-value pairs include the following: + % 'neighbors' - . + % 'radius' - . + % 'counts' - . + % 'steps' - . + % + % Returns: + % userParams - . + + userParams = struct('maketimebar', false, 'test', []); + + % Parse the Property/Value pairs +% if rem(length(varargin), 2) ~= 0 +% error('Damkjer:PropertyValueNotPair', ... +% ['Additional arguments must take the form of Property/Value '... +% 'pairs']); +% end + + propertyNames = {'timebar', 'test'}; + + while ~isempty(varargin) + + if islogical(varargin{1}) + userParams.maketimebar = varargin{1}; + varargin(1) = []; + elseif iscell(varargin{1}) + userParams.test = varargin{1}; + varargin(1) = []; + elseif (ischar(varargin{1}) && ... + (length(varargin) < 2 || ischar(varargin{2})) && ... + strcmp('timebar',varargin{1})) + userParams.maketimebar = true; + varargin(1) = []; + else + property = varargin{1}; + value = varargin{2}; + + % If the property has been supplied in a shortened form, lengthen it + iProperty = find(strncmpi(property, propertyNames, length(property))); + + if isempty(iProperty) + error('Damkjer:InvalidProperty', 'Invalid Property'); + elseif length(iProperty) > 1 + error('Damkjer:AmbiguousProperty', ... + 'Supplied shortened property name is ambiguous'); + end + + property = propertyNames{iProperty}; + + switch property + case 'timebar' + if (islogical(value)) + userParams.maketimebar = value; + else + error('Damkjer:InvalidCount', ... + 'Timebar must be called with a logical parameter'); + end + case 'test' + if (isnumeric(value)) + userParams.test = fix(value); + else + error('Damkjer:InvalidRadius', ... + 'Test must be called with a set of indices'); + end + end + + varargin(1:2) = []; + end + end +end Index: Damkjer/PointProcessing/SpatialAnalyzer/spanalyze.m =================================================================== --- Damkjer/PointProcessing/SpatialAnalyzer/spanalyze.m (revision 13) +++ Damkjer/PointProcessing/SpatialAnalyzer/spanalyze.m (revision 14) @@ -79,4 +79,5 @@ % Sizing metrics dimensions=size(queries, 1); +% dimensions=3; elements=size(queries, 2); @@ -84,9 +85,12 @@ norms=zeros(size(queries)); feats=zeros(size(queries)); +% norms=zeros(3,length(queries)); +% feats=zeros(3, length(queries)); biases=zeros(1,elements); ints=zeros(1,elements); % Build database - database=VpTree(points); +% database=VpTree(points); + database=VpTree(points(1:3,:)); % This process may take a while. Display time bar while processing. @@ -108,27 +112,27 @@ if (userParams.radius > 0 && userParams.neighbors <= 0) % Perform a fixed radius search - NN=database.rnn(queries(:,elem:last),... + NN=database.rnn(queries(1:3,elem:last),... userParams.radius); elseif (userParams.radius <= 0 && userParams.neighbors > 0) % Search unconstrained neighbors - NN=database.knn(queries(:,elem:last),... + NN=database.knn(queries(1:3,elem:last),... userParams.neighbors); elseif (userParams.radius > 0 && userParams.neighbors > 0) % Search constrained to radius - NN=database.knn(queries(:,elem:last),... + NN=database.knn(queries(1:3,elem:last),... userParams.neighbors,... 'lim',userParams.radius); elseif (~isempty(userParams.counts) && userParams.radius <= 0) % Search unconstrained neighbors - NN=database.knn(queries(:,elem:last),... + NN=database.knn(queries(1:3,elem:last),... max(userParams.counts)); elseif (~isempty(userParams.counts) && userParams.radius > 0) % Search constrained to radius - NN=database.knn(queries(:,elem:last),... + NN=database.knn(queries(1:3,elem:last),... max(userParams.counts),... 'lim',userParams.radius); elseif (~isempty(userParams.steps)) % Perform a fixed radius search - [NN,DISTS]=database.rnn(queries(:,elem:last),... + [NN,DISTS]=database.rnn(queries(1:3,elem:last),... max(userParams.steps)); end @@ -141,8 +145,18 @@ for n=1:length(NN) cells{n}=points(:,NN{n})'; +% cells{n}=points(4:6,NN{n})'; end [covs,bias,inty]=fastcov(cells); +% points(:,1:10) +% NN{1:10} +% cells{1:10} +% covs{1:10} +% bias{1:10} +% inty{1:10} +% +% error('stop'); + % elseif (~isempty(userParams.counts)) % @@ -284,5 +298,5 @@ % Option 2: Define features as singular values. Eigenvalues are not % guaranteed to be in any order. We want them in descending order. - feats(:,elem+nbr-1)=sqrt(sort(D{nbr},'descend')); + feats(:,elem+nbr-1)=sqrt(max(sort(D{nbr},'descend'),0)); biases(elem+nbr-1)=bias{nbr}; @@ -304,6 +318,13 @@ % Compute Normalized Features energy=sum(feats,1); - nfeats=bsxfun(@times,feats',1./energy')'; + + nfeats=bsxfun(@times,feats.',1./energy.').'; nfeats(isnan(nfeats))=1; + nfeats(nfeats==0)=1; + +% max(energy) +% gnfeats=feats./max(energy); + +% error('Stop'); % Compute Omnivariance @@ -314,5 +335,8 @@ % Compute Eigen Entropy entropy=-sum(nfeats.*log(nfeats))./log(dimensions); - entropy(entropy==1)=0; +% entropy(entropy==1)=0; + entropy(isnan(entropy))=0; + +% global_entropy=-sum(gnfeats.*log(gnfeats))./log(dimensions); % Compute Eigen Structure @@ -344,5 +368,18 @@ % alpha=[dims;iso]; alpha=dims; + alpha(isnan(alpha))=1; + alpha(alpha==0)=1; de=-sum(alpha.*log(alpha))./log(dimensions); + + % Compute dimensional degree +% gdims=zeros(size(feats)); + +% max(feats(1,:)) + +% gdims(1:dimensions-1,:)=(gnfeats(1:dimensions-1,:)-gnfeats(2:dimensions,:))./... +% repmat(gnfeats(1,:),dimensions-1,1); +% gdims(dimensions,:)=gnfeats(dimensions,:)./(gnfeats(1,:)); + +% global_de=-sum(gdims.*log(gdims))./log(dimensions); % Populate feature classes @@ -353,4 +390,5 @@ classes.anisotropy=ani; classes.FA=fa; + classes.energy=energy; classes.entropy=entropy; classes.structure=structure; @@ -360,4 +398,7 @@ classes.intensity=ints; classes.de=de; +% classes.global_dimensionality=gdims; +% classes.global_entropy=global_entropy; +% classes.global_de=global_de; % TODO: Add curvature analysis back into spatial analyzer. Will require Index: Damkjer/PointProcessing/Thinning_Damkjer/thin_Damkjer.m =================================================================== --- Damkjer/PointProcessing/Thinning_Damkjer/thin_Damkjer.m (revision 13) +++ Damkjer/PointProcessing/Thinning_Damkjer/thin_Damkjer.m (revision 14) @@ -75,17 +75,32 @@ userParams = parseInputs(varargin{:}); % Parse Inputs - % Make sure we have unique points... we barf otherwise... + %*** + % CONFIRM UNIQUENESS + %*** + + % Make sure that all points in the set are unique... disp('Culling duplicate points...'); tstart=tic; - points = unique(points.','rows').'; + check = unique(points.','rows').'; disp(['...done in ' num2str(toc(tstart)) 's']); - + + if (any(size(check) ~= size(points))) + error('Points must be unique. Try: unique(X.'',''rows'',''stable'').'''); + end + + %*** + % INDEX POINT SET + %*** disp('Indexing points...'); tstart=tic; - database = VpTree(points); + database = VpTree(points(1:3,:)); queries = points; disp(['...done in ' num2str(toc(tstart)) 's']); - database.excludeWithin(0.001); % Set mm precision on searches. + %*** + % CREATE NEIGHBORHOODS + %*** + + database.excludeWithin(0.000000001); % Set um precision on searches. % Sizing metrics @@ -94,8 +109,6 @@ nbrs = cell(1,elements); - - %*** - % CREATE NEIGHBORHOODS - %*** + dists = cell(1,elements); + test = num2cell(1:elements, 1); % This process may take a while. Display time bar while processing. @@ -119,28 +132,34 @@ if (userParams.radius > 0 && userParams.neighbors <= 0) % Perform a fixed radius search - nbrs(elem:last)=database.rnn(queries(:,elem:last),... - userParams.radius); + [nbrs(elem:last), ... + dists(elem:last)] = database.rnn(queries(1:3,elem:last), ... + userParams.radius); elseif (userParams.radius <= 0 && userParams.neighbors > 0) % Search unconstrained neighbors - nbrs(elem:last)=database.knn(queries(:,elem:last),... - userParams.neighbors); + [nbrs(elem:last), ... + dists(elem:last)] = database.knn(queries(1:3,elem:last), ... + userParams.neighbors); elseif (userParams.radius > 0 && userParams.neighbors > 0) % Search constrained to radius - nbrs(elem:last)=database.knn(queries(:,elem:last),... - userParams.neighbors,... - 'lim',userParams.radius); + [nbrs(elem:last), ... + dists(elem:last)] = database.knn(queries(1:3,elem:last), ... + userParams.neighbors, ... + 'lim', userParams.radius); elseif (~isempty(userParams.counts) && userParams.radius <= 0) % Search unconstrained neighbors - nbrs(elem:last)=database.knn(queries(:,elem:last),... - max(userParams.counts)); + [nbrs(elem:last), ... + dists(elem:last)] = database.knn(queries(1:3,elem:last), ... + max(userParams.counts)); elseif (~isempty(userParams.counts) && userParams.radius > 0) % Search constrained to radius - nbrs(elem:last)=database.knn(queries(:,elem:last),... - max(userParams.counts),... - 'lim',userParams.radius); + [nbrs(elem:last), ... + dists(elem:last)] = database.knn(queries(1:3,elem:last), ... + max(userParams.counts), ... + 'lim', userParams.radius); elseif (~isempty(userParams.steps)) % Perform a fixed radius search - [nbrs(elem:last),DISTS]=database.rnn(queries(:,elem:last),... - max(userParams.steps)); + [nbrs(elem:last), ... + dists(elem:last)] = database.rnn(queries(1:3,elem:last), ... + max(userParams.steps)); end @@ -159,41 +178,6 @@ disp(['...done in ' num2str(toc(tstart)) 's']); -%% Obsolete code - %*** - % COMPUTE DUAL NEIGHBORHOODS (THIS IS SLOW. LET'S MEX IT) - %*** - -% % This process may take a while. Display time bar while processing. -% msg='Computing Dual Neigborhoods...'; -% tstart=tic; -% h = timebar(1, elements, msg, tstart); -% -% disp(msg); -% -% tic; -% for elem=1:elements -% -% % Build up dual set -% for nbr=1:length(nbrs{elem}) -% n = nbrs{elem}(nbr); -% duals{n}=[duals{n}(:); elem]; -% end -% -% % Update the time bar. -% if (toc > 1) -% tic; -% h = timebar(elem, elements, msg, tstart, h); -% end -% end -% -% % Close the time bar, if still open. -% if (all(ishghandle(h, 'figure'))) -% close(h); -% end -% -% disp(['...done in ' num2str(toc(tstart)) 's']); -%% - %*** - % COMPUTE DUAL NEIGHBORHOODS -- MEX'D + %*** + % COMPUTE DUAL NEIGHBORHOODS %*** @@ -204,16 +188,79 @@ %*** + % COMPUTE NATIVE ATTRIBUTES + %*** + disp('Computing Native Attributes...'); + tstart=tic; + + natives = nbrspanalyze(points, nbrs, test, 'timebar'); +% natives = nbrspanalyze(points, nbrs, 'timebar'); + +% if (any(isnan(natives.entropy))) +% disp('Ugh... seriously?'); +% end +% + disp(['...done in ' num2str(toc(tstart)) 's']); + + %*** + % COMPUTE NEIGHBORHOOD ATTRIBUTES + %*** + disp('Computing Neighborhood Attributes...'); + tstart=tic; + + feats = nbrspanalyze(points, nbrs, 'timebar'); + +% if (any(isnan(feats.entropy))) +% disp('Ugh... seriously?'); +% end +% + disp(['...done in ' num2str(toc(tstart)) 's']); + + %*** % COMPUTE SIGNIFICANCES %*** - disp('Computing Significances...'); + sigs=zeros(1,elements); + + msg='Computing Significances...'; + + disp(msg); + tstart=tic; - - feats = nbrspanalyze(points, nbrs, 'timebar'); - + h = timebar(1, elements, msg, tstart); + + tic; + + for elem=1:elements + sigs(elem)=abs(feats.entropy(elem)-natives.entropy(elem)); + + % Update the time bar. + if (toc > 1) + tic; + h = timebar(elem, elements, msg, tstart, h); + end + end + + % Close the time bar, if still open. + if (all(ishghandle(h, 'figure'))) + close(h); + end + +% if (any(isnan(sigs))) +% disp('Ugh... seriously?'); +% end +% disp(['...done in ' num2str(toc(tstart)) 's']); - sigs = feats.de; -% sigs = sqrt(feats.biases.*feats.de); -% sigs = 0.5 * (feats.biases + feats.de); +% origFormat = get(0, 'format'); +% format('longG'); +% min(sigs) +% max(sigs) +% set(0, 'format', origFormat); + +% results = natives; +% end +% +% function cont +% error('Stop'); + %*** @@ -221,17 +268,17 @@ %*** - disp('Partition by Label...'); - tstart=tic; - - parts = cell(1,dimensions); - - for dim=1:dimensions - parts{dim}=uint64(find(feats.labeling==dim)); - end - - [~,idx]=sort(cellfun('length',parts),'ascend'); - parts=parts(idx); - - disp(['...done in ' num2str(toc(tstart)) 's']); +% disp('Partition by Label...'); +% tstart=tic; +% +% parts = cell(1,dimensions); +% +% for dim=1:dimensions +% parts{dim}=sort(uint64(find(feats.labeling==dim)),'ascend'); +% end +% +% [~,idx]=sort(cellfun('length',parts),'descend'); +% parts=parts(idx); +% +% disp(['...done in ' num2str(toc(tstart)) 's']); %*** @@ -239,5 +286,5 @@ %*** - fraction=0; +% fraction=0; removed=zeros(1,elements-threshold); @@ -245,24 +292,24 @@ current=1; - removed_wgts=cell(1,3); +% removed_wgts=cell(1,3); fellback = 0; - for dim=1:dimensions - desired = threshold * length(parts{dim}) / elements + fraction; - tau = floor(desired); - fraction = desired-tau; - heap = SplayTree(sigs(parts{dim}),parts{dim}); - n = length(parts{dim})-tau; - - removed_wgts{dim}=zeros(1,n); - -% tau=threshold; -% heap = SplayTree(feats.de); -% n = length(feats.de) - tau; -% removed_wgts=zeros(1,n); - - msg=['Thin Dimension ' num2str(dim) '(' num2str(length(parts{dim})) ... - '->' num2str(tau) ')...']; -% msg='Thin...'; +% for dim=1:dimensions +% desired = threshold * length(parts{dim}) / elements + fraction; +% tau = floor(desired); +% fraction = desired-tau; +% heap = SplayTree(sigs(parts{dim}),parts{dim}); +% n = length(parts{dim})-tau; + +% removed_wgts{dim}=zeros(1,n); + + tau=threshold; + heap = SplayTree(sigs,1:length(sigs)); + n = length(sigs) - tau; + removed_wgts=zeros(1,n); + +% msg=['Thin Dimension ' num2str(dim) '(' num2str(length(parts{dim})) ... +% '->' num2str(tau) ')...']; + msg='Thin...'; tstart=tic; h = timebar(1, n, msg, tstart); @@ -271,30 +318,119 @@ tic; +% disp(nbrs{318935}) +% disp(test{318935}) +% disp(duals{318935}) +% +% disp(nbrs{458702}) +% disp(test{458702}) +% disp(duals{458702}) + for x=1:n -% if false - - % Remove Element + +% disp('Pop'); + +% heap_size = heap.size(); + + % Locate the least significant point to remove [wgt, idx]=heap.pop_head(); +% disp('Popped'); + +% if (heap.size() ~= (heap_size - 1)) +% disp(['...and there it is: ' num2str(heap_size) ' - 1 vs. ' num2str(heap.size())]); +% disp(idx); +% error('oops'); +% end + +% if (isempty(nbrs{idx}) || isempty(duals{idx})) +% disp(idx); +% end + +% if (isempty(nbrs{idx}) && isempty(duals{idx})) +% disp('uh oh'); +% disp(idx); +% end +% idx +% wgt + +% if (any(removed==idx)) +% disp('uh oh'); +% disp(['Index already removed ' num2str(idx)]); +% end + + % Mark the point as removed removed(current)=idx; - removed_wgts{dim}(x)=wgt; - - % Find neighborhoods that contain removed element - neighborhoods=sort(duals{idx}, 'ascend'); - dual_size(current)=length(neighborhoods); +% removed_wgts{dim}(x)=wgt; + removed_wgts(x)=wgt; + + + % Distribute the points in the test set to the neighbor's test sets. + neighbors = nbrs{idx}; + samples = test{idx}; + +% neighbors +% samples + + for samp=1:length(samples) + % Find the closest +% [dsq,new_idx]=min(sum((bsxfun(@minus, points(:,neighbors), ... +% points(:,samples(samp)))).^2, 1)); + + [dsq,new_idx]=min(sum((bsxfun(@minus, points(1:3,neighbors), ... + points(1:3,samples(samp)))).^2, 1)); + +% samples(samp) +% points(1:3,neighbors) +% points(1:3,samples(samp)) +% dsq +% new_idx +% neighbors(new_idx) + + % Update the test set + test{neighbors(new_idx)}(end + 1,1) = samples(samp); + +% if (samples(samp) == 458702) +% disp(['relocating 458702 from test{' num2str(idx) '} to test{'... +% num2str(neighbors(new_idx)) '}']); +% +% disp(test{neighbors(new_idx)}); +% end +% +% if (neighbors(new_idx) == 458702) +% disp(['relocating ' num2str(samples(samp)) ' from test{' num2str(idx) '} to test{458702}']); +% +% disp(test{neighbors(new_idx)}); +% end + +% tdist(neighbors(new_idx)) = max(tdist(neighbors(new_idx)),... +% sqrt(dsq)); + end % The removed element neighborhood is eliminated, so it no longer % contains the neighbors. Update the neighbor duals appropriately. - neighbors=nbrs{idx}; - for nbr=1:length(neighbors) +% if (neighbors(nbr) == 458702) +% disp([num2str(idx) '''s neighborhood no longer contains 458702']); +% disp(duals{neighbors(nbr)}); +% end + duals{neighbors(nbr)}=duals{neighbors(nbr)}(duals{neighbors(nbr)}~=idx); - end - + +% if (neighbors(nbr) == 458702) +% disp(duals{neighbors(nbr)}); +% end + end + + % Find neighborhoods that contain removed element + neighborhoods=duals{idx}; + dual_size(current)=length(neighborhoods); + neighborhoods=sort(neighborhoods(neighborhoods~=idx), 'ascend'); + % The fallback neighborhood is identical for all neighborhoods, but % may be expensive to compute. Allow it to be lazily instantiated. fallback_pool = []; - % Replace removed element in each neighborhood with nearest neighbor + % Replace removed element in each neighborhood with nearest neighbor's + % neighbor for hood=1:length(neighborhoods) @@ -302,9 +438,4 @@ % disp('building candidates...'); -%% Obsolete code -% nbr_pool=unique(vertcat(nbrs{nbrs{neighborhoods(hood)}})); -% nbr_pool=setdiff(nbr_pool, idx); -% nbr_pool=setdiff(nbr_pool, nbrs{neighborhoods(hood)}); -%% nbr_pool=fastsetunion({nbrs{nbrs{neighborhoods(hood)}}}); non_nbrs=fastsetunion({nbrs{neighborhoods(hood)},idx,neighborhoods(hood)}); @@ -319,5 +450,4 @@ nbr_pool=nbr_pool(~ismembc(nbr_pool(:), non_nbrs(:))); end - % Find the closest @@ -331,72 +461,159 @@ % new_idx=nn{1} %% - [~,new_idx]=min(sum((bsxfun(@minus, points(:,nbr_pool), ... - points(:,idx))).^2, 1)); - +% points(1:3,nbr_pool) +% points(1:3,neighborhoods(hood)) +% + [dsq,new_idx]=min(sum((bsxfun(@minus, points(1:3,nbr_pool), ... + points(1:3,neighborhoods(hood)))).^2, 1)); + +% points(1:3,nbr_pool) +% points(1:3,neighborhoods(hood)) +% bsxfun(@minus, points(1:3,nbr_pool), points(1:3,neighborhoods(hood))) +% (bsxfun(@minus, points(1:3,nbr_pool), points(1:3,neighborhoods(hood)))).^2 +% sum((bsxfun(@minus, points(1:3,nbr_pool), points(1:3,neighborhoods(hood)))).^2,1) +% [dsq,new_idx]=min(sum((bsxfun(@minus, points(1:3,nbr_pool), points(1:3,neighborhoods(hood)))).^2,1)) + % Update the neighborhood % disp('updating neighborhood...'); - if (length(nbr_pool(new_idx)) ~= length(nbrs{neighborhoods(hood)}(nbrs{neighborhoods(hood)}==idx))) - neighborhoods - size(fallback_pool) - non_nbrs - nbr_pool - nbr_pool(new_idx) + new_dist = sqrt(dsq); + +% if (neighborhoods(hood) == 458702) +% disp(['replace ' num2str(idx) ' in 458702''s neighborhood']); +% disp('Pool...'); +% disp(nbr_pool(new_idx)); +% disp('Before...'); +% disp(nbrs{neighborhoods(hood)}); +% end +% +% if (nbr_pool(new_idx) == 458702) +% disp(['458702 is replacing ' idx ' in ' neighborhoods(hood) '''s neighborhood']); +% disp(nbrs{neighborhoods(hood)}); +% end +% + % Maintain a partially sorted list of indices by max distance. + if (new_dist > dists{neighborhoods(hood)}(end)) + dists{neighborhoods(hood)}(nbrs{neighborhoods(hood)}==idx) = ... + dists{neighborhoods(hood)}(end); - nbrs{neighborhoods(hood)}(nbrs{neighborhoods(hood)}==idx) + nbrs{neighborhoods(hood)}(nbrs{neighborhoods(hood)}==idx) = ... + nbrs{neighborhoods(hood)}(end); + + dists{neighborhoods(hood)}(end) = sqrt(dsq); - error('What happened?'); + nbrs{neighborhoods(hood)}(end) = nbr_pool(new_idx); + else + dists{neighborhoods(hood)}(nbrs{neighborhoods(hood)}==idx) = ... + sqrt(dsq); + + nbrs{neighborhoods(hood)}(nbrs{neighborhoods(hood)}==idx) = ... + nbr_pool(new_idx); end - nbrs{neighborhoods(hood)}(nbrs{neighborhoods(hood)}==idx) = ... - nbr_pool(new_idx); - + if (isempty(nbrs{neighborhoods(hood)})) + disp(neighborhoods(hood)); + end +% nbrs{neighborhoods(hood)} +% dists{neighborhoods(hood)} + % Update the dual % disp('updating dual...'); + + +% duals{nbr_pool(new_idx)} duals{nbr_pool(new_idx)}(end + 1,1) = neighborhoods(hood); - +% duals{nbr_pool(new_idx)} % disp(nbrs{neighborhoods(hood)}); % disp(duals{nbr_pool(new_idx)}); - end - + +% error('stop'); +% if (neighborhoods(hood) == 458702) +% disp('After...'); +% disp(nbrs{neighborhoods(hood)}); +% disp([num2str(nbr_pool(new_idx)) ' is now in 458702''s neighborhood']); +% disp(duals{nbr_pool(new_idx)}); +% end +% +% if (nbr_pool(new_idx) == 458702) +% disp(nbrs{neighborhoods(hood)}); +% disp(['458702 is now in ' num2str(nbr_pool(new_idx)) '''s neighborhood']); +% disp(duals{nbr_pool(new_idx)}); +% end + + end + + nbrs{idx}=[]; duals{idx}=[]; - + test{idx}=[]; + % Update analysis for elements % disp('updating analysis...'); - locfeats = nbrspanalyze(points, nbrs(neighborhoods)); - - locsigs = locfeats.de; -% locsigs = sqrt(locfeats.biases.*locfeats.de); -% locsigs = 0.5 * (locfeats.biases + locfeats.de); + +% if (isempty(neighborhoods)) +% warning('Damkjer:strange','empty neighborhoods'); +% end +% +% if (isempty(neighbors)) +% warning('Damkjer:stranger','empty neighbors'); +% end +% + if (isempty(neighborhoods)) + neighborhoods=neighbors; + end + +% if (~isempty(neighborhoods) && ~isempty(neighbors)) + neighborhoods=fastsetunion({neighborhoods, neighbors}); +% elseif (isempty(neighborhoods) && ~isempty(neighbors)) +% neighborhoods=neighbors; +% end + +% nbrs{neighborhoods} +% test{neighborhoods} + +% error('Stop'); + + if (~isempty(neighborhoods)) + locfeats = nbrspanalyze(points, ... + nbrs(neighborhoods)); +% locfeats = nbrspanalyze(points, ... +% nbrs(neighborhoods), ... +% test(neighborhoods)); + + locsigs=zeros(1,length(neighborhoods)); + + for hood=1:length(neighborhoods) + locsigs(hood)=max(abs(natives.entropy(test{neighborhoods(hood)})-locfeats.entropy(hood))); + end % Update heap, where required -% size(neighborhoods') -% -% if (size(sigs(neighborhoods)) ~= size(neighborhoods')) -% sigsize=size(sigs(neighborhoods)) -% nbrsize=size(neighborhoods') -% error('Stop'); -% end - -% sigs(neighborhoods) -% neighborhoods' -% locsigs - - updates=parts{dim}(ismember(parts{dim},neighborhoods)); - sig_upd=locsigs(ismember(neighborhoods,updates)); - - if (length(updates) ~= length(sig_upd)) - fellback - neighborhoods - updates - sig_upd - error('OK, seriously!'); - end - - heap.erase(sigs(updates),updates); - heap.insert(sig_upd,updates); - sigs(neighborhoods)=locsigs; - feats.de(neighborhoods)=locfeats.de; - feats.biases(neighborhoods)=locfeats.biases; +% updates=neighborhoods(ismembc(neighborhoods,parts{dim})).'; +% sig_upd=locsigs(ismembc(neighborhoods,updates)); + updates=neighborhoods'; + sig_upd=locsigs; + +% for pt=1:length(updates) +% if (any(removed==updates(pt))) +% disp('uh oh'); +% disp(['Index already removed ' num2str(updates(pt))]); +% end +% end + + +% heap_size = heap.size(); +% +% heap.erase(sigs(updates),updates); +% +% if (heap.size() ~= (heap_size - length(updates))) +% disp(['...and there it is: ' num2str(heap_size) ' - ' num2str(length(updates)) ' vs. ' num2str(heap.size())]); +% disp(sigs(updates)); +% disp(updates); +% error('oops'); +% end +% +% heap.insert(sig_upd,updates); + + heap.update(sig_upd,updates); + sigs(neighborhoods)=locsigs; + end current = current + 1; @@ -415,24 +632,23 @@ disp(['...done in ' num2str(toc(tstart)) 's']); - end +% end disp(['Generated fallback neighborhood pool ' num2str(fellback) ' time(s).']); - x_max=max(cellfun(@length, removed_wgts)); - y=nan(length(removed_wgts),x_max); - - for dim=1:dimensions - y(dim,1:length(removed_wgts{dim}))=removed_wgts{dim}; - end - -% x_max=length(removed_wgts); -% y = removed_wgts; +% x_max=max(cellfun(@length, removed_wgts)); +% y=nan(length(removed_wgts),x_max); + +% for dim=1:dimensions +% y(dim,1:length(removed_wgts{dim}))=removed_wgts{dim}; +% end + + x_max=length(removed_wgts); + y = removed_wgts; figure, plot(1:x_max,y); -% figure, semilogx(1:x_max,y); - figure, plot(1:length(dual_size),dual_size); - results = points(:,~ismember(1:size(points,2),removed)); +% results = points(:,~ismember(1:size(points,2),removed)); + results = removed; end Index: Damkjer/PointProcessing/Thinning_Damkjer/thin_Damkjer_old.m =================================================================== --- Damkjer/PointProcessing/Thinning_Damkjer/thin_Damkjer_old.m (revision 14) +++ Damkjer/PointProcessing/Thinning_Damkjer/thin_Damkjer_old.m (revision 14) @@ -0,0 +1,812 @@ +% Thin_Dyn Point Cloud Thinning (Dyn et al) +% +% File: thin_Dyn.m +% +% Description: +% Perform thinning according to Dyn 2004. +% +% Limitations: +% No known limitations. +% +% Synopsis: +% [results] = spanalyse(points, threshold, ...) +% +% Inputs: +% points - . +% threshold - must be positive. between 0-1, treated as percentage, > 1 +% points. +% +% Option strings may be provided as abbreviations as long as they resolve to +% a unique option. +% +% 'neighbors' - . +% 'radius' - . +% 'counts' - . +% 'steps' - . +% +% Outputs: +% results - Thinned point cloud. +% +% Toolbox requirements: +% None. +% +% Code requirements: +% None. +% +% Data requirements: +% None. +% +% References: +% Dyn reference. +% +% See Also: +% None. +% + +% Copyright (C) 2013 Kristian L. Damkjer. +% +% Software History: +% 2013-DEC-28 Initial coding. +% + +%****************************************************************************** +% thin_Dyn +%****************************************************************************** +function [ results ] = thin_Damkjer( points, threshold, varargin ) + % Perform local feature attribution via eigenspace analysis. + % + % Parameters: + % points - . + % + % threshold - . + % + % varargin - Variable-length input argument list. Option strings may be + % provided as abbreviations as long as they resolve to a + % unique option.Defined key-value pairs include the following: + % 'neighbors' - . + % 'radius' - . + % 'counts' - . + % 'steps' - . + % + % Returns: + % classes - Structure containing the spatial analysis classifications. + % + + userParams = parseInputs(varargin{:}); % Parse Inputs + + %*** + % CONFIRM UNIQUENESS + %*** + + % Make sure that all points in the set are unique... + disp('Culling duplicate points...'); + tstart=tic; + points = unique(points.','rows').'; + disp(['...done in ' num2str(toc(tstart)) 's']); + + %*** + % INDEX POINT SET + %*** + disp('Indexing points...'); + tstart=tic; + database = VpTree(points(1:3,:)); + queries = points; + disp(['...done in ' num2str(toc(tstart)) 's']); + + %*** + % CREATE NEIGHBORHOODS + %*** + + database.excludeWithin(0.000001); % Set um precision on searches. + + % Sizing metrics + dimensions=size(queries, 1); + elements=size(queries, 2); + + nbrs = cell(1,elements); + dists = cell(1,elements); + test = num2cell(1:elements, 1); + +% sigfeats=zeros(dimensions, elements); + + % This process may take a while. Display time bar while processing. + msg='Creating Neigborhoods...'; + tstart=tic; + h = timebar(1, elements, msg, tstart); + + disp(msg); + + % Step through the queries in chunks. + step=1000; + + tic; + + for elem=1:step:elements + + % The last available element may be closer than elem + step. + last=min(elem+step-1,elements); + + % Get the nearest neighbors of elem + if (userParams.radius > 0 && userParams.neighbors <= 0) + % Perform a fixed radius search + [nbrs(elem:last), ... + dists(elem:last)] = database.rnn(queries(1:3,elem:last), ... + userParams.radius); + elseif (userParams.radius <= 0 && userParams.neighbors > 0) + % Search unconstrained neighbors + [nbrs(elem:last), ... + dists(elem:last)] = database.knn(queries(1:3,elem:last), ... + userParams.neighbors); + elseif (userParams.radius > 0 && userParams.neighbors > 0) + % Search constrained to radius + [nbrs(elem:last), ... + dists(elem:last)] = database.knn(queries(1:3,elem:last), ... + userParams.neighbors, ... + 'lim', userParams.radius); + elseif (~isempty(userParams.counts) && userParams.radius <= 0) + % Search unconstrained neighbors + [nbrs(elem:last), ... + dists(elem:last)] = database.knn(queries(1:3,elem:last), ... + max(userParams.counts)); + elseif (~isempty(userParams.counts) && userParams.radius > 0) + % Search constrained to radius + [nbrs(elem:last), ... + dists(elem:last)] = database.knn(queries(1:3,elem:last), ... + max(userParams.counts), ... + 'lim', userParams.radius); + elseif (~isempty(userParams.steps)) + % Perform a fixed radius search + [nbrs(elem:last), ... + dists(elem:last)] = database.rnn(queries(1:3,elem:last), ... + max(userParams.steps)); + end + + % Update the time bar. + if (toc > 1) + tic; + h = timebar(elem, elements, msg, tstart, h); + end + end + + % Close the time bar, if still open. + if (all(ishghandle(h, 'figure'))) + close(h); + end + + disp(['...done in ' num2str(toc(tstart)) 's']); + +%% Obsolete code +% %*** +% % COMPUTE DUAL NEIGHBORHOODS (THIS IS SLOW. LET'S MEX IT) +% %*** +% +% % This process may take a while. Display time bar while processing. +% msg='Computing Dual Neigborhoods...'; +% tstart=tic; +% h = timebar(1, elements, msg, tstart); +% +% disp(msg); +% +% tic; +% for elem=1:elements +% +% % Build up dual set +% for nbr=1:length(nbrs{elem}) +% n = nbrs{elem}(nbr); +% duals{n}=[duals{n}(:); elem]; +% end +% +% % Update the time bar. +% if (toc > 1) +% tic; +% h = timebar(elem, elements, msg, tstart, h); +% end +% end +% +% % Close the time bar, if still open. +% if (all(ishghandle(h, 'figure'))) +% close(h); +% end +% +% disp(['...done in ' num2str(toc(tstart)) 's']); +%% + + %*** + % COMPUTE DUAL NEIGHBORHOODS + %*** + + disp('Computing Dual Neigborhoods...'); + tstart=tic; + duals=fastsetdual(nbrs); + disp(['...done in ' num2str(toc(tstart)) 's']); + + %*** + % COMPUTE NATIVE ATTRIBUTES + %*** + disp('Computing Native Attributes...'); + tstart=tic; + + natives = nbrspanalyze(points, nbrs, test, 'timebar'); + + disp(['...done in ' num2str(toc(tstart)) 's']); + +% error('Stop'); + +% radii=cellfun(@max,dists); +% biases=feats.dists2means./radii; +% tdist = zeros(1,elements); + + %*** + % COMPUTE NEIGHBORHOOD ATTRIBUTES + %*** + disp('Computing Neighborhood Attributes...'); + tstart=tic; + + feats = nbrspanalyze(points, nbrs, 'timebar'); + + disp(['...done in ' num2str(toc(tstart)) 's']); + + %*** + % COMPUTE SIGNIFICANCES + %*** + sigs=zeros(1,elements); + + msg='Computing Significances...'; + + disp(msg); + + tstart=tic; + h = timebar(1, elements, msg, tstart); + + tic; + +% nbrdims=zeros(dimensions,1); + + for elem=1:elements + +% % last=min(elem+step-1,elements); +% +% % feats.features(:,nbrs{elem}) +% % feats.features(:,nbrs{elem}).^2 +% % feats.features(:,nbrs{elem}).^2./length(nbrs{elem}) +% % sum(feats.features(:,nbrs{elem}).^2./length(nbrs{elem}), 2) +% % rss=sqrt(sum(feats.features(:,nbrs{elem}).^2./length(nbrs{elem}), 2)); +% +% % nbrdims(1:dimensions-1)=(rss(1:dimensions-1)-rss(2:dimensions))./... +% % repmat(rss(1),dimensions-1,1); +% % nbrdims(dimensions)=rss(dimensions)./rss(1); +% +% % de=-sum(nbrdims.*log(nbrdims))./log(dimensions); +% +% natives.de(elem) +% feats.de(elem) +% feats.de(elem)-natives.de(elem) +% abs(feats.de(elem)-natives.de(elem)) + +% sigs(elem)=abs(feats.de(elem)-natives.de(elem)); + sigs(elem)=abs(feats.entropy(elem)-natives.entropy(elem)); + +% % % Hausdorff DE is no-go +% % A=abs(bsxfun(@minus, feats.de(test{elem}), feats.de(nbrs{elem}))); +% % sigs(elem)=max(max(min(A,[],1)),max(min(A,[],2))); + +% error('Stop'); + + % Update the time bar. + if (toc > 1) + tic; + h = timebar(elem, elements, msg, tstart, h); + end + end + + % Close the time bar, if still open. + if (all(ishghandle(h, 'figure'))) + close(h); + end + + disp(['...done in ' num2str(toc(tstart)) 's']); + + min(sigs) + max(sigs) + +% results = sigs; +%end + +%function cont +% error('Stop'); + +% sigs = feats.de; +% sigs = sqrt(biases.*feats.de); +% sigs = 0.5 * (biases + feats.de); +% sigs = max(biases, feats.de); +% sigs = min(biases, feats.de); +% sigs = biases; + +% RMS DE for sigs in nbrhood and test? + + %*** + % PARTITION THE DATA SET + %*** + +% disp('Partition by Label...'); +% tstart=tic; +% +% parts = cell(1,dimensions); +% +% for dim=1:dimensions +% parts{dim}=sort(uint64(find(feats.labeling==dim)),'ascend'); +% end +% +% [~,idx]=sort(cellfun('length',parts),'descend'); +% parts=parts(idx); +% +% disp(['...done in ' num2str(toc(tstart)) 's']); + + %*** + % THIN BY LABEL + %*** + +% fraction=0; + + removed=zeros(1,elements-threshold); + dual_size=zeros(1,elements-threshold); + current=1; + +% removed_wgts=cell(1,3); + fellback = 0; + +% for dim=1:dimensions +% desired = threshold * length(parts{dim}) / elements + fraction; +% tau = floor(desired); +% fraction = desired-tau; +% heap = SplayTree(sigs(parts{dim}),parts{dim}); +% n = length(parts{dim})-tau; + +% removed_wgts{dim}=zeros(1,n); + + tau=threshold; + heap = SplayTree(sigs,1:length(sigs)); + n = length(sigs) - tau; + removed_wgts=zeros(1,n); + +% msg=['Thin Dimension ' num2str(dim) '(' num2str(length(parts{dim})) ... +% '->' num2str(tau) ')...']; + msg='Thin...'; + tstart=tic; + h = timebar(1, n, msg, tstart); + + disp(msg); + tic; + + for x=1:n + + % Locate the least significant point to remove + [wgt, idx]=heap.pop_head(); + +% idx +% wgt + + % Mark the point as removed + removed(current)=idx; +% removed_wgts{dim}(x)=wgt; + removed_wgts(x)=wgt; + + + % Distribute the points in the test set to the neighbor's test sets. + neighbors = nbrs{idx}; + samples = test{idx}; + +% neighbors +% samples + + for samp=1:length(samples) + % Find the closest +% [dsq,new_idx]=min(sum((bsxfun(@minus, points(:,neighbors), ... +% points(:,samples(samp)))).^2, 1)); + [dsq,new_idx]=min(sum((bsxfun(@minus, points(1:3,neighbors), ... + points(1:3,samples(samp)))).^2, 1)); + +% samples(samp) +% points(1:3,neighbors) +% points(1:3,samples(samp)) +% dsq +% new_idx +% neighbors(new_idx) + + % Update the test set + test{neighbors(new_idx)}(end + 1,1) = samples(samp); + +% tdist(neighbors(new_idx)) = max(tdist(neighbors(new_idx)),... +% sqrt(dsq)); + end + +% duals{neighbors} + + % The removed element neighborhood is eliminated, so it no longer + % contains the neighbors. Update the neighbor duals appropriately. + for nbr=1:length(neighbors) + duals{neighbors(nbr)}=duals{neighbors(nbr)}(duals{neighbors(nbr)}~=idx); + end + +% duals{neighbors} + + % Find neighborhoods that contain removed element + neighborhoods=duals{idx}; + dual_size(current)=length(neighborhoods); + neighborhoods=sort(neighborhoods(neighborhoods~=idx), 'ascend'); + + % The fallback neighborhood is identical for all neighborhoods, but + % may be expensive to compute. Allow it to be lazily instantiated. + fallback_pool = []; + + % Replace removed element in each neighborhood with nearest neighbor's + % neighbor + for hood=1:length(neighborhoods) + + % Build up the list of potential new neighbors +% disp('building candidates...'); + +%% Obsolete code +% nbr_pool=unique(vertcat(nbrs{nbrs{neighborhoods(hood)}})); +% nbr_pool=setdiff(nbr_pool, idx); +% nbr_pool=setdiff(nbr_pool, nbrs{neighborhoods(hood)}); +%% + nbr_pool=fastsetunion({nbrs{nbrs{neighborhoods(hood)}}}); + non_nbrs=fastsetunion({nbrs{neighborhoods(hood)},idx,neighborhoods(hood)}); + nbr_pool=nbr_pool(~ismembc(nbr_pool(:), non_nbrs(:))); + + if (isempty(nbr_pool)) + if (isempty(fallback_pool)) + fallback_pool = uint64(find(~ismember(1:size(points,2),removed))); + end + fellback = fellback + 1; + nbr_pool = fallback_pool; + nbr_pool=nbr_pool(~ismembc(nbr_pool(:), non_nbrs(:))); + end + + % Find the closest +% disp('finding closest...'); + +%% Obsolete code +% % Don't use VP Tree for single nearest neighbor... overhead not +% % worth it. +% locdb = VpTree(points(:,nbr_pool)); +% nn=locdb.knn(points(:,idx),1); +% new_idx=nn{1} +%% +% points(1:3,nbr_pool) +% points(1:3,neighborhoods(hood)) +% + [dsq,new_idx]=min(sum((bsxfun(@minus, points(1:3,nbr_pool), ... + points(1:3,neighborhoods(hood)))).^2, 1)); + +% points(1:3,nbr_pool) +% points(1:3,neighborhoods(hood)) +% bsxfun(@minus, points(1:3,nbr_pool), points(1:3,neighborhoods(hood))) +% (bsxfun(@minus, points(1:3,nbr_pool), points(1:3,neighborhoods(hood)))).^2 +% sum((bsxfun(@minus, points(1:3,nbr_pool), points(1:3,neighborhoods(hood)))).^2,1) +% [dsq,new_idx]=min(sum((bsxfun(@minus, points(1:3,nbr_pool), points(1:3,neighborhoods(hood)))).^2,1)) + + % Update the neighborhood +% disp('updating neighborhood...'); + + new_dist = sqrt(dsq); + +% new_idx +% new_dist +% nbr_pool(new_idx) +% neighborhoods(hood) +% nbrs{neighborhoods(hood)} +% dists{neighborhoods(hood)} + + % Maintain a partially sorted list of indices by max distance. + if (new_dist > dists{neighborhoods(hood)}(end)) + dists{neighborhoods(hood)}(nbrs{neighborhoods(hood)}==idx) = ... + dists{neighborhoods(hood)}(end); + + nbrs{neighborhoods(hood)}(nbrs{neighborhoods(hood)}==idx) = ... + nbrs{neighborhoods(hood)}(end); + + dists{neighborhoods(hood)}(end) = sqrt(dsq); + + nbrs{neighborhoods(hood)}(end) = nbr_pool(new_idx); + else + dists{neighborhoods(hood)}(nbrs{neighborhoods(hood)}==idx) = ... + sqrt(dsq); + + nbrs{neighborhoods(hood)}(nbrs{neighborhoods(hood)}==idx) = ... + nbr_pool(new_idx); + end + +% nbrs{neighborhoods(hood)} +% dists{neighborhoods(hood)} + + % Update the dual +% disp('updating dual...'); + + +% duals{nbr_pool(new_idx)} + duals{nbr_pool(new_idx)}(end + 1,1) = neighborhoods(hood); +% duals{nbr_pool(new_idx)} +% disp(nbrs{neighborhoods(hood)}); +% disp(duals{nbr_pool(new_idx)}); + +% error('stop'); + + end + + nbrs{idx}=[]; + duals{idx}=[]; + test{idx}=[]; + + % Update analysis for elements +% disp('updating analysis...'); + + if (isempty(neighborhoods)) + warning('id1','empty neighborhoods'); + end + + if (isempty(neighbors)) + warning('id2','empty neighbors'); + end + +% if (~isempty(neighborhoods) && ~isempty(neighbors)) + neighborhoods=fastsetunion({neighborhoods, neighbors}); +% elseif (isempty(neighborhoods) && ~isempty(neighbors)) +% neighborhoods=neighbors; +% end + +% nbrs{neighborhoods} +% test{neighborhoods} + +% error('Stop'); + + if (~isempty(neighborhoods)) + locfeats = nbrspanalyze(points, ... + nbrs(neighborhoods)); +% locfeats = nbrspanalyze(points, ... +% nbrs(neighborhoods), ... +% test(neighborhoods)); + + locsigs=zeros(1,length(neighborhoods)); + + % Hausdorff DE is no-go + for hood=1:length(neighborhoods) +% rss=sqrt(sum(feats.features(:,nbrs{neighborhoods(hood)}).^2./length(nbrs{neighborhoods(hood)}), 2)); + +% nbrdims(1:dimensions-1)=(rss(1:dimensions-1)-rss(2:dimensions))./... +% repmat(rss(1),dimensions-1,1); +% nbrdims(dimensions)=rss(dimensions)./rss(1); + +% de=-sum(nbrdims.*log(nbrdims))./log(dimensions); + +% feats.de(test{elem}) +% feats.de(test{elem})-de +% abs(feats.de(test{elem})-de) + +% locsigs(hood)=max(abs(feats.de(test{neighborhoods(hood)})-de)); + +% hood +% neighborhoods(hood) +% sigs(neighborhoods(hood)) +% locfeats.entropy(hood) +% feats.entropy(test{neighborhoods(hood)}) +% natives.entropy(test{neighborhoods(hood)}) +% natives.entropy(test{neighborhoods(hood)})-locfeats.entropy(hood) +% abs(natives.entropy(test{neighborhoods(hood)})-locfeats.entropy(hood)) +% max(abs(natives.entropy(test{neighborhoods(hood)})-locfeats.entropy(hood))) + + + locsigs(hood)=max(abs(natives.entropy(test{neighborhoods(hood)})-locfeats.entropy(hood))); + +% error('Stop'); + +% A=abs(bsxfun(@minus, feats.de(test{neighborhoods(hood)}).', feats.de(nbrs{neighborhoods(hood)}))); +% locsigs(hood)=max(max(min(A,[],1)),max(min(A,[],2))); + end + +% radii(neighborhoods)=cellfun(@max,dists(neighborhoods)); +% radii(neighborhoods)=max(radii(neighborhoods), tdist(neighborhoods)); +% biases(neighborhoods)=(locfeats.dists2means)./(radii(neighborhoods)); + +% locsigs = locfeats.de; +% locsigs = sqrt(biases(neighborhoods).*locfeats.de); +% locsigs = 0.5 * (biases(neighborhoods) + locfeats.de); +% locsigs = max(biases(neighborhoods), locfeats.de); +% locsigs = min(biases(neighborhoods), locfeats.de); +% locsigs = biases(neighborhoods); + +% locsigs +% sigs(neighborhoods) +% error('Stop'); + + % Update heap, where required +% updates=neighborhoods(ismembc(neighborhoods,parts{dim})).'; +% sig_upd=locsigs(ismembc(neighborhoods,updates)); + updates=neighborhoods'; + sig_upd=locsigs; + + heap.erase(sigs(updates),updates); + heap.insert(sig_upd,updates); + sigs(neighborhoods)=locsigs; +% feats.de(neighborhoods)=locfeats.de; +% feats.dists2means(neighborhoods)=locfeats.dists2means; +% feats.biases(neighborhoods)=locfeats.biases; + end + + current = current + 1; + + % Update the time bar. + if (toc > 1) + tic; + h = timebar(x, n, msg, tstart, h); + end + end + + % Close the time bar, if still open. + if (all(ishghandle(h, 'figure'))) + close(h); + end + + disp(['...done in ' num2str(toc(tstart)) 's']); +% end + + disp(['Generated fallback neighborhood pool ' num2str(fellback) ' time(s).']); + +% x_max=max(cellfun(@length, removed_wgts)); +% y=nan(length(removed_wgts),x_max); + +% for dim=1:dimensions +% y(dim,1:length(removed_wgts{dim}))=removed_wgts{dim}; +% end + + x_max=length(removed_wgts); + y = removed_wgts; + + figure, plot(1:x_max,y); + figure, plot(1:length(dual_size),dual_size); + +% results = points(:,~ismember(1:size(points,2),removed)); + results = removed; + +end + +%****************************************************************************** +% parseInputs +%****************************************************************************** +function [userParams] = parseInputs(varargin) + % Support function to parse inputs into userParams structure + % + % Parameters: + % varargin - Variable-length input argument list. Option strings may be + % provided as abbreviations as long as they resolve to a + % unique option.Defined key-value pairs include the following: + % 'neighbors' - . + % 'radius' - . + % 'counts' - . + % 'steps' - . + % + % Returns: + % userParams - . + + userParams = struct('radius', 0, 'neighbors', 0, 'counts', [], 'steps', []); + + if isempty(varargin) + error('Damkjer:InvalidOptions', ... + ['A neighborhood size must be specified, either directly or '... + 'via optional parameters.']); + end + + if length(varargin) == 1 || ~isnumeric(varargin{2}) + value = varargin{1}; + + if (isscalar(value) && ... + isreal(value) && ... + value >= 5) + userParams.neighbors = fix(value); + else + error('Damkjer:InvalidCount', ... + ['Number of Neighbors must be a real valued positive '... + 'integer greater or equal to 5: ' num2str(value)]); + end + + varargin(1) = []; + end + + % Parse the Property/Value pairs + if rem(length(varargin), 2) ~= 0 + error('Damkjer:PropertyValueNotPair', ... + ['Additional arguments must take the form of Property/Value '... + 'pairs']); + end + + propertyNames = {'neighbors', 'radius', 'counts', 'steps'}; + + while ~isempty(varargin) + property = varargin{1}; + value = varargin{2}; + + % If the property has been supplied in a shortened form, lengthen it + iProperty = find(strncmpi(property, propertyNames, length(property))); + + if isempty(iProperty) + error('Damkjer:InvalidProperty', 'Invalid Property'); + elseif length(iProperty) > 1 + error('Damkjer:AmbiguousProperty', ... + 'Supplied shortened property name is ambiguous'); + end + + property = propertyNames{iProperty}; + + switch property + case 'neighbors' + if (isscalar(value) && ... + isreal(value) && ... + value >= 5) + userParams.neighbors = fix(value); + else + error('Damkjer:InvalidCount', ... + ['Number of Neighbors must be a real valued positive '... + 'integer greater or equal to 5']); + end + case 'radius' + if (isscalar(value) && ... + isnumeric(value) && ... + isreal(value) && ... + value > 0) + userParams.radius = value; + else + error('Damkjer:InvalidRadius', ... + 'Radius must be a real valued positive scalar'); + end + case 'counts' + if (isvector(value) && ... + isreal(value) && ... + issorted(value) && ... + all(value >= 5)) + userParams.counts = fix(value); + else + error('Damkjer:InvalidCount', ... + ['Counts must be a sorted vector of real valued positive '... + 'integers greater or equal to 5']); + end + case 'steps' + if (isvector(value) && ... + isreal(value) && ... + issorted(value) && ... + all(value > 0)) + userParams.steps = value; + else + error('Damkjer:InvalidSteps', ... + ['Steps must be a sorted vector of real valued positive '... + 'values']); + end + end + + varargin(1:2) = []; + end + + % Check for mutually exclusive options. + if (~isempty(userParams.counts) && userParams.neighbors >= 5) + error('Damkjer:MutExOpts', ... + '''neighbors'' and ''counts'' options are mutually exclusive'); + end + + if (~isempty(userParams.steps) && userParams.radius > 0) + error('Damkjer:MutExOpts', ... + '''steps'' and ''radius'' options are mutually exclusive'); + end + + if (~isempty(userParams.counts) && ~isempty(userParams.steps)) + error('Damkjer:MutExOpts', ... + '''steps'' and ''counts'' options are mutually exclusive'); + end + + % Default, if necessary. + if (userParams.neighbors <= 0 && ... + userParams.radius <= 0 && ... + isempty(userParams.counts) && ... + isempty(userParams.steps)) + userParams.radius = 1; + end +end Index: Damkjer/PointProcessing/Thinning_Dyn/thin_Dyn.m =================================================================== --- Damkjer/PointProcessing/Thinning_Dyn/thin_Dyn.m (revision 13) +++ Damkjer/PointProcessing/Thinning_Dyn/thin_Dyn.m (revision 14) @@ -75,47 +75,44 @@ userParams = parseInputs(varargin{:}); % Parse Inputs - % Make sure we have unique points... we barf otherwise... + %*** + % CONFIRM UNIQUENESS + %*** + + % Make sure that all points in the set are unique... disp('Culling duplicate points...'); tstart=tic; points = unique(points.','rows').'; disp(['...done in ' num2str(toc(tstart)) 's']); - + + %*** + % INDEX POINT SET + %*** disp('Indexing points...'); + tstart=tic; database = VpTree(points); queries = points; disp(['...done in ' num2str(toc(tstart)) 's']); - database.excludeWithin(0.001); % Set mm precision on searches. + %*** + % CREATE NEIGHBORHOODS + %*** + + database.excludeWithin(0.00001); % Set um precision on searches. % Sizing metrics dimensions=size(queries, 1); elements=size(queries, 2); - - % Pre-allocate classes -% norms=zeros(size(queries)); -% feats=zeros(size(queries)); -% biases=zeros(1,elements); -% ints=zeros(1,elements); - - + nbrs = cell(1,elements); - duals = cell(1,elements); - tests = cell(1,elements); - - sigs=zeros(1,elements); - - for elem=1:elements - tests{elem}=elem; - end - - %*** - % COMPUTE NEIGHBORHOODS - %*** + dists = cell(1,elements); + test = num2cell(1:elements, 1); % This process may take a while. Display time bar while processing. - msg='Computing Neigborhoods...'; + msg='Creating Neigborhoods...'; tstart=tic; h = timebar(1, elements, msg, tstart); - + + disp(msg); + % Step through the queries in chunks. step=1000; @@ -131,28 +128,34 @@ if (userParams.radius > 0 && userParams.neighbors <= 0) % Perform a fixed radius search - nbrs(elem:last)=database.rnn(queries(:,elem:last),... - userParams.radius); + [nbrs(elem:last), ... + dists(elem:last)] = database.rnn(queries(:,elem:last), ... + userParams.radius); elseif (userParams.radius <= 0 && userParams.neighbors > 0) % Search unconstrained neighbors - nbrs(elem:last)=database.knn(queries(:,elem:last),... - userParams.neighbors); + [nbrs(elem:last), ... + dists(elem:last)] = database.knn(queries(:,elem:last), ... + userParams.neighbors); elseif (userParams.radius > 0 && userParams.neighbors > 0) % Search constrained to radius - nbrs(elem:last)=database.knn(queries(:,elem:last),... - userParams.neighbors,... - 'lim',userParams.radius); + [nbrs(elem:last), ... + dists(elem:last)] = database.knn(queries(:,elem:last), ... + userParams.neighbors, ... + 'lim', userParams.radius); elseif (~isempty(userParams.counts) && userParams.radius <= 0) % Search unconstrained neighbors - nbrs(elem:last)=database.knn(queries(:,elem:last),... - max(userParams.counts)); + [nbrs(elem:last), ... + dists(elem:last)] = database.knn(queries(:,elem:last), ... + max(userParams.counts)); elseif (~isempty(userParams.counts) && userParams.radius > 0) % Search constrained to radius - nbrs(elem:last)=database.knn(queries(:,elem:last),... - max(userParams.counts),... - 'lim',userParams.radius); + [nbrs(elem:last), ... + dists(elem:last)] = database.knn(queries(:,elem:last), ... + max(userParams.counts), ... + 'lim', userParams.radius); elseif (~isempty(userParams.steps)) % Perform a fixed radius search - [nbrs(elem:last),DISTS]=database.rnn(queries(:,elem:last),... - max(userParams.steps)); + [nbrs(elem:last), ... + dists(elem:last)] = database.rnn(queries(:,elem:last), ... + max(userParams.steps)); end @@ -169,9 +172,80 @@ end - %*** - % COMPUTE SIGNIFICANCES (THIS IS SLOW AND SHOULD BE EASY TO PARALLELIZE) - %*** - - % This process may take a while. Display time bar while processing. + disp(['...done in ' num2str(toc(tstart)) 's']); + + min(cellfun(@min, dists)) + +% error('Stop'); + +%% Obsolete code +% %*** +% % COMPUTE DUAL NEIGHBORHOODS (THIS IS SLOW. LET'S MEX IT) +% %*** +% +% % This process may take a while. Display time bar while processing. +% msg='Computing Dual Neigborhoods...'; +% tstart=tic; +% h = timebar(1, elements, msg, tstart); +% +% disp(msg); +% +% tic; +% for elem=1:elements +% +% % Build up dual set +% for nbr=1:length(nbrs{elem}) +% n = nbrs{elem}(nbr); +% duals{n}=[duals{n}(:); elem]; +% end +% +% % Update the time bar. +% if (toc > 1) +% tic; +% h = timebar(elem, elements, msg, tstart, h); +% end +% end +% +% % Close the time bar, if still open. +% if (all(ishghandle(h, 'figure'))) +% close(h); +% end +% +% disp(['...done in ' num2str(toc(tstart)) 's']); +%% + + %*** + % COMPUTE DUAL NEIGHBORHOODS + %*** + + disp('Computing Dual Neigborhoods...'); + tstart=tic; + duals=fastsetdual(nbrs); + disp(['...done in ' num2str(toc(tstart)) 's']); + + %*** + % COMPUTE SIGNIFICANCES + %*** + disp('Computing Initial Significances...'); +% tstart=tic; +% +% % feats = nbrspanalyze(points, nbrs, 'timebar'); +% +% feats = nbrspanalyze(points, nbrs, test, 'timebar'); +% +% disp(['...done in ' num2str(toc(tstart)) 's']); +% +% % error('Stop'); +% +% radii=cellfun(@max,dists); +% biases=feats.dists2means./radii; +% tdist = zeros(1,elements); +% +% % sigs = feats.de; +% % sigs = sqrt(biases.*feats.de); +% % sigs = 0.5 * (biases + feats.de); +% sigs = max(biases, feats.de); +% % sigs = min(biases, feats.de); +% % sigs = biases; + msg='Computing Significances...'; tstart=tic; @@ -183,4 +257,6 @@ tic; + sigs=zeros(1,elements); + for elem=1:step:elements @@ -189,5 +265,5 @@ % Compute significance - sigs(elem:last)=sig_Dyn(points,tests(elem:last),nbrs(elem:last)); + sigs(elem:last)=sig_Dyn(points,test(elem:last),nbrs(elem:last)); % Update the time bar. @@ -202,47 +278,261 @@ close(h); end - - results = [points;sigs]; - -% %*** -% % COMPUTE DUAL NEIGHBORHOODS (THIS IS SLOW AND SHOULD BE EASY TO PARALLELIZE!) -% %*** -% -% % This process may take a while. Display time bar while processing. -% msg='Computing Dual Neigborhoods...'; -% tstart=tic; -% h = timebar(1, elements, msg, tstart); -% -% tic; -% -% for elem=1:elements -% % Build up dual set -% for nbr=1:length(nbrs{elem}) -% n = nbrs{elem}(nbr); -% duals{n}=[duals{n}(:); elem]; -% end -% -% % if rem(elem-1,100000)==0 -% % for insp=1:10 -% % disp(insp); -% % disp(duals{insp}'); -% % end -% % end -% -% % end -% -% % Update the time bar. -% if (toc > 1) -% tic; -% h = timebar(elem, elements, msg, tstart, h); -% end -% end -% -% % Close the time bar, if still open. -% if (all(ishghandle(h, 'figure'))) -% close(h); -% end -% -% results = points; + + results=sigs; + +end + +function cont +% error('Stop') + + %*** + % PARTITION THE DATA SET + %*** + +% disp('Partition by Label...'); +% tstart=tic; +% +% parts = cell(1,dimensions); +% +% for dim=1:dimensions +% parts{dim}=sort(uint64(find(feats.labeling==dim)),'ascend'); +% end +% +% [~,idx]=sort(cellfun('length',parts),'descend'); +% parts=parts(idx); +% +% disp(['...done in ' num2str(toc(tstart)) 's']); + + %*** + % THIN BY LABEL + %*** + +% fraction=0; + + removed=zeros(1,elements-threshold); + dual_size=zeros(1,elements-threshold); + current=1; + +% removed_wgts=cell(1,3); + fellback = 0; + +% for dim=1:dimensions +% desired = threshold * length(parts{dim}) / elements + fraction; +% tau = floor(desired); +% fraction = desired-tau; +% heap = SplayTree(sigs(parts{dim}),parts{dim}); +% n = length(parts{dim})-tau; + +% removed_wgts{dim}=zeros(1,n); + + tau=threshold; + heap = SplayTree(sigs,1:length(sigs)); + n = length(sigs) - tau; + removed_wgts=zeros(1,n); + +% msg=['Thin Dimension ' num2str(dim) '(' num2str(length(parts{dim})) ... +% '->' num2str(tau) ')...']; + msg='Thin...'; + tstart=tic; + h = timebar(1, n, msg, tstart); + + disp(msg); + tic; + + for x=1:n + + % Locate the least significant point to remove + [wgt, idx]=heap.pop_head(); + + % Mark the point as removed + removed(current)=idx; +% removed_wgts{dim}(x)=wgt; + removed_wgts(x)=wgt; + + + % Distribute the points in the test set to the neighbor's test sets. + neighbors = nbrs{idx}; + samples = test{idx}; + + for samp=1:length(samples) + % Find the closest + [~,new_idx]=min(sum((bsxfun(@minus, points(:,neighbors), ... + points(:,samples(samp)))).^2, 1)); + + % Update the test set + test{neighbors(new_idx)}(end + 1,1) = samples(samp); + +% tdist(neighbors(new_idx)) = max(tdist(neighbors(new_idx)),... +% sqrt(dsq)); + end + + % The removed element neighborhood is eliminated, so it no longer + % contains the neighbors. Update the neighbor duals appropriately. + for nbr=1:length(neighbors) + duals{neighbors(nbr)}=duals{neighbors(nbr)}(duals{neighbors(nbr)}~=idx); + end + + % Find neighborhoods that contain removed element + neighborhoods=duals{idx}; + dual_size(current)=length(neighborhoods); + neighborhoods=sort(neighborhoods(neighborhoods~=idx), 'ascend'); + + % The fallback neighborhood is identical for all neighborhoods, but + % may be expensive to compute. Allow it to be lazily instantiated. + fallback_pool = []; + + % Replace removed element in each neighborhood with nearest neighbor's + % neighbor + for hood=1:length(neighborhoods) + + % Build up the list of potential new neighbors +% disp('building candidates...'); + +%% Obsolete code +% nbr_pool=unique(vertcat(nbrs{nbrs{neighborhoods(hood)}})); +% nbr_pool=setdiff(nbr_pool, idx); +% nbr_pool=setdiff(nbr_pool, nbrs{neighborhoods(hood)}); +%% + nbr_pool=fastsetunion({nbrs{nbrs{neighborhoods(hood)}}}); + non_nbrs=fastsetunion({nbrs{neighborhoods(hood)},idx,neighborhoods(hood)}); + nbr_pool=nbr_pool(~ismembc(nbr_pool(:), non_nbrs(:))); + + if (isempty(nbr_pool)) + if (isempty(fallback_pool)) + fallback_pool = uint64(find(~ismember(1:size(points,2),removed))); + end + fellback = fellback + 1; + nbr_pool = fallback_pool; + nbr_pool=nbr_pool(~ismembc(nbr_pool(:), non_nbrs(:))); + end + + % Find the closest +% disp('finding closest...'); + +%% Obsolete code +% % Don't use VP Tree for single nearest neighbor... overhead not +% % worth it. +% locdb = VpTree(points(:,nbr_pool)); +% nn=locdb.knn(points(:,idx),1); +% new_idx=nn{1} +%% + [dsq,new_idx]=min(sum((bsxfun(@minus, points(:,nbr_pool), ... + points(:,neighborhoods(hood)))).^2, 1)); + + % Update the neighborhood +% disp('updating neighborhood...'); + + new_dist = sqrt(dsq); + + % Maintain a partially sorted list of indices by max distance. + if (new_dist > dists{neighborhoods(hood)}(end)) + dists{neighborhoods(hood)}(nbrs{neighborhoods(hood)}==idx) = ... + dists{neighborhoods(hood)}(end); + + nbrs{neighborhoods(hood)}(nbrs{neighborhoods(hood)}==idx) = ... + nbrs{neighborhoods(hood)}(end); + + dists{neighborhoods(hood)}(end) = sqrt(dsq); + + nbrs{neighborhoods(hood)}(end) = nbr_pool(new_idx); + else + dists{neighborhoods(hood)}(nbrs{neighborhoods(hood)}==idx) = ... + sqrt(dsq); + + nbrs{neighborhoods(hood)}(nbrs{neighborhoods(hood)}==idx) = ... + nbr_pool(new_idx); + end + + % Update the dual +% disp('updating dual...'); + duals{nbr_pool(new_idx)}(end + 1,1) = neighborhoods(hood); + +% disp(nbrs{neighborhoods(hood)}); +% disp(duals{nbr_pool(new_idx)}); + end + + nbrs{idx}=[]; + duals{idx}=[]; + test{idx}=[]; + + % Update analysis for elements +% disp('updating analysis...'); + + neighborhoods=fastsetunion({neighborhoods, neighbors}); + +% nbrs{neighborhoods} +% test{neighborhoods} + +% error('Stop'); + + if (~isempty(neighborhoods)) + + locsigs=sig_Dyn(points,test(neighborhoods),nbrs(neighborhoods)); + + +% locfeats = nbrspanalyze(points, ... +% nbrs(neighborhoods), ... +% test(neighborhoods)); + +% radii(neighborhoods)=cellfun(@max,dists(neighborhoods)); +% radii(neighborhoods)=max(radii(neighborhoods), tdist(neighborhoods)); +% biases(neighborhoods)=(locfeats.dists2means)./(radii(neighborhoods)); + +% locsigs = locfeats.de; +% locsigs = sqrt(biases(neighborhoods).*locfeats.de); +% locsigs = 0.5 * (biases(neighborhoods) + locfeats.de); +% locsigs = max(biases(neighborhoods), locfeats.de); +% locsigs = min(biases(neighborhoods), locfeats.de); +% locsigs = biases(neighborhoods); + + % Update heap, where required +% updates=neighborhoods(ismembc(neighborhoods,parts{dim})).'; +% sig_upd=locsigs(ismembc(neighborhoods,updates)); + updates=neighborhoods'; + sig_upd=locsigs; + + heap.erase(sigs(updates),updates); + heap.insert(sig_upd,updates); + sigs(neighborhoods)=locsigs; +% feats.de(neighborhoods)=locfeats.de; +% feats.dists2means(neighborhoods)=locfeats.dists2means; +% feats.biases(neighborhoods)=locfeats.biases; + end + + current = current + 1; + + % Update the time bar. + if (toc > 1) + tic; + h = timebar(x, n, msg, tstart, h); + end + end + + % Close the time bar, if still open. + if (all(ishghandle(h, 'figure'))) + close(h); + end + + disp(['...done in ' num2str(toc(tstart)) 's']); +% end + + disp(['Generated fallback neighborhood pool ' num2str(fellback) ' time(s).']); + +% x_max=max(cellfun(@length, removed_wgts)); +% y=nan(length(removed_wgts),x_max); + +% for dim=1:dimensions +% y(dim,1:length(removed_wgts{dim}))=removed_wgts{dim}; +% end + + x_max=length(removed_wgts); + y = removed_wgts; + + figure, plot(1:x_max,y); + figure, plot(1:length(dual_size),dual_size); + +% results = points(:,~ismember(1:size(points,2),removed)); + results = removed; + end Index: Damkjer/PointProcessing/Thinning_Dyn/thin_Dyn_old.m =================================================================== --- Damkjer/PointProcessing/Thinning_Dyn/thin_Dyn_old.m (revision 14) +++ Damkjer/PointProcessing/Thinning_Dyn/thin_Dyn_old.m (revision 14) @@ -0,0 +1,387 @@ +% Thin_Dyn Point Cloud Thinning (Dyn et al) +% +% File: thin_Dyn.m +% +% Description: +% Perform thinning according to Dyn 2004. +% +% Limitations: +% No known limitations. +% +% Synopsis: +% [results] = spanalyse(points, threshold, ...) +% +% Inputs: +% points - . +% threshold - must be positive. between 0-1, treated as percentage, > 1 +% points. +% +% Option strings may be provided as abbreviations as long as they resolve to +% a unique option. +% +% 'neighbors' - . +% 'radius' - . +% 'counts' - . +% 'steps' - . +% +% Outputs: +% results - Thinned point cloud. +% +% Toolbox requirements: +% None. +% +% Code requirements: +% None. +% +% Data requirements: +% None. +% +% References: +% Dyn reference. +% +% See Also: +% None. +% + +% Copyright (C) 2013 Kristian L. Damkjer. +% +% Software History: +% 2013-DEC-28 Initial coding. +% + +%****************************************************************************** +% thin_Dyn +%****************************************************************************** +function [ results ] = thin_Dyn( points, threshold, varargin ) + % Perform local feature attribution via eigenspace analysis. + % + % Parameters: + % points - . + % + % threshold - . + % + % varargin - Variable-length input argument list. Option strings may be + % provided as abbreviations as long as they resolve to a + % unique option.Defined key-value pairs include the following: + % 'neighbors' - . + % 'radius' - . + % 'counts' - . + % 'steps' - . + % + % Returns: + % classes - Structure containing the spatial analysis classifications. + % + + userParams = parseInputs(varargin{:}); % Parse Inputs + + % Make sure we have unique points... we barf otherwise... + disp('Culling duplicate points...'); + tstart=tic; + points = unique(points.','rows').'; + disp(['...done in ' num2str(toc(tstart)) 's']); + + disp('Indexing points...'); + database = VpTree(points); + queries = points; + disp(['...done in ' num2str(toc(tstart)) 's']); + + database.excludeWithin(0.001); % Set mm precision on searches. + + % Sizing metrics + dimensions=size(queries, 1); + elements=size(queries, 2); + + % Pre-allocate classes +% norms=zeros(size(queries)); +% feats=zeros(size(queries)); +% biases=zeros(1,elements); +% ints=zeros(1,elements); + + + nbrs = cell(1,elements); + duals = cell(1,elements); + tests = cell(1,elements); + + sigs=zeros(1,elements); + + for elem=1:elements + tests{elem}=elem; + end + + %*** + % COMPUTE NEIGHBORHOODS + %*** + + % This process may take a while. Display time bar while processing. + msg='Computing Neigborhoods...'; + tstart=tic; + h = timebar(1, elements, msg, tstart); + + % Step through the queries in chunks. + step=1000; + + tic; + + for elem=1:step:elements + + % The last available element may be closer than elem + step. + last=min(elem+step-1,elements); + + % Get the nearest neighbors of elem + if (userParams.radius > 0 && userParams.neighbors <= 0) + % Perform a fixed radius search + nbrs(elem:last)=database.rnn(queries(:,elem:last),... + userParams.radius); + elseif (userParams.radius <= 0 && userParams.neighbors > 0) + % Search unconstrained neighbors + nbrs(elem:last)=database.knn(queries(:,elem:last),... + userParams.neighbors); + elseif (userParams.radius > 0 && userParams.neighbors > 0) + % Search constrained to radius + nbrs(elem:last)=database.knn(queries(:,elem:last),... + userParams.neighbors,... + 'lim',userParams.radius); + elseif (~isempty(userParams.counts) && userParams.radius <= 0) + % Search unconstrained neighbors + nbrs(elem:last)=database.knn(queries(:,elem:last),... + max(userParams.counts)); + elseif (~isempty(userParams.counts) && userParams.radius > 0) + % Search constrained to radius + nbrs(elem:last)=database.knn(queries(:,elem:last),... + max(userParams.counts),... + 'lim',userParams.radius); + elseif (~isempty(userParams.steps)) + % Perform a fixed radius search + [nbrs(elem:last),DISTS]=database.rnn(queries(:,elem:last),... + max(userParams.steps)); + end + + % Update the time bar. + if (toc > 1) + tic; + h = timebar(elem, elements, msg, tstart, h); + end + end + + % Close the time bar, if still open. + if (all(ishghandle(h, 'figure'))) + close(h); + end + + %*** + % COMPUTE SIGNIFICANCES (THIS IS SLOW AND SHOULD BE EASY TO PARALLELIZE) + %*** + + % This process may take a while. Display time bar while processing. + msg='Computing Significances...'; + tstart=tic; + h = timebar(1, elements, msg, tstart); + + % Step through the queries in chunks. + step=1000; + + tic; + + for elem=1:step:elements + + % The last available element may be closer than elem + step. + last=min(elem+step-1,elements); + + % Compute significance + sigs(elem:last)=sig_Dyn(points,tests(elem:last),nbrs(elem:last)); + + % Update the time bar. + if (toc > 1) + tic; + h = timebar(elem, elements, msg, tstart, h); + end + end + + % Close the time bar, if still open. + if (all(ishghandle(h, 'figure'))) + close(h); + end + + results = [points;sigs]; + +% %*** +% % COMPUTE DUAL NEIGHBORHOODS (THIS IS SLOW AND SHOULD BE EASY TO PARALLELIZE!) +% %*** +% +% % This process may take a while. Display time bar while processing. +% msg='Computing Dual Neigborhoods...'; +% tstart=tic; +% h = timebar(1, elements, msg, tstart); +% +% tic; +% +% for elem=1:elements +% % Build up dual set +% for nbr=1:length(nbrs{elem}) +% n = nbrs{elem}(nbr); +% duals{n}=[duals{n}(:); elem]; +% end +% +% % if rem(elem-1,100000)==0 +% % for insp=1:10 +% % disp(insp); +% % disp(duals{insp}'); +% % end +% % end +% +% % end +% +% % Update the time bar. +% if (toc > 1) +% tic; +% h = timebar(elem, elements, msg, tstart, h); +% end +% end +% +% % Close the time bar, if still open. +% if (all(ishghandle(h, 'figure'))) +% close(h); +% end +% +% results = points; +end + +%****************************************************************************** +% parseInputs +%****************************************************************************** +function [userParams] = parseInputs(varargin) + % Support function to parse inputs into userParams structure + % + % Parameters: + % varargin - Variable-length input argument list. Option strings may be + % provided as abbreviations as long as they resolve to a + % unique option.Defined key-value pairs include the following: + % 'neighbors' - . + % 'radius' - . + % 'counts' - . + % 'steps' - . + % + % Returns: + % userParams - . + + userParams = struct('radius', 0, 'neighbors', 0, 'counts', [], 'steps', []); + + if isempty(varargin) + error('Damkjer:InvalidOptions', ... + ['A neighborhood size must be specified, either directly or '... + 'via optional parameters.']); + end + + if length(varargin) == 1 || ~isnumeric(varargin{2}) + value = varargin{1}; + + if (isscalar(value) && ... + isreal(value) && ... + value >= 5) + userParams.neighbors = fix(value); + else + error('Damkjer:InvalidCount', ... + ['Number of Neighbors must be a real valued positive '... + 'integer greater or equal to 5: ' num2str(value)]); + end + + varargin(1) = []; + end + + % Parse the Property/Value pairs + if rem(length(varargin), 2) ~= 0 + error('Damkjer:PropertyValueNotPair', ... + ['Additional arguments must take the form of Property/Value '... + 'pairs']); + end + + propertyNames = {'neighbors', 'radius', 'counts', 'steps'}; + + while ~isempty(varargin) + property = varargin{1}; + value = varargin{2}; + + % If the property has been supplied in a shortened form, lengthen it + iProperty = find(strncmpi(property, propertyNames, length(property))); + + if isempty(iProperty) + error('Damkjer:InvalidProperty', 'Invalid Property'); + elseif length(iProperty) > 1 + error('Damkjer:AmbiguousProperty', ... + 'Supplied shortened property name is ambiguous'); + end + + property = propertyNames{iProperty}; + + switch property + case 'neighbors' + if (isscalar(value) && ... + isreal(value) && ... + value >= 5) + userParams.neighbors = fix(value); + else + error('Damkjer:InvalidCount', ... + ['Number of Neighbors must be a real valued positive '... + 'integer greater or equal to 5']); + end + case 'radius' + if (isscalar(value) && ... + isnumeric(value) && ... + isreal(value) && ... + value > 0) + userParams.radius = value; + else + error('Damkjer:InvalidRadius', ... + 'Radius must be a real valued positive scalar'); + end + case 'counts' + if (isvector(value) && ... + isreal(value) && ... + issorted(value) && ... + all(value >= 5)) + userParams.counts = fix(value); + else + error('Damkjer:InvalidCount', ... + ['Counts must be a sorted vector of real valued positive '... + 'integers greater or equal to 5']); + end + case 'steps' + if (isvector(value) && ... + isreal(value) && ... + issorted(value) && ... + all(value > 0)) + userParams.steps = value; + else + error('Damkjer:InvalidSteps', ... + ['Steps must be a sorted vector of real valued positive '... + 'values']); + end + end + + varargin(1:2) = []; + end + + % Check for mutually exclusive options. + if (~isempty(userParams.counts) && userParams.neighbors >= 5) + error('Damkjer:MutExOpts', ... + '''neighbors'' and ''counts'' options are mutually exclusive'); + end + + if (~isempty(userParams.steps) && userParams.radius > 0) + error('Damkjer:MutExOpts', ... + '''steps'' and ''radius'' options are mutually exclusive'); + end + + if (~isempty(userParams.counts) && ~isempty(userParams.steps)) + error('Damkjer:MutExOpts', ... + '''steps'' and ''counts'' options are mutually exclusive'); + end + + % Default, if necessary. + if (userParams.neighbors <= 0 && ... + userParams.radius <= 0 && ... + isempty(userParams.counts) && ... + isempty(userParams.steps)) + userParams.radius = 1; + end +end Index: Damkjer/PointProcessing/sig_Dyn.m =================================================================== --- Damkjer/PointProcessing/sig_Dyn.m (revision 13) +++ Damkjer/PointProcessing/sig_Dyn.m (revision 14) @@ -128,19 +128,25 @@ PhiM = (1-rM).^4.*(4*rM+1); - if rcond(PhiM) < 1e-12 - fmt=get(0, 'format'); - format('hex'); - disp(rM); - disp(Z); - disp(f); - disp(PhiM); - hold on - scatter3(cells{elem}(:,1),cells{elem}(:,2),cells{elem}(:,3), 5, 'blue', 'o'); - scatter3(Z(:,1),Z(:,2),Z(:,3), 5, 'green', 'o'); - scatter3(x(1),x(2),x(3), 5, 'red', 'x'); - hold off - set(0, 'format', fmt); - error('oops'); - end +% if rcond(PhiM) < 1e-12 +% % fmt=get(0, 'format'); +% % format('hex'); +% disp(Z); +% disp(rM); +% disp(f); +% disp(PhiM); +% 1-rM +% (1-rM).^4 +% 4*rM +% 4*rM+1 +% ((1-rM).^4).*(4*rM+1) +% 1-(((1-rM).^4).*(4*rM+1)) +% hold on +% scatter3(cells{elem}(:,1),cells{elem}(:,2),cells{elem}(:,3), 5, 'blue', 'o'); +% scatter3(Z(:,1),Z(:,2),Z(:,3), 5, 'green', 'o'); +% scatter3(x(1),x(2),x(3), 5, 'red', 'x'); +% hold off +% % set(0, 'format', fmt); +% error('oops'); +% end as = PhiM\f; Index: Damkjer/Util/Dictionaries/SplayTree.h =================================================================== --- Damkjer/Util/Dictionaries/SplayTree.h (revision 13) +++ Damkjer/Util/Dictionaries/SplayTree.h (revision 14) @@ -47,4 +47,5 @@ , theRBranch(0) , theElement(element) + , theIndex(0) { } @@ -81,4 +82,5 @@ Node* theRBranch; T theElement; + std::size_t theIndex; }; //> A node interface for a splay tree. @@ -97,5 +99,13 @@ //< + void update(const T&, const std::size_t&); + //> Remove an element from the tree. + //< + void erase(const T&); + //> Remove an element from the tree. + //< + + void eraseIndex(const std::size_t&); //> Remove an element from the tree. //< @@ -148,4 +158,6 @@ //> The total number of elements contained in the tree. //< + + std::deque theIndexedList; }; Index: Damkjer/Util/Dictionaries/SplayTree.hpp =================================================================== --- Damkjer/Util/Dictionaries/SplayTree.hpp (revision 13) +++ Damkjer/Util/Dictionaries/SplayTree.hpp (revision 14) @@ -37,4 +37,5 @@ : theRoot(0) , theSize(0) + , theIndexedList(1) { } @@ -62,4 +63,6 @@ { Node* node = new Node(element); + node->theIndex = theIndexedList.size(); + theIndexedList.push_back(node); if (!theRoot) @@ -134,4 +137,51 @@ //***************************************************************************** +// SplayTree::update(const T&) +//> Add a new element to the tree. +//< +//***************************************************************************** +template +void +SplayTree::update(const T& element, const std::size_t& index) +{ + // This can probably be made more efficient, but this will suffice for now. + Node* found = theIndexedList[index]; + + if (!found) return; // What? Who's mismanaging the structure? + + eraseIndex(index); + + Node* node = new Node(element); + node->theIndex = index; + theIndexedList[index] = node; + + //if (!theRoot) // How can we be found, but not have a root? + //{ + // theRoot = node; + // ++theSize; + // return; + //} + + splay(element); + + if (theOrderIs(element, theRoot->theElement)) + { + node->theLBranch = theRoot->theLBranch; + node->theRBranch = theRoot; + theRoot->theLBranch = 0; + theRoot = node; + ++theSize; + } + else // if (theOrderIs(theRoot->theElement, element)) + { + node->theRBranch = theRoot->theRBranch; + node->theLBranch = theRoot; + theRoot->theRBranch = 0; + theRoot = node; + ++theSize; + } +} + +//***************************************************************************** // SplayTree::erase(const T&) //> Remove an element from the tree. @@ -161,4 +211,5 @@ found->theLBranch = 0; found->theRBranch = 0; + theIndexedList[found->theIndex] = 0; delete found; @@ -205,4 +256,38 @@ //***************************************************************************** +// SplayTree::erase_index(const std::size_t&) +//> Remove an element from the tree. +//< +//***************************************************************************** +template +void +SplayTree::eraseIndex(const std::size_t& index) +{ + Node* found = theIndexedList.at(index); + + if (!found) return; + + splay(found->theElement); + + if (!found->theLBranch) + { + theRoot = found->theRBranch; + } + else + { + theRoot = found->theLBranch; + splay(found->theElement); + theRoot->theRBranch = found->theRBranch; + } + + --theSize; + + found->theLBranch = 0; + found->theRBranch = 0; + theIndexedList[index] = 0; + delete found; +} + +//***************************************************************************** // SplayTree::find(const T&) //> Find an element in the tree. Index: Damkjer/Util/MATLAB/SplayTree.m =================================================================== --- Damkjer/Util/MATLAB/SplayTree.m (revision 13) +++ Damkjer/Util/MATLAB/SplayTree.m (revision 14) @@ -131,4 +131,19 @@ %************************************************************************ + % VpTree/update + %************************************************************************ + function update(this, weights, elems) + % Perform a k-nearest neighbor search on the database with the + % set of queries. + % + % Parameters: + % + % Returns: + % + + SplayTreeAPI('update', this.theTree, weights, elems); + end + + %************************************************************************ % VpTree/erase %************************************************************************ @@ -144,4 +159,34 @@ SplayTreeAPI('erase', this.theTree, weights, elems); end + + %************************************************************************ + % VpTree/erase_index + %************************************************************************ + function erase_index(this, elems) + % Perform a k-nearest neighbor search on the database with the + % set of queries. + % + % Parameters: + % + % Returns: + % + + SplayTreeAPI('erase_index', this.theTree, elems); + end + + %************************************************************************ + % VpTree/size + %************************************************************************ + function [tree_size] = size(this) + % Perform a k-nearest neighbor search on the database with the + % set of queries. + % + % Parameters: + % + % Returns: + % + + [tree_size] = SplayTreeAPI('size', this.theTree); + end end end Index: Damkjer/Util/MATLAB/SplayTreeAPI/SplayTreeAPI.cpp =================================================================== --- Damkjer/Util/MATLAB/SplayTreeAPI/SplayTreeAPI.cpp (revision 13) +++ Damkjer/Util/MATLAB/SplayTreeAPI/SplayTreeAPI.cpp (revision 14) @@ -167,5 +167,6 @@ TreeT* returnTree = new TreeT(); - for (mwSize elem = iElems; elem --> 0;) +// for (mwSize elem = iElems; elem --> 0;) + for (mwSize elem = 0; elem < iElems; ++elem) { returnTree->insert((WeightedPointT)std::make_pair(wData[elem], @@ -307,11 +308,12 @@ index[0] = head.second; - tree.erase(head); - - if (tree.find(head) != 0) - { - mexErrMsgIdAndTxt("Damkjer:splayTreePopHead:nargout", - "Pop failed."); - } +// tree.erase(head); + tree.eraseIndex(head.second); + +// if (tree.find(head) != 0) +// { +// mexErrMsgIdAndTxt("Damkjer:splayTreePopHead:nargout", +// "Pop failed."); +// } } @@ -419,6 +421,7 @@ } -//***************************************************************************** -// FUNCTION: splayTreeErase + +//***************************************************************************** +// FUNCTION: splayTreeUpdate //> Find minumum element in splay tree. // @@ -428,5 +431,5 @@ //***************************************************************************** void -splayTreeErase(int nrhs, const mxArray** prhs) +splayTreeUpdate(int nrhs, const mxArray** prhs) { //*** @@ -437,6 +440,6 @@ if (nrhs != 4 || !mxIsNumeric(prhs[1])) { - mexErrMsgIdAndTxt("Damkjer:splayTreeErase:nargin", - "The splayTreeErase function requires three inputs."); + mexErrMsgIdAndTxt("Damkjer:splayTreeInsert:nargin", + "The splayTreeInsert function requires three inputs."); } @@ -457,7 +460,13 @@ !mxIsNumeric(indices)) { - mexErrMsgIdAndTxt("Damkjer:splayTreeErase:prhs", - "Index input to splayTreeErase must be a full vector " - "of real-valued data representing point significance."); + std::stringstream msg; + msg << ((mxIsSparse(indices)) ? "" : "not ") << "sparse " + << (mxGetNumberOfDimensions(indices)) << " " + << ((mxIsComplex(indices)) ? "" : "not ") << "complex " + << ((mxIsNumeric(indices)) ? "" : "not ") << "numeric " + << "Index input to splayTreeInsert must be a full vector " + "of real-valued data representing point indices."; + mexErrMsgIdAndTxt("Damkjer:splayTreeInsert:prhs", + msg.str().c_str()); } @@ -468,7 +477,13 @@ !mxIsNumeric(weights)) { - mexErrMsgIdAndTxt("Damkjer:splayTreeErase:prhs", - "Weight input to splayTreeErase must be a full vector " - "of real-valued data representing point significance."); + std::stringstream msg; + msg << ((mxIsSparse(weights)) ? "" : "not ") << "sparse " + << (mxGetNumberOfDimensions(weights)) << " " + << ((mxIsComplex(weights)) ? "" : "not ") << "complex " + << ((mxIsNumeric(weights)) ? "" : "not ") << "numeric " + << "Weight input to splayTreeInsert must be a full vector " + "of real-valued data representing point significance."; + mexErrMsgIdAndTxt("Damkjer:splayTreeInsert:prhs", + msg.str().c_str()); } @@ -480,6 +495,10 @@ if (iElems != wElems) { - mexErrMsgIdAndTxt("Damkjer:splayTreeErase:prhs", - "Must provide exactly one weight for each index."); + std::stringstream msg; + msg << "Index elements : " << mxGetM(indices) << "x" << iElems << "\n" + << "Weight elements: " << mxGetM(weights) << "x" << wElems << "\n" + << "Must provide exactly one weight for each index."; + mexErrMsgIdAndTxt("Damkjer:splayTreeInsert:prhs", + msg.str().c_str()); } @@ -502,5 +521,159 @@ for (mwSize elem = iElems; elem --> 0;) { + tree.update((WeightedPointT)std::make_pair(wData[elem], iData[elem]), + iData[elem]); + } +} + +//***************************************************************************** +// FUNCTION: splayTreeErase +//> Find minumum element in splay tree. +// +// @param nrhs the number of right-hand side parameters. +// @param prhs the array of right-hand side parameters. +//< +//***************************************************************************** +void +splayTreeErase(int nrhs, const mxArray** prhs) +{ + //*** + // Check parameters. Remember that we always have one "extra" input + // parameter to handle function dispatching through the MEX gateway. It + // always occupies the first input parameter position. + //*** + if (nrhs != 4 || !mxIsNumeric(prhs[1])) + { + mexErrMsgIdAndTxt("Damkjer:splayTreeErase:nargin", + "The splayTreeErase function requires three inputs."); + } + + // Retrieve the tree object through the ClassHandle helper method. + TreeT& tree = HandleT::handleReference(prhs[1]); + + //*** + // The two inputs required by splayTreeErase is the collection of indices + // and the collection of weights. + //*** + const mxArray* weights=prhs[2]; + const mxArray* indices=prhs[3]; + + // Check to make sure that weights are real-valued numerics + if (mxIsSparse(indices) || + mxGetNumberOfDimensions(indices) != 2 || + mxIsComplex(indices) || + !mxIsNumeric(indices)) + { + mexErrMsgIdAndTxt("Damkjer:splayTreeErase:prhs", + "Index input to splayTreeErase must be a full vector " + "of real-valued data representing point significance."); + } + + // Check to make sure that weights are real-valued numerics + if (mxIsSparse(weights) || + mxGetNumberOfDimensions(weights) != 2 || + mxIsComplex(weights) || + !mxIsNumeric(weights)) + { + mexErrMsgIdAndTxt("Damkjer:splayTreeErase:prhs", + "Weight input to splayTreeErase must be a full vector " + "of real-valued data representing point significance."); + } + + // Attempt to make the Splay Tree. + const mwSize iElems = mxGetN(indices); + const mwSize wElems = mxGetN(weights); + + // Check to make sure that weights are real-valued numerics + if (iElems != wElems) + { + mexErrMsgIdAndTxt("Damkjer:splayTreeErase:prhs", + "Must provide exactly one weight for each index."); + } + + mwIndex* iData = (mwIndex*)mxGetData(indices); + double* wData = mxGetPr(weights); + + //*** + // Selection of data structure for the point database is not arbitrary. + // + // 1. Point coordinates are passed in as double-precision values. Preserve + // the data fidelity. + // 2. Individual points are often low-dimensional. Keep the coordinates for + // a single point in contiguous memory by using a std::vector container. + // 3. The point database can be quite large. Since VP Tree is an online data + // structure, allow references to points to occupy non-contiguous memory. + // This also speeds construction and destruction of the point database. + //*** +// WeightedPointSetT pointData(elems); + + for (mwSize elem = iElems; elem --> 0;) + { tree.erase((WeightedPointT)std::make_pair(wData[elem], iData[elem])); + } +} + +//***************************************************************************** +// FUNCTION: splayTreeEraseIndex +//> Find minumum element in splay tree. +// +// @param nrhs the number of right-hand side parameters. +// @param prhs the array of right-hand side parameters. +//< +//***************************************************************************** +void +splayTreeEraseIndex(int nrhs, const mxArray** prhs) +{ + //*** + // Check parameters. Remember that we always have one "extra" input + // parameter to handle function dispatching through the MEX gateway. It + // always occupies the first input parameter position. + //*** + if (nrhs != 3 || !mxIsNumeric(prhs[1])) + { + mexErrMsgIdAndTxt("Damkjer:splayTreeErase:nargin", + "The splayTreeErase function requires two inputs."); + } + + // Retrieve the tree object through the ClassHandle helper method. + TreeT& tree = HandleT::handleReference(prhs[1]); + + //*** + // The two inputs required by splayTreeErase is the collection of indices + // and the collection of weights. + //*** + const mxArray* indices=prhs[2]; + + // Check to make sure that weights are real-valued numerics + if (mxIsSparse(indices) || + mxGetNumberOfDimensions(indices) != 2 || + mxIsComplex(indices) || + !mxIsNumeric(indices)) + { + mexErrMsgIdAndTxt("Damkjer:splayTreeErase:prhs", + "Index input to splayTreeErase must be a full vector " + "of real-valued data representing point significance."); + } + + // Attempt to make the Splay Tree. + const mwSize iElems = mxGetN(indices); + + mwIndex* iData = (mwIndex*)mxGetData(indices); + + //*** + // Selection of data structure for the point database is not arbitrary. + // + // 1. Point coordinates are passed in as double-precision values. Preserve + // the data fidelity. + // 2. Individual points are often low-dimensional. Keep the coordinates for + // a single point in contiguous memory by using a std::vector container. + // 3. The point database can be quite large. Since VP Tree is an online data + // structure, allow references to points to occupy non-contiguous memory. + // This also speeds construction and destruction of the point database. + //*** +// WeightedPointSetT pointData(elems); + + for (mwSize elem = iElems; elem --> 0;) + { + tree.eraseIndex(iData[elem]); } } @@ -573,7 +746,25 @@ splayTreeInsert(nrhs, prhs); } + else if (!strcmp("update", operation)) + { + splayTreeUpdate(nrhs, prhs); + } else if (!strcmp("erase", operation)) { splayTreeErase(nrhs, prhs); + } + else if (!strcmp("erase_index", operation)) + { + splayTreeEraseIndex(nrhs, prhs); + } + else if (!strcmp("size", operation)) + { + TreeT& tree = HandleT::handleReference(prhs[1]); + + const std::size_t& size=tree.size(); + + plhs[0]=mxCreateNumericMatrix(1, 1, mxINDEX_CLASS, mxREAL); + mwIndex* return_array = (mwIndex*) mxGetData(plhs[0]); + return_array[0] = size; } else Index: Damkjer/Util/MATLAB/VpTreeAPI/VpTreeAPI.cpp =================================================================== --- Damkjer/Util/MATLAB/VpTreeAPI/VpTreeAPI.cpp (revision 13) +++ Damkjer/Util/MATLAB/VpTreeAPI/VpTreeAPI.cpp (revision 14) @@ -757,5 +757,5 @@ mwSize neighbors = results[q].first.size(); - nbr_dists[q] = mxCreateDoubleMatrix(0, 0, mxREAL); + nbr_dists.push_back(mxCreateDoubleMatrix(0, 0, mxREAL)); mxSetM(nbr_dists[q], neighbors); mxSetN(nbr_dists[q], 1); Index: Damkjer/Util/Math/fastcov.cpp =================================================================== --- Damkjer/Util/Math/fastcov.cpp (revision 12) +++ Damkjer/Util/Math/fastcov.cpp (revision 14) @@ -58,5 +58,5 @@ "Missing or invalid input argument."); } - + if (nlhs > 3) { @@ -65,6 +65,13 @@ } - mwSize cells = mxGetNumberOfElements (prhs[0]); - + mwSize cells = mxGetNumberOfElements(prhs[0]); +// mwSize ctrs = mxGetNumberOfElements(prhs[1]); + +// if (useCenters && cells != ctrs) +// { +// mexErrMsgIdAndTxt("Damkjer:fastcov:varargout", +// "Size of neighborhoods and centers does not agree."); +// } + plhs[0] = mxCreateCellMatrix(cells, 1); @@ -89,6 +96,6 @@ //HACK: Probably a more efficient way to do this - std::vector bias(cells,0); - std::vector bias_vals(cells,0); + std::vector dist2mean(cells,0); + std::vector dist2mean_vals(cells,0); //HACK: Probably a more efficient way to do this @@ -111,6 +118,6 @@ //HACK: Probably a more efficient way to do this - bias[cell] = mxCreateDoubleMatrix(1, 1, mxREAL); - bias_vals[cell] = mxGetPr(bias[cell]); + dist2mean[cell] = mxCreateDoubleMatrix(1, 1, mxREAL); + dist2mean_vals[cell] = mxGetPr(dist2mean[cell]); //HACK: Probably a more efficient way to do this @@ -128,4 +135,5 @@ { std::vector mean(Ns[cellp], 0.); + std::vector skewmean(Ns[cellp], 0.); // Comment out to remove mean calculation when centering to first point. @@ -134,8 +142,12 @@ for (mwSize n = Ns[cellp]; n --> 0;) { - for (mwSize m = Ms[cellp]; m --> 0;) + for (mwSize m = Ms[cellp]; m --> 1;) { mean[n] += vals[cellp][m + Ms[cellp] * n] * w1; } + + skewmean[n] = mean[n] / ((Ms[cellp] - 1) * w1); + + mean[n] += vals[cellp][Ms[cellp] * n] * w1; } // End comment. @@ -173,5 +185,5 @@ if (nlhs > 1) { - // Experiment: Estimate bias, i.e., how "centered" is our mass? + // Experiment: Estimate distance from first item to mean //std::vector diff_loc(Ns[cellp], 0.); @@ -184,5 +196,6 @@ for (mwSize n = Ns[cellp]; n --> 0;) { - temp = vals[cellp][Ms[cellp] * n] - mean[n]; + // We can't assume anything about point ordering + temp = vals[cellp][Ms[cellp] * n] - skewmean[n]; diff_loc += temp * temp; @@ -192,5 +205,6 @@ } - bias_vals[cellp][0] = std::sqrt(diff_loc)/std::sqrt(diff_ext); +// dist2mean_vals[cellp][0] = std::sqrt(diff_loc)/std::sqrt(diff_ext); + dist2mean_vals[cellp][0] = std::sqrt(diff_loc); // End comment. @@ -198,4 +212,5 @@ { // Experiment: Estimate intensity, i.e., how "dense" is our mass? + // Assumes last point is farthest from first. inty_vals[cellp][0] = Ms[cellp]/(4./3.*M_PI*diff_ext*sqrt(diff_ext)); // End comment. @@ -213,5 +228,5 @@ for (int cell = 0; cell < cells; ++cell) { - mxSetCell(plhs[1], cell, bias[cell]); + mxSetCell(plhs[1], cell, dist2mean[cell]); } } Index: Damkjer/Util/SpatialIndexing/VpTree/VpTree_m.cc =================================================================== --- Damkjer/Util/SpatialIndexing/VpTree/VpTree_m.cc (revision 14) +++ Damkjer/Util/SpatialIndexing/VpTree/VpTree_m.cc (revision 14) @@ -0,0 +1,1 @@ +Util\SpatialIndexing\VpTree\VpTree.m