pgRouting  2.2
pgRouting extends the PostGIS / PostgreSQL geospatial database to provide geospatial routing functionality.
 All Classes Functions Variables Pages
alpha.c
1 /*PGR-GNU*****************************************************************
2 FILE: alpha.c
3 
4 Copyright (c) 2006 Anton A. Patrushev, Orkney, Inc.
5 Mail: project@pgrouting.org
6 
7 ------
8 
9 This program is free software; you can redistribute it and/or modify
10 it under the terms of the GNU General Public License as published by
11 the Free Software Foundation; either version 2 of the License, or
12 (at your option) any later version.
13 
14 This program is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU General Public License for more details.
18 
19 You should have received a copy of the GNU General Public License
20 along with this program; if not, write to the Free Software
21 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
22 
23 ********************************************************************PGR-GNU*/
24 
25 #include "../../common/src/pgr_types.h"
26 #include "postgres.h"
27 #include "executor/spi.h"
28 #include "funcapi.h"
29 #include "catalog/pg_type.h"
30 #if PGSQL_VERSION > 92
31 #include "access/htup_details.h"
32 #endif
33 
34 #include "alpha.h"
35 
36 #include "fmgr.h"
37 
38 
39 #ifndef PG_MODULE_MAGIC
40 PG_MODULE_MAGIC;
41 #endif
42 
43 /*
44  * Define this to have profiling enabled
45  */
46 //#define PROFILE
47 
48 #ifdef PROFILE
49 #include <sys/time.h>
50 
54 #define MAX_NODES 1000000
55 
56 struct timeval prof_alpha, prof_store, prof_extract, prof_total;
57 long proftime[5];
58 long profipts1, profipts2, profopts;
59 #define profstart(x) do { gettimeofday(&x, NULL); } while (0);
60 #define profstop(n, x) do { struct timeval _profstop; \
61  long _proftime; \
62  gettimeofday(&_profstop, NULL); \
63  _proftime = ( _profstop.tv_sec*1000000+_profstop.tv_usec) - \
64  ( x.tv_sec*1000000+x.tv_usec); \
65  elog(NOTICE, \
66  "PRF(%s) %lu (%f ms)", \
67  (n), \
68  _proftime, _proftime / 1000.0); \
69  } while (0);
70 
71 #else
72 
73 #define profstart(x) do { } while (0);
74 #define profstop(n, x) do { } while (0);
75 
76 
77 #endif // PROFILE
78 
79 
80 Datum alphashape(PG_FUNCTION_ARGS);
81 
82 #undef DEBUG
83 //#define DEBUG 1
84 #include "../../common/src/debug_macro.h"
85 
86 // The number of tuples to fetch from the SPI cursor at each iteration
87 #define TUPLIMIT 1000
88 
89 static char *
90 text2char(text *in)
91 {
92  char *out = palloc(VARSIZE(in));
93 
94  memcpy(out, VARDATA(in), VARSIZE(in) - VARHDRSZ);
95  out[VARSIZE(in) - VARHDRSZ] = '\0';
96  return out;
97 }
98 
99 static int
100 finish(int code, int ret)
101 {
102  code = SPI_finish();
103  if (code != SPI_OK_FINISH )
104  {
105  elog(ERROR,"couldn't disconnect from SPI");
106  return -1 ;
107  }
108  return ret;
109 }
110 
111 
112 typedef struct vertex_columns
113 {
114  int id;
115  int x;
116  int y;
117 
119 
120 
121 
122 static int
123 fetch_vertices_columns(SPITupleTable *tuptable,
125 {
126  vertex_columns->id = SPI_fnumber(SPI_tuptable->tupdesc, "id");
127  vertex_columns->x = SPI_fnumber(SPI_tuptable->tupdesc, "x");
128  vertex_columns->y = SPI_fnumber(SPI_tuptable->tupdesc, "y");
129 
130  if (vertex_columns->id == SPI_ERROR_NOATTRIBUTE ||
131  vertex_columns->x == SPI_ERROR_NOATTRIBUTE ||
132  vertex_columns->y == SPI_ERROR_NOATTRIBUTE)
133  {
134  elog(ERROR, "Error, query must return columns "
135  "'id', 'x' and 'y'");
136  return -1;
137  }
138 
139  if (SPI_gettypeid(SPI_tuptable->tupdesc, vertex_columns->id) != INT4OID ||
140  SPI_gettypeid(SPI_tuptable->tupdesc, vertex_columns->x) != FLOAT8OID ||
141  SPI_gettypeid(SPI_tuptable->tupdesc, vertex_columns->y) != FLOAT8OID)
142  {
143  elog(ERROR,
144  "Error, column 'id' must be of type int4, 'x' and 'y' must be of type float8");
145  return -1;
146  }
147 
148  return 0;
149 }
150 
151 static void
152 fetch_vertex(HeapTuple *tuple, TupleDesc *tupdesc,
153  vertex_columns_t *vertex_columns, vertex_t *target_vertex)
154 {
155  Datum binval;
156  bool isnull;
157 
158  binval = SPI_getbinval(*tuple, *tupdesc, vertex_columns->x, &isnull);
159  if (isnull)
160  elog(ERROR, "x contains a null value");
161  target_vertex->x = DatumGetFloat8(binval);
162 
163  binval = SPI_getbinval(*tuple, *tupdesc, vertex_columns->y, &isnull);
164  if (isnull)
165  elog(ERROR, "y contains a null value");
166  target_vertex->y = DatumGetFloat8(binval);
167 }
168 
169 static int compute_alpha_shape(char* sql, float8 alpha, vertex_t **res, size_t *res_count)
170 {
171 
172  int SPIcode;
173  void *SPIplan;
174  Portal SPIportal;
175  bool moredata = TRUE;
176  uint32_t ntuples;
177  vertex_t *vertices = NULL;
178  uint32_t total_tuples = 0;
179  vertex_columns_t vertex_columns = {.id= -1, .x= -1, .y= -1};
180  char *err_msg;
181  int ret = -1;
182 
183  PGR_DBG("start alpha_shape\n");
184 
185  SPIcode = SPI_connect();
186  if (SPIcode != SPI_OK_CONNECT)
187  {
188  elog(ERROR, "alpha_shape: couldn't open a connection to SPI");
189  return -1;
190  }
191 
192  SPIplan = SPI_prepare(sql, 0, NULL);
193  if (SPIplan == NULL)
194  {
195  elog(ERROR, "alpha_shape: couldn't create query plan via SPI");
196  return -1;
197  }
198 
199  if ((SPIportal = SPI_cursor_open(NULL, SPIplan, NULL, NULL, true)) == NULL)
200  {
201  elog(ERROR, "alpha_shape: SPI_cursor_open('%s') returns NULL", sql);
202  return -1;
203  }
204 
205  while (moredata == TRUE)
206  {
207  SPI_cursor_fetch(SPIportal, TRUE, TUPLIMIT);
208 
209  if (vertex_columns.id == -1)
210  {
211  if (fetch_vertices_columns(SPI_tuptable, &vertex_columns) == -1)
212  return finish(SPIcode, ret);
213  }
214 
215  ntuples = SPI_processed;
216  total_tuples += ntuples;
217  if (!vertices)
218  vertices = palloc(total_tuples * sizeof(vertex_t));
219  else
220  vertices = repalloc(vertices, total_tuples * sizeof(vertex_t));
221 
222  if (vertices == NULL)
223  {
224  elog(ERROR, "Out of memory");
225  return finish(SPIcode, ret);
226  }
227 
228  if (ntuples > 0)
229  {
230  uint32_t t;
231  SPITupleTable *tuptable = SPI_tuptable;
232  TupleDesc tupdesc = SPI_tuptable->tupdesc;
233 
234  for (t = 0; t < ntuples; t++)
235  {
236  HeapTuple tuple = tuptable->vals[t];
237  fetch_vertex(&tuple, &tupdesc, &vertex_columns,
238  &vertices[total_tuples - ntuples + t]);
239  }
240  SPI_freetuptable(tuptable);
241  }
242  else
243  {
244  moredata = FALSE;
245  }
246  }
247 
248 
249  // if (total_tuples < 2) //this was the buggy code of the pgrouting project.
250  // TODO: report this as a bug to the pgrouting project
251  // the CGAL alpha-shape function crashes if called with less than three points!!!
252 
253  if (total_tuples == 0) {
254  elog(ERROR, "Distance is too short. no vertex for alpha shape calculation. alpha shape calculation needs at least 3 vertices.");
255  }
256  if (total_tuples == 1) {
257  elog(ERROR, "Distance is too short. only 1 vertex for alpha shape calculation. alpha shape calculation needs at least 3 vertices.");
258  }
259  if (total_tuples == 2) {
260  elog(ERROR, "Distance is too short. only 2 vertices for alpha shape calculation. alpha shape calculation needs at least 3 vertices.");
261  }
262  if (total_tuples < 3)
263  {
264  // elog(ERROR, "Distance is too short ....");
265  return finish(SPIcode, ret);
266  }
267 
268  PGR_DBG("Calling CGAL alpha-shape\n");
269 
270  profstop("extract", prof_extract);
271  profstart(prof_alpha);
272 
273  ret = alpha_shape(vertices, total_tuples, alpha, res, res_count, &err_msg);
274 
275  profstop("alpha", prof_alpha);
276  profstart(prof_store);
277 
278  if (ret < 0)
279  {
280  //elog(ERROR, "Error computing shape: %s", err_msg);
281  ereport(ERROR, (errcode(ERRCODE_E_R_E_CONTAINING_SQL_NOT_PERMITTED), errmsg("Error computing shape: %s", err_msg)));
282  }
283 
284  return finish(SPIcode, ret);
285 }
286 
287 PG_FUNCTION_INFO_V1(alphashape);
288 
289 Datum alphashape(PG_FUNCTION_ARGS)
290 {
291  FuncCallContext *funcctx;
292  uint32_t call_cntr;
293  uint32_t max_calls;
294  TupleDesc tuple_desc;
295  vertex_t *res = 0;
296 
297  /* stuff done only on the first call of the function */
298  if (SRF_IS_FIRSTCALL())
299  {
300  MemoryContext oldcontext;
301  size_t res_count;
302 
303  // XXX profiling messages are not thread safe
304  profstart(prof_total);
305  profstart(prof_extract);
306 
307  /* create a function context for cross-call persistence */
308  funcctx = SRF_FIRSTCALL_INIT();
309 
310  /* switch to memory context appropriate for multiple function calls */
311  oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
312 
313  compute_alpha_shape(text2char(PG_GETARG_TEXT_P(0)),
314  PG_GETARG_FLOAT8(1), &res, &res_count);
315 
316  /* total number of tuples to be returned */
317  PGR_DBG("Conting tuples number\n");
318  funcctx->max_calls = (uint32_t)res_count;
319  funcctx->user_fctx = res;
320 
321  PGR_DBG("Total count %i", res_count);
322 
323  if (get_call_result_type(fcinfo, NULL, &tuple_desc) != TYPEFUNC_COMPOSITE)
324  ereport(ERROR,
325  (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
326  errmsg("function returning record called in context "
327  "that cannot accept type record")));
328 
329  funcctx->tuple_desc = BlessTupleDesc(tuple_desc);
330 
331  MemoryContextSwitchTo(oldcontext);
332  }
333 
334  /* stuff done on every call of the function */
335  PGR_DBG("Strange stuff doing\n");
336  funcctx = SRF_PERCALL_SETUP();
337 
338  call_cntr = funcctx->call_cntr;
339  max_calls = funcctx->max_calls;
340  tuple_desc = funcctx->tuple_desc;
341  res = (vertex_t*) funcctx->user_fctx;
342 
343  PGR_DBG("Trying to allocate some memory\n");
344 
345  if (call_cntr < max_calls) /* do when there is more left to send */
346  {
347  HeapTuple tuple;
348  Datum result;
349  Datum *values;
350  char* nulls;
351  double x;
352  double y;
353 
354  /* This will work for some compilers. If it crashes with segfault, try to change the following block with this one
355 
356  values = palloc(3 * sizeof(Datum));
357  nulls = palloc(3 * sizeof(char));
358 
359  values[0] = call_cntr;
360  nulls[0] = ' ';
361  values[1] = Float8GetDatum(res[call_cntr].x);
362  nulls[1] = ' ';
363  values[2] = Float8GetDatum(res[call_cntr].y);
364  nulls[2] = ' ';
365  */
366 
367  values = palloc(2 * sizeof(Datum));
368  nulls = palloc(2 * sizeof(char));
369 
370  x = res[call_cntr].x;
371  y = res[call_cntr].y;
372  if (x == DBL_MAX && y == DBL_MAX)
373  {
374  values[0] = 0;
375  values[1] = 0;
376  nulls[0] = 'n';
377  nulls[1] = 'n';
378  }
379  else
380  {
381  values[0] = Float8GetDatum(x);
382  values[1] = Float8GetDatum(y);
383  nulls[0] = ' ';
384  nulls[1] = ' ';
385  }
386 
387  PGR_DBG("Heap making\n");
388 
389  tuple = heap_formtuple(tuple_desc, values, nulls);
390 
391  PGR_DBG("Datum making\n");
392 
393  /* make the tuple into a datum */
394  result = HeapTupleGetDatum(tuple);
395 
396  PGR_DBG("Trying to free some memory\n");
397 
398  /* clean up (this is not really necessary) */
399  pfree(values);
400  pfree(nulls);
401 
402  SRF_RETURN_NEXT(funcctx, result);
403  }
404  else /* do when there is no more left */
405  {
406  if (res) free(res);
407  profstop("store", prof_store);
408  profstop("total", prof_total);
409 #ifdef PROFILE
410  elog(NOTICE, "_________");
411 #endif
412  SRF_RETURN_DONE(funcctx);
413  }
414 }
Definition: alpha.h:33