1 /*
2
3 Copyright (C) 2017 Andreas Stahel
4
5 This program is free software; you can redistribute it and/or modify it
6 under the terms of the GNU General Public License as published by the
7 Free Software Foundation; either version 3 of the License, or (at your
8 option) any later version.
9
10 This program is distributed in the hope that it will be useful, but WITHOUT
11 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
13 for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program. If not, see
17 <http://www.gnu.org/licenses/>.
18
19 */
20
21 #include <R.h>
22 #include <Rdefines.h>
23 #include <R_ext/Rdynload.h>
24 #include <Rinternals.h>
25 #include "Rgeometry.h"
26 #include "qhull_ra.h"
27
C_tsearchn(const SEXP dt,const SEXP p)28 SEXP C_tsearchn(const SEXP dt, const SEXP p)
29 {
30 int debug = 0;
31 /* Get the qh object from the delaunayn object */
32 SEXP ptr, tag;
33 qhT *qh;
34 tag = PROTECT(allocVector(STRSXP, 1));
35 SET_STRING_ELT(tag, 0, mkChar("delaunayn"));
36 ptr = PROTECT(getAttrib(dt, tag));
37 if (ptr == R_NilValue) {
38 error("Delaunay triangulation has no delaunayn attribute");
39 }
40 qh = R_ExternalPtrAddr(ptr);
41 UNPROTECT(2);
42
43 /* Check input matrix */
44 if(!isMatrix(p) || !isReal(p)){
45 error("Second argument should be a real matrix.");
46 }
47 unsigned int dim, n;
48 dim = ncols(p) + 1;
49 n = nrows(p);
50 if(dim <= 0 || n <= 0){
51 error("Invalid input matrix.");
52 }
53 if (dim != qh->hull_dim)
54 error("Invalid input matrix.");
55
56 /* Construct map from facet id to index */
57 facetT *facet;
58
59 /* Count the number of facets so we know how much space to
60 allocate in R */
61 int nf = 0; /* Number of facets */
62 int max_facet_id = 0;
63 int exitcode = 0;
64 FORALLfacets {
65 if (!facet->upperdelaunay) {
66 nf++;
67 if (facet->id > max_facet_id)
68 max_facet_id = facet->id;
69
70 /* Double check. Non-simplicial facets will cause segfault
71 below */
72 if (!facet->simplicial) {
73 Rprintf("Qhull returned non-simplicial facets -- try delaunayn with different options");
74 exitcode = 1;
75 break;
76 }
77 }
78 }
79
80 int *idmap = (int *) R_alloc(max_facet_id + 1, sizeof(int));
81 int i = 0;
82 FORALLfacets {
83 if (!facet->upperdelaunay) {
84 i++;
85 if (debug & 1) Rprintf("Facet id %d; index %d\n;", facet->id, i);
86 if (facet->id < 1 || facet->id > max_facet_id) {
87 Rf_error("facet_id %d (at index %d) is not in {1,...,%d}", facet->id, i, max_facet_id);
88 }
89 idmap[facet->id] = i;
90 }
91 }
92
93 /* Make space for output */
94 SEXP retlist, retnames; /* Return list and names */
95 int retlen = 2; /* Length of return list */
96 SEXP idx, points;
97 idx = PROTECT(allocVector(INTSXP, n));
98 int *iidx = INTEGER(idx);
99 points = PROTECT(allocMatrix(REALSXP, qh->num_points, dim - 1));
100
101 int j, k;
102
103 /* Output points */
104 pointT *point;
105 pointT *pointtemp;
106 if (debug & 2) Rprintf("%d POINTS\n", qh->num_points);
107 i = 0;
108 FORALLpoints {
109 for (k=0; k<(dim - 1); k++) {
110 REAL(points)[i+k*qh->num_points] = point[k];
111 if (debug & 2) Rprintf("%f ", point[k]);
112 }
113 i++;
114 if (debug & 2) Rprintf("\n");
115 }
116
117 /* Run through the matrix using qh_findbestfacet to determine
118 whether in hull or not */
119 boolT isoutside;
120 realT bestdist;
121 vertexT *vertex, **vertexp;
122
123 /* The name point is reserved for use with FORALLpoints */
124 coordT *testpoint;
125 testpoint = (coordT *) R_alloc(dim, sizeof(coordT));
126
127 for(i=0; i < n; i++) {
128 if (debug) Rprintf("\nTestpoint\n");
129 for(k=0; k < (dim - 1); k++) {
130 testpoint[k] = REAL(p)[i+n*k]; /* could have been pt_array = REAL(p) if p had been transposed */
131 if (debug) Rprintf(" %f", testpoint[k]);
132 }
133 if (debug) Rprintf("\n");
134 qh_setdelaunay(qh, dim, 1, testpoint);
135 facet = qh_findbestfacet(qh, testpoint, qh_ALL, &bestdist, &isoutside);
136 if (facet->tricoplanar) {
137 exitcode = 1;
138 break;
139 }
140 if (debug) Rprintf("Facet id %d; index %d\n", facet->id, idmap[facet->id]);
141 /* Convert facet id to id of triangle */
142 iidx[i] = idmap[facet->id];
143 /* /\* Return vertices of triangle *\/ */
144 j = 0;
145 FOREACHvertex_ (facet->vertices) {
146 for (j=0; j<dim - 1; j++) {
147 if (debug) Rprintf("%f ", vertex->point[j]);
148 }
149 if (debug) Rprintf("\n");
150 }
151
152 }
153
154 retlist = PROTECT(allocVector(VECSXP, retlen));
155 retnames = PROTECT(allocVector(VECSXP, retlen));
156 SET_VECTOR_ELT(retlist, 0, idx);
157 SET_VECTOR_ELT(retnames, 0, mkChar("idx"));
158 SET_VECTOR_ELT(retlist, 1, points);
159 SET_VECTOR_ELT(retnames, 1, mkChar("P"));
160 setAttrib(retlist, R_NamesSymbol, retnames);
161 UNPROTECT(4);
162
163 if (exitcode)
164 error("findDelaunay: not implemented for triangulated, non-simplicial Delaunay regions (tricoplanar facet, f%d).", facet->id);
165
166 return retlist;
167 }
168