########################################################################
##
## Copyright (C) 2011-2021 The Octave Project Developers
##
## See the file COPYRIGHT.md in the top-level directory of this
## distribution or .
##
## This file is part of Octave.
##
## Octave is free software: you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING. If not, see
## .
##
########################################################################
classdef scatteredInterpolant
## -*- texinfo -*-
## @deftypefn {} {@var{f} =} scatteredInterpolant
## @deftypefnx {} {@var{f} =} scatteredInterpolant (@var{x}, @var{y}, @var{q})
## @deftypefnx {} {@var{f} =} scatteredInterpolant (@var{x}, @var{y}, @var{z}, @var{q})
## @deftypefnx {} {@var{f} =} scatteredInterpolant (@var{P}, @var{q})
## @deftypefnx {} {@var{f} =} scatteredInterpolant (@dots{}, @var{method})
## @deftypefnx {} {@var{f} =} scatteredInterpolant (@dots{}, @var{method}, @var{extrapolationmethod})
##
## Generate a 2D or 3D interpolating function @var{f} defined by an input set
## of scattered datapoints @var{x}, @var{y}, @var{z}, and function values,
## @var{q} at those points.
##
## @var{x}, @var{y}, @var{z} must be vectors of the same length as @var{q}.
## The input points may also be specified as a single matrix, @var{P} where
## each row of @var{P} contains the coordinates of a single point
## corresponding to function values in @var{q}.
##
## The returned function can then be evaluated as any of:
##
## @example
## @group
## @var{Si} = @var{f} (@var{xi}, @var{yi})
## @var{Si} = @var{f} (@var{xi}, @var{yi}, @var{zi})
## @var{Si} = @var{f} (@var{Pi})
## @var{Si} = @var{f} (@{@var{Xg}, @var{Yg}@})
## @var{Si} = @var{f} (@{@var{Xg}, @var{Yg}, @var{Zg}@})
## @end group
## @end example
##
## where @var{xi}, @var{yi}, (and @var{zi} for 3D points) are equal sized
## arrays containing coordinates of every query point, and @var{Pi} refers to
## a single matrix of query points, with each row consisting of the @var{x},
## @var{y} (and @var{z}) coordinates of a single point. @{@var{Xg}, @var{Yg},
## @var{Zg}@} refers to a cell array of grid vectors defining each axis of
## a grid space to be interpolated.
##
## The output Si will contain the interpolated values with the same array
## shape as the inputs, except for grid vector inputs, where the output will
## be the shape of the defined grid. Note that the grid vector ordering uses
## @code{ndgrid}, where the first vector is expanded along dim-1, and not
## @code{meshgrid}, where the first vector is expanded along the 'horizontal'
## or dim-2 dimension. @xref{Three-Dimensional Plots} for more detail.
##
## Interpolation is based on Delaunay triangulation and can be controlled by
## supplying options for @var{method}, which can take the values of:
##
## @table @code
## @item nearest
## Discontinuous interpolation assigning value of nearest sample data point.
##
## @item linear
## (Default) Linear interpolation between nearest grid points according to
## local triangulation.
##
## @item natural
## Natural neighbor interpoltaion
## @end table
##
## Extrapolation behavior can be defined by providing by supplying an optional
## @var{ExtrapolationMethod} option. Valid extrapolation methods are:
##
## @table @code
## @item none
## No extrapolation is performed and any points outside of the convex hull of
## the points used to define @var{f} will return a value of NaN.
##
## @item linear
## Linear extrapolation will be performed based on the gradient at nearby
## triangulation boundaries. This is the default behavior when interpolation
## method is either @option{linear} or @option{natural}.
##
## @item nearest
## Nearest neighbor extrapolation will return the value of the nearest
## neigboring point on the boundary. This is the default when interpolation
## method is set to @option{nearest}.
## @end table
##
## @sc{Matlab} compatibility note: Octave's @code{scatteredInterpolant} class
## uses @code{delaunayn} to calculate the triangulation for interpolation,
## whereas @sc{Matlab} likely uses the newer delaunayTriangulation objects.
## Certain point inputs fail for delaunayn in both programs that can be
## handled by delaunayTriangulation (for example the 8 corner points of a 3D
## cube). Additionally, the existance of multiple valid Delaunay
## triangulations for a given point set can produce valid, but differing
## differing interpolation results between the two programs. This will be an
## unavoidable scatteredInterpolant output compatibilty limitation at least
## until a compatible delaunayTriangulation is implemented.
##
## @seealso{TriScatteredInterp, griddata, delaunay, delaunayn}
## @end deftypefn
properties (Access = public)
Points = [] # 2D or 3D array of sampling Points used to define interpolation
Values = [] # function values at sampling Points
Method = "linear" # interpolation method: linear, nearest, or natural
ExtrapolationMethod = "linear" # extrapolation method: linear, nearest, or none
%% endproperties
%%
%% properties (Access = private, Hidden = true)
dimension = 0 #values - 0 (empty), 2D, or 3D, from columns(Points)
tri = [] #stored delaunay triangulation
## state variables for warnings/errors. start all as 'valid' for empty
## constructor
valid_tri = true; #is current triangulation valid
enough_points = true; # are there enough points (>3 for 2D, >4 for 3D)
valid_points_vals = true; # do points and values match in number/dimensions
endproperties
methods (Access = public)
function this = scatteredInterpolant (varargin) ##object constructor
if (nargin == 0)
## do nothing, defaults already set, avoid other error checks.
elseif (nargin == 1 || nargin > 6)
## if not empty, varargin must contain 2-6 elements
print_usage (this);
#error ("scatteredInterpolant: invalid number of inputs");
elseif (isempty (varargin{1}) || isempty (varargin{2}))
error ("scatteredInterpolant: input points and values cannot be empty");
else
## must have at least two array inputs at the front. can have one or two
## char inputs, but must be last.
## array inputs must either be 3 or 4 equal length vectors, or one
## 2 or 3 column array and one vector, with number of rows in the array
## equal to vector length. The last numeric input is the value vector q.
## anything that follows must be 1 or 2 char inputs.
##
## Compatibility note - Matlab 2021a errors if not column vector inputs, but
## docs just say vectors. Allowing vectors as a superset of compatible
## function and in case later behavior changes to allow more general inputs.
numer_input_loc = cellfun(@isnumeric, varargin);
numer_vector_loc = cellfun(@isvector, varargin) & numer_input_loc;
char_input_loc = cellfun(@ischar, varargin);
## make sure no input that isn't numeric or a char
if !all (numer_input_loc + char_input_loc)
error ("scatteredInterpolant: inputs must be numeric arrays or option strings ");
endif
## make sure all numeric input comes before char input
## no char inputs produce empty comparison result which if counts as false
if (find (numer_input_loc, 1, "last") >= find (char_input_loc, 1, "first"))
print_usage (this);
endif
num_input_count = sum (numer_input_loc);
char_input_count = sum (char_input_loc);
## check that number of data and method inputs don't fall out of bounds
if ((num_input_count < 2 || num_input_count > 4) || (char_input_count > 2))
print_usage (this);
endif
## make sure last numeric input, should be q, is a vector
if !(numer_vector_loc(num_input_count) == 1)
error ("scatteredInterpolant: Value input must be a numeric vector.");
endif
## make all vectors are column vectors for easier processing later.
varargin(numer_vector_loc) = cellfun (@(x) x(:), ...
varargin(numer_vector_loc), "UniformOutput", false);
## make sure all numeric inputs have equal number of rows as q
numpoints = numel (varargin{num_input_count});
if any (cellfun(@rows, varargin(numer_input_loc)) != numpoints);
error ("scatteredInterpolant: Point and Value inputs must have the same number of rows.");
endif
## process numeric inputs, assign to this.Points and this.Values.
## check that numeric inputs are either 3-4 columns or
## first one is 2-3 wide array and second is a column
switch num_input_count
case 2
if ((!ismatrix (varargin{1})) || (!any (size (varargin{1}, 2) == [2,3])))
error ("scatteredInterpolant: Point input must be a 2 or 3 column array.");
endif
## already verified second input is column vector numeric input.
this.Points = varargin{1};
this.Values = varargin{2};
this = check_points (this);
this = setTriangulation (this);
case {3, 4}
## already verified numeric inputs are column vectors.
this.Points = [varargin{1:num_input_count-1}];
this.Values = varargin{num_input_count};
this = check_points (this);
this = setTriangulation (this);
otherwise
print_usage (this);
endswitch
switch char_input_count
case 0
# Do nothing, Method defaults already set when this created
# this.Method = "linear"; this.ExtrapolationMethod = "linear";
case 1
this.Method = tolower (varargin{num_input_count + 1});
switch this.Method
case {"linear", "natural"}
# Do nothing. extrap method default 'linear' set when this created.
case {"nearest"}
this.ExtrapolationMethod = "nearest";
otherwise
error ("scatteredInterpolant: Invalid METHOD '%s'", this.Method);
endswitch
case 2
this.Method = tolower (varargin{num_input_count + 1});
this.ExtrapolationMethod = tolower (varargin{num_input_count + 2});
#verify both methods are valid.
if !any (strcmp (this.Method, {"linear", "nearest", "natural"}))
error ("scatteredInterpolant: Invalid METHOD '%s'", this.Method);
elseif !any (strcmp (this.ExtrapolationMethod, {"linear", "nearest", "none"}))
error ("scatteredInterpolant: Invalid EXTRAPOLATIONMETHOD '%s'", ...
this.ExtrapolationMethod);
endif
otherwise
print_usage ();
endswitch
endif
endfunction
function v = subsref (this, S)
## subsref either returns a property value, or if () does the actual interpolation
## performing interpolation should return errors or warnings if in a bad state
#issue warning if Point status flags set
this.dimension = columns (this.Points); ## 0 (empty) or 2/3 for 2D/3D
if (! this.valid_tri)
if (! this.enough_points)
warning("scatteredInterpolant: not enough points to create a triangulation\n");
else
warning("scatteredInterpolant: cannot calculate triangulation from given points\n");
endif
endif
for S_idx = 1: numel(S)
if (S_idx == 1)
switch S(1).type(1)
case "("
if (isempty (this.Points) || (! this.valid_tri)) ...
|| (any (cellfun (@isempty, S(1).subs))) ...
|| isempty (S(1).subs)
v = [];
elseif (! this.valid_points_vals)
## unequal number of Points and values, throw error
error ("scatteredInterpolant: unequal number of points and values, cannot interpolate")
else
## Query points input validation
## interp output will be a column vector. If input has any other
## size, flag must be changed to true and sz_output must be set
## for final reshape
reshape_flag = false;
num_query_elements = numel (S(1).subs);
switch num_query_elements
case 1
if isnumeric (S(1).subs{1})
## single numeric input array. must be 2D array.
## columns must match dim.
## can allow a col vector to be accepted, handled as a
## single point by switching to row vector
qpts = S(1).subs{1};
if (ndims (qpts) > 2)
error ("scatteredInterpolant: query points must be 2D vectors or arrays");
endif
## if vector, ensure row vector
if isvector (qpts)
qpts = qpts(:).';
endif
#check for correct dimensionality
if (columns (qpts) != this.dimension)
error ("scatteredInterpolant: query points dimension must match interpolant");
endif
elseif iscell (S(1).subs{1})
## cell inputs must the same number of vectors as
## dimension. process as grid vectors to build query point
## array.
keyboard
if (! all (cellfun (@isnumeric, S(1).subs{1})))
error ("scatteredInterpolant: query grid vectors must be numeric");
elseif numel (S(1).subs{1}) != this.dimension
error ("scatteredInterpolant: query grid vector count must match interpolant dimension");
elseif (! all (cellfun (@isvector, S(1).subs{1})))
error ("scatteredInterpolant: query grid vectors must be row or column vectors");
endif
## extract grid vectors and produce point array
## vector orientation doesn't matter for ndgrid
switch this.dimension
case 2
[qpts_x, qpts_y] = ndgrid (S(1).subs{:}{:});
qpts = [qpts_x(:), qpts_y(:)];
case 3
[qpts_x, qpts_y, qpts_z] = ndgrid (S(1).subs{:}{:});
qpts = [qpts_x(:), qpts_y(:), qpts_z(:)];
endswitch
sz_output = size(qpts_x);
reshape_flag = true;
else
print_query_points_usage (this);
endif
case {2,3}
## all query inputs need to be numeric vectors or arrays
if ! all (cellfun (@isnumeric, S(1).subs))
print_query_points_usage (this);
elseif (num_query_elements != this.dimension)
error ("scatteredInterpolant: query points dimension must match interpolant");
endif
## check for vectors to be equal length
if ((all (cellfun (@isvector, S(1).subs))) ...
&& (! all (isequal (cellfun (@numel, S(1).subs, ...
"UniformOutput", false){:}))))
error ("scatteredInterpolant: query point vectors must have equal length");
## and nd arrays to be equal size
elseif ! isequal (cellfun (@size, S(1).subs, "UniformOutput", false){:})
error ("scatteredInterpolant: query point inputs must have equal size");
endif
## set output size based on first input element
if (! iscolumn (S(1).subs{1}))
sz_output = size(S(1).subs{1});
reshape_flag = true;
endif
switch this.dimension
case 2
qpts = [S(1).subs{1}(:), S(1).subs{2}(:)];
case 3
qpts = [S(1).subs{1}(:), S(1).subs{2}(:), S(1).subs{3}(:)];
endswitch
otherwise
## must be 1,2, or 3 input elements. call query usage error
print_query_points_usage (this);
endswitch
## tsearchn outputs vector of containing simplex, or NaN for Outside point
## plus barycentric coordinates of point within the simplex
## vi = [barycentric_coord].[local.Values]
v = NaN (rows(qpts), 1);
##Perform interpolation and extrapolation according to method
## note - Matlab provides little indication of the underlying
## linear extrapolation method, compatibility unsure.
switch this.Method
case "nearest"
## add an IF to even check extrap
switch this.ExtrapolationMethod
case "none"
nearest_pt_idx = dsearchn (this.Points, this.tri, qpts, NaN);
v(! isnan (nearest_pt_idx)) = ...
this.Values(nearest_pt_idx(! isnan (nearest_pt_idx)));
case "nearest"
nearest_pt_idx = dsearchn (this.Points, this.tri, qpts);
v = this.Values(nearest_pt_idx);
case "linear"
outside_qpt_idx = isnan (tsearchn (this.Points, this.tri, qpts));
nearest_pt_idx = dsearchn (this.Points, this.tri, qpts);
v(! outside_qpt_idx) = ...
this.Values(nearest_pt_idx(! outside_qpt_idx));
if any (outside_qpt_idx)
error ('not finished')
endif
otherwise
error ("scatteredInterpolant: invalid EXTRAPOLATIONMETHOD %s", this.ExtrapolationMethod);
endswitch
case "linear"
[tri_with_qpts, qpt_bary_coords] = tsearchn (this.Points, this.tri, qpts);
inside_qpt_idx = ! isnan (tri_with_qpts);
qpt_vertices = this.tri(tri_with_qpts(inside_qpt_idx), :);
if (isvector (qpt_vertices))
#ensure vertex values in row vector
qpt_vertex_values = this.Values(qpt_vertices)(:)';
else
qpt_vertex_values = this.Values(qpt_vertices);
endif
v(inside_qpt_idx) = dot (qpt_vertex_values, qpt_bary_coords(inside_qpt_idx,:), 2);
switch this.ExtrapolationMethod
case "none"
#nothing to do. NaN already pre-filled for outside points.
case "nearest"
outside_qpt_idx = ! inside_qpt_idx;
nearest_pt_idx = dsearchn (this.Points, this.tri, qpts(outside_qpt_idx,:));
v (outside_qpt_idx) = this.Values(nearest_pt_idx);
case "linear"
outside_qpt_idx = isnan (tsearchn (this.Points, this.tri, qpts));
##check for external pts
if any (outside_qpt_idx)
## find external triangulation points (repeats last pt)
boundary_pts = convhull (this.Points);
switch this.dimension
case 2
## find visible boundary points for each qpt
## return cell array of col vectors for each qpt
boundary_pts_vis = find_convex_vis_pts (this, this.Points(boundary_pts(1:end-1),:), qpts(outside_qpt_idx,:));
boundary_pts_count = cellfun ('numel', boundary_pts_vis);
all_boundary_vis_pts = sort (boundary_pts(unique (cell2mat (boundary_pts_vis'))'));
[boundary_pts_vis_idx_to_all, ~] = arrayfun(@(A) find(all_boundary_pts_vis==boundary_pts_vis{A}),[1:numel(boundary_pts_vis)]', "UniformOutput",false);
## compute gradients at visible boundary pts
## find points connecting to each boundary pt
[boundary_tris ~] = arrayfun (@(A) find (this.tri==A), all_boundary_vis_pts, "UniformOutput", false);
boundary_connections = cellfun (@(A) unique (this.tri(A,:)), boundary_tris, "UniformOutput", false);
## for each boundary point, make gradient system of equations
boundary_gradeqns = arrayfun(@(A) [this.Points(boundary_connections{A}(boundary_connections{A}~=A),:) - this.Points(all_boundary_vis_pts(A),:), ...
this.Values(boundary_connections{A}(boundary_connections{A}~=A)) - this.Values(all_boundary_vis_pts(A))], ...
[1:numel(all_boundary_vis_pts)]', "UniformOutput", false);
## solve overdefined eqn system for [dQdx, dQdy] at each boundary pt (each row n is for boundary pt n)
boundary_pt_grads = cell2mat(cellfun (@(A) mldivide(A(:,1:end-1), A(:,end)), boundary_gradeqns', "UniformOutput",false))';
## use visible boundary point gradients to extrapolate value at each point
##do 2-point visibility
v(outside_qpt_idx)(boundary_pts_count==2)
#do > 2-pt visibility
## split boundary into 'visibility regions' projected from boundary edges.
## perform linear extrapolation based on that region and visible gradients
case 3 ##3D
endswitch
endif
otherwise
error ("scatteredInterpolant: invalid EXTRAPOLATIONMETHOD %s", this.ExtrapolationMethod);
endswitch
case "natural"
## error ("scatteredInterpolant: 'natural' interpolation method not yet implemented")
##currently requires matgeom package
if (! pkg('list','matgeom'){1}.loaded)
error ("scatteredInterpolant: 'NATURAL' method requires 'matgeom' package be installed and loaded");
endif
#identify point/tri associations and index of points inside convex hull
[tri_with_qpts, qpt_bary_coords] = tsearchn (this.Points, this.tri, qpts);
inside_qpt_idx = !isnan (tri_with_qpts);
switch this.dimension
case 2
## find circumcircles of delaunay tri (center & rad)
## which are also the voronoi facet intersection points.
## [v_facet_pts, v_facet_list] = voronoin (this.Points);
## tri_circ_radii = distancePoints (v_facet_pts(2:end,:), this.Points(this.tri(:,1),:), "diag");
## ##cirs/spheres containing qp are it's nn points
#find centers
## use circumcenters instead of voronoi to keep data
## matched to delaunay tri's.
tri_circ_centers = circumCenter (this.Points(this.tri(:,1),:),this.Points(this.tri(:,2),:),this.Points(this.tri(:,3),:));
tri_circ_radii = distancePoints (tri_circ_centers, this.Points(this.tri(:,1),:), "diag");
## find which circles cover which qpts. those qpts
## natural neighbors will be all of the points of
## the triangles in those circles.
in_which_circle = distancePoints (tri_circ_centers, qpts(inside_qpt_idx,:)) <= tri_circ_radii; #rows are tri_circs, cols are for each qpt
## Get the interpolation points that are nat neighbors
## for each qpt. likely unequal lengths, so store in
## cells.
qpt_nat_neighbors_pts = arrayfun(...
@(x) unique(this.tri(in_which_circle(:,x),:))(:),...
[1:columns(in_which_circle)],"UniformOutput",false);
switch this.ExtrapolationMethod
case "none"
case "nearest"
case "linear"
otherwise
error ("scatteredInterpolant: invalid EXTRAPOLATIONMETHOD %s", this.ExtrapolationMethod);
endswitch
case 3
switch this.ExtrapolationMethod
case "none"
case "nearest"
case "linear"
otherwise
error ("scatteredInterpolant: invalid EXTRAPOLATIONMETHOD %s", this.ExtrapolationMethod);
endswitch
endswitch
endswitch
if (reshape_flag)
v = reshape (v, sz_output);
endif
endif
case {".", "{"}
v = builtin ("subsref", this, S(1));
otherwise
error ("scatteredInterpolant: invalid scatteredInterpolant index type '%s'", S(1).type);
endswitch
else
v = builtin ("subsref", v, S(S_idx));
endif
endfor
endfunction
function this = subsasgn (this, S, val)
# there are only four public properties (even if empty) at object
# creation. They can be updated via subasgn, but not created, and new
# properties cannot be added. Points and values may be individually updated
# to have a count mismatch without error. Error will be produced on attempted
# interpolation.
#
# subsasgn should error on attempt to add new property, add invalid method, or
# add poorly formed points/values. Points/values can be updated individually.
# so numrows doesn't have to match on assignment. But the object should
# track whether or not the points/values are in a matched state.
#
# Points - can be empty, must be a 2 x n (2D) or 3 x n (3D) double array
# Values - can be empty, must be a 1 x m double vector. if not column, store as column.
# Method - must be char/string: 'nearest', 'linear', or 'natural' (can't be empty)
# ExtrapolationMethod - must be char/string: 'none', 'linear', 'nearest' (can't be empty)
#
# Updating either points or values should recompute the triangulation.
update_triangulation = false;
update_point_value_state = false;
switch S(1).type(1)
case "."
if (numel(S) == 1)
#if only 1 level to S, validate input values/forms before assignment
switch (S(1).subs)
case "Points"
if (!isnumeric (val))
error ("scatteredInterpolant: Points must be numeric");
elseif ((!isempty(val)) && ((!ismatrix (val)) || (!any (size (val, 2) == [2,3]))))
error ("scatteredInterpolant: Points input must be a 2 or 3 column array.");
endif
update_triangulation = true;
update_point_value_state = true;
case "Values"
if (!isnumeric (val))
error ("scatteredInterpolant: Values must be numeric");
elseif (!isempty (val))
if (!isvector (val))
error ("scatteredInterpolant: Values must be in vector form");
elseif (!iscolumn (val))
val = val(:);
endif
endif
update_point_value_state = true;
case "Method"
if (!ischar (val))
error ("scatteredInterpolant: METHOD input must be a string");
elseif (!any (strcmp (val, {"linear", "nearest", "natural"})))
error ("scatteredInterpolant: invalid METHOD '%s'", ...
val);
endif
case "ExtrapolationMethod"
if (!ischar (val))
error ("scatteredInterpolant: EXTRAPOLATIONMETHOD input must be a string");
elseif (!any (strcmp (val, {"linear","nearest","none"})))
error ("scatteredInterpolant: invalid EXTRAPOLATIONMETHOD '%s'", ...
val);
endif
otherwise
error ("scatteredInterpolant: invalid property '%s'", ...
S(1).subs);
endswitch
this = builtin ("subsasgn", this, S, val);
elseif (numel(S) == 2)
## if 2nd index level, must be (), check value before calling builtin
switch S(2).type(1)
case "("
switch (S(1).subs)
case "Points"
if (!isnumeric (val))
error ("scatteredInterpolant: Points must be numeric");
elseif (isempty (val) || isempty (this.Points))
# empty values can cause invalid array shape
testval = builtin ("subsasgn", this.Points, S(2), val);
if (!any (size (testval, 2) == [2,3]))
error ("scatteredInterpolant: invalid Points shape");
endif
endif
update_triangulation = true;
update_point_value_state = true;
case "Values"
if (!isnumeric (val))
error ("scatteredInterpolant: Values must be numeric");
elseif ((isempty (val) || isempty (this.Values)))
# empty values can cause invalid array shape
testval = builtin ("subsasgn", this.Values, S(2), val);
if (!isvector (testval))
error ("scatteredInterpolant: invalid Values shape");
elseif (isempty (this.Values) && !isempty (val) && !iscolumn (val))
val = val(:);
endif
endif
update_point_value_state = true;
case "Method"
if (!ischar (val))
error ("scatteredInterpolant: METHOD input must be a string");
endif
#test substitution, error if it results in a bad value
testval = builtin ("subsasgn", this.Method, S(2), val);
if (!any (strcmp (testval, {"linear", "nearest", "natural"})))
error ("scatteredInterpolant: invalid METHOD '%s'", ...
testval);
endif
case "ExtrapolationMethod"
if (!ischar (val))
error ("scatteredInterpolant: EXTRAPOLATIONMETHOD input must be a string");
endif
#test substitution, error if it results in a bad value
testval = builtin ("subsasgn", this.ExtrapolationMethod, S(2), val);
if (!any (strcmp (val, {"linear","nearest","none"})))
error ("scatteredInterpolant: invalid EXTRAPOLATIONMETHOD '%s'", ...
testval);
endif
otherwise
error ("scatteredInterpolant: invalid property '%s'", ...
S(1).subs);
endswitch
this = builtin ("subsasgn", this, S, val);
case "."
error ("scatteredInterpolant: scatteredInterpolant properties cannot be subindexed by '.'");
case "{"
error ("scatteredInterpolant: scatteredInterpolant cannot be indexed by {}");
otherwise
endswitch
else
## no valid function for levels greater than 2
error ("scatteredInterpolant: scatteredInterpolant assignment depth > 2 undefined");
endif
if (update_triangulation);
## update triangulation only if Points changed
this = setTriangulation (this);
endif
if (update_point_value_state)
##update p/v state for points or values change
this = check_points (this);
endif
case "("
# () indexing reserved for interpolation. No object array definition,
# cannot be used for array assignment.
error ("scatteredInterpolant: scatteredInterpolant array value assignment undefined");
case "{"
error ("scatteredInterpolant: scatteredInterpolant cannot be indexed by {}");
otherwise
error ("scatteredInterpolant: invalid scatteredInterpolant index type '%s'", S(1).type);
endswitch
endfunction
function disp (this)
if (nargout > 0)
error ("scatteredInterpolant: output assignment not defined");
elseif (nargin != 1)
error ("scatteredInterpolant: only defined for one input.");
elseif !strcmp (class (this), "scatteredInterpolant")
error ("scatteredInterpolant: only defined for scatteredInterpolant objects.")
endif
printf('%s object with properties:\n\n', class (this));
printf(' Points = %dx%d matrix\n', rows(this.Points), columns(this.Points));
printf(' Values = %dx%d matrix\n', rows(this.Values), columns(this.Values));
printf(' Method = %s\n', this.Method);
printf(' ExtrapolationMethod = %s\n', this.ExtrapolationMethod);
endfunction
endmethods
methods (Access = public, Hidden = false)
function this = setTriangulation (this)
##called by constructor and subsasgn whenever tri needs a recalc
if (isempty (this.Points))
this.tri = [];
else
#compute delaunay triangulation. if invalid, catch error but continue
try
##disable qhull warnings??
this.tri = delaunayn (this.Points);
catch
this.tri = [];
end_try_catch
endif
if (isempty (this.tri) || isequal(this.tri, 0))
this.valid_tri = false;
else
this.valid_tri = true;
endif
endfunction
function this = check_points (this)
this.dimension = columns (this.Points); ## 0 (empty) or 2/3 for 2D/3D
this.valid_points_vals = isequal (rows (this.Points), rows (this.Values));
this.enough_points = this.dimension && (rows (this.Points) > this.dimension);
endfunction
function visible_points = find_convex_vis_pts (this, pts, viewpts)
%% pts = pointset given as 2d array, 1 pt per row. [x1, y1, (z1); x2, y2 (z2); etc]
%% pts already assumed to be a convex hull, in CCW order, without repetition.
%% viewpt = viewpoint origin (campera location) as row vector [x y (z)]
#recast multiple viewpts along dim3
viewpts = permute (viewpts, [3 2 1]);
% shift all points so that viewpt is at the origin (simplifies vector math)
pts = pts - viewpts;
[num_pts, dim, num_vpts] = size(pts);
% find point center of mass to set reference vector
pts_com = mean (pts, 1);
if dim == 2
% 2D calculate sin(theta) for each angle relative to C.O.M. vector
% as defined, in 2D, neg angles are CCW w.r.t. com vector, pos angle are CW
% ( sin(theta) = A x B / |A||B| )
pt_sin_thetas = (pts(:,1,:).*pts_com(:,2,:) - pts(:,2,:).*pts_com(:,1,:)) ...
./ (sqrt (sumsq (pts, 2) .* sumsq (pts_com, 2)));
[minangle, min_theta_idx] = min (pt_sin_thetas, [], 1);
[maxangle, max_theta_idx] = max (pt_sin_thetas, [], 1);
## check for numerically colinear points. for a convex hull, max/min
## angle points can only have same angle as another point if
## viewpoint-point vector is colinear with pt-pt edge. in that case, need
## to select closer point, with shorter viewpoint-point distance.
[same_min_idx, same_min_vpt] = find (pt_sin_thetas == minangle);
[same_max_idx, same_max_vpt] = find (pt_sin_thetas == maxangle);
for idx = 1 : num_vpts
repeated = same_min_idx(same_min_vpt == idx);
if (numel (repeated) > 1)
[~, closer_min] = min (sumsq (pts(repeated,:,idx), 2));
min_theta_idx(:,:,idx) = repeated(closer_min);
minangle(:,:,idx) = pt_sin_thetas(min_theta_idx(:,:,idx));
endif
repeated = same_max_idx(same_max_vpt == idx);
if (numel (repeated) > 1)
[~, closer_max] = max (sumsq (pts(repeated,:,idx), 2));
max_theta_idx(:,:,idx) = repeated(closer_max);
maxangle(:,:,idx) = pt_sin_thetas(max_theta_idx(:,:,idx));
endif
endfor
## visible points are all those between min&max index. points from
## convex hull are in CCW order (cyclical), min-max angle always CW
## around viewpoint, which for external viewpoint is always CCW around
## convex hull. so points go from [min_theta_idx : max_theta_idx],
## including around the limits.
if any (min_theta_idx == max_theta_idx)
error ("only 1 visible point, shouldn't be possible for a 2D convex hull");
endif
simple_order = max_theta_idx > min_theta_idx;
visible_points = cell (num_vpts, 1);
visible_points (simple_order) = arrayfun (@colon, min_theta_idx(simple_order), max_theta_idx(simple_order), "UniformOutput", false);
visible_points (!simple_order) = arrayfun(@(A,B) sort([A:num_pts, 1:B]), min_theta_idx(!simple_order), max_theta_idx(!simple_order),"UniformOutput",false);
else
# 3D
error ("no 3D yet");
endif
endfunction
function print_usage (this)
## FIXME: overloading print_usage until bug ___ is fixed for classdefs.
## Can be removed and errors reverted back to 'print_usage ()' once that
## is fixed.
ident = "Octave:invalid-fun-call";
msg = sprintf(["Invalid call to scatteredInterpolant. Correct usage is:\n\n", ...
" -- F = scatteredInterpolant\n"...
" -- F = scatteredInterpolant (X, Y, Q)\n"...
" -- F = scatteredInterpolant (X, Y, Z, Q)\n"...
" -- F = scatteredInterpolant (P, Q)\n"...
" -- F = scatteredInterpolant (..., METHOD)\n"...
" -- F = scatteredInterpolant (..., METHOD, EXTRAPOLATIONMETHOD)\n\n"...
"Additional help for built-in functions and operators is\n"...
"available in the online version of the manual. Use the command\n"...
"'doc ' to search the manual index.\n\n"...
"Help and information about Octave is also available on the WWW\n"...
"at https://www.octave.org and via the help@octave.org\n"...
"mailing list.\n"]);
error (struct ("message", msg, "identifier", ident, "stack", dbstack (1)));
endfunction
function print_query_points_usage (this)
msg = sprintf(["scatteredInterpolant: invalid query points form. Correct usage is:\n\n", ...
" -- Si = f(xi, yi)\n", ...
" -- Si = f(xi, yi, zi)\n", ...
" -- Si = f(Pi)\n", ...
" -- Si = f({Xg, Yg})\n", ...
" -- Si = f({Xg, Yg, Zg})\n\n", ...
" See 'help scatteredInterpolant' for more information."]);
error (struct ("message", msg, "identifier", "", "stack", dbstack (1)));
endfunction
endmethods
endclassdef
##TEST CONSTRUCTOR
%!test # empty object with defaults
%! A = scatteredInterpolant ();
%! assert (class (A), "scatteredInterpolant");
%! assert ({A.Points, A.Values, A.Method, A.ExtrapolationMethod}, {[], [], 'linear', 'linear'});
%!test # simple object, class, verify value orientation independence
%! pts = [magic(3); 2*magic(3)];
%! A = scatteredInterpolant (pts, [1:6]');
%! B = scatteredInterpolant (pts, [1:6]);
%! assert ({class(A), class(B)}, {"scatteredInterpolant", "scatteredInterpolant"});
%! assert ({A.Points, A.Values, A.Method, A.ExtrapolationMethod}, {pts, [1:6]', 'linear', 'linear'});
%! assert ({B.Points, B.Values, B.Method, B.ExtrapolationMethod}, {pts, [1:6]', 'linear', 'linear'});
%!test # check that objects unable to produce triangulations still produce object
%! pts = magic (3);
%! A = scatteredInterpolant (pts, [1:3]');
%! warning ("off");
%! assert (class(A), "scatteredInterpolant");
%! assert ({A.Points, A.Values, A.Method, A.ExtrapolationMethod}, {pts, [1:3]', 'linear', 'linear'});
%! warning ("on");
%!test # 2D vector orientation checks
%! A = scatteredInterpolant ([1 2 3 4 5]', [5 1 3 2 4]', [1:5]');
%! B = scatteredInterpolant ([1 2 3 4 5], [5 1 3 2 4]', [1:5]');
%! C = scatteredInterpolant ([1 2 3 4 5]', [5 1 3 2 4], [1:5]');
%! D = scatteredInterpolant ([1 2 3 4 5]', [5 1 3 2 4]', [1:5]);
%! E = scatteredInterpolant ([1 2 3 4 5]', [5 1 3 2 4], [1:5]);
%! F = scatteredInterpolant ([1 2 3 4 5], [5 1 3 2 4]', [1:5]);
%! G = scatteredInterpolant ([1 5; 2 1; 3 3; 4 2; 5 4], [1:5]');
%! H = scatteredInterpolant ([1 5; 2 1; 3 3; 4 2; 5 4], [1:5]);
%! assert (class (A), "scatteredInterpolant");
%! assert ({A.Points, A.Values, A.Method, A.ExtrapolationMethod}, {[1 2 3 4 5; 5 1 3 2 4]', [1:5]', 'linear', 'linear'});
%! assert (isequal (A.Points, B.Points, C.Points, D.Points, E.Points, F.Points, G.Points, H.Points));
%! assert (isequal (A.Values, B.Values, C.Values, D.Values, E.Values, F.Values, G.Values, H.Values));
%!test # 3D vector orientation checks
%! A = scatteredInterpolant ([1 2 3 4 5]', [5 1 3 2 4]', [2 1 5 4 3]', [1:5]');
%! B = scatteredInterpolant ([1 2 3 4 5], [5 1 3 2 4]', [2 1 5 4 3]', [1:5]');
%! C = scatteredInterpolant ([1 2 3 4 5]', [5 1 3 2 4], [2 1 5 4 3]', [1:5]');
%! D = scatteredInterpolant ([1 2 3 4 5]', [5 1 3 2 4]', [2 1 5 4 3], [1:5]');
%! E = scatteredInterpolant ([1 2 3 4 5]', [5 1 3 2 4]', [2 1 5 4 3]', [1:5]);
%! F = scatteredInterpolant ([1 2 3 4 5]', [5 1 3 2 4], [2 1 5 4 3], [1:5]);
%! G = scatteredInterpolant ([1 2 3 4 5]', [5 1 3 2 4]', [2 1 5 4 3], [1:5]);
%! H = scatteredInterpolant ([1 5 2; 2 1 1; 3 3 5; 4 2 4; 5 4 3], [1:5]');
%! I = scatteredInterpolant ([1 5 2; 2 1 1; 3 3 5; 4 2 4; 5 4 3], [1:5]);
%! assert (class (A), "scatteredInterpolant");
%! assert ({A.Points, A.Values, A.Method, A.ExtrapolationMethod}, {[1 5 2; 2 1 1; 3 3 5; 4 2 4; 5 4 3], [1:5]', 'linear', 'linear'});
%! assert (isequal (A.Points, B.Points, C.Points, D.Points, E.Points, F.Points, G.Points, H.Points, I.Points));
%! assert (isequal (A.Values, B.Values, C.Values, D.Values, E.Values, F.Values, G.Values, H.Values, I.Values));
%!test # method input checks
%! A = scatteredInterpolant ([1 5; 2 1; 3 3; 4 2; 5 4], [1:5]', "nearest");
%! B = scatteredInterpolant ([1 5; 2 1; 3 3; 4 2; 5 4], [1:5]', "nearest", "none");
%! C = scatteredInterpolant ([1 5; 2 1; 3 3; 4 2; 5 4], [1:5]', "linear", "none");
%! D = scatteredInterpolant ([1 2 3 4 5]', [5 1 3 2 4]', [2 1 5 4 3], [1:5]', "nearest");
%! E = scatteredInterpolant ([1 2 3 4 5]', [5 1 3 2 4]', [2 1 5 4 3], [1:5]', "nearest", "none");
%! F = scatteredInterpolant ([1 2 3 4 5]', [5 1 3 2 4]', [2 1 5 4 3], [1:5]', "linear", "none");
%! assert ({A.Method, A.ExtrapolationMethod}, {"nearest", "nearest"});
%! assert ({B.Method, B.ExtrapolationMethod}, {"nearest", "none"});
%! assert ({C.Method, C.ExtrapolationMethod}, {"linear", "none"});
%! assert ({D.Method, D.ExtrapolationMethod}, {"nearest", "nearest"});
%! assert ({E.Method, E.ExtrapolationMethod}, {"nearest", "none"});
%! assert ({F.Method, F.ExtrapolationMethod}, {"linear", "none"});
## Test input validation
%!error scatteredInterpolant (1)
%!error scatteredInterpolant (1, 2, 3, 4, 5, 6, 7)
%!error scatteredInterpolant (1, 2, 3, 4, 5, 6)
%!error scatteredInterpolant (1, 2, "abc", "def", "ghi")
%!error scatteredInterpolant ([], [])
%!error scatteredInterpolant ([], [1 2 3]')
%!error scatteredInterpolant (magic(3), [])
%!error scatteredInterpolant (1, 1, false, 1, "abc" , "def")
%!error scatteredInterpolant (1, 1, 1, "abc", 1, "def")
%!error scatteredInterpolant (magic (3), magic (3))
%!error scatteredInterpolant (magic (3), [1:4])
%!error scatteredInterpolant (magic (3), [1:2])
%!error scatteredInterpolant ([1:8]', [1:10]', [1:10]', [1:10]')
%!error scatteredInterpolant ([1:10]', [1:8]', [1:10]', [1:10]')
%!error scatteredInterpolant ([1:10]', [1:10]', [1:8]', [1:10]')
%!error scatteredInterpolant ([1:10]', [1:10]', [1:10]', [1:8]')
%!error scatteredInterpolant ([1:3]', [1:3]')
%!error scatteredInterpolant ([[1:4];[3:6];[4:7]], [1:3]')
%!error scatteredInterpolant (cat(3,magic(3),magic(3)), [1:3]')
%!error scatteredInterpolant (magic (3), [1:3], "foo")
%!error scatteredInterpolant (magic (3), [1:3], "foo", "none")
%!error scatteredInterpolant (magic (3), [1:3], "foo", "bar")
%!error scatteredInterpolant (magic (3), [1:3], "linear", "foo")
##TEST SUBSREF AND INTERPOLATION FUNCTION
%!assert (subsref (scatteredInterpolant (), substruct (".", "Points")), [])
%!assert (subsref (scatteredInterpolant (), substruct (".", "Values")), [])
%!assert (subsref (scatteredInterpolant (), substruct (".", "Method")), "linear")
%!assert (subsref (scatteredInterpolant (), substruct (".", "ExtrapolationMethod")), "linear")
%!assert (subsref (scatteredInterpolant (), substruct (".", "Method", "()", {2})), "i")
%!error subsref (scatteredInterpolant (), substruct (".", "blah"))
%!error subsref (scatteredInterpolant (), substruct ("{}", {1, 1}))
%!error subsref (scatteredInterpolant (), struct ("type", "/", "subs", 1))
%!test ## check point & triangulation errors
%! A = scatteredInterpolant ([magic(3); 2*magic(3)],[1:6]');
%! fail ("A(1,2)", "query points dimension");
%! A.Values = [1:5]';
%! fail ("A(1,2,3)", "unequal number");
%!test ##check triangulation warnings
%! A = scatteredInterpolant ([[1:5]',[1:5]'], [1:5]');
%! fail ("A(1,2)", "warning", "cannot calculate triangulation");
%!test ##check point sufficiency warnings
%! A = scatteredInterpolant (magic(3), [1:3]');
%! fail ("A(1,2,3)", "warning", "not enough points");
%! A = scatteredInterpolant (magic(2), [1:2]');
%! fail ("A(1,2)", "warning", "not enough points");
##interpolation input handling
%!test
%! A = scatteredInterpolant ([magic(3); 2*magic(3)], [1:6]);
%! qpts = cat (3, [1 2 3], [4 5 6]);
%! fail ("A(qpts)", "query points must be 2D vectors or arrays");
%! fail ("A([1,2])", "query points dimension must match");
%! fail ("A({'abc'})", "must be numeric");
%! fail ("A({[1,2]})", "query grid vector count must match");
%! fail ("A({[1,2,3]})", "query grid vector count must match");
%! fail ("A({[1,2],[3,4]})", "query grid vector count must match");
%! fail ("A({[1,2;3,4],1,1})", "query grid vectors must be row or column");
%! fail ("A(1,2)", "query points dimension must match interpolant");
%! fail ("A(1,2,[3 4])", "query point vectors must have equal length");
%! fail ("A([1 2 3],[1 2],[3 4 5])", "query point vectors must have equal length");
%! fail ("A([1 2 3]',[1 2]',[3 4 5]')", "query point vectors must have equal length");
%! fail ("A([1 2 3],[1 2],[3 4 5]')", "query point vectors must have equal length");
%! fail ("A([1,2;3,4],[1,2;3,4],[1,2])", "query point inputs must have equal size");
%! fail ("A(qpts,qpts,[1 2 3])", "query point inputs must have equal size");
%! fail ("A('foo')", "invalid query points form");
%! fail ("A({1 2 3},{1 2 3})", "invalid query points form");
%! fail ("A(1,2,3,4)", "invalid query points form");
##interpolation output tests
%!test
%! pts = [0 0 0; 0 0 1; 0 1 0; 1 0 0; 0 1 1; 1 0 1; 1 1 0; 1 1 1; 0.5 0.5 0.5];
%! A = scatteredInterpolant (pts, [1:9]', 'nearest', 'none');
%! qpts = magic(3)/10 -0.2;
%! assert (A(qpts), [NaN, 9, 3]')
%! A.ExtrapolationMethod = "nearest";
%! assert (A(qpts), [4, 9, 3]')
##TEST DISP
%!error disp (scatteredInterpolant (), 1)
%!error