branches/2d_boundary/Tests/svgfig/_curve.c

   1 #include <Python.h>
   2
   3 #ifndef PyMODINIT_FUNC
   4 #define PyMODINIT_FUNC void
   5 #endif
   6
   7 struct sample {
   8   double t, x, y;
   9   int discontinuity;
  10   struct sample *left;
  11   struct sample *right;
  12 };
  13
  14 struct common_block {
  15   PyObject *parametric;
  16   PyObject *listoftrans;
  17   int counter;
  18   PyObject *random_sampling;
  19   int recursion_limit;
  20   double linearity_limit;
  21   double discontinuity_limit;
  22   unsigned int MT[624];
  23   int count624;
  24 };
  25
  26 /* An implementation of the Mersenne Twistor random algorithm */
  27 /* Copyright (C) 1997 Makoto Matsumoto and Takuji Nishimura. */
  28 /* Copied from ROOT's TRandom3 */
  29 /* (wanted to avoid compile-time dependencies for such a small part of the program) */
  30 static void _curve_setseed(int seed, struct common_block *block) {
  31   block->count624 = 624;
  32   if (seed > 0) {
  33     block->MT[0] = seed;
  34   }
  35   int i;
  36   for (i = 1;  i < 624;  i++) {
  37     block->MT[i] = (1812433253 * (block->MT[i-1] ^ (block->MT[i-1] >> 30)) + i);
  38   }
  39 }
  40
  41 static double _curve_random(struct common_block *block) {
  42    unsigned int y;
  43    int kM = 397;
  44    int kN = 624;
  45    unsigned int kTemperingMaskB = 0x9d2c5680;
  46    unsigned int kTemperingMaskC = 0xefc60000;
  47    unsigned int kUpperMask = 0x80000000;
  48    unsigned int kLowerMask = 0x7fffffff;
  49    unsigned int kMatrixA = 0x9908b0df;
  50
  51    if (block->count624 >= kN) {
  52       int i;
  53       for (i = 0;  i < kN-kM;  i++) {
  54          y = (block->MT[i] & kUpperMask) | (block->MT[i+1] & kLowerMask);
  55          block->MT[i] = block->MT[i+kM] ^ (y >> 1) ^ ((y & 0x1) ? kMatrixA : 0x0);
  56       }
  57
  58       for (;  i < kN-1;  i++) {
  59          y = (block->MT[i] & kUpperMask) | (block->MT[i+1] & kLowerMask);
  60          block->MT[i] = block->MT[i+kM-kN] ^ (y >> 1) ^ ((y & 0x1) ? kMatrixA : 0x0);
  61       }
  62
  63       y = (block->MT[kN-1] & kUpperMask) | (block->MT[0] & kLowerMask);
  64       block->MT[kN-1] = block->MT[kM-1] ^ (y >> 1) ^ ((y & 0x1) ? kMatrixA : 0x0);
  65       block->count624 = 0;
  66    }
  67
  68    y = block->MT[block->count624++];
  69    y ^= (y >> 11);
  70    y ^= ((y << 7 ) & kTemperingMaskB);
  71    y ^= ((y << 15) & kTemperingMaskC);
  72    y ^= (y >> 18);
  73
  74    if (y != 0) return ((double) y * 2.3283064365386963e-10); // * pow(2, -32)
  75    return _curve_random(block);
  76 }
  77
  78 /* evaluate parametric function at a point, passing it through a list of coordinate transformations */
  79 static int _curve_eval(PyObject *parametric, PyObject *listoftrans, double t, double *fx, double *fy) {
  80   PyObject *args = Py_BuildValue("(d)", t);
  81   PyObject *result = PyObject_CallObject(parametric, args);
  82   if (result == NULL) {
  83     Py_DECREF(args);
  84     return 0;
  85   }
  86   Py_DECREF(args);
  87
  88   if (!PySequence_Check(result)  ||  PySequence_Size(result) != 2) {
  89     PyErr_SetString(PyExc_TypeError, "The parametric function must return two real values.");
  90     Py_DECREF(result);
  91     return 0;
  92   }
  93
  94   PyObject *x = PySequence_GetItem(result, 0);
  95   PyObject *y = PySequence_GetItem(result, 1);
  96   Py_DECREF(result);
  97
  98   if (!PyNumber_Check(x)  ||  !PyNumber_Check(y)) {
  99     PyErr_SetString(PyExc_TypeError, "The parametric function must return two real values.");
 100     Py_DECREF(x);
 101     Py_DECREF(y);
 102     return 0;
 103   }
 104
 105   int i;
 106   int lenlistoftrans = PySequence_Size(listoftrans);
 107   for (i = 0;  i < lenlistoftrans;  i++) {
 108     PyObject *trans = PySequence_GetItem(listoftrans, i);
 109
 110     args = Py_BuildValue("(OO)", x, y);
 111     Py_DECREF(x);
 112     Py_DECREF(y);
 113     result = PyObject_CallObject(trans, args);
 114     if (result == NULL) {
 115       Py_DECREF(trans);
 116       Py_DECREF(args);
 117       return 0;
 118     }
 119     Py_DECREF(trans);
 120     Py_DECREF(args);
 121
 122     if (!PySequence_Check(result)  ||  PySequence_Size(result) != 2) {
 123       PyErr_SetString(PyExc_TypeError, "The transformation functions must return two real values.");
 124       Py_DECREF(result);
 125       return 0;
 126     }
 127
 128     x = PySequence_GetItem(result, 0);  // only borrowed
 129     y = PySequence_GetItem(result, 1);  // only borrowed
 130     Py_DECREF(result);
 131
 132     if (!PyNumber_Check(x)  ||  !PyNumber_Check(y)) {
 133       PyErr_SetString(PyExc_TypeError, "The transformation functions must return two real values.");
 134       Py_DECREF(x);
 135       Py_DECREF(y);
 136       return 0;
 137     }
 138   }
 139
 140   *fx = PyFloat_AsDouble(x);
 141   *fy = PyFloat_AsDouble(y);
 142   Py_DECREF(x);
 143   Py_DECREF(y);
 144
 145   return 1;
 146 }
 147
 148 /* recursively called to fill a (doubly-linked) list of sample points where it needs it most */
 149 /* with the default parameters, it computes a few more points than are typically needed */
 150 /* after pruning the extras, that guarantees a nice smooth curve */
 151 int _curve_subsample(struct sample *left, struct sample *right, int depth, struct common_block *block) {
 152   /* make new mid node and link it up */
 153   struct sample *mid = (struct sample*)malloc(sizeof(struct sample));
 154   if (mid == NULL) {
 155     PyErr_Format(PyExc_MemoryError, "Ran out of memory while sampling function (%d nodes created)", block->counter);
 156     return 0;
 157   }
 158   block->counter++;
 159
 160   if (block->random_sampling == Py_True) {
 161     mid->t = left->t + (0.3 + 0.4*_curve_random(block))*(right->t - left->t);
 162   }
 163   else {
 164     mid->t = left->t + 0.5*(right->t - left->t);
 165   }
 166
 167   if (!_curve_eval(block->parametric, block->listoftrans, mid->t, &(mid->x), &(mid->y))) {
 168     free(mid);
 169     return 0;
 170   }
 171   mid->discontinuity = 0;
 172
 173   left->right = mid;
 174   mid->left = left;
 175   mid->right = right;
 176   right->left = mid;
 177
 178   /* calculate the distance of closest approach of mid to the line between left and right */
 179   double numer = (left->x)*(right->y - mid->y) + (mid->x)*(left->y - right->y) + (right->x)*(mid->y - left->y);
 180   double denom = sqrt((left->x - right->x)*(left->x - right->x) + (left->y - right->y)*(left->y - right->y));
 181
 182   if (depth < 3  ||
 183       (denom == 0.  &&  left->t != right->t)  ||
 184       (denom > block->discontinuity_limit)  ||
 185       (denom != 0.  &&  fabs(numer/denom) > block->linearity_limit)) {
 186
 187     if (depth < block->recursion_limit) {
 188       if (!_curve_subsample(left, mid, depth+1, block)) return 0;
 189       if (!_curve_subsample(mid, right, depth+1, block)) return 0;
 190     }
 191
 192     else {
 193       /* we've sampled many points and yet it's still not a small linear gap */
 194       /* break the line: we've found a discontinuity */
 195       mid->discontinuity = 1;
 196     }
 197   }
 198
 199   return 1;
 200 }
 201
 202 /* the only function which is called from the outside: the interface to Python */
 203 static PyObject *_curve_curve(PyObject *self, PyObject *args, PyObject *kwds) {
 204   const char *errstring = "arguments are: parametric function to plot, list of transformations to apply to each point, low endpoint, high endpoint.  \nkeyword arguments are: random_sampling (True), random_seed (12345), recursion_limit (15), linearity_limit (0.05), discontinuity_limit (5.)";
 205
 206   PyObject *parametric;
 207   PyObject *listoftrans;
 208   double low, high;
 209   PyObject *random_sampling = Py_True;
 210   int random_seed = 12345;
 211   int recursion_limit = 15;
 212   double linearity_limit = 0.05;
 213   double discontinuity_limit = 5.;
 214
 215   static char *kwlist[] = {"parametric", "listoftrans", "low", "high", "random_sampling", "random_seed", "recursion_limit", "linearity_limit", "discontinuity_limit", NULL};
 216   if (!PyArg_ParseTupleAndKeywords(args, kwds, "OOdd|Oiidd", kwlist, &parametric, &listoftrans, &low, &high, &random_sampling, &random_seed, &recursion_limit, &linearity_limit, &discontinuity_limit)) {
 217     PyErr_SetString(PyExc_TypeError, errstring);
 218     return NULL;
 219   }
 220
 221   if (!PyCallable_Check(parametric)) {
 222     PyErr_SetString(PyExc_TypeError, errstring);
 223     return NULL;
 224   }
 225
 226   if (!PySequence_Check(listoftrans)) {
 227     PyErr_SetString(PyExc_TypeError, errstring);
 228     return NULL;
 229   }
 230
 231   int i;
 232   int lenlistoftrans = PySequence_Size(listoftrans);
 233   for (i = 0;  i < lenlistoftrans;  i++) {
 234     PyObject *trans = PySequence_GetItem(listoftrans, i);
 235     if (!PyCallable_Check(trans)) {
 236       PyErr_SetString(PyExc_TypeError, errstring);
 237       Py_DECREF(trans);
 238       return NULL;
 239     }
 240     Py_DECREF(trans);
 241   }
 242
 243   if (random_sampling != Py_True  &&  random_sampling != Py_False) {
 244     PyErr_SetString(PyExc_TypeError, errstring);
 245     return NULL;
 246   }
 247
 248   /* build doubly-linked list of samples */
 249   struct sample *samplelow = (struct sample*)malloc(sizeof(struct sample));
 250   if (samplelow == NULL) {
 251     PyErr_SetString(PyExc_MemoryError, "Ran out of memory while sampling function (1 node created)");
 252     return NULL;
 253   }
 254   samplelow->t = low;
 255   if (!_curve_eval(parametric, listoftrans, samplelow->t, &(samplelow->x), &(samplelow->y))) {
 256     free(samplelow);
 257     return NULL;
 258   }
 259   samplelow->discontinuity = 0;
 260
 261   struct sample *samplehigh = (struct sample*)malloc(sizeof(struct sample));
 262   if (samplelow == NULL) {
 263     PyErr_SetString(PyExc_MemoryError, "Ran out of memory while sampling function (2 nodes created)");
 264     free(samplelow);
 265     return NULL;
 266   }
 267   samplehigh->t = high;
 268   if (!_curve_eval(parametric, listoftrans, samplehigh->t, &(samplehigh->x), &(samplehigh->y))) {
 269     free(samplelow);
 270     free(samplehigh);
 271     return NULL;
 272   }
 273   samplehigh->discontinuity = 0;
 274
 275   samplelow->left = NULL;
 276   samplelow->right = samplehigh;
 277   samplehigh->left = samplelow;
 278   samplehigh->right = NULL;
 279
 280   struct common_block block;
 281   block.parametric = parametric;
 282   block.listoftrans = listoftrans;
 283   block.counter = 2;
 284   block.random_sampling = random_sampling;
 285   block.recursion_limit = recursion_limit;
 286   block.linearity_limit = linearity_limit;
 287   block.discontinuity_limit = discontinuity_limit;
 288   _curve_setseed(random_seed, &block);
 289
 290   /* recursively find most of the points */
 291   if (!_curve_subsample(samplelow, samplehigh, 0, &block)) {
 292     /* free nodes in case of error */
 293     struct sample *p = samplelow;
 294     while (p != NULL) {
 295       struct sample *next = p->right;
 296       free(p);
 297       p = next;
 298     }
 299     return NULL;
 300   }
 301
 302   /* prune excess points that are within the linearity bounds */
 303   struct sample *left = samplelow;
 304   int length = 1; /* get the post-pruning length */
 305   while (left->right != NULL) {
 306     struct sample *mid = left->right;
 307     struct sample *right = mid->right;
 308
 309     if (right != NULL  &&  !left->discontinuity  &&  !mid->discontinuity  &&  !right->discontinuity) {
 310       double numer = (left->x)*(right->y - mid->y) + (mid->x)*(left->y - right->y) + (right->x)*(mid->y - left->y);
 311       double denom = sqrt((left->x - right->x)*(left->x - right->x) + (left->y - right->y)*(left->y - right->y));
 312
 313       if (denom != 0.  &&  fabs(numer/denom) < linearity_limit) {
 314         free(mid); /* drop this point; it doesn't contribute to the smoothness of the curve */
 315         left->right = right;
 316         right->left = left;
 317       }
 318       else {
 319         length++;
 320         left = left->right; /* increment left */
 321       }
 322     }
 323     else {
 324       length++;
 325       left = left->right; /* increment left */
 326     }
 327   }
 328
 329   /* return a Python tuple of numbers and free all nodes */
 330   PyObject *output = PyTuple_New(length);
 331   int failure = 0;
 332
 333   struct sample *p = samplelow;
 334   for (i = 0;  i < length;  i++) {
 335     if (p->discontinuity) {
 336       Py_INCREF(Py_None);
 337       if (PyTuple_SetItem(output, i, Py_None) != 0) failure = 1;
 338     }
 339     else {
 340       if (PyTuple_SetItem(output, i, Py_BuildValue("dd", p->x, p->y)) != 0) failure = 1;
 341     }
 342
 343     struct sample *next = p->right;
 344     free(p);
 345     p = next;
 346   }
 347
 348   /* it's important to finish freeing nodes first */
 349   if (failure) {
 350     Py_DECREF(output);
 351     return NULL;
 352   }
 353
 354   return output;
 355 }
 356
 357 static PyMethodDef _curve_methods[] = {
 358   {"curve", ((PyCFunction)(_curve_curve)), METH_VARARGS | METH_KEYWORDS, ""},
 359   {NULL}
 360 };
 361
 362 PyMODINIT_FUNC init_curve() {
 363   Py_InitModule3("_curve", _curve_methods, "");
 364 }