Is it possible numerical problem ?

1 view (last 30 days)
AI-CHI Chang
AI-CHI Chang on 7 Mar 2023
Commented: Askic V on 7 Mar 2023
I'm sorry to post long code here. But I just can't find the problem...
I tried to boost my code efficiency by removing for loop, but it came out two different results from them...
I check the first point's variable d_d, d_n, w_d, w_n in first loop. All them are the same but became different in variable sum_den
61.330229160510050 and 61.330229160510020.....That's so weired. Is it possible the numerical errors?
The input of the former one is 3D points location, neighbors for each point (num_pt*num_neighbor*3), and normal vector for each points.
for i = 1:num_pt
p_nn(i,:,:) = p(knn(i,2:end),:);
end
function [p] = BilateralFilter(p,p_nn,normal)
min_z = min(p(:,3));
maxdepth = max(p(:,3)) - min_z;
sigma_n = 0.5;
sigma_d = 0.01 + 0.5*(p(:,3) - min_z)/maxdepth;
p_permute = permute(p, [1 3 2]);
n_permute = permute(normal, [1 3 2]);
d_d = sqrt(sum((p_permute - p_nn).^2,3));
d_n = abs(sum((p_nn - p_permute).*n_permute,3));
w_d = exp(-d_d.^2./2./sigma_d.^2);
w_n = exp(-d_n.^2./2./sigma_n.^2);
residual = sum((p_nn - p_permute).*n_permute,3);
sum_num = sum(w_d.*w_n.*residaul,2);
sum_den = sum(w_d.*w_n,2);
del_p = sum_num./sum_den;
p = p + del_p.*normal;
The latter one is 3D points location, index of neighbors for each point (num_pt*num_neighbor), and normal vector for each points.
function [p] = BilateralFilter(p,knn,normal)
num_pt = size(p,1);
del_p = zeros(num_pt,1);
min_z = min(p(:,3));
maxdepth = max(p(:,3)) - min_z;
for i = 1:num_pt
idx_nn = knn(i,2:end);
p_nn = p(idx_nn,:);
sigma_n = 0.5;
sigma_d = 0.01 + 0.5*(p(i,3)-min_z)/maxdepth;
d_d = abs(sqrt(sum((p(i,:)-p_nn).*(p(i,:)-p_nn),2)));
d_n = abs(sum((p_nn-p(i,:)).*normal(i,:),2));
w_d = exp(-d_d.^2./2./sigma_d^2);
w_n = exp(-d_n.^2./2./sigma_n^2);
residual = sum((p_nn-p(i,:)).*normal(i,:),2);
sum_num = sum(w_d.*w_n.*residaul);
sum_den = sum(w_d.*w_n);
del_p(i) = sum_num/sum_den;
end
p = p + del_p.*normal;

Answers (1)

Askic V
Askic V on 7 Mar 2023
This is called a round-off error caused by floating point arithmetic.
Please read this (it is for Python, but equally applicable):
https://pythonnumericalmethods.berkeley.edu/notebooks/chapter09.03-Roundoff-Errors.html
  2 Comments
AI-CHI Chang
AI-CHI Chang on 7 Mar 2023
Edited: AI-CHI Chang on 7 Mar 2023
@Askic V Thank you for this information! However, can round-off errors cause significant differences ?( The website run 1000 times ) I ran the iterations only five times, and they all looked completely different.
Askic V
Askic V on 7 Mar 2023
Please check whther this might be helpful:
https://www.mathworks.com/help/symbolic/vpa.html?searchHighlight=vpa&s_tid=srchtitle_vpa_1

Sign in to comment.

Categories

Find more on Interpolation in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!