Code owners
Assign users and groups as approvers for specific file changes. Learn more.
delaunay.c 7.11 KiB
/*
============================================================================
Name : constant-workspace-algorithms-c.c
Author : Mika Delor
Version :
Copyright : Your copyright notice
Description : Hello World in C, Ansi-style
============================================================================
*/
#include "delaunay.h"
int comparePointsByX (const void * a, const void * b) {
vertex* pointA = (vertex*) a;
vertex* pointB = (vertex*) b;
if(pointA->x > pointB->x)
return 1;
if(pointA->x < pointB->x)
return -1;
return 0;
}
/*int reportLine(vertex* p, int size, FILE *out) {
line l = {p[0], p[1]};
svgDrawLine(l, out);
l.a = p[0];
l.b = p[2];
svgDrawLine(l, out);
return 1;
}*/
void reportTriangle(vertex a, vertex b, vertex c) {
// qsort (p, size, sizeof(*p), comparePointsByX);
//
// char *color = "255,0,0";
#ifdef OUTPUTFILE
svgDrawTriangle(a, b, c);
#endif
// edge l = {p[0], p[1]};
// svgDrawLine(l, out, color);
// l.b = p[2];
// svgDrawLine(l, out, color);
// l.a = p[1];
// svgDrawLine(l, out, color);
}
//http://mathforum.org/library/drmath/view/54323.html
circle circle_vvv(vertex v1, vertex v2, vertex v3)
{
circle c;
double bx = v1.x; double by = v1.y;
double cx = v2.x; double cy = v2.y;
double dx = v3.x; double dy = v3.y; double temp = cx*cx+cy*cy;
double bc = (bx*bx + by*by - temp)/2.0;
double cd = (temp - dx*dx - dy*dy)/2.0;
double det = (bx-cx)*(cy-dy)-(cx-dx)*(by-cy);
if (fabs(det) < 1.0e-6) {
c.center.x = c.center.y = 1.0;
return c;
}
det = 1/det;
c.center.x = (bc*(cy-dy)-cd*(by-cy))*det;
c.center.y = ((bx-cx)*cd-(cx-dx)*bc)*det;
cx = c.center.x; cy = c.center.y;
c.radius = distance(v1, c.center);
return c;
}
double det2(double a, double b, double c, double d) {
return a*d-b*c;
}
void circleCenter(const vertex* v1, const vertex* v2, const vertex* v3,
vertex* center) {
double ax = 2 * (v1->x - v2->x);
double bx = 2 * (v1->y - v2->y);
double cx = v1->x * v1->x - v2->x * v2->x + v1->y * v1->y - v2->y * v2->y;
double ay = 2 * (v1->x - v3->x);
double by = 2 * (v1->y - v3->y);
double cy = v1->x * v1->x - v3->x * v3->x + v1->y * v1->y - v3->y * v3->y;
double det = det2(ax, bx, ay, by);
center->x = det2(cx, bx, cy, by) / det;
center->y = det2(ax, cx, ay, cy) / det;
}
/*
* https://stackoverflow.com/questions/22791951/algorithm-to-find-an-arc-its-center-radius-and-angles-given-3-points
* https://de.wikipedia.org/wiki/Cramersche_Regel
* https://de.wikipedia.org/wiki/Determinante
*/
//TODO seperate radius from this function
circle circleFrom3Points(vertex v1, vertex v2, vertex v3) {
circle c;
circleCenter(&v1, &v2, &v3, &c.center);
c.radius = sqrt((v1.x-c.center.x)*(v1.x-c.center.x)+(v1.y-c.center.y)*(v1.y-c.center.y));
return c;
}
/*void reportTriangle(struct Triangle triangle) {
printf("%f, %f - %f, %f - %f, %f\n", triangle.a.x, triangle.a.y, triangle.b.x, triangle.b.y, triangle.c.x, triangle.c.y);
}*/
void printPoint(vertex point) {
//printf("%f, %f\n", point.x, point.y);
}
double calcArea(vertex p, vertex q, vertex r) {
//return 0.5*(-q.x*p.y+r.x*p.y+p.x*q.y-r.x*q.y-p.x*r.y+q.x*r.y);
return 0.5*(p.x*(q.y-r.y)-p.y*(q.x-r.x)+(q.x*r.y-q.y*r.x));
}
/**
*
* @param p1
* @param p2
* @return the distance between p1 and p2
*/
double distance(vertex p1, vertex p2) { return hypot(p1.x - p2.x, p1.y-p2.y);
}
//https://math.stackexchange.com/questions/1232773/is-the-point-on-the-left-or-the-right-of-the-vector-in-2d-space
//https://betterexplained.com/articles/cross-product/
/**
*
* @param p
* @param q
* @param r
* @return true, if r lies right of the vector pointing from p to q, false otherwise
*/
int isRight(vertex p, vertex q, vertex r) {
vertex a;
vertex b;
a = vector(p, q);
b = vector(p, r);
return a.x*b.y-a.y*b.x > 0;
}
/**
*
* @param a Point
* @param b Point
* @return the vector pointing from a to b
*/
vertex vector(vertex a, vertex b) {
double x = a.x - b.x;
double y = a.y - b.y;
vertex p = {x, y};
return p;
}
/**
* law of cosines https://www.mathsisfun.com/algebra/trig-cosine-law.html
* other solutions https://stackoverflow.com/questions/1211212/how-to-calculate-an-angle-from-three-points
* @param q
* @param p
* @param r
* @return angle qpr (angle in p)
*/
double calcAngle(vertex q, vertex p, vertex r) {
double pq = distance(p, q);
double pr = distance(p, r);
double qr = distance(q, r);
return acos((pq*pq + pr*pr - qr*qr)/(2*pq*pr));
}
int firstDelaunayEdgeIncidentTo(int j, vertex* S, int size) {
int k;
double minDist, dist;
if(j != 0) {
k = 0;
minDist = distance(S[j], S[0]);
} else /*if(j != 1)*/ {
k = 1;
minDist = distance(S[j], S[1]);
}
//TODO erst bei j anfangen? problem: es ist nicht unbedingt die kleinste seite, wenn wir spaeter anfangen
for(int i = 0; i < size; i++) {
if(j != i) {
dist = distance(S[j], S[i]);
if(dist < minDist) {
k = i;
minDist = dist;
}
}
}
return k;
}
int arrayContains(int *jAlreadyUsed, int numberOfJ, int j){
for(int i = 0; i < numberOfJ; i++) {
if(jAlreadyUsed[i] == j) {
return 1;
}
}
return 0;
}
void clockwiseNextDelaunayEdge(const int *p, const int *q, int *r, const vertex *S, const int *size, int *right) {
//0: r not initialized, 1: r is left, 2: r is right
*right = 0;
double maxAngle;
circle c;
//area and angle into 1 var
double angle;
for(int i = 0; i < *size; i++) {
if(i == *p || i == *q) {
} else if(isRight(S[*p], S[*q], S[i])) {
//maximize signed area
//http://mathworld.wolfram.com/TriangleArea.html
if(*right < 2) {
//calc circle
*right = 2;
*r = i;
c = circleFrom3Points(S[*p], S[*q], S[*r]);
} else {
//if inside circle: calculate new circle, else: do nothing
if(distance(c.center, S[i]) < c.radius) {
*r = i;
c = circleFrom3Points(S[*p], S[*q], S[*r]);
}
}
} else if(*right < 2) {
//maximize angle qpr
angle = calcAngle(S[*q], S[*p], S[i]);
//TODO is c 'or' lazy, if yes combine if and if else
//TODO only return edge and not entire triangle, if left
if(*right < 1) {
*right = 1;
maxAngle = angle;
*r = i;
} else if(angle > maxAngle){
maxAngle = angle;
*r = i;
}
}
}
//if you want to draw the circle comment this line out
//svgDrawCircle(c, out);
}
void delaunayCws(vertex* S, int size) {
//struct Triangle delaunayTriangles[size-2];
for(int i = 0; i < size; i++){
//TODO eventuell nur j zurueckgeben (i haben wir ja schon)
int j = firstDelaunayEdgeIncidentTo(i, S, size);
// Find the point p[j] e S that is nearest to p[i]
int j0 = j;
do {
//TODO eventuell nur k zurueckgeben (i haben wir ja schon)
int k;
int isTriangle;
clockwiseNextDelaunayEdge(&i, &j, &k, S, &size, &isTriangle);
if(i < j && i < k) { if(isTriangle == 2){
#ifdef OUTPUTFILE
reportTriangle(S[i], S[j], S[k]);
#endif
}
}
j = k;
} while(j != j0);
}
}
//https://stackoverflow.com/questions/29762048/c-structure-to-string
char* vertexToString(vertex v) {
int len = 0;
len = snprintf(NULL, len, "%f,%f", v.x, v.y);
char* str = malloc(sizeof(char) * len + 1);
snprintf(str, sizeof(str), "%f,%f", v.x, v.y);
return str;
}