Publish the code.

Author: gyshgx868

.DS_Store

@@ -1 +1,30 @@
-# feature-graph-learning
+# Feature Graph Learning for 3D Point Cloud Denoising
+
+This repository is the official MATLAB implementation of the following paper:
+
+Wei Hu, Xiang Gao, Gene Cheung, Zongming Guo, "Feature Graph Learning for 3D Point Cloud Denoising," IEEE Transactions on Signal Processing (TSP), February, 2020.
+
+# Usage
+
+This code is tested on `MATLAB R2018a` under Windows 10 x64 platform. You can also run the following command in the MATLAB command line on Linux/macOS:
+
+```
+MATLAB> main
+```
+
+**Note:** you can directly run the script from the MATLAB command window.
+
+To support more datasets and different noise levels, you can run the `add_gaussian_noise` script to generate your own data and place them in the `models` folder.
+
+# Reference
+
+Please cite our paper if you use any part of the code from this repository:
+
+```
+@article{hu2019feature,
+  title={Feature Graph Learning for 3D Point Cloud Denoising},
+  author={Hu, Wei and Gao, Xiang and Cheung, Gene and Guo, Zongming},
+  journal={IEEE Transactions on Signal Processing (TSP)},
+  year={2020}
+}
+```

README.md

@@ -0,0 +1,19 @@
+function add_gaussian_noise(cloud, shapename, variance)
+  pts = cloud.Location;
+  SAMPLING_SET = 200;
+  srf = struct('X', pts(:, 1), 'Y', pts(:, 2), 'Z', pts(:, 3));
+  % estimate diameter
+  ifps = fps_euc(srf, SAMPLING_SET);
+  Dfps = pdist2(pts(ifps, :), pts(ifps, :));
+  diam = sqrt(max(Dfps(:)));
+
+  sig = diam * variance;
+  new_x = pts + randn(size(pts)) * sig;
+  new_pc = pointCloud(new_x);
+  new_pc.Color = cloud.Color;
+  new_pc.Normal = cloud.Normal;
+
+  ply_filename = [shapename, '_gaussian_noise_', num2str(variance)];        
+  pcwrite(new_pc, ['models/noise/', ply_filename, '.ply']);
+  disp(['models/noise/', ply_filename, '.ply saved.']);
+end

add_gaussian_noise.m

@@ -0,0 +1,82 @@
+%% Create currupt version of the model
+% Takes a shape filename (assumed to be of the name XXX_GT.mat), and
+% creates noisy versions of it.
+% options is a structure with the following optional fields:
+% 'GaussianNoise' - a list of standard deviations for Gaussian additive corrupted versions of the surface. The standard deviation is expressed as a fraction of the surface' diameter
+% 'ShotNoiseProbability','ShotNoiseNumber','ShotNoiseMagnitude' - parameter lists for shot noise (i.e point are replaced by arbitrary 3D values. Parameter lists should be of equal size
+% 'DepthShotNoiseProbability','DepthShotNoiseMagnitude','DepthShotNoiseNumber' - for depth shot noise (i.e point are perturbed in the direction of the ray leading from the range scanner to the object
+% 'CameraCenter' - assumed camera location, for depth noise
+
+function X = addnoise(X_gt, shapename, options)
+if (exist('options','var')==0)
+    options=[];
+end
+% defaults=struct('GaussianNoise',[0.00125 0.0025 0.005 0.01 ],...
+%     'ShotNoiseProbability',[0 0],...
+%     'ShotNoiseNumber',[10 100],...
+%     'ShotNoiseMagnitude',[1 1]*0.2,...
+%     'DepthShotNoiseProbability',0,...
+%     'DepthShotNoiseMagnitude',0.2,...
+%     'DepthShotNoiseNumber',100,...
+%     'CameraCenter',[ 0 0 -6]); % optical center for GIP camera
+% options=incorporate_defaults(options,defaults);
+
+SAMPLING_SET=200;
+srf=struct('X',X_gt(:,1),'Y',X_gt(:,2),'Z',X_gt(:,3));
+% estimate diameter
+ifps=fps_euc(srf,SAMPLING_SET);
+Dfps=pdist2(X_gt(ifps,:),X_gt(ifps,:));
+diam=sqrt(max(Dfps(:)))
+
+% % Create surface with holes
+% for i=1:length(options.DepthShotNoiseProbability)
+%     prb=options.DepthShotNoiseProbability(i);
+%     X=X_gt;
+%     if (prb>0)
+%         idx=find(rand(size(X_gt(:,1)))<prb);
+%         prb_name=num2str(prb);
+%     else
+%         idx=randperm(size(X_gt,1),options.DepthShotNoiseNumber(i));
+%         prb_name=num2str(options.DepthShotNoiseNumber(i));
+%     end
+%     dr=(X_gt(idx,:)-repmat(options.CameraCenter(:)',[numel(idx),1]));
+%     ndr=sqrt(sum(dr.^2,2)+1e-6);
+%     X(idx,:)=X_gt(idx,:)+(repmat(randn(size(X_gt(idx,1))),[1 3]).*bsxfun(@rdivide,dr,ndr)*diam*options.DepthShotNoiseMagnitude(i));
+%     mat_filename=[shapename,'_depth_shot_noise_',prb_name,'.mat'];
+%     save(mat_filename,'X');
+%     ply_filename=[shapename,'_depth_shot_noise_',prb_name,'.ply'];
+%     write_ply_only_points(X,ply_filename);
+% end
+
+% Create a Gaussian noised version of the surface
+%for i=1:length(options.GaussianNoise)
+    sig=diam*options.GaussianNoise
+    X=X_gt+randn(size(X_gt))*sig;    
+    %mat_filename=[shapename,'_gaussian_noise_',num2str(sig),'.mat'];
+    %save(mat_filename,'X');    
+    ply_filename=[shapename,'_gaussian_noise_',num2str(options.GaussianNoise)];        
+    %ply_filename=[shapename,'_gaussian_noise_',num2str(sig),'.ply'];
+    %#####################################
+    pcwrite(pointCloud(X), ['models/noise/' ply_filename '.ply']);
+    %####################################
+%end
+
+% % Create a 3D shot noise version of the surface
+% for i=1:length(options.ShotNoiseProbability)
+%     prb=options.ShotNoiseProbability(i);
+%     X=X_gt;
+%     if (prb>0)
+%         idx=find(rand(size(X_gt(:,1)))<prb);
+%         prb_name=num2str(prb);
+%     else
+%         idx=randperm(size(X_gt,1),options.ShotNoiseNumber(i));
+%         prb_name=num2str(options.ShotNoiseNumber(i));
+%     end
+%     X(idx,:)=X_gt(idx,:)+randn(size(X_gt(idx,:)))*diam*options.ShotNoiseMagnitude(i);
+%     mat_filename=[shapename,'_shot_noise_',prb_name,'.mat'];
+%     save(mat_filename,'X');
+%     ply_filename=[shapename,'_shot_noise_',prb_name,'.ply'];
+%     write_ply_only_points(X,ply_filename);
+% end
+
+end

addnoise.m

@@ -0,0 +1,121 @@
+function [M, weights, cur_val, min_iter] = block_gradient_descent( ...
+  patch_s, patch_t, use_normal, use_cosine, use_color, C, log_file)
+
+if ~exist('use_normal', 'var')
+  use_normal = true;
+end
+if ~exist('use_cosine', 'var')
+  use_cosine = false;
+end
+if ~exist('use_color', 'var')
+  use_color = false;
+end
+if ~exist('C', 'var')
+  C = 5.0;
+end
+if ~exist('log_file', 'var')
+  log_file = './log.txt';
+end
+
+indices = 1:size(patch_s, 1);
+
+% initialize variables
+diff = patch_s(indices, 4:6) - patch_t(indices, 4:6);  % 9
+if use_normal && use_cosine
+  norm_s = sum(patch_s(indices, 7:9).^2, 2).^0.5;
+  norm_t = sum(patch_t(indices, 7:9).^2, 2).^0.5;
+  cos_diff = sum(patch_s(indices, 7:9) .* patch_t(indices, 7:9), 2);
+  cos_diff = 1.0 - abs(cos_diff ./ norm_s ./ norm_t);
+  diff = [diff, cos_diff];
+elseif use_normal
+  diff = [diff, patch_s(indices, 7:9) - patch_t(indices, 7:9)];
+end
+if use_color
+  diff = [diff, patch_s(indices, 10:12) - patch_t(indices, 10:12)];
+end
+dist = sum(abs(patch_s(indices, 1:3) - patch_t(indices, 1:3)).^2, 2);
+
+[node_dim, feat_dim] = size(diff);
+
+max_algo_iter = 10;
+max_diag_iter = 10;
+max_block_iter = 10;
+step = 1e-6;
+
+M = eye(feat_dim) * C / feat_dim;
+min_val = Inf;
+min_M = M;
+min_iter = 1;
+cur_val = zeros(max_algo_iter, 1);
+diag_val = zeros(max_diag_iter, 1);
+
+for iter = 1 : max_algo_iter
+  % optimize block
+  for f = 1 : feat_dim
+    M21 = [M(1:f-1, f); M(f+1:end, f)];
+    M22 = [M(1:f-1, 1:f-1), M(1:f-1, f+1:end); ...
+      M(f+1:end, 1:f-1), M(f+1:end, f+1:end)];
+    g1 = diff(:, f);
+    g2 = [diff(:, 1:f-1), diff(:, f+1:end)];
+    theta_min = eigs(M22, 1, 'smallestabs');
+    dist_hat = exp(-sum(g2*M22.*g2, 2)) .* dist;
+    m11 = M(f, f);
+    for block_iter = 1 : max_block_iter
+      exp_obj = exp(-(g1.^2*m11 + 2.0*sum(g1*M21'.*g2,2)));
+      grad = -2.0 * sum(g1.*g2.*exp_obj.*dist_hat);
+      M21 = M21 - step * grad';
+      if norm(M21, 2) > sqrt(theta_min * m11)
+        M21 = M21 * sqrt(theta_min * m11) / norm(M21, 2);
+      end
+    end
+    M(1:f-1, f) = M21(1:f-1);
+    M(f+1:end, f) = M21(f:end);
+    M(f, 1:f-1) = M21(1:f-1)';
+    M(f, f+1:end) = M21(f:end)';
+  end
+  
+  % optimize diagonal
+  for diag_iter = 1 : max_diag_iter
+    diag_M = diag(M);
+    dist_tilde = exp(-sum(diff * (M - diag(diag_M)) .* diff, 2)) .* dist;
+    row_M = sum(abs(M), 2) - abs(diag_M);
+
+    grad_M = zeros(feat_dim, 1);
+    temp = 0;
+    for f = 1 : feat_dim
+      temp = temp + diff(:, f).^2 * diag_M(f);
+    end
+    temp = exp(-temp);
+    for f = 1 : feat_dim
+      grad_M(f) = -sum(temp .* diff(:, f).^2 .* dist_tilde);
+    end
+    diag_M = diag_M - step * grad_M;
+    temp_diag = min(max(diag_M, row_M), C);
+    if sum(temp_diag) <= C
+      diag_M = temp_diag;
+    else
+      diag_M = diag_M - (sum(diag_M) - C) / feat_dim;
+      diag_M = min(max(diag_M, row_M), C);
+    end
+    M(1:feat_dim+1:end) = diag_M;
+    
+    weights = exp(-sum((diff * M) .* diff, 2));
+    diag_val(diag_iter) = sum(weights .* dist);
+  end
+
+  weights = exp(-sum((diff * M) .* diff, 2));
+  iter_val = sum(weights .* dist);
+  cur_val(iter) = iter_val;
+  if iter_val < min_val
+    min_val = iter_val;
+    min_iter = iter;
+    min_M = M;
+  end
+
+  write_log(log_file, ['  ', int2str(iter), '/', int2str(max_algo_iter), ' done'], false);
+end
+
+M = min_M;
+weights = exp(-sum((diff * M) .* diff, 2));
+
+end

block_gradient_descent.m

@@ -0,0 +1,7 @@
+function result = bool2str(value)
+  if value
+    result = 'true';
+  else
+    result = 'false';
+  end
+end

bool2str.m

@@ -0,0 +1,7 @@
+function [result] = build_patches(cloud, centers, size)
+  result = zeros(length(centers), size);
+  for i = 1 : length(centers)
+    [indices, ~] = findNearestNeighbors(cloud, centers(i, :), size);
+    result(i, :) = indices';
+  end
+end

build_patches.m

@@ -0,0 +1,22 @@
+function [result] = comput_total_variations(cloud_normal, patches)
+  [rows, ~] = size(patches);
+  result = zeros(rows, 1);
+  for index = 1 : rows
+    patch = patches(index, :);
+    sum_angle = 0.0;
+    pair_count = 0;
+    for i = 1 : length(patch)
+      na = cloud_normal(patch(i), :);
+      for j = 1 : length(patch)
+        if patch(i) == patch(j)
+          continue;
+        end
+        nb = cloud_normal(patch(j), :);
+        sum_angle = sum_angle + (1.0 - abs(na * nb' / (norm(na) * norm(nb))));
+        pair_count = pair_count + 1;
+      end
+    end
+    sum_angle = pair_count;
+    result(index) = sum_angle / pair_count;
+  end
+end

comput_total_variations.m

@@ -0,0 +1,13 @@
+function [result] = compute_principal_normals(cloud_normal, patches)
+  [rows, ~] = size(patches);
+  result = zeros(rows, 3);
+  for i = 1 : rows
+    patch_normal = cloud_normal(patches(i, :), :);
+    sum_normal = sum(patch_normal);
+    normalized = norm(sum_normal);
+    if normalized < 1e-8
+      normalized = 1.0;
+    end
+    result(i, :) = sum_normal / normalized;
+  end
+end

compute_principal_normals.m

@@ -0,0 +1,4 @@
+function [snr] = compute_snr(cloud, mse)
+  s = mean(sum(cloud.Location .^ 2, 2) .^ 0.5);
+  snr = 10.0 * log(s / mse);
+end