matplotlib/src/qhull_wrap.c at v3.2.x · matplotlib/matplotlib · GitHub

Repository navigation

History

354 lines (304 loc) · 10.6 KB

Raw

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

319

320

321

322

323

324

325

326

327

328

329

330

331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

* Wrapper module for libqhull, providing Delaunay triangulation.

* This module's methods should not be accessed directly. To obtain a Delaunay

* triangulation, construct an instance of the matplotlib.tri.Triangulation

* class without specifying a triangles array.

#define PY_SSIZE_T_CLEAN

#include "Python.h"

#include "numpy/ndarrayobject.h"

#include "libqhull/qhull_a.h"

#include <stdio.h>

#ifndef MPL_DEVNULL

#error "MPL_DEVNULL must be defined as the OS-equivalent of /dev/null"

#endif

#define STRINGIFY(x) STR(x)

#define STR(x) #x

static const char* qhull_error_msg[6] = {

"", /* 0 = qh_ERRnone */

"input inconsistency", /* 1 = qh_ERRinput */

"singular input data", /* 2 = qh_ERRsingular */

"precision error", /* 3 = qh_ERRprec */

"insufficient memory", /* 4 = qh_ERRmem */

"internal error"}; /* 5 = qh_ERRqhull */

/* Return the indices of the 3 vertices that comprise the specified facet (i.e.

* triangle). */

static void

get_facet_vertices(const facetT* facet, int indices[3])

{

vertexT *vertex, **vertexp;

FOREACHvertex_(facet->vertices)

*indices++ = qh_pointid(vertex->point);

}

/* Return the indices of the 3 triangles that are neighbors of the specified

* facet (triangle). */

static void

get_facet_neighbours(const facetT* facet, const int* tri_indices,

int indices[3])

{

facetT *neighbor, **neighborp;

FOREACHneighbor_(facet)

*indices++ = (neighbor->upperdelaunay ? -1 : tri_indices[neighbor->id]);

}

/* Return 1 if the specified points arrays contain at least 3 unique points,

* or 0 otherwise. */

static int

at_least_3_unique_points(int npoints, const double* x, const double* y)

{

int i;

const int unique1 = 0; /* First unique point has index 0. */

int unique2 = 0; /* Second unique point index is 0 until set. */

if (npoints < 3)

return 0;

for (i = 1; i < npoints; ++i) {

if (unique2 == 0) {

/* Looking for second unique point. */

if (x[i] != x[unique1] || y[i] != y[unique1])

unique2 = i;

}

else {

/* Looking for third unique point. */

if ( (x[i] != x[unique1] || y[i] != y[unique1]) &&

(x[i] != x[unique2] || y[i] != y[unique2]) ) {

/* 3 unique points found, with indices 0, unique2 and i. */

return 1;

}

/* Run out of points before 3 unique points found. */

return 0;

}

/* Delaunay implementation methyod. If hide_qhull_errors is 1 then qhull error

* messages are discarded; if it is 0 then they are written to stderr. */

static PyObject*

delaunay_impl(int npoints, const double* x, const double* y,

int hide_qhull_errors)

{

coordT* points = NULL;

facetT* facet;

int i, ntri, max_facet_id;

FILE* error_file = NULL; /* qhull expects a FILE* to write errors to. */

int exitcode; /* Value returned from qh_new_qhull(). */

int* tri_indices = NULL; /* Maps qhull facet id to triangle index. */

int indices[3];

int curlong, totlong; /* Memory remaining after qh_memfreeshort. */

PyObject* tuple; /* Return tuple (triangles, neighbors). */

const int ndim = 2;

npy_intp dims[2];

PyArrayObject* triangles = NULL;

PyArrayObject* neighbors = NULL;

int* triangles_ptr;

int* neighbors_ptr;

double x_mean = 0.0;

double y_mean = 0.0;

QHULL_LIB_CHECK

/* Allocate points. */

points = (coordT*)malloc(npoints*ndim*sizeof(coordT));

if (points == NULL) {

PyErr_SetString(PyExc_MemoryError,

"Could not allocate points array in qhull.delaunay");

goto error_before_qhull;

}

/* Determine mean x, y coordinates. */

for (i = 0; i < npoints; ++i) {

x_mean += x[i];

y_mean += y[i];

}

x_mean /= npoints;

y_mean /= npoints;

/* Prepare points array to pass to qhull. */

for (i = 0; i < npoints; ++i) {

points[2*i ] = x[i] - x_mean;

points[2*i+1] = y[i] - y_mean;

}

/* qhull expects a FILE* to write errors to. */

if (hide_qhull_errors) {

/* qhull errors are ignored by writing to OS-equivalent of /dev/null.

* Rather than have OS-specific code here, instead it is determined by

* setupext.py and passed in via the macro MPL_DEVNULL. */

error_file = fopen(STRINGIFY(MPL_DEVNULL), "w");

if (error_file == NULL) {

PyErr_SetString(PyExc_RuntimeError,

"Could not open devnull in qhull.delaunay");

goto error_before_qhull;

}

else {

/* qhull errors written to stderr. */

error_file = stderr;

}

/* Perform Delaunay triangulation. */

exitcode = qh_new_qhull(ndim, npoints, points, False,

"qhull d Qt Qbb Qc Qz", NULL, error_file);

if (exitcode != qh_ERRnone) {

PyErr_Format(PyExc_RuntimeError,

"Error in qhull Delaunay triangulation calculation: %s (exitcode=%d)%s",

qhull_error_msg[exitcode], exitcode,

hide_qhull_errors ? "; use python verbose option (-v) to see original qhull error." : "");

goto error;

}

/* Split facets so that they only have 3 points each. */

qh_triangulate();

/* Determine ntri and max_facet_id.

Note that libqhull uses macros to iterate through collections. */

ntri = 0;

FORALLfacets {

if (!facet->upperdelaunay)

++ntri;

}

max_facet_id = qh facet_id - 1;

/* Create array to map facet id to triangle index. */

tri_indices = (int*)malloc((max_facet_id+1)*sizeof(int));

if (tri_indices == NULL) {

PyErr_SetString(PyExc_MemoryError,

"Could not allocate triangle map in qhull.delaunay");

goto error;

}

/* Allocate python arrays to return. */

dims[0] = ntri;

dims[1] = 3;

triangles = (PyArrayObject*)PyArray_SimpleNew(ndim, dims, NPY_INT);

if (triangles == NULL) {

PyErr_SetString(PyExc_MemoryError,

"Could not allocate triangles array in qhull.delaunay");

goto error;

}

neighbors = (PyArrayObject*)PyArray_SimpleNew(ndim, dims, NPY_INT);

if (neighbors == NULL) {

PyErr_SetString(PyExc_MemoryError,

"Could not allocate neighbors array in qhull.delaunay");

goto error;

}

triangles_ptr = (int*)PyArray_DATA(triangles);

neighbors_ptr = (int*)PyArray_DATA(neighbors);

/* Determine triangles array and set tri_indices array. */

i = 0;

FORALLfacets {

if (!facet->upperdelaunay) {

tri_indices[facet->id] = i++;

get_facet_vertices(facet, indices);

*triangles_ptr++ = (facet->toporient ? indices[0] : indices[2]);

*triangles_ptr++ = indices[1];

*triangles_ptr++ = (facet->toporient ? indices[2] : indices[0]);

}

else

tri_indices[facet->id] = -1;

}

/* Determine neighbors array. */

FORALLfacets {

if (!facet->upperdelaunay) {

get_facet_neighbours(facet, tri_indices, indices);

*neighbors_ptr++ = (facet->toporient ? indices[2] : indices[0]);

*neighbors_ptr++ = (facet->toporient ? indices[0] : indices[2]);

*neighbors_ptr++ = indices[1];

}

/* Clean up. */

qh_freeqhull(!qh_ALL);

qh_memfreeshort(&curlong, &totlong);

if (curlong || totlong)

PyErr_WarnEx(PyExc_RuntimeWarning,

"Qhull could not free all allocated memory", 1);

if (hide_qhull_errors)

fclose(error_file);

free(tri_indices);

free(points);

tuple = PyTuple_New(2);

PyTuple_SetItem(tuple, 0, (PyObject*)triangles);

PyTuple_SetItem(tuple, 1, (PyObject*)neighbors);

return tuple;

error:

/* Clean up. */

Py_XDECREF(triangles);

Py_XDECREF(neighbors);

qh_freeqhull(!qh_ALL);

qh_memfreeshort(&curlong, &totlong);

/* Don't bother checking curlong and totlong as raising error anyway. */

if (hide_qhull_errors)

fclose(error_file);

free(tri_indices);

error_before_qhull:

free(points);

return NULL;

}

/* Process python arguments and call Delaunay implementation method. */

static PyObject*

delaunay(PyObject *self, PyObject *args)

{

PyObject* xarg;

PyObject* yarg;

PyArrayObject* xarray;

PyArrayObject* yarray;

PyObject* ret;

int npoints;

const double* x;

const double* y;

if (!PyArg_ParseTuple(args, "OO", &xarg, &yarg)) {

PyErr_SetString(PyExc_ValueError, "expecting x and y arrays");

return NULL;

}

xarray = (PyArrayObject*)PyArray_ContiguousFromObject(xarg, NPY_DOUBLE,

1, 1);

yarray = (PyArrayObject*)PyArray_ContiguousFromObject(yarg, NPY_DOUBLE,

1, 1);

if (xarray == 0 || yarray == 0 ||

PyArray_DIM(xarray,0) != PyArray_DIM(yarray, 0)) {

Py_XDECREF(xarray);

Py_XDECREF(yarray);

PyErr_SetString(PyExc_ValueError,

"x and y must be 1D arrays of the same length");

return NULL;

}

npoints = PyArray_DIM(xarray, 0);

if (npoints < 3) {

Py_XDECREF(xarray);

Py_XDECREF(yarray);

PyErr_SetString(PyExc_ValueError,

"x and y arrays must have a length of at least 3");

return NULL;

}

x = (const double*)PyArray_DATA(xarray);

y = (const double*)PyArray_DATA(yarray);

if (!at_least_3_unique_points(npoints, x, y)) {

Py_XDECREF(xarray);

Py_XDECREF(yarray);

PyErr_SetString(PyExc_ValueError,

"x and y arrays must consist of at least 3 unique points");

return NULL;

}

ret = delaunay_impl(npoints, x, y, Py_VerboseFlag == 0);

Py_XDECREF(xarray);

Py_XDECREF(yarray);

return ret;

}

/* Return qhull version string for assistance in debugging. */

static PyObject*

version(void)

{

return PyBytes_FromString(qh_version);

}

static PyMethodDef qhull_methods[] = {

{"delaunay", (PyCFunction)delaunay, METH_VARARGS, ""},

{"version", (PyCFunction)version, METH_NOARGS, ""},

{NULL, NULL, 0, NULL}

};

static struct PyModuleDef qhull_module = {

PyModuleDef_HEAD_INIT,

"qhull",

"Computing Delaunay triangulations.\n",

-1,

qhull_methods,

NULL, NULL, NULL, NULL

};

PyMODINIT_FUNC

PyInit__qhull(void)

{

PyObject* m;

m = PyModule_Create(&qhull_module);

if (m == NULL) {

return NULL;

}

import_array();

return m;

}

matplotlib/src/qhull_wrap.c at v3.2.x · matplotlib/matplotlib · GitHub

Repository navigation

Expand file tree

Breadcrumbs

qhull_wrap.c

Latest commit

History

Breadcrumbs

qhull_wrap.c

File metadata and controls