Index: Damkjer/PointProcessing/Colorize/normcolor.m =================================================================== --- Damkjer/PointProcessing/Colorize/normcolor.m (revision 0) +++ Damkjer/PointProcessing/Colorize/normcolor.m (revision 3) @@ -3,4 +3,9 @@ % % Color points by surface normal estimate orientation. +% +% Software History: +% 2013-JAN-29 Damkjer, K. +% Updated Error Statement to reflect appropriate package. +% function colors=normcolor(QUERIES, DATABASE, varargin) @@ -58,7 +63,5 @@ [V,D]=eig(covs{nbr}); - [~,index]=min(abs(diag(D))); - norms(:,elem+nbr-1)=abs(V(:,index)); % norms(:,elem+nbr-1)=abs(V(:,dimensions)); @@ -118,7 +121,7 @@ if isempty(iProperty) - error('psfConv:InvalidProperty', 'Invalid Property'); + error('Damkjer:InvalidProperty', 'Invalid Property'); elseif length(iProperty) > 1 - error('psfConv:AmbiguousProperty', ... + error('Damkjer:AmbiguousProperty', ... 'Supplied shortened property name is ambiguous'); end @@ -144,5 +147,5 @@ userParams.radius = value; else - error('psfConv:InvalidPointWeights', ... + error('Damkjer:InvalidPointWeights', ... 'Radius must be a real valued positive scalar'); end Index: Damkjer/PointProcessing/SpatialAnalyzer/spanalyze.m =================================================================== --- Damkjer/PointProcessing/SpatialAnalyzer/spanalyze.m (revision 3) +++ Damkjer/PointProcessing/SpatialAnalyzer/spanalyze.m (revision 3) @@ -0,0 +1,267 @@ +function [ CLASSES ] = spanalyze( QUERIES, DATABASE, varargin ) +%% +% [CLASSES] = SPANALYZE(QUERIES, DATABASE, ...) +% +% Attribute points by eigen analysis. + +% Parse Inputs +userParams = parseInputs(varargin{:}); + +dimensions=size(QUERIES, 1); +elements=size(QUERIES, 2); + +norms=zeros(size(QUERIES)); +feats=zeros(size(QUERIES)); + +biases=zeros(1,elements); + +ints=zeros(1,elements); + +msg='Computing Structure Tensors...'; +tstart=tic; +h = timebar(1, elements, msg, tstart); + +step=1000; +tic; +for elem=1:step:elements + + last=min(elem+step-1,elements); + + % Get the nearest neighbors of elem + if (userParams.neighbors <= 0) + % Perform a fixed radius search + NN=DATABASE.rnn(QUERIES(:,elem:last),... + userParams.radius); + else + if (userParams.radius <= 0) + % Search unconstrained + NN=DATABASE.knn(QUERIES(:,elem:last),... + userParams.neighbors); + else + % Search constrained to radius + NN=DATABASE.knn(QUERIES(:,elem:last),... + userParams.neighbors,... + 'lim',userParams.radius); + end + end + + [covs,bias,inty]=fastcov(cellfun(@(x) QUERIES(:,x)',NN,'UniformOutput',false)); + +% Same solution via SVD appears to be slower via emperical testing +% recentered=fastcenter(cellfun(@(x) QUERIES(:,x)',NN,'UniformOutput',false)); + + for nbr=1:size(NN,1) + % skip underconstrained (possible with radius searches) + if (length(NN{nbr})<5) + continue; + end + + [V,D]=eig(covs{nbr}); + [~,index]=min(diag(D)); +% norms(:,elem+nbr-1)=abs(V(:,index)) .* sqrt(D(index, index)); + norms(:,elem+nbr-1)=abs(V(:,index)); + feats(:,elem+nbr-1)=sort(diag(D),'descend'); + + biases(elem+nbr-1)=bias{nbr}; + ints(elem+nbr-1)=inty{nbr}; + +% Same solution via SVD appears to be slower via emperical testing +% [~,S,V]=svd(recentered{nbr},0); +% S = S ./ sqrt(length(NN{nbr})-1); +% norms(:,elem+nbr-1)=abs(V(:,dimensions) .* S(dimensions, dimensions)); +% feats(:,elem+nbr-1)=diag(S).^2; + + end + + if (toc > 1) + tic; + h = timebar(elem, elements, msg, tstart, h); + end +end + +if (all(ishghandle(h, 'figure'))) + close(h); +end + +% Compute Omnivariance +omnivar=prod(feats,1).^(1/dimensions); + +% Compute Eigen Entropy +energy=sum(feats,1); +nfeats=bsxfun(@times,feats',1./energy')'; +nfeats(isnan(nfeats))=1; +entropy=-sum(nfeats.*log(nfeats))./log(dimensions); +entropy(entropy==1)=0; + +% Compute Eigen Structure +structure=1-entropy; + +% Compute Fractional Anisotropy +evmeans=mean(feats,1); +evdevs=feats-repmat(evmeans,dimensions,1); +numer=dimensions.*sum(evdevs.^2,1); +denom=(dimensions-1).*sum(feats.^2,1); +fa=(numer./denom).^(1/2); + +% Compute Anisotropy +ani=(feats(1,:)-feats(dimensions,:))./(feats(1,:)); + +% Compute Isotropy +iso=1-ani; + +% Compute dimensional degree +dims=(feats(1:dimensions-1,:)-feats(2:dimensions,:))./... + repmat(feats(1,:),dimensions-1,1); + +% Compute Dimensional Embedding +[~,embed]=max([dims;iso],[],1); + +% Populate feature classes +CLASSES=[feats;norms;dims;iso;ani;fa;entropy;structure;omnivar;embed;biases;ints]; + +% st=zeros(size(QUERIES)); +% +% msg='Computing Normal Structure Tensors...'; +% tstart=tic; +% h = timebar(1, elements, msg, tstart); +% +% for elem=1:step:elements +% +% last=min(elem+step-1,elements); +% +% % Get the nearest neighbors of elem +% if (userParams.neighbors <= 0) +% % Perform a fixed radius search +% NN=DATABASE.rnn(QUERIES(:,elem:last),... +% userParams.radius); +% else +% if (userParams.radius <= 0) +% % Search unconstrained +% NN=DATABASE.knn(QUERIES(:,elem:last),... +% userParams.neighbors); +% else +% % Search constrained to radius +% NN=DATABASE.knn(QUERIES(:,elem:last),... +% userParams.neighbors,... +% 'lim',userParams.radius); +% end +% end +% +% for nbr=1:size(NN,1) +% normals=norms(:,NN{nbr}); +% +% S=svd(normals',0); +% st(:,elem+nbr-1)=S.*S./(length(NN{nbr})-1); +% +% % D=eig(1/(length(NN{nbr})-1)*(normals*normals')); +% % st(:,elem+nbr-1)=sort(D,'descend'); +% +% % ST=normals*normals'; +% +% % if (nbr == 1) +% % D=eig(1/(length(NN{nbr})-1)*(normals*normals')); +% % disp('EIG:'); +% % disp(sqrt(D)); +% % +% % S=svd(normals',0); +% % S = S ./ sqrt(length(NN{nbr})-1); +% % disp('SVD:'); +% % disp(abs(S)); +% % +% % end +% +% end +% +% if(toc > 1) +% tic; +% h = timebar(elem, elements, msg, tstart, h); +% end +% end +% +% if (all(ishghandle(h, 'figure'))) +% close(h); +% end + +% CLASSES = QUERIES; + +end + +%% +% PARSEINPUTS Support function to parse inputs into userParams structure +function [userParams] = parseInputs(varargin) + +userParams = struct('radius', 0, 'neighbors', 0); + +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'}; + +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:InvalidPointWeights', ... + 'Radius must be a real valued positive scalar'); + end + end + + varargin(1:2) = []; +end + +if (userParams.neighbors <= 0 && userParams.radius <= 0) + userParams.radius = 1; +end + +end + Index: Damkjer/Util/Math/fastcenter.cpp =================================================================== --- Damkjer/Util/Math/fastcenter.cpp (revision 0) +++ Damkjer/Util/Math/fastcenter.cpp (revision 3) @@ -4,6 +4,6 @@ // Copyright (C) 2012 Kristian Damkjer. // -// DESCRIPTION: This MEX source file provides a fast implementation of cov -// method for cell-arrays of real valued matrices. +// DESCRIPTION: This MEX source file provides a fast implementation of +// centering method for cell-arrays of real valued matrices. // // LIMITATIONS: Does not work for cell-arrays of complex matrices. Index: Damkjer/Util/Math/fastcoeffvar.cpp =================================================================== --- Damkjer/Util/Math/fastcoeffvar.cpp (revision 3) +++ Damkjer/Util/Math/fastcoeffvar.cpp (revision 3) @@ -0,0 +1,114 @@ +//========================================================================= +// FILE: fastcoeffvar.cpp +// +// Copyright (C) 2012 Kristian Damkjer. +// +// DESCRIPTION: This MEX source file provides a fast implementation of cov +// method for cell-arrays of real valued matrices. +// +// LIMITATIONS: Does not work for cell-arrays of complex matrices. +// +// SOFTWARE HISTORY: +//> 2012-SEP-11 K. Damkjer +// Initial Coding. +//< +//========================================================================= + +#ifdef _OPENMP +#include +#endif + +#include + +#ifdef _CHAR16T +#define CHAR16_T +#endif + +#include "mex.h" + +void mexFunction( + int nlhs, mxArray* plhs[], + int nrhs, const mxArray* prhs[]) +{ + if (nrhs != 1 || !mxIsCell(prhs[0])) + { + mexErrMsgIdAndTxt("Damkjer:fastcoeffvar:varargin", + "Missing or invalid input argument."); + } + + if (nlhs > 1) + { + mexErrMsgIdAndTxt("Damkjer:fastcoeffvar:varargout", + "Too many output arguments."); + } + + mwSize cells = mxGetNumberOfElements (prhs[0]); + + plhs[0] = mxCreateCellMatrix(cells, 1); + + std::vector vals(cells,0); + std::vector Ms(cells,0); + std::vector Ns(cells,0); + + std::vector covs(cells,0); + std::vector cov_vals(cells,0); + + for (int cell = 0; cell < cells; ++cell) + { + vals[cell]=mxGetPr(mxGetCell(prhs[0], cell)); + Ms[cell]=mxGetM(mxGetCell(prhs[0], cell)); + Ns[cell]=mxGetN(mxGetCell(prhs[0], cell)); + + covs[cell] = mxCreateDoubleMatrix(0, 0, mxREAL); + mxSetM(covs[cell], Ns[cell]); + mxSetN(covs[cell], Ns[cell]); + mxSetData(covs[cell], mxMalloc(sizeof(double))); + cov_vals[cell] = mxGetPr(covs[cell]); + } + +#ifdef _OPENMP + omp_set_dynamic(1); + omp_set_num_threads(omp_get_num_procs()); +#endif + + #pragma omp parallel for + for (int cellp = 0; cellp < cells; ++cellp) + { + std::vector dist(Ms[cellp], 0.); + double mean=0; + + // calculate the average distance to this point + double w1 = 1./Ms[cellp]; + + for (mwSize m = Ms[cellp]; m --> 0;) + { + for (mwSize n = Ns[cellp]; n --> 0;) + { + double diff = vals[cellp][m + Ms[cellp] * n]-vals[cellp][Ms[cellp] * n]; + dist[n]+=diff*diff; + } + mean+=sqrt(dist[m])*w1; + } + + double w2 = 1./(Ms[cellp]-1); + + for (mwSize n = Ns[cellp]; n --> 0;) + { + cov_vals[cellp] = 0; + + for (mwSize mc = Ms[cellp]; mc --> 0;) + { + cov_vals[cellp] += + w2 + * (vals[cellp][mc + Ms[cellp] * n]-mean[n] + * (vals[cellp][mc + Ms[cellp] * n]-vals[cellp][Ms[cellp] * n2]); + } + } + } + } + + for (int cell = 0; cell < cells; ++cell) + { + mxSetCell(plhs[0], cell, covs[cell]); + } +} Index: Damkjer/Util/Math/fastcov.cpp =================================================================== --- Damkjer/Util/Math/fastcov.cpp (revision 0) +++ Damkjer/Util/Math/fastcov.cpp (revision 3) @@ -20,4 +20,5 @@ #include +#include #ifdef _CHAR16T @@ -37,5 +38,5 @@ } - if (nlhs > 1) + if (nlhs > 3) { mexErrMsgIdAndTxt("Damkjer:fastcov:varargout", @@ -47,4 +48,16 @@ plhs[0] = mxCreateCellMatrix(cells, 1); + // Better way? + if (nlhs > 1) + { + plhs[1] = mxCreateCellMatrix(cells, 1); + } + + // Better way? + if (nlhs > 2) + { + plhs[2] = mxCreateCellMatrix(cells, 1); + } + std::vector vals(cells,0); std::vector Ms(cells,0); @@ -54,4 +67,13 @@ std::vector cov_vals(cells,0); + //HACK: Probably a more efficient way to do this + std::vector bias(cells,0); + std::vector bias_vals(cells,0); + + //HACK: Probably a more efficient way to do this + std::vector inty(cells,0); + std::vector inty_vals(cells,0); + + // Note for future: Ms - points, Ns - dimensions for (int cell = 0; cell < cells; ++cell) { @@ -59,5 +81,6 @@ Ms[cell]=mxGetM(mxGetCell(prhs[0], cell)); Ns[cell]=mxGetN(mxGetCell(prhs[0], cell)); - + + // We will be setting each value, so don't bother to initialize to zero. covs[cell] = mxCreateDoubleMatrix(0, 0, mxREAL); mxSetM(covs[cell], Ns[cell]); @@ -65,4 +88,12 @@ mxSetData(covs[cell], mxMalloc(sizeof(double)*Ns[cell]*Ns[cell])); cov_vals[cell] = mxGetPr(covs[cell]); + + //HACK: Probably a more efficient way to do this + bias[cell] = mxCreateDoubleMatrix(1, 1, mxREAL); + bias_vals[cell] = mxGetPr(bias[cell]); + + //HACK: Probably a more efficient way to do this + inty[cell] = mxCreateDoubleMatrix(1, 1, mxREAL); + inty_vals[cell] = mxGetPr(inty[cell]); } @@ -77,4 +108,5 @@ std::vector mean(Ns[cellp], 0.); + // Comment out to remove mean calculation when centering to first point. double w1 = 1./Ms[cellp]; @@ -86,4 +118,5 @@ } } + // End comment. double w2 = 1./(Ms[cellp]-1); @@ -97,11 +130,51 @@ for (mwSize mc = Ms[cellp]; mc --> 0;) { + // Center to mean. +// cov_vals[cellp][n2 + Ns[cellp] * n1] += +// w2 +// * (vals[cellp][mc + Ms[cellp] * n1]-mean[n1]) +// * (vals[cellp][mc + Ms[cellp] * n2]-mean[n2]); + + // Center to first point. cov_vals[cellp][n2 + Ns[cellp] * n1] += w2 - * (vals[cellp][mc + Ms[cellp] * n1]-mean[n1]) - * (vals[cellp][mc + Ms[cellp] * n2]-mean[n2]); + * (vals[cellp][mc + Ms[cellp] * n1]-vals[cellp][Ms[cellp] * n1]) + * (vals[cellp][mc + Ms[cellp] * n2]-vals[cellp][Ms[cellp] * n2]); + // End comment. } } } + + if (nlhs > 1) + { + // Experiment: Estimate bias, i.e., how "centered" is our mass? + + //std::vector diff_loc(Ns[cellp], 0.); + //std::vector diff_ext(Ns[cellp], 0.); + + double diff_loc = 0; + double diff_ext = 0; + double temp; + + for (mwSize n = Ns[cellp]; n --> 0;) + { + temp = vals[cellp][Ms[cellp] * n] - mean[n]; + diff_loc += temp * temp; + + temp = vals[cellp][Ms[cellp] * n] - + vals[cellp][Ms[cellp] * (n + 1) - 1]; + diff_ext += temp * temp; + } + + bias_vals[cellp][0] = std::sqrt(diff_loc)/std::sqrt(diff_ext); + // End comment. + + if (nlhs > 2) + { + // Experiment: Estimate intensity, i.e., how "dense" is our mass? + inty_vals[cellp][0] = Ms[cellp]/(4./3.*M_PI*diff_ext*sqrt(diff_ext)); + // End comment. + } + } } @@ -110,3 +183,19 @@ mxSetCell(plhs[0], cell, covs[cell]); } + + if (nlhs > 1) + { + for (int cell = 0; cell < cells; ++cell) + { + mxSetCell(plhs[1], cell, bias[cell]); + } + } + + if (nlhs > 2) + { + for (int cell = 0; cell < cells; ++cell) + { + mxSetCell(plhs[2], cell, inty[cell]); + } + } } Index: Damkjer/Util/Math/makeMath.m =================================================================== --- Damkjer/Util/Math/makeMath.m (revision 0) +++ Damkjer/Util/Math/makeMath.m (revision 3) @@ -22,5 +22,5 @@ function makeMath(varargin) -command = 'mex -largeArrayDims'; +command = 'mex -largeArrayDims -D_USE_MATH_DEFINES'; flags = {'verbose', 'debug', 'warnings', 'parallel'};