scpolytope.c
Go to the documentation of this file.
1 #include "scpolytope.h"
2 
3 #include "defines.h"
4 #include "random.h"
5 #include "algebra.h"
6 #include "geom.h"
7 
8 #include "math.h"
9 #include "string.h"
10 
21 void InitEmptySPolytope(double delta,unsigned int k,double r,double sr,Tscpolytope *mp)
22 {
23  mp->r=r;
24 
25  mp->sr=sr;
26  mp->lsr=mp->r+2*delta;
27  if (mp->sr<mp->lsr)
28  mp->sr=mp->lsr;
29  mp->msr=mp->sr;
30 
31  mp->k=k;
32 
33  mp->nFaces=0;
34  mp->maxFaces=0;
35 
36  /* Empty space for faces*/
37  mp->face=NULL;
38  mp->ID=NULL;
39 
40  mp->v=SPolytopeMaxVolume(mp);
41 }
42 
44 {
45  if (mp->maxFaces==0)
46  {
47  /* Define space for the faces */
49  NEW(mp->face,mp->maxFaces,double *);
50  NEW(mp->ID,mp->maxFaces,unsigned int);
51  }
52 }
53 
54 void CopySPolytope(Tscpolytope *mp_dst,Tscpolytope *mp_src)
55 {
56  mp_dst->k=mp_src->k;
57  mp_dst->r=mp_src->r;
58  mp_dst->sr=mp_src->sr;
59  mp_dst->lsr=mp_src->lsr;
60  mp_dst->msr=mp_src->msr;
61  mp_dst->v=mp_src->v;
62 
63  mp_dst->nFaces=mp_src->nFaces;
64  mp_dst->maxFaces=mp_src->maxFaces;
65 
66  if (mp_src->maxFaces>0)
67  {
68  unsigned int i,j;
69 
70  NEW(mp_dst->face,mp_dst->maxFaces,double *);
71  NEW(mp_dst->ID,mp_dst->maxFaces,unsigned int);
72  for(i=0;i<mp_dst->nFaces;i++)
73  {
74  NEW(mp_dst->face[i],mp_dst->k+1,double);
75  for(j=0;j<mp_dst->k+1;j++)
76  mp_dst->face[i][j]=mp_src->face[i][j];
77 
78  mp_dst->ID[i]=mp_src->ID[i];
79  }
80  }
81  else
82  {
83  mp_dst->face=NULL;
84  mp_dst->ID=NULL;
85  }
86 }
87 
88 boolean InsideSPolytope(double *t,Tscpolytope *mp)
89 {
90  unsigned int i;
91  boolean inside;
92 
93  inside=(Norm(mp->k,t)<=mp->msr);
94 
95  i=0;
96  while((i<mp->nFaces)&&(inside))
97  {
98  inside=(GeneralDotProduct(mp->k,mp->face[i],t)+mp->face[i][mp->k]<=0);
99  i++;
100  }
101 
102  return(inside);
103 }
104 
105 unsigned int DetermineSPolytopeNeighbour(double epsilon,double *t,Tscpolytope *mp)
106 {
107  if ((Norm(mp->k,t)<=mp->r)&&(!InsideSPolytope(t,mp)))
108  {
109  unsigned int nv,s;
110  unsigned int *v;
111  double a1,a2,ac;
112  double *p;
113  unsigned int i,id;
114  boolean found;
115 
116  NEW(v,mp->nFaces,unsigned int);
117  nv=0;
118  for(i=0;i<mp->nFaces;i++)
119  {
120  if (GeneralDotProduct(mp->k,mp->face[i],t)+mp->face[i][mp->k]>0)
121  {
122  v[nv]=i;
123  nv++;
124  }
125  }
126 
127  if (nv>1)
128  {
129  NEW(p,mp->k,double);
130  a1=0.0;
131  a2=1.0;
132  ac=0.9;
133  found=FALSE;
134  while(!found)
135  {
136  memcpy(p,t,sizeof(double)*mp->k);
137  ScaleVector(ac,mp->k,p);
138 
139  s=0;
140  for(i=0;i<nv;i++)
141  {
142  id=v[i];
143  if (GeneralDotProduct(mp->k,mp->face[id],p)+mp->face[id][mp->k]>0)
144  s++;
145  }
146 
147  found=((s==1)||((a2-a1)<epsilon));
148  if (!found)
149  {
150  if (nv==0)
151  a1=ac;
152  else
153  a2=ac;
154  ac=0.1*a1+0.9*a2;
155  }
156  }
157  free(p);
158  }
159 
160  id=mp->ID[v[0]];
161  free(v);
162 
163  return(id);
164  }
165  else
166  return(NO_UINT);
167 }
168 
169 
170 void EnlargeSPolytope(double *t,Tscpolytope *mp)
171 {
172  unsigned int i;
173  double v;
174 
175  for(i=0;i<mp->nFaces;i++)
176  {
177  v=(GeneralDotProduct(mp->k,mp->face[i],t)+mp->face[i][mp->k]);
178  if (v>0)
179  {
180  mp->face[i][mp->k]-=v;
181  mp->v=-1.0; /* the volume is invalidated */
182  }
183  }
184 }
185 
186 void CutSPolytope(double *t,double r,unsigned int id,Tscpolytope *mp)
187 {
188  double n;
189 
190  /*The plane separating the two spheres has parameters 't'
191  and the offset is computed here */
192  n=Norm(mp->k,t);
193 
194  /* Note that typically mp->r==r and thus
195  -(mp->r*mp->r-r*r+n*n)/2.0 = -n*n/2.0*/
196  CutSPolytopeWithFace(t,-(mp->r*mp->r-r*r+n*n)/2.0,id,mp);
197 }
198 
199 void CutSPolytopeWithFace(double *t,double offset,unsigned int id,Tscpolytope *mp)
200 {
201  unsigned int i;
202 
203  /* empty identifiers are used in normal polytopes to bound the
204  initial polytope (a box). Here the box is explicitly defined
205  (no need of this extra faces). */
206  if (id!=NO_UINT)
207  {
208  if (mp->maxFaces==0)
209  DefineSPolytope(mp); /*make initial room for the faces*/
210 
211  /*Store the new face*/
212  if (mp->nFaces==mp->maxFaces)
213  {
214  MEM_DUP(mp->face,mp->maxFaces,double *);
215  MEM_EXPAND(mp->ID,mp->maxFaces,unsigned int);
216  }
217 
218  NEW(mp->face[mp->nFaces],mp->k+1,double);
219  for(i=0;i<mp->k;i++)
220  mp->face[mp->nFaces][i]=t[i];
221  mp->face[mp->nFaces][mp->k]=offset;
222  mp->ID[mp->nFaces]=id;
223  mp->nFaces++;
224 
225  mp->v=-1.0; /* the volume is invalidated */
226  }
227 }
228 
230 {
231  return(mp->r);
232 }
233 
235 {
236  return(2*mp->sr);
237 }
238 
239 unsigned int GetSPolytopeDim(Tscpolytope *mp)
240 {
241  return(mp->k);
242 }
243 
245 {
246  return(mp->nFaces);
247 }
248 
249 void GetSPolytopeFace(unsigned int n,double *f,Tscpolytope *mp)
250 {
251  if (n<mp->nFaces)
252  memcpy(f,mp->face[n],(mp->k+1)*sizeof(double));
253 }
254 
255 boolean SPolytopeRandomPointOnBoundary(double rSample,double *t,Tscpolytope *mp)
256 {
257  randomOnBall(rSample,mp->k,t);
258  return(InsideSPolytope(t,mp));
259 }
260 
261 
262 boolean RandomPointInSPolytope(double scale,double *t,Tscpolytope *mp)
263 {
264  boolean havePoint;
265  double sr;
266 
267  sr=scale*mp->sr;
268  if (sr>mp->msr) sr=mp->msr;
269  else if (sr<mp->lsr) sr=mp->lsr;
270 
271  randomInBall(sr,mp->k,t);
272 
273  if (mp->nFaces==0)
274  havePoint=TRUE;
275  else
276  havePoint=InsideSPolytope(t,mp);
277 
278  return(havePoint);
279 }
280 
282 {
283  return(mp->sr);
284 }
285 
287 {
288  mp->sr*=(1.0+MOV_AVG_UP);
289  if (mp->sr>mp->msr)
290  mp->sr=mp->msr;
291 }
292 
294 {
295  mp->sr*=(1.0-MOV_AVG_DOWN);
296  if (mp->sr<mp->lsr)
297  mp->sr=mp->lsr;
298 }
299 
301 {
302  return(BallVolume(mp->k,mp->sr));
303 }
304 
306 {
307  if (mp->v<0)
308  {
309  unsigned int i,s,n;
310  double *t;
311 
312  NEW(t,mp->k,double);
313  s=0;
314  n=(unsigned int)floor(pow(10,mp->k));
315  if (n>10000) n=10000;
316  for(i=0;i<n;i++)
317  {
318  randomInBall(mp->sr,mp->k,t);
319  if (InsideSPolytope(t,mp))
320  s++;
321  }
322  free(t);
323  mp->v=BallVolume(mp->k,mp->sr)*((double)s)/((double)n);
324  }
325  return(mp->v);
326 }
327 
329 {
330  return(mp->nFaces);
331 }
332 
333 unsigned int SPolytopeNeighbourID(unsigned int n,Tscpolytope *mp)
334 {
335  if (n<mp->nFaces)
336  return(mp->ID[n]);
337  else
338  return(NO_UINT);
339 }
340 
341 unsigned int SPolytopeMemSize(Tscpolytope *mp)
342 {
343  unsigned int b;
344 
345  b=sizeof(double)*mp->nFaces*(mp->k+1); /* face[][] */
346  b+=sizeof(unsigned int)*mp->nFaces; /* ID[] */
347 
348  return(b);
349 }
350 
351 void SaveSPolytope(FILE *f,Tscpolytope *mp)
352 {
353  unsigned int i,j;
354 
355  fprintf(f,"%u\n",mp->k);
356  fprintf(f,"%f\n",mp->r);
357  fprintf(f,"%f\n",mp->sr);
358  fprintf(f,"%f\n",mp->lsr);
359  fprintf(f,"%f\n",mp->msr);
360  fprintf(f,"%f\n",mp->v);
361 
362  fprintf(f,"%u\n",mp->nFaces);
363  fprintf(f,"%u\n",mp->maxFaces);
364 
365  if (mp->nFaces>0)
366  {
367  for(i=0;i<mp->nFaces;i++)
368  {
369  for(j=0;j<mp->k+1;j++)
370  fprintf(f,"%.12f ",mp->face[i][j]);
371  fprintf(f,"\n");
372  fprintf(f,"%u\n",mp->ID[i]);
373  }
374  }
375 }
376 
377 void LoadSPolytope(FILE *f,Tscpolytope *mp)
378 {
379  unsigned int i,j;
380 
381  fscanf(f,"%u",&(mp->k));
382  fscanf(f,"%lf",&(mp->r));
383  fscanf(f,"%lf",&(mp->sr));
384  fscanf(f,"%lf",&(mp->lsr));
385  fscanf(f,"%lf",&(mp->msr));
386  fscanf(f,"%lf",&(mp->v));
387 
388  fscanf(f,"%u",&(mp->nFaces));
389  fscanf(f,"%u",&(mp->maxFaces));
390 
391  if (mp->maxFaces>0)
392  {
393  NEW(mp->face,mp->maxFaces,double *);
394  NEW(mp->ID,mp->maxFaces,unsigned int);
395  for(i=0;i<mp->nFaces;i++)
396  {
397  NEW(mp->face[i],mp->k+1,double);
398  for(j=0;j<mp->k+1;j++)
399  fscanf(f,"%lf",&(mp->face[i][j]));
400 
401  fscanf(f,"%u",&(mp->ID[i]));
402  }
403  }
404  else
405  {
406  mp->face=NULL;
407  mp->ID=NULL;
408  }
409 }
410 
412 {
413  if (mp->maxFaces>0)
414  {
415  unsigned int i;
416 
417  for(i=0;i<mp->nFaces;i++)
418  {
419  if (mp->face[i]!=NULL)
420  free(mp->face[i]);
421  }
422  free(mp->face);
423  mp->face=NULL;
424 
425  free(mp->ID);
426  mp->ID=NULL;
427 
428  mp->nFaces=0;
429  mp->maxFaces=0;
430  }
431 }
CBLAS_INLINE void ScaleVector(double f, unsigned int s, double *v)
Scales a vector.
Definition: basic_algebra.c:30
double v
Definition: scpolytope.h:54
void SPolytopeDecreaseSamplingRadius(Tscpolytope *mp)
Decreases the sampling radious.
Definition: scpolytope.c:293
Definition of basic functions.
double ** face
Definition: scpolytope.h:58
#define FALSE
FALSE.
Definition: boolean.h:30
unsigned int GetSPolytopeDim(Tscpolytope *mp)
Returns the simple polytope dimensionality.
Definition: scpolytope.c:239
double lsr
Definition: scpolytope.h:50
#define NEW(_var, _n, _type)
Allocates memory space.
Definition: defines.h:385
void CutSPolytopeWithFace(double *t, double offset, unsigned int id, Tscpolytope *mp)
Cuts a simple polytope with a given plane.
Definition: scpolytope.c:199
unsigned int SPolytopeNumNeighbours(Tscpolytope *mp)
Number of neighbours of the simple polytope.
Definition: scpolytope.c:328
boolean RandomPointInSPolytope(double scale, double *t, Tscpolytope *mp)
Random point on the polytope with uniform distribution.
Definition: scpolytope.c:262
void CutSPolytope(double *t, double r, unsigned int id, Tscpolytope *mp)
Crops the polytope bounding chart with a plane.
Definition: scpolytope.c:186
CBLAS_INLINE double Norm(unsigned int s, double *v)
Computes the norm of a vector.
void EnlargeSPolytope(double *t, Tscpolytope *mp)
Ensures that a polytompe includes a given point.
Definition: scpolytope.c:170
unsigned int DetermineSPolytopeNeighbour(double epsilon, double *t, Tscpolytope *mp)
Identifes the neighbour containing a given point.
Definition: scpolytope.c:105
#define TRUE
TRUE.
Definition: boolean.h:21
boolean InsideSPolytope(double *t, Tscpolytope *mp)
Identifies points inside a chart simple polytope.
Definition: scpolytope.c:88
double SPolytopeVolume(Tscpolytope *mp)
Volume of the simple polytope.
Definition: scpolytope.c:305
void SaveSPolytope(FILE *f, Tscpolytope *mp)
Saves the chart polytope to a file.
Definition: scpolytope.c:351
void DefineSPolytope(Tscpolytope *mp)
Initial definition of the simple polytope bounding the local chart.
Definition: scpolytope.c:43
double SPolytopeGetSamplingRadius(Tscpolytope *mp)
Returns the current sampling radius.
Definition: scpolytope.c:281
A simpleifed polytope associated with chart on a manifold.
Definition: scpolytope.h:45
unsigned int GetSPolytopeNFaces(Tscpolytope *mp)
Number of faces of a simple chart polytope.
Definition: scpolytope.c:244
unsigned int * ID
Definition: scpolytope.h:60
void LoadSPolytope(FILE *f, Tscpolytope *mp)
Reads the chart polytope from a file.
Definition: scpolytope.c:377
CBLAS_INLINE double GeneralDotProduct(unsigned int s, double *v1, double *v2)
Computes the dot product of two general vectors.
Definition: basic_algebra.c:15
unsigned int nFaces
Definition: scpolytope.h:56
Definitions of constants and macros used in several parts of the cuik library.
double SPolytopeMaxVolume(Tscpolytope *mp)
Maximum volume of the simple polytope.
Definition: scpolytope.c:300
boolean SPolytopeRandomPointOnBoundary(double rSample, double *t, Tscpolytope *mp)
Random point on the boundary of the chart.
Definition: scpolytope.c:255
#define MOV_AVG_UP
Weight of new data when computing moving averages.
Definition: defines.h:451
double GetSPolytopeBoxSide(Tscpolytope *mp)
Returns the size of the box side.
Definition: scpolytope.c:234
double GetSPolytopeRadius(Tscpolytope *mp)
Returns the simple polytope radius.
Definition: scpolytope.c:229
unsigned int k
Definition: scpolytope.h:46
double msr
Definition: scpolytope.h:53
Definition of a smple polytope associated to a chart.
#define MEM_DUP(_var, _n, _type)
Duplicates a previously allocated memory space.
Definition: defines.h:414
#define NO_UINT
Used to denote an identifier that has not been initialized.
Definition: defines.h:435
void InitEmptySPolytope(double delta, unsigned int k, double r, double sr, Tscpolytope *mp)
Defines an empty chart simplieifed polytope.
Definition: scpolytope.c:21
void SPolytopeIncreaseSamplingRadius(Tscpolytope *mp)
Increases the sampling radius.
Definition: scpolytope.c:286
unsigned int maxFaces
Definition: scpolytope.h:57
void randomOnBall(double r, unsigned int k, double *p)
Random number on a k dimensional ball.
Definition: random.c:110
unsigned int SPolytopeNeighbourID(unsigned int n, Tscpolytope *mp)
Returns the identifier of one of the neighbours of a polytope.
Definition: scpolytope.c:333
#define MEM_EXPAND(_var, _n, _type)
Expands a previously allocated memory space.
Definition: defines.h:404
void GetSPolytopeFace(unsigned int n, double *f, Tscpolytope *mp)
Gets a face.
Definition: scpolytope.c:249
void CopySPolytope(Tscpolytope *mp_dst, Tscpolytope *mp_src)
Copies the simplified polytope from one chart to another.
Definition: scpolytope.c:54
Definition of basic randomization functions.
#define INIT_NUM_FACES
Initial space for faces.
Definition: scpolytope.h:22
void randomInBall(double r, unsigned int k, double *p)
Random number in a k dimensional ball.
Definition: random.c:126
unsigned int SPolytopeMemSize(Tscpolytope *mp)
Computes the memory used by the polytope.
Definition: scpolytope.c:341
double sr
Definition: scpolytope.h:49
double r
Definition: scpolytope.h:47
double BallVolume(unsigned int n, double r)
Volume of of a n-ball.
Definition: geom.c:682
void DeleteSPolytope(Tscpolytope *mp)
Deletes the structure allocated by DefineSPolytope.
Definition: scpolytope.c:411
#define MOV_AVG_DOWN
Weight of new data when computing moving averages.
Definition: defines.h:467