Mercurial > octave-nkf

comparison src/DLD-FUNCTIONS/__voronoi__.cc @ 13871:adf60d6dc1dd

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
more compatibility fixes for __voronoi__ * __voronoi__.cc (F__voronoi__): Use v option for qhull, not Fv. Improve argument handling. Call qh_findgood_all to obtain number of Voronoi vertices directly. Correctly size output values. Use qh_pointid to place cell indices in the same order as the given points.
author John W. Eaton <jwe@octave.org>
date Wed, 16 Nov 2011 16:17:13 -0500
parents 6d7e133a4bed
children 440d7914cf01
comparison
equal deleted inserted replaced
13868:87f78c11d725 13871:adf60d6dc1dd
74 { 74 {
75 print_usage (); 75 print_usage ();
76 return retval; 76 return retval;
77 } 77 }
78 78
79 std::string options = ""; 79 std::string cmd = "qhull v";
80 80
81 if (nargin == 2) 81 if (nargin == 2 && ! args(1).is_empty ())
82 { 82 {
83 if (args(1).is_string ()) 83 if (args(1).is_string ())
84 options = args(1).string_value (); 84 cmd += " " + args(1).string_value ();
85 else if (args(1).is_empty ())
86 ; // Use default options
87 else if (args(1).is_cellstr ()) 85 else if (args(1).is_cellstr ())
88 { 86 {
89 options = "";
90 Array<std::string> tmp = args(1).cellstr_value (); 87 Array<std::string> tmp = args(1).cellstr_value ();
91 88
92 for (octave_idx_type i = 0; i < tmp.numel (); i++) 89 for (octave_idx_type i = 0; i < tmp.numel (); i++)
93 options += tmp(i) + " "; 90 cmd += " " + tmp(i);
94 } 91 }
95 else 92 else
96 { 93 {
97 error ("__voronoi__: OPTIONS argument must be a string, cell array of strings, or empty"); 94 error ("__voronoi__: OPTIONS argument must be a string, cell array of strings, or empty");
98 return retval; 95 return retval;
99 } 96 }
100 } 97 }
101 98
102 Matrix p (args(0).matrix_value ()); 99 Matrix points (args(0).matrix_value ());
103 const octave_idx_type dim = p.columns (); 100 const octave_idx_type dim = points.columns ();
104 const octave_idx_type np = p.rows (); 101 const octave_idx_type num_points = points.rows ();
105 102
106 p = p.transpose (); 103 points = points.transpose ();
107 double *pt_array = p.fortran_vec ();
108 104
109 boolT ismalloc = false; 105 boolT ismalloc = false;
110 106
111 // Replace the 0 pointer with stdout for debugging information 107 // Replace the 0 pointer with stdout for debugging information
112 FILE *outfile = 0; 108 FILE *outfile = 0;
113 FILE *errfile = stderr; 109 FILE *errfile = stderr;
114 110
115 // Qhull flags argument is not const char* 111 // Qhull flags and points arguments are not const...
116 OCTAVE_LOCAL_BUFFER (char, flags, 12 + options.length ()); 112
117 113 OCTAVE_LOCAL_BUFFER (char, cmd_str, cmd.length () + 1);
118 sprintf (flags, "qhull v Fv %s", options.c_str ()); 114
119 115 strcpy (cmd_str, cmd.c_str ());
120 int exitcode = qh_new_qhull (dim, np, pt_array, 116
121 ismalloc, flags, outfile, errfile); 117 int exitcode = qh_new_qhull (dim, num_points, points.fortran_vec (),
118 ismalloc, cmd_str, outfile, errfile);
122 if (! exitcode) 119 if (! exitcode)
123 { 120 {
121 // Calling findgood_all provides the number of Voronoi vertices
122 // (sets qh num_good).
123
124 qh_findgood_all (qh facet_list);
125
126 octave_idx_type num_voronoi_regions
127 = qh num_vertices - qh_setsize (qh del_vertices);
128
129 octave_idx_type num_voronoi_vertices = qh num_good;
130
131 // Find the voronoi centers for all facets.
132
133 qh_setvoronoi_all ();
134
124 facetT *facet; 135 facetT *facet;
125 vertexT *vertex; 136 vertexT *vertex;
126 137 octave_idx_type k;
127 octave_idx_type i = 0, n = 1, k = 0, m = 0, j = 0, r = 0; 138
128 139 // Find the number of Voronoi vertices for each Voronoi cell and
129 // Count number of vertices for size of NI. FIXME -- does Qhull 140 // store them in NI so we can use them later to set the dimensions
130 // have a way to query this value directly? 141 // of the RowVector objects used to collect them.
131 octave_idx_type nv = 0; 142
132 FORALLvertices 143 FORALLfacets
133 { 144 {
134 nv++; 145 facet->seen = false;
135 } 146 }
136 147
137 OCTAVE_LOCAL_BUFFER (octave_idx_type, ni, nv); 148 OCTAVE_LOCAL_BUFFER (octave_idx_type, ni, num_voronoi_regions);
138 149 for (octave_idx_type i = 0; i < num_voronoi_regions; i++)
139 for (i = 0; i < nv; i++)
140 ni[i] = 0; 150 ni[i] = 0;
141 151
142 qh_setvoronoi_all (); 152 k = 0;
143
144 bool infinity_seen = false;
145 facetT *neighbor, **neighborp;
146 coordT *voronoi_vertex;
147
148 FORALLfacets
149 {
150 facet->seen = false;
151 }
152 153
153 FORALLvertices 154 FORALLvertices
154 { 155 {
155 if (qh hull_dim == 3) 156 if (qh hull_dim == 3)
156 qh_order_vertexneighbors (vertex); 157 qh_order_vertexneighbors (vertex);
157 158
158 infinity_seen = false; 159 bool infinity_seen = false;
160
161 facetT *neighbor, **neighborp;
159 162
160 FOREACHneighbor_ (vertex) 163 FOREACHneighbor_ (vertex)
161 { 164 {
162 if (! neighbor->upperdelaunay) 165 if (neighbor->upperdelaunay)
163 { 166 {
164 if (! neighbor->seen) 167 if (! infinity_seen)
165 { 168 {
166 n++; 169 infinity_seen = true;
167 neighbor->seen = true; 170 ni[k]++;
168 } 171 }
169 ni[k]++;
170 } 172 }
171 else if (! infinity_seen) 173 else
172 { 174 {
173 infinity_seen = true; 175 neighbor->seen = true;
174 ni[k]++; 176 ni[k]++;
175 } 177 }
176 } 178 }
177 179
178 k++; 180 k++;
179 } 181 }
180 182
181 Matrix v (n, dim); 183 // If Qhull finds fewer regions than points, we will pad the end
184 // of the at_inf and C arrays so that they always contain at least
185 // as many elements as the given points array.
186
187 // FIXME -- is it possible (or does it make sense) for
188 // num_voronoi_regions to ever be larger than num_points?
189
190 octave_idx_type nr = (num_points > num_voronoi_regions
191 ? num_points : num_voronoi_regions);
192
193 boolMatrix at_inf (nr, 1, false);
194
195 // The list of Voronoi vertices. The first element is always
196 // Inf.
197 Matrix F (num_voronoi_vertices+1, dim);
198
182 for (octave_idx_type d = 0; d < dim; d++) 199 for (octave_idx_type d = 0; d < dim; d++)
183 v(0,d) = octave_Inf; 200 F(0,d) = octave_Inf;
184 201
185 boolMatrix AtInf (nv, 1, false); 202 // The cell array of vectors of indices into F that represent the
186 std::list<RowVector> F; 203 // vertices of the Voronoi regions (cells).
204
205 Cell C (nr, 1);
206
207 // Now loop through the list of vertices again and store the
208 // coordinates of the Voronoi vertices and the lists of indices
209 // for the cells.
210
211 FORALLfacets
212 {
213 facet->seen = false;
214 }
215
216 octave_idx_type i = 0;
187 k = 0; 217 k = 0;
188 i = 0;
189
190 FORALLfacets
191 {
192 facet->seen = false;
193 }
194 218
195 FORALLvertices 219 FORALLvertices
196 { 220 {
197 if (qh hull_dim == 3) 221 if (qh hull_dim == 3)
198 qh_order_vertexneighbors (vertex); 222 qh_order_vertexneighbors (vertex);
199 223
200 infinity_seen = false; 224 bool infinity_seen = false;
201 225
202 octave_idx_type n_vertices = ni[k++]; 226 octave_idx_type idx = qh_pointid (vertex->point);
203 227
204 // Qhull seems to sometimes produce "facets" with a single 228 octave_idx_type num_vertices = ni[k++];
205 // vertex. Is that a bug? How can a facet have just one 229
230 // Qhull seems to sometimes produces regions with a single
231 // vertex. Is that a bug? How can a region have just one
206 // vertex? Let's skip it. 232 // vertex? Let's skip it.
207 233
208 if (n_vertices == 1) 234 if (num_vertices == 1)
209 continue; 235 continue;
210 236
211 RowVector facet_list (n_vertices); 237 RowVector facet_list (num_vertices);
212 m = 0; 238
239 octave_idx_type m = 0;
240
241 facetT *neighbor, **neighborp;
213 242
214 FOREACHneighbor_(vertex) 243 FOREACHneighbor_(vertex)
215 { 244 {
216 if (neighbor->upperdelaunay) 245 if (neighbor->upperdelaunay)
217 { 246 {
218 if (! infinity_seen) 247 if (! infinity_seen)
219 { 248 {
220 infinity_seen = true; 249 infinity_seen = true;
221 facet_list(m++) = 1; 250 facet_list(m++) = 1;
222 AtInf(j) = true; 251 at_inf(idx) = true;
223 } 252 }
224 } 253 }
225 else 254 else
226 { 255 {
227 if (! neighbor->seen) 256 if (! neighbor->seen)
228 { 257 {
229 voronoi_vertex = neighbor->center;
230 i++; 258 i++;
231 for (octave_idx_type d = 0; d < dim; d++) 259 for (octave_idx_type d = 0; d < dim; d++)
232 { 260 F(i,d) = neighbor->center[d];
233 v(i,d) = *(voronoi_vertex++); 261
234 }
235 neighbor->seen = true; 262 neighbor->seen = true;
236 neighbor->visitid = i; 263 neighbor->visitid = i;
237 } 264 }
265
238 facet_list(m++) = neighbor->visitid + 1; 266 facet_list(m++) = neighbor->visitid + 1;
239 } 267 }
240 } 268 }
241 F.push_back (facet_list); 269
242 j++; 270 C(idx) = facet_list;
243 } 271 }
244 272
245 // For compatibility with Matlab, pad the cell array of vertex 273 retval(2) = at_inf;
246 // lists with empty matrices if there are fewer facets than
247 // points.
248 octave_idx_type f_len = F.size ();
249 Cell C(np > f_len ? np : f_len, 1);
250
251 i = 0;
252 for (std::list<RowVector>::const_iterator it = F.begin ();
253 it != F.end (); it++)
254 C.elem (i++) = *it;
255
256 F.clear ();
257
258 AtInf.resize (f_len, 1);
259 retval(2) = AtInf;
260 retval(1) = C; 274 retval(1) = C;
261 retval(0) = v; 275 retval(0) = F;
262 } 276 }
263 else 277 else
264 error ("__voronoi__: qhull failed"); 278 error ("__voronoi__: qhull failed");
265 279
266 // free memory from Qhull 280 // Free memory from Qhull
267 qh_freeqhull (! qh_ALL); 281 qh_freeqhull (! qh_ALL);
268 282
269 int curlong, totlong; 283 int curlong, totlong;
270 qh_memfreeshort (&curlong, &totlong); 284 qh_memfreeshort (&curlong, &totlong);
271 285