From b944c1d7e1fd1033082b1d9ae81265362787b0b4 Mon Sep 17 00:00:00 2001 From: Miguel Date: Tue, 28 Mar 2023 17:54:05 +0100 Subject: [PATCH 1/4] debugged --- src/tools/define_input.m | 20 ++------------------ src/treeqsm.m | 2 ++ 2 files changed, 4 insertions(+), 18 deletions(-) diff --git a/src/tools/define_input.m b/src/tools/define_input.m index 76602ac..0257675 100644 --- a/src/tools/define_input.m +++ b/src/tools/define_input.m @@ -49,24 +49,8 @@ create_input Inputs = inputs; -% If given multiple clouds, extract the names -if ischar(Clouds) || isstring(Clouds) - matobj = matfile([Clouds,'.mat']); - names = fieldnames(matobj); - i = 1; - n = max(size(names)); - while i <= n && ~strcmp(names{i,:},'Properties') - i = i+1; - end - I = (1:1:n); - I = setdiff(I,i); - names = names(I,1); - names = sort(names); - nt = max(size(names)); % number of trees/point clouds -else - P = Clouds; - nt = 1; -end +P = transpose(Clouds); +nt = 1; inputs(nt).PatchDiam1 = 0; diff --git a/src/treeqsm.m b/src/treeqsm.m index 7f37142..2493108 100644 --- a/src/treeqsm.m +++ b/src/treeqsm.m @@ -265,6 +265,8 @@ disp('Progress:') end +P=transpose(P) + %% Make the point cloud into proper form % only 3-dimensional data if size(P,2) > 3 From 223ca8e6f0a9b68f5f39bc81aeaea12b057a9b11 Mon Sep 17 00:00:00 2001 From: Miguel Nogueira Date: Fri, 14 Apr 2023 15:12:56 +0100 Subject: [PATCH 2/4] fixed bugs --- src/tools/define_input.m | 6 +- src/tools/pdist2.m | 183 +++++++++++++++++++++++++++++++++++++++ src/treeqsm.m | 2 - 3 files changed, 187 insertions(+), 4 deletions(-) create mode 100644 src/tools/pdist2.m diff --git a/src/tools/define_input.m b/src/tools/define_input.m index 0257675..5a8e4fb 100644 --- a/src/tools/define_input.m +++ b/src/tools/define_input.m @@ -49,7 +49,9 @@ create_input Inputs = inputs; -P = transpose(Clouds); +% P = transpose(Clouds); +P = Clouds; + nt = 1; inputs(nt).PatchDiam1 = 0; @@ -81,7 +83,7 @@ hSec = hSecTop-hSecBot; Sec = Hei > hSecBot & Hei < hSecTop; StemBot = P(Sec,1:3); - + % Estimate stem axis (point and direction) AxisPoint = mean(StemBot); V = StemBot-AxisPoint; diff --git a/src/tools/pdist2.m b/src/tools/pdist2.m new file mode 100644 index 0000000..4585ae2 --- /dev/null +++ b/src/tools/pdist2.m @@ -0,0 +1,183 @@ +function D = pdist2( X, Y, metric ) + % Calculates the distance between sets of vectors. + % + % Let X be an m-by-p matrix representing m points in p-dimensional space + % and Y be an n-by-p matrix representing another set of points in the same + % space. This function computes the m-by-n distance matrix D where D(i,j) + % is the distance between X(i,:) and Y(j,:). This function has been + % optimized where possible, with most of the distance computations + % requiring few or no loops. + % + % The metric can be one of the following: + % + % 'euclidean' / 'sqeuclidean': + % Euclidean / SQUARED Euclidean distance. Note that 'sqeuclidean' + % is significantly faster. + % + % 'chisq' + % The chi-squared distance between two vectors is defined as: + % d(x,y) = sum( (xi-yi)^2 / (xi+yi) ) / 2; + % The chi-squared distance is useful when comparing histograms. + % + % 'cosine' + % Distance is defined as the cosine of the angle between two vectors. + % + % 'emd' + % Earth Mover's Distance (EMD) between positive vectors (histograms). + % Note for 1D, with all histograms having equal weight, there is a simple + % closed form for the calculation of the EMD. The EMD between histograms + % x and y is given by the sum(abs(cdf(x)-cdf(y))), where cdf is the + % cumulative distribution function (computed simply by cumsum). + % + % 'L1' + % The L1 distance between two vectors is defined as: sum(abs(x-y)); + % + % + % USAGE + % D = pdist2( X, Y, [metric] ) + % + % INPUTS + % X - [m x p] matrix of m p-dimensional vectors + % Y - [n x p] matrix of n p-dimensional vectors + % metric - ['sqeuclidean'], 'chisq', 'cosine', 'emd', 'euclidean', 'L1' + % + % OUTPUTS + % D - [m x n] distance matrix + % + % EXAMPLE + % % simple example where points cluster well + % [X,IDX] = demoGenData(100,0,5,4,10,2,0); + % D = pdist2( X, X, 'sqeuclidean' ); + % distMatrixShow( D, IDX ); + % % comparison to pdist + % n=500; d=200; r=100; X=rand(n,d); + % tic, for i=1:r, D1 = pdist( X, 'euclidean' ); end, toc + % tic, for i=1:r, D2 = pdist2( X, X, 'euclidean' ); end, toc + % D1=squareform(D1); del=D1-D2; sum(abs(del(:))) + % + % See also pdist, distMatrixShow + % + % Piotr's Computer Vision Matlab Toolbox Version 2.52 + % Copyright 2014 Piotr Dollar. [pdollar-at-gmail.com] + % Licensed under the Simplified BSD License [see external/bsd.txt] + + if( nargin<3 || isempty(metric) ); metric=0; end; + + switch metric + case {0,'sqeuclidean'} + D = distEucSq( X, Y ); + case 'euclidean' + D = sqrt(distEucSq( X, Y )); + case 'L1' + D = distL1( X, Y ); + case 'cosine' + D = distCosine( X, Y ); + case 'emd' + D = distEmd( X, Y ); + case 'chisq' + D = distChiSq( X, Y ); + otherwise + error(['pdist2 - unknown metric: ' metric]); + end + D = max(0,D); + end + + function D = distL1( X, Y ) + m = size(X,1); n = size(Y,1); + mOnes = ones(1,m); D = zeros(m,n); + for i=1:n + yi = Y(i,:); yi = yi( mOnes, : ); + D(:,i) = sum( abs( X-yi),2 ); + end + end + + function D = distCosine( X, Y ) + p=size(X,2); + XX = sqrt(sum(X.*X,2)); X = X ./ XX(:,ones(1,p)); + YY = sqrt(sum(Y.*Y,2)); Y = Y ./ YY(:,ones(1,p)); + D = 1 - X*Y'; + end + + function D = distEmd( X, Y ) + Xcdf = cumsum(X,2); + Ycdf = cumsum(Y,2); + m = size(X,1); n = size(Y,1); + mOnes = ones(1,m); D = zeros(m,n); + for i=1:n + ycdf = Ycdf(i,:); + ycdfRep = ycdf( mOnes, : ); + D(:,i) = sum(abs(Xcdf - ycdfRep),2); + end + end + + function D = distChiSq( X, Y ) + % note: supposedly it's possible to implement this without a loop! + m = size(X,1); n = size(Y,1); + mOnes = ones(1,m); D = zeros(m,n); + for i=1:n + yi = Y(i,:); yiRep = yi( mOnes, : ); + s = yiRep + X; d = yiRep - X; + D(:,i) = sum( d.^2 ./ (s+eps), 2 ); + end + D = D/2; + end + + function D = distEucSq( X, Y ) + Yt = Y'; + XX = sum(X.*X,2); + YY = sum(Yt.*Yt,1); + D = bsxfun(@plus,XX,YY)-2*X*Yt; + end + + %%%% code from Charles Elkan with variables renamed + % function D = distEucSq( X, Y ) + % m = size(X,1); n = size(Y,1); + % D = sum(X.^2, 2) * ones(1,n) + ones(m,1) * sum(Y.^2, 2)' - 2.*X*Y'; + % end + + %%% LOOP METHOD - SLOW + % [m p] = size(X); + % [n p] = size(Y); + % D = zeros(m,n); + % onesM = ones(m,1); + % for i=1:n + % y = Y(i,:); + % d = X - y(onesM,:); + % D(:,i) = sum( d.*d, 2 ); + % end + + %%% PARALLEL METHOD THAT IS SUPER SLOW (slower than loop)! + % % From "MATLAB array manipulation tips and tricks" by Peter J. Acklam + % Xb = permute(X, [1 3 2]); + % Yb = permute(Y, [3 1 2]); + % D = sum( (Xb(:,ones(1,n),:) - Yb(ones(1,m),:,:)).^2, 3); + + %%% USELESS FOR EVEN VERY LARGE ARRAYS X=16000x1000!! and Y=100x1000 + % call recursively to save memory + % if( (m+n)*p > 10^5 && (m>1 || n>1)) + % if( m>n ) + % X1 = X(1:floor(end/2),:); + % X2 = X((floor(end/2)+1):end,:); + % D1 = distEucSq( X1, Y ); + % D2 = distEucSq( X2, Y ); + % D = cat( 1, D1, D2 ); + % else + % Y1 = Y(1:floor(end/2),:); + % Y2 = Y((floor(end/2)+1):end,:); + % D1 = distEucSq( X, Y1 ); + % D2 = distEucSq( X, Y2 ); + % D = cat( 2, D1, D2 ); + % end + % return; + % end + + %%% L1 COMPUTATION WITH LOOP OVER p, FAST FOR SMALL p. + % function D = distL1( X, Y ) + % + % m = size(X,1); n = size(Y,1); p = size(X,2); + % mOnes = ones(1,m); nOnes = ones(1,n); D = zeros(m,n); + % for i=1:p + % yi = Y(:,i); yi = yi( :, mOnes ); + % xi = X(:,i); xi = xi( :, nOnes ); + % D = D + abs( xi-yi' ); + % end \ No newline at end of file diff --git a/src/treeqsm.m b/src/treeqsm.m index 2493108..7f37142 100644 --- a/src/treeqsm.m +++ b/src/treeqsm.m @@ -265,8 +265,6 @@ disp('Progress:') end -P=transpose(P) - %% Make the point cloud into proper form % only 3-dimensional data if size(P,2) > 3 From b4b02d2e55d2ae320be9e92ba9fe4068b182081f Mon Sep 17 00:00:00 2001 From: Miguel Nogueira Date: Fri, 14 Apr 2023 18:40:16 +0100 Subject: [PATCH 3/4] removed alphaShape depedency, removed outputs --- src/create_input.m | 8 ++++---- src/main_steps/tree_data.m | 12 ++++++++---- 2 files changed, 12 insertions(+), 8 deletions(-) diff --git a/src/create_input.m b/src/create_input.m index bc74e42..7513c9b 100644 --- a/src/create_input.m +++ b/src/create_input.m @@ -105,14 +105,14 @@ % Save the output struct QSM as a matlab-file into \result folder. % If name = 'pine', tree = 2, model = 5, the name of the saved file is % 'QSM_pine_t2_m5.mat': -inputs.savemat = 1; +inputs.savemat = 0; % Save the models in .txt-files (check "save_model_text.m"): -inputs.savetxt = 1; +inputs.savetxt = 0; % What are plotted during reconstruction process: % 2 = plots the QSM, the segmentated point cloud and distributions, % 1 = plots the QSM and the segmentated point cloud % 0 = plots nothing -inputs.plot = 2; +inputs.plot = 0; % What are displayed during the reconstruction: 2 = display all; % 1 = display name, parameters and distances; 0 = display only the name: -inputs.disp = 2; +inputs.disp = 0; diff --git a/src/main_steps/tree_data.m b/src/main_steps/tree_data.m index aaadbaf..7db634a 100644 --- a/src/main_steps/tree_data.m +++ b/src/main_steps/tree_data.m @@ -593,8 +593,10 @@ %% Crown areas from convex hull and alpha shape: treedata.CrownAreaConv = A; alp = max(0.5,treedata.CrownDiamAve/10); -shp = alphaShape(X(:,1),X(:,2),alp); -treedata.CrownAreaAlpha = shp.area; +% Miguel: alphaShape is NOT implemented in Octave +%shp = alphaShape(X(:,1),X(:,2),alp); +%treedata.CrownAreaAlpha = shp.area; +treedata.CrownAreaAlpha = -9999.; %% Crown base % Define first major branch as the branch whose diameter > min(0.05*dbh,5cm) @@ -682,8 +684,10 @@ [K,V] = convhull(X(:,1),X(:,2),X(:,3)); treedata.CrownVolumeConv = V; alp = max(0.5,treedata.CrownDiamAve/5); - shp = alphaShape(X(:,1),X(:,2),X(:,3),alp,'HoleThreshold',10000); - treedata.CrownVolumeAlpha = shp.volume; + % Miguel: alphaShape is NOT implemented in Octave + %shp = alphaShape(X(:,1),X(:,2),X(:,3),alp,'HoleThreshold',10000); + %treedata.CrownVolumeAlpha = shp.volume; + treedata.CrownVolumeAlpha = -9999.; else % No branches From 918c5b446f4131d00b818d3cffdcacbe025a079f Mon Sep 17 00:00:00 2001 From: Miguel Nogueira Date: Thu, 11 May 2023 07:52:31 +0100 Subject: [PATCH 4/4] fix to plotting --- src/create_input.m | 2 +- src/plotting/plot_cylinder_model.m | 2 +- src/plotting/plot_models_segmentations.m | 2 +- src/treeqsm.m | 2 ++ 4 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/create_input.m b/src/create_input.m index 7513c9b..e2b391e 100644 --- a/src/create_input.m +++ b/src/create_input.m @@ -91,7 +91,7 @@ inputs.filter.EdgeLength = 0.004; % Plot the filtering results automatically after the filtering if % filter.plot > 0 -inputs.filter.plot = 1; +inputs.filter.plot = 0; %% Other inputs % These parameters don't affect the QSM-reconstruction but define what is diff --git a/src/plotting/plot_cylinder_model.m b/src/plotting/plot_cylinder_model.m index a12154c..2f111b2 100644 --- a/src/plotting/plot_cylinder_model.m +++ b/src/plotting/plot_cylinder_model.m @@ -167,7 +167,7 @@ function plot_cylinder_model(cylinder,Color,fig,nf,alp,Ind) figure(fig) plot3(Vert(1,1),Vert(1,2),Vert(1,3)) patch('Vertices',Vert,'Faces',Facets,'FaceVertexCData',fvd,'FaceColor','flat') -alpha(alp) +% alpha(alp) axis equal grid on view(-37.5,30) diff --git a/src/plotting/plot_models_segmentations.m b/src/plotting/plot_models_segmentations.m index c94b900..de5a718 100644 --- a/src/plotting/plot_models_segmentations.m +++ b/src/plotting/plot_models_segmentations.m @@ -62,7 +62,7 @@ function plot_models_segmentations(P,cover,segment,cylinder,trunk,triangulation) plot_cylinder_model(cylinder,'order',2,10) subplot(1,2,2) plot_cylinder_model(cylinder,'branch',2,10) - +saveas(gcf,'/mnt/c/contenda/modelo.png') %% figure 3, segmented point cloud and cylinder model plot_branch_segmentation(P,cover,segment,'order',3,1) hold on diff --git a/src/treeqsm.m b/src/treeqsm.m index 7f37142..23b6658 100644 --- a/src/treeqsm.m +++ b/src/treeqsm.m @@ -492,10 +492,12 @@ %% Plot models and segmentations if inputs.plot >= 1 if inputs.Tria + disp('plotting tria') plot_models_segmentations(P,cover2,segment2,cylinder,trunk,... triangulation) else plot_models_segmentations(P,cover2,segment2,cylinder) + disp('plotting ') end if nd > 1 || na > 1 || ni > 1 pause