Index: Damkjer/PointProcessing/Thinning_Damkjer/stochastic_thin.m =================================================================== --- Damkjer/PointProcessing/Thinning_Damkjer/stochastic_thin.m (revision 16) +++ Damkjer/PointProcessing/Thinning_Damkjer/stochastic_thin.m (revision 16) @@ -0,0 +1,633 @@ +% stochastic_thin 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 ] = stochastic_thin( 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); + 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); + + % 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(:,elem:last), ... + userParams.radius); + elseif (userParams.radius <= 0 && userParams.neighbors > 0) + % Search unconstrained 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), ... + 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), ... + 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), ... + 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(elem:last)] = 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 + + 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 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; + + %*** + % 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; + + removed=sigs; +% 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 +% [dsq,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)) +% 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; + results = sigs; + +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/Util/FileIO/FXYZ/readFXYZ.m =================================================================== --- Damkjer/Util/FileIO/FXYZ/readFXYZ.m (revision 16) +++ Damkjer/Util/FileIO/FXYZ/readFXYZ.m (revision 16) @@ -0,0 +1,57 @@ +%%READXYZ Read XYZ formatted file. +% [POINTS,FRAMES,PIXELS] = READXYZ(FILENAME) parses the file FILENAME for +% binary point tuples. POINTS is a 3-by-N collection of world coordinate +% points, FRAMES is a 1-by-N collection of IFOV frame indices for the +% associated points, and PIXELS is a 1-by-N collection of IFOV pixel +% indices for the associated points. +% +% FILENAME is a string containing the name of the file to be opened. +% +% PIXELS indices follow the convention that negative indices capture +% information about IFOV. Index -5 is the effective center of the IFOV +% while indices -1 through -4 capture the IFOV corners at the top of the +% range gate. +% +% See also FOPEN. + +% Copyright 2012 Kristian Linn Damkjer +function [X,F,P] = readXYZ(filename) + +% Attempt to open the file for reading. +[fid, msg] = fopen(filename, 'r'); + +% If unsuccessful, let the user know why and return empy tuples. +if (fid < 0) + disp(msg); + X = zeros(3,0); + F = zeros(1,0,'int32'); + P = zeros(1,0,'int32'); + return; +end + +% Otherwise, parse the data. +dirStruct = dir(filename); + + +numPts = (dirStruct.bytes-24)/20; + +tic; + +offsets = fread(fid, 3, 'double=>double'); + +temp = reshape(fread(fid, inf, 'single=>single'),5,numPts); + +size(offsets) +size(temp) + +X = bsxfun(@plus, temp(1:3,:), offsets); +F = typecast(temp(4,:),'int32'); +P = typecast(temp(5,:),'int32'); + +disp(['Time to parse file: ' num2str(toc) 's']); + +fclose(fid); + +return; + +end Index: Damkjer/Util/LPt.cc =================================================================== --- Damkjer/Util/LPt.cc (revision 16) +++ Damkjer/Util/LPt.cc (revision 16) @@ -0,0 +1,11 @@ +#include "LPt.h" + + +LPt::LPt(void) +{ +} + + +LPt::~LPt(void) +{ +} Index: Damkjer/Util/LPt.h =================================================================== --- Damkjer/Util/LPt.h (revision 16) +++ Damkjer/Util/LPt.h (revision 16) @@ -0,0 +1,57 @@ +//**************************************************************************** +// FILE: SplayTree.h +// +// Copyright (C) 2013 Kristian Damkjer. +// +// DESCRIPTION: +//> The interface definition for splay trees. +//< +// +// LIMITATIONS: +//> This class template file follows the template inclusion pattern. This +// header file should be the only file included by clients wishing to +// instantiate a VpTree specialization. +//< +// +// SOFTWARE HISTORY: +// +//> 2013-DEC-02 K. Damkjer +// Initial Coding. +//< +//**************************************************************************** +#ifndef Damkjer_SplayTree_HEADER +#define Damkjer_SplayTree_HEADER + +namespace Damkjer +{ + +//***************************************************************************** +// CLASS: LPt +//> A LiDAR point type. +// +// @tparam MetricT The metric search space. +//< +//***************************************************************************** +template +class LPt +{ +public: + typedef T value_type; + typedef T& reference; + typedef const T& const_reference; + typedef T* iterator; + // const_iterator; + // differe + + LPt(); + + LPt(const LPt&); + + virtual ~LPt(void); + + +}; + +} + +#endif