# https://savannah.gnu.org/bugs/?55751#comment15
# Test coplanar of points (x,y,z) by least square fit the data points to a planar.
# Eddy

function [copl, normal_vector, eig_vals] = is_coplanar_lsm (x, y, z, tol_dist, lsm_mode)

if ~exist('tol_dist', 'var')
  lsm_mode = 'eig';      % 'eig' or 'svd' or 'lsm'
end

vecs = [x(:) y(:) z(:)];
vecs -= mean(vecs, 1);  % So that the least square planar pass through origin.

switch lsm_mode
case 'svd'
  [~, s, eig_vecs] = svd(vecs / sqrt(rows(vecs)), 'econ');
  eig_vals = diag(s) .^ 2;

case 'eig'
  [eig_vecs, eig_vals] = eig(vecs' * vecs / rows(vecs));
  eig_vals = diag(eig_vals);

case 'lsm'
  var_each = var(vecs, 1, 1);
  [~, id_var_each] = sort(var_each);  % So that maximize var(z)
  if max(var_each) == 0
    copl = true;
    normal_vector = [1 0 0];
    eig_vals = 0;
    return
  end

  % Least square fit: [x y] * [a; b] = z + error, minimize var(error).
  vecs2 = vecs(:, id_var_each(2:3));
  coef = (vecs2' * vecs2) \ (vecs2' * vecs(:, id_var_each(1)))
%  coef = vecs2 \ vecs(:, id_var_each(3));

  eig_vals = [var(vecs(:, id_var_each(1)) - vecs2 * coef, 1), ...
              var_each(id_var_each(3))];
  
  coef = [coef; -1];
  coef = coef / norm(coef);
  eig_vecs = zeros(3,1);
  eig_vecs(id_var_each) = coef;
end

% a normal vector of the planar, if any.
[~, id_min_eig] = min(eig_vals);
normal_vector = eig_vecs(:, id_min_eig);

% default tolerance of distance deviate from the planar.
if ~exist('tol_dist', 'var') || isempty(tol_dist)
  % Note: the factor sqrt(eps) comes from vecs' * vecs, which squared 
  % the conditional number (shadows smaller eigenvalue).
  % To improve, use lsm_mode = 'svd'.
  tol_dist = sqrt(max(eig_vals) * 10*eps());
end

criteria_norm = 2;  % L2 norm or L0 norm.

if criteria_norm == 2
  % min(eig_vals) is the "variance" of the points deviate from the planar.
  % NOTE: seems a bug in eig(): sometimes return negative eigenvalue for 
  %       symmetric positive semi-defined matrix. e.g. avoid sqrt(min(eig_vals)).
  copl = min(eig_vals) < tol_dist^2;
else
  % alternatively, check the max deviation from planar.
  max_dist = max(abs(vecs * normal_vector));
  copl = max_dist < tol_dist;
end

endfunction
