summaryrefslogtreecommitdiff
path: root/inst/private
diff options
context:
space:
mode:
Diffstat (limited to 'inst/private')
-rw-r--r--inst/private/__dbl2int64__.m47
-rw-r--r--inst/private/clipplg.m207
-rw-r--r--inst/private/clippln.m178
-rw-r--r--inst/private/spheres_radius.m2
4 files changed, 70 insertions, 364 deletions
diff --git a/inst/private/__dbl2int64__.m b/inst/private/__dbl2int64__.m
new file mode 100644
index 0000000..0fa8478
--- /dev/null
+++ b/inst/private/__dbl2int64__.m
@@ -0,0 +1,47 @@
+## Copyright (C) 2017-2020 Philip Nienhuis
+##
+## This program 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.
+##
+## This program 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 this program. If not, see <http://www.gnu.org/licenses/>.
+
+function [opol, xy_mean, xy_magn] = __dbl2int64__ (inpoly, clippoly=[], xy_mean=[], xy_magn=[])
+
+ if (isempty (clippoly))
+ clippoly = zeros (0, size (inpoly, 2));
+ endif
+ if (isempty (xy_mean))
+ ## Convert & scale to int64 (that's what Clipper works with)
+ ## Find (X, Y) translation
+ xy_mean = mean ([inpoly; clippoly] (isfinite ([inpoly; clippoly](:, 1)), :));
+ ## Find (X, Y) magnitude
+ xy_magn = max ([inpoly; clippoly] (isfinite ([inpoly; clippoly](:, 1)), :)) ...
+ - min ([inpoly; clippoly] (isfinite ([inpoly; clippoly](:, 1)), :));
+ ## Apply (X,Y) multiplication (floor (log10 (intmax ("int64"))) = 18)
+ xy_magn = 10^(17 - ceil (max (log10 (xy_magn))));
+ endif
+
+ ## Scale inpoly coordinates to optimally use int64
+ inpoly(:, 1:2) -= xy_mean;
+ inpoly *= xy_magn;
+
+ idin = [ 0 find(isnan (inpoly (:, 1)))' numel(inpoly (:, 1))+1 ];
+ ## Provisional preallocation. npolx is average nr. of vertices per polygon
+ npolx = fix (size (inpoly, 1) / numel (idin)-1);
+ npoly = size (inpoly, 2);
+ opol = repmat (struct ("x", zeros (npolx, npoly), "y", zeros (npolx, npoly)), ...
+ 1, numel(idin) - 1);
+ for ii=1:numel (idin) - 1
+ opol(ii).x = int64 (inpoly(idin(ii)+1:idin(ii+1)-1, 1));
+ opol(ii).y = int64 (inpoly(idin(ii)+1:idin(ii+1)-1, 2));
+ endfor
+
+endfunction
diff --git a/inst/private/clipplg.m b/inst/private/clipplg.m
deleted file mode 100644
index 7c523aa..0000000
--- a/inst/private/clipplg.m
+++ /dev/null
@@ -1,207 +0,0 @@
-## Copyright (C) 2014 Philip Nienhuis
-##
-## This program 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.
-##
-## This program 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 this program. If not, see <http://www.gnu.org/licenses/>.
-
-## -*- texinfo -*-
-## @deftypefn {Function File} [@var{val}, @var{npts}, @var{pptr}] = clipplg (@var{val}, @var{npts}, @var{pptr}, @var{sbox}[, @var{styp]})
-## Undocumented internal function for clipping polygons & interpolating Z & M
-## values.
-##
-## @seealso{}
-## @end deftypefn
-
-## Author: Philip Nienhuis <prnienhuis@users.sf.net>
-## Created: 2014-11-18
-## Updates:
-## 2014-11-22 Create stanza for finding back vertices and interpolating Z & M
-## 2014-11-25 Fix interpolating M & Z values, properly treat new corner points
-## of clipped polygons
-## 2014-11-26 Update tnpt and tnpr arrays
-## 2015-07-10 Provision for segments crossed twice by clipping polygon
-
-function [val, tnpt, tnpr] = clipplg (val, tnpt, tnpr, sbox, styp=5)
-
- ## Indices to start of each subfeature, plus end+1
- tnprt = [(tnpr + 1) tnpt+1];
-
- ## Initialize total number of clipped vertices
- tnptclp = 0;
- tnprclp = cell (numel (tnpr, 1));
-
- for kk=1:numel (tnpr)
- ## Work from end back to start of subfeatures to avoid mixing up index arrays
- jj = numel (tnpr) - kk + 1;
- ## Select rows belonging to this partial feature. First save non-selected rows
- b_val = e_val = [];
- if (jj > 1)
- ## There's one or more subfeatures lower down
- b_val = val(1:tnprt(jj)-1, :);
- endif
- if (jj < numel (tnpr))
- ## There's one or more subfeatures higher up
- e_val = val(tnprt(jj+1):end, :);
- endif
- tval = val(tnprt(jj):tnprt(jj+1)-1, :);
- ## oc_polybool is in OF geometry package
- [X, Y, npol, b, c] = oc_polybool (tval(:, 1:2), sbox, 'AND');
-
- ## Initialize new number of points & new part pointers in clipped polygon(s)
- nptclp = 0;
- nprclp = 0;
- valn = [];
- if (npol)
- ## Make an XY matrix, remove NaNs on upper and lower row
- valc = [X Y](2:end-1, :);
- ## Augment NaNs for Z and M, and augment type + shape record index columns
- ncl = size (valc, 1);
- valc = [ valc NaN(ncl, 2) tval(1, 5)*ones(ncl, 1) tval(1, 6)*ones(ncl, 1) ];
- ## Pointers to subpolygons resulting from clipping
- ipt = find (isnan (valc(:, 1)))';
- ipt = [ 0 ipt (size (valc, 1) + 1) ];
- ## For each new polygon...
- for ipol=1:npol
- valn = valc(ipt(ipol)+1:ipt(ipol+1)-1, :);
- ## Update total number of points in clipped polygon(s)
- nptclp += size (valn, 1);
- tnptclp += size (valn, 1);
- ## Add a new 0-based pointer to next part
- nprclp = [ nprclp nptclp ];
- ## Compute all interdistances. distancePoints is in OF geometry package
- ## Avoid polygon end point ( = start point)
- dsts = distancePoints (valn(1:end-1, 1:2), tval(1:end-1, 1:2));
- ## Find matching points in sub and out polygon (row, col)
- [rw, cl] = ind2sub (size (dsts), find (abs (dsts) < eps));
- ## Transfer known Z and M-values
- valn(rw, 3:4) = tval(cl, 3:4);
- ## cl indices refer to original shape, rw indices to clipped shape
- if (numel (cl) >= 1)
- ## Separate polygon segments clipped, or vertex on bounding box side
- ## For each valn row coords not in tval, interpolate Z and M values
- im = setdiff ([1:size(valn, 1)-1], rw);
- ## mi equals cl filled with zeros for non-matches, to easen indexing
- mi = zeros (1, size (valn, 1) - 1);
- mi(rw) = cl;
- ## Find direction of polyline
- pdir = find (abs (diff (rw)) - 1 < eps);
- if (isempty (pdir))
- ## Single point within bounding box. Direction doesn't matter then
- drctn = 1;
- else
- drctn = sign (diff (rw([pdir pdir+1])))(1);
- endif
- for ii=1:numel (im)
- ## Get matching outer vertex. Below IF-ELSEIF order = critical to
- ## avoid index out-of-range errors
- if (im(ii) == 1)
- ## Clipped off outer vertex = previous in tval. diff(cl) = direction
- intpl = true;
- idx = mi(im(ii)+1) - drctn;
- ovtx = tval(idx, :);
- cvtx = tval(mi(im(ii)+1), :);
- elseif ((! ismember (im(ii)-1, rw)) && (! ismember (im(ii)+1, rw)))
- ## Probably a corner point. Just retain NaN values
- intpl = false;
- elseif (! ismember (im(ii)-1, rw))
- ## Clipped off outer vertex = previous in tval. diff(cl) = direction
- intpl = true;
- idx = mi(im(ii)+1) - drctn;
- ovtx = tval(idx, :);
- cvtx = tval(mi(im(ii)+1), :);
- elseif (! ismember (im(ii)+1, rw))
- ## Clipped off outer vertex = next in tval. diff(cl) = direction
- intpl = true;
- idx = mi(im(ii)-1) + drctn;
- ovtx = tval(idx, :);
- cvtx = tval(mi(im(ii)-1), :);
- endif
- ## Parent points found, now interpolate M and Z (if appropriate)
- if (intpl && styp > 5)
- ## Compute missing M and Z values. Invoke largest diff of X/Y coordinates
- difx = abs (cvtx(1) - ovtx(1));
- dify = abs (cvtx(2) - ovtx(2));
- if (difx > dify)
- ## X distance is greater
- fac = (valn(im(ii), 1) - cvtx(1)) / difx;
- else
- ## Y distance is greater
- fac = (valn(im(ii), 2) - cvtx(2)) / dify;
- endif
- fac = abs(fac);
- ## FIXME a debug stmt to detect wrong interpolation => wrong vertices
- if (fac > 1.0)
- printf ("Oops - fac > 1..\n");
- % keyboard
- endif
- if (isfinite (ovtx(3)))
- valn(im(ii), 3) = fac * (ovtx(3) - cvtx(3)) + cvtx(3); ## Z-value
- endif
- if (isfinite (ovtx(4)))
- valn(im(ii), 4) = fac * (ovtx(4) - cvtx(4)) + cvtx(4); ## M-value
- endif
- endif
- endfor
- ## Remove last nprclp entry and temporarily store it in a cell arr
- tnprclp(jj) = nprclp;
-
- elseif (numel (cl) == 0)
- ## One polygon segment clipped twice. Simply assign nearest Z & M values
- ## FIXME proper interpolation required
- ## Find points interpolated on segment(s); they're not in sbox
- [im, ix] = min (distancePoints (sbox(1:end-1, :), valn(1:end-1, 1:2)));
- im = find (im > 0);
- ix = ix(im);
- ## Find nearest polygon points (could be on another polygon segment !)
- [~, ix] = min (distancePoints (tval(1:end-1, 1:2), valn(im, 1:2)));
- ## Assign Z and M values
- valn(im, 3:4) = tval(ix, 3:4);
- ## Remove last nprclp entry and temporarily store it in a cell arr
- tnprclp(jj) = nprclp;
-
- endif
- ## Last row of polygon equals first
- valn(end, :) = valn(1, :);
- ## Augment new polygon after (yet untouched) previous polygons
- b_val = [ b_val; valn ];
-
- endfor ## clipped subpolygons
-
- else
- ## No intersection at all. Just drop tval
- % tnprclp = {};
- endif
-
- val = [b_val ; e_val];
- tnprt(jj+1:end) -= tnprt(jj+1) - tnprt(jj) - size (valn, 1);
- if (isempty (valn))
- ## This subfeature has no points in boundingbox +> drop from list
- tnprt(jj+1) = [];
- endif
-
- endfor
-
- ## Adapt & clean up npt
- tnpt = tnptclp;
- ## Adapt & clean up npr. Concatenate all pointers created by oc_polybool
- tnpr = [ tnprclp{1} ];
- for ii=2:numel (tnprclp)
- ## Skip empty entries
- if (! isempty (tnprclp{ii}))
- tnpr = [tnpr(1:end-1) (tnprclp{ii} + tnpr(end)) ];
- endif
- endfor
- if (! isempty (tnpr))
- tnpr(end) = [];
- endif
-
-endfunction
diff --git a/inst/private/clippln.m b/inst/private/clippln.m
index 7d8a575..cfcf32c 100644
--- a/inst/private/clippln.m
+++ b/inst/private/clippln.m
@@ -1,4 +1,4 @@
-## Copyright (C) 2014 Philip Nienhuis
+## Copyright (C) 2014-2020 Philip Nienhuis
##
## This program is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
@@ -16,7 +16,7 @@
## -*- texinfo -*-
## @deftypefn {Function File} [@var{val}, @var{npts}, @var{pprt}] = clipplg (@var{val}, @var{npts}, @var{pprt}, @var{sbox}, @var{styp})
## Undocumented internal function for clipping (poly)lines within a bounding
-## box and interpolating M and Z-values.
+## box and copying M and Z-values from nearest old vertex.
##
## @seealso{}
## @end deftypefn
@@ -24,162 +24,28 @@
## Author: Philip Nienhuis <prnienhuis@users.sf.net>
## Created: 2014-12-01
-function [val, tnpt, tnpr] = clippln (val, tnpt, tnpr, sbox, styp=3)
+function [valn, tnpt, tnpr] = clippln (val, tnpt, tnpr, sbox, styp=3)
- ## Backup
- btnpr = tnpr;
- ## Prepare augmented npr array for easier indexing into subfeatures
- ttnpr = [tnpr tnpt];
- ctnpr = {};
- ## Initialize output array
- valt = [];
-
- ## First step is to check all polyline vertices whether they lie in sbox
- [aa, bb] = inpolygon (val(:, 1), val(:, 2), sbox(:, 1), sbox(:, 2));
- if (any (aa) || any (bb))
- ## So there are. Morph bounding box into something that OF geometry likes:
- ## Parametric representation of bounding box. "createLine" is in OF geometry
- ssbox = createLine (sbox(1:4, :), sbox (2:5, :));
-## Create sbx from sbox ([minx max miny maxy])
-sbx = [min(sbox(:, 1)), max(sbox(:, 1)), min(sbox(:, 2)), max(sbox(:, 2))];
- ## For those points ON the boundary no interpolation is needed.
- ## First treat those points requiring interpolation (1st inside,
- ## 2nd outside or vice versa)
- for ll=1:numel (tnpr)
- ## Start at end of polyline-part to avoid mixing up tnpr pointers
- ii = numel (tnpr) - ll + 1;
- nprt = [];
- ## Select polyline part
- valn = val(ttnpr(ii)+1:ttnpr(ii+1), :);
- a = aa(ttnpr(ii)+1:ttnpr(ii+1));
- if (any (a))
- b = bb(ttnpr(ii)+1:ttnpr(ii+1));
- da = diff (a);
- idx = find (da);
- ## Find out which bounding box segments have been crossed where
- ## 1. Parametric representation of line
- lines = createLine (valn(idx, 1:2), valn(idx+1, 1:2));
- ## 2. Compute intersections. clipLine & distancePoints are in OF geometry
- for jj=1:numel (idx)
- ## Start at end of polyline-part to avoid mixing up tnpr pointers
- kk = numel (idx) - jj + 1;
- ## Skip points on border
- if (! b(idx(kk)+1))
- intsc = reshape (clipLine (lines(kk, :), sbx), 2, 2)';
- dst = distancePoints (intsc, valn(idx(kk):idx(kk)+1, 1:2));
- ## If segment goes out, take nearest point to valn(idx(kk)+1)
- if (da(idx(kk)) < 1)
- [~, ix] = min (dst(:, 2));
- else
- [~, ix] = min (dst(:, 1));
- endif
- fac = dst(ix, 2) / sum (dst(2, :));
- intsc = intsc (ix, :);
- ## Insert new intersection point
- valn(idx(kk)+2:end+1, :) = valn(idx(kk)+1:end, :);
- valn(idx(kk)+1, 1:2) = intsc(1:2);
- if (styp > 3)
- valn(idx(kk)+1, 3:4) = valn(idx(kk)+2, 3:4) + fac * ...
- (valn(idx(kk), 3:4) - valn(idx(kk)+2, 3:4));
- else
- valn(idx(kk)+1, 3:4) = [NaN NaN];
- endif
- ## Also update a, b, da and idx. Also mark first point outside sbox
- aa(ttnpr(ii)+idx(kk):ttnpr(ii)+idx(kk)+1) = 1;
- a = [ a(1:idx(kk)); 1; a(idx(kk)+1:end) ];
- b = [ b(1:idx(kk)); 1; b(idx(kk)+1:end) ];
- endif
- if (da(idx(kk)) < 0)
- ## Segment leaves bbox. Replace first outer point with NaN row
- valn(idx(kk)+2, :) = NaN (1, 6);
- ## Be sure to include NaN row
- a(idx(kk)+2) = 1;
- endif
- endfor
- ## Update valt, ttnpr
- valn = valn(find(a), :);
- if (! isempty (valn))
- ## Remove last NaN row if present
- if (all (isnan (valn(end, 1:2))))
- valn(end,:) = [];
- endif
- isn = 1;
- while (! isempty (isn))
- isn = find (isnan (valn(:, 1)));
- if (! isempty (isn))
- ## Insert new subfeature pointer. npart pointers are 0-based
- nprt = [ nprt isn(1) ];
- valn(isn(1), :) = [];
- endif
- endwhile
- ## Build up new points from top down
- valt = [ valn; valt ];
- ## nprt only gets non-empty if a polyline part was split up
- ## FIXME moet net als bij polygons
- ttnpr = [ ttnpr(1:ii) (nprt+ttnpr(ii)-1) ttnpr(ii+1:end) ];
- endif
- endif
- endfor
- tnpt = size (valt, 1);
- tnpr = ttnpr(1:end-1);
+ ## Clip (intersection)
+ if (mod (styp, 10) == 3)
+ [valn, nparts] = clipPolyline_clipper (val(:, 1:2), sbox, 1);
+ else
+ [valn, nparts] = clipPolygon_clipper (val(:, 1:2), sbox, 1);
endif
- ## There's still the possibility of segments crossing through Bounding Box
- ## w/o having points inside. Search these segments one by one
- cc = find (! aa);
- if (numel (cc) > 1)
- ## Find discontinuities; and merge wit tnpr (multipart) discontinuities
- ci = find (diff(cc) > 1)';
- ## Make sure to only include multipart discontinuities within range of cc
- ctnpr = unique ([max(btnpr, cc(1)-1) cc(ci)' min(cc(end), size (val, 1))]);
- ## Create sbx from sbox ([minx max miny maxy])
- sbx = [min(sbox(:, 1)), max(sbox(:, 1)), min(sbox(:, 2)), max(sbox(:, 2))];
- for ic=1:numel (ctnpr) - 1
- valc = val(ctnpr(ic)+1:ctnpr(ic+1), :);
- ## Eliminate segments that lie to one outside of the bounding box.
- ## Step 1: assess which side of the bounding box the segments lie
- f1 = valc(:, 1) < min (sbox(:, 1));
- f2 = valc(:, 1) > max (sbox(:, 1));
- f3 = valc(:, 2) < min (sbox(:, 2));
- f4 = valc(:, 2) > max (sbox(:, 2));
- ff = [f1'; f2' ; f3' ; f4']';
- ## Step 2: find those segments that cross > 2 extended sbox boundary line
- dd = find (sum ([abs(diff(f1)'); abs(diff(f2)'); ...
- abs(diff(f3)'); abs(diff(f4)')]) > 1);
- ## If there are any vertices changing orientation w.r.t. bounding box:
- if (! isempty (dd))
- ## Create array of lines potentially crossing sbox
- lines = [];
- for ii=1:numel (dd)
- lines = [ lines; createLine(valc(dd(ii), 1:2), valc(dd(ii)+1, 1:2))];
- endfor
- ## For each of that segment, try to compute intersections with sbox
- for ii=1:numel (dd)
- ## clipLine is in OF geometry
- jntsc = reshape (clipLine (lines, sbx), 2, 2)';
- ## Check if there are any intersections at all
- if (! any (isnan (jntsc)))
- valn = [ jntsc NaN(2, 4) ];
- ## Interpolate M and Z values. valc(dd(ii)+1) also required ...
- dist = distancePoints (valc(dd(ii):dd(ii)+1, 1:2), valn(:, 1:2));
- ## ... for length of original segment:
- sl = sum (dist(:, 1));
- ## Find closest point of original segment
- valn(1, 3:4) = valc(dd(ii), 3:4) + dist(1, 1) / sl * ...
- (valc(dd(ii)+1, 3:4) - valc(dd(ii), 3:4));
- valn(2, 3:4) = valc(dd(ii), 3:4) + dist(1, 2) / sl * ...
- (valc(dd(ii)+1, 3:4) - valc(dd(ii), 3:4));
- ## Adapt tnpr
- tnpr = [ tnpr size(valt, 1) ];
- valt = [ valt; valn ];
- endif
- endfor
- tnpr = unique (tnpr);
- endif
- endfor
- endif
-
- valt(:, 5:6) = repmat (val(1, 5:6), size (valt, 1), 1);
- val = valt;
+ ## Set up pointers to subpolygons
+ idn = [0 find(isnan (valn(:, 1)))' size(valn, 1)+1];
+ tnpr = idn(1:end-1);
+ ## Setup Z/M/nr/type columns
+ for ii=1:nparts
+ dists = distancePoints (val(:, 1:2), valn(idn(ii)+1:idn(ii+1)-1, 1:2));
+ [mind, idm] = min (dists);
+ ## Simply copy over Z & M values & no. and type of nearest "old" vertex
+ valn(idn(ii)+1:idn(ii+1)-1, 3:6) = val(idm, 3:6);
+ endfor
+ ## Fix up pointers and vertices
+ valn(idn(2:end-1), :) = [];
+ tnpr = tnpr - [0 : numel(tnpr)-1];
+ tnpt = size (valn, 1);
endfunction
diff --git a/inst/private/spheres_radius.m b/inst/private/spheres_radius.m
index 2e22da0..10f1778 100644
--- a/inst/private/spheres_radius.m
+++ b/inst/private/spheres_radius.m
@@ -1,4 +1,4 @@
-## Copyright (C) 2013 Carnë Draug <carandraug@octave.org>
+## Copyright (C) 2013-2020 Carnë Draug <carandraug@octave.org>
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-27 21:44:03 +0000