Ticket #198: alpha.c

File alpha.c, 10.5 KB (added by stefanzweig, 9 months ago)

bugfix for alpha.c

Line 
1/*
2 * Alpha-Shapes for PostgreSQL
3 *
4 * Copyright (c) 2006 Anton A. Patrushev, Orkney, Inc.
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
19 *
20 */
21
22
23#include "postgres.h"
24#include "executor/spi.h"
25#include "funcapi.h"
26#include "catalog/pg_type.h"
27
28#include "alpha.h"
29
30#include "fmgr.h"
31
32
33#ifdef PG_MODULE_MAGIC
34PG_MODULE_MAGIC;
35#endif
36
37/*
38 * Define this to have profiling enabled
39 */
40//#define PROFILE
41
42#ifdef PROFILE
43#include <sys/time.h>
44
45/*! \def MAX_NODES
46  \brief Maximal number of nodes in the path (to avoid infinite loops)
47*/
48#define MAX_NODES 1000000
49
50struct timeval prof_alpha, prof_store, prof_extract, prof_total;
51long proftime[5];
52long profipts1, profipts2, profopts;
53#define profstart(x) do { gettimeofday(&x, NULL); } while (0);
54#define profstop(n, x) do { struct timeval _profstop;   \
55        long _proftime;                         \
56        gettimeofday(&_profstop, NULL);                         \
57        _proftime = ( _profstop.tv_sec*1000000+_profstop.tv_usec) -     \
58                ( x.tv_sec*1000000+x.tv_usec); \
59        elog(NOTICE, \
60                "PRF(%s) %lu (%f ms)", \
61                (n), \
62             _proftime, _proftime / 1000.0);    \
63        } while (0);
64
65#else
66
67#define profstart(x) do { } while (0);
68#define profstop(n, x) do { } while (0);
69
70
71#endif // PROFILE
72
73
74Datum alphashape(PG_FUNCTION_ARGS);
75
76#undef DEBUG
77//#define DEBUG 1
78
79#ifdef DEBUG
80#define DBG(format, arg...)                     \
81    elog(NOTICE, format , ## arg)
82#else
83#define DBG(format, arg...) do { ; } while (0)
84#endif
85   
86// The number of tuples to fetch from the SPI cursor at each iteration
87#define TUPLIMIT 1000
88
89static char *
90text2char(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
99static int
100finish(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
112typedef struct vertex_columns
113{
114  int id;
115  int x;
116  int y;
117
118} vertex_columns_t;
119
120
121
122static int
123fetch_vertices_columns(SPITupleTable *tuptable,
124                       vertex_columns_t *vertex_columns)
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
151static void
152fetch_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
169static int compute_alpha_shape(char* sql, vertex_t **res, int *res_count)
170{
171
172  int SPIcode;
173  void *SPIplan;
174  Portal SPIportal;
175  bool moredata = TRUE;
176  int ntuples;
177  vertex_t *vertices = NULL;
178  int 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  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          int 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  DBG("Calling CGAL alpha-shape\n");
269       
270  profstop("extract", prof_extract);
271  profstart(prof_alpha);
272
273  ret = alpha_shape(vertices, total_tuples, 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
287PG_FUNCTION_INFO_V1(alphashape);
288
289Datum alphashape(PG_FUNCTION_ARGS)
290{
291  FuncCallContext      *funcctx;
292  int                  call_cntr;
293  int                  max_calls;
294  TupleDesc            tuple_desc;
295  vertex_t     *res;
296                   
297  /* stuff done only on the first call of the function */
298  if (SRF_IS_FIRSTCALL())
299    {
300      MemoryContext   oldcontext;
301      int res_count;
302      int ret;
303                           
304      // XXX profiling messages are not thread safe
305      profstart(prof_total);
306      profstart(prof_extract);
307                                           
308      /* create a function context for cross-call persistence */
309      funcctx = SRF_FIRSTCALL_INIT();
310                                                                   
311      /* switch to memory context appropriate for multiple function calls */
312      oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
313
314      ret = compute_alpha_shape(text2char(PG_GETARG_TEXT_P(0)),
315                                &res, &res_count);
316
317      /* total number of tuples to be returned */
318      DBG("Conting tuples number\n");
319      funcctx->max_calls = res_count;
320      funcctx->user_fctx = res;
321
322      DBG("Total count %i", res_count);
323
324      funcctx->tuple_desc = BlessTupleDesc(
325                           RelationNameGetTupleDesc("vertex_result"));
326
327      MemoryContextSwitchTo(oldcontext);
328    }
329
330  /* stuff done on every call of the function */
331  DBG("Strange stuff doing\n");
332  funcctx = SRF_PERCALL_SETUP();
333
334  call_cntr = funcctx->call_cntr;
335  max_calls = funcctx->max_calls;
336  tuple_desc = funcctx->tuple_desc;
337  res = (vertex_t*) funcctx->user_fctx;
338
339  DBG("Trying to allocate some memory\n");
340
341  if (call_cntr < max_calls)    /* do when there is more left to send */
342    {
343      HeapTuple    tuple;
344      Datum        result;
345      Datum *values;
346      char* nulls;
347
348      /* This will work for some compilers. If it crashes with segfault, try to change the following block with this one   
349
350      values = palloc(3 * sizeof(Datum));
351      nulls = palloc(3 * sizeof(char));
352
353      values[0] = call_cntr;
354      nulls[0] = ' ';
355      values[1] = Float8GetDatum(res[call_cntr].x);
356      nulls[1] = ' ';
357      values[2] = Float8GetDatum(res[call_cntr].y);
358      nulls[2] = ' ';
359      */
360   
361      values = palloc(2 * sizeof(Datum));
362      nulls = palloc(2 * sizeof(char));
363
364      values[0] = Float8GetDatum(res[call_cntr].x);
365      nulls[0] = ' ';
366      values[1] = Float8GetDatum(res[call_cntr].y);
367      nulls[1] = ' ';
368       
369      DBG("Heap making\n");
370
371      tuple = heap_formtuple(tuple_desc, values, nulls);
372
373      DBG("Datum making\n");
374
375      /* make the tuple into a datum */
376      result = HeapTupleGetDatum(tuple);
377
378      DBG("Trying to free some memory\n");
379
380      /* clean up (this is not really necessary) */
381      pfree(values);
382      pfree(nulls);
383
384      SRF_RETURN_NEXT(funcctx, result);
385    }
386  else    /* do when there is no more left */
387    {
388      profstop("store", prof_store);
389      profstop("total", prof_total);
390#ifdef PROFILE
391      elog(NOTICE, "_________");
392#endif
393      SRF_RETURN_DONE(funcctx);
394    }
395}