struct FPOINT { double x; double y; }; struct LINESEGMENT { struct FPOINT a; struct FPOINT b; }; #include #define LineTest(one, two) (-((one.b.y - two.b.y)*(two.a.x - two.b.x) - (one.b.x - two.b.x)*(two.a.y - two.b.y))/((one.a.y - one.b.y)*(two.a.x - two.b.x) - (one.a.x - one.b.x)*(two.a.y - two.b.y))) #define Parallel(one, two) (((one.a.y - one.b.y)*(two.a.x - two.b.x) - (one.a.x - one.b.x)*(two.a.y - two.b.y)) == 0.0) #define GetX(line, A) (line.a.x * A + line.b.x * (1.0 - A)) #define GetY(line, A) (line.a.y * A + line.b.y * (1.0 - A)) #define DoubleEquals(one, two) (fabs(one - two) < 0.0000000001) #define LineEquals(one, two) ((DoubleEquals(one.a.x, two.a.x) && DoubleEquals(one.a.y, two.a.y) && DoubleEquals(one.b.x, two.b.x) && DoubleEquals(one.b.y, two.b.y)) || (DoubleEquals(one.a.x, two.b.x) && DoubleEquals(one.a.y, two.b.y) && DoubleEquals(one.b.x, two.a.x) && DoubleEquals(one.b.y, two.a.y))) #define PointEquals(one, two) (DoubleEquals(one.x, two.x) && DoubleEquals(one.y, two.y)) #define LineConnects(one, two) (PointEquals(one.a, two.a) || PointEquals(one.a, two.b) || PointEquals(one.b, two.a) || PointEquals(one.b, two.b)) #define Length(one, two) (sqrt((one.x - two.x) * (one.x - two.x) + (one.y - two.y) * (one.y - two.y))) void set_alphabeta(double *Quad_x, double *Quad_y, double *alpha, double *beta) { alpha[0] = Quad_x[0]; alpha[1] = Quad_x[1] - Quad_x[0]; alpha[2] = Quad_x[3] - Quad_x[0]; alpha[3] = (Quad_x[0] + Quad_x[2]) - (Quad_x[1] + Quad_x[3]); beta[0] = Quad_y[0]; beta[1] = Quad_y[1] - Quad_y[0]; beta[2] = Quad_y[3] - Quad_y[0]; beta[3] = (Quad_y[0] + Quad_y[2]) - (Quad_y[1] + Quad_y[3]); } int xy_n_quad(double *Quad_x, double *Quad_y, double *xy) { double aa, bb, cc; double l, m; double alpha[4], beta[4]; int unity_map; unity_map = (Quad_x[0] == 0. && Quad_x[1] == 1. && Quad_x[2] == 1. && Quad_x[3] == 0.); if( unity_map ) unity_map = (Quad_y[0] == 0. && Quad_y[1] == 0. && Quad_y[2] == 1. && Quad_y[3] == 1.); if( !unity_map ) { set_alphabeta(Quad_x, Quad_y, alpha, beta); aa = alpha[3]*beta[2] - alpha[2]*beta[3]; bb = alpha[3]*beta[0] + alpha[1]*beta[2] + xy[0]*beta[3]; bb -= (alpha[0]*beta[3] + alpha[2]*beta[1] + xy[1]*alpha[3]); cc = alpha[1]*beta[0] + xy[0]*beta[1] - (alpha[0]*beta[1] + xy[1]*alpha[1]); m = .5*(sqrt(bb*bb - 4.*aa*cc) -bb)/aa; if( m < 0.0 || m > 1.0 ) return 0; else { l = (xy[0] - (alpha[0] + alpha[2]*m))/(alpha[1] + alpha[3]*m); if( l < 0.0 || l > 1.0 ) return 0; else return 1; } } else { if( xy[1] < 0.0 || xy[1] > 1.0 ) return 0; else if( xy[0] < 0.0 || xy[0] > 1.0 ) return 0; else return 1; } } double poly_area(int nv, double *x, double *y) { int i,im1,ip1; double accum; accum = 0.0; for(i = 0; i < nv; i++ ) { ip1 = (i+1)%nv; im1 = i - 1; if( im1 < 0 ) im1 = nv - 1; accum += (x[ip1] - x[im1])*y[i]; } return( -.5*accum ); } double quadintersectarea_(double *Quad1_x, double *Quad1_y, double *Quad2_x, double *Quad2_y) { //----------------------------------------------------------- // Convert from points to line segments... //----------------------------------------------------------- struct LINESEGMENT quad[2][4]; struct LINESEGMENT shapes[16]; int shapeIndex = 0; int i, n, p, r; int cnt; int mask[16]; //----------------------------------------------------------- // first a few diagnostics //----------------------------------------------------------- { double area; area = poly_area(4, Quad1_x, Quad1_y); area = poly_area(4, Quad2_x, Quad2_y); } //----------------------------------------------------------- // check for complete inclusion //----------------------------------------------------------- p = 0; r = 0; for (int i = 0; i < 4; i++) // Quad one in Quad two { double vertex[2]; vertex[0] = Quad1_x[i]; vertex[1] = Quad1_y[i]; n = xy_n_quad(Quad2_x, Quad2_y, vertex); if( n == 0 ) break; else p += n; } if( p != 4 ) { p = 0; for (int i = 0; i < 4; i++) // Quad two in Quad one { double vertex[2]; vertex[0] = Quad2_x[i]; vertex[1] = Quad2_y[i]; n = xy_n_quad(Quad1_x, Quad1_y, vertex); if( n == 0 ) break; else p += n; } if( p == 4 ) r = 1; } if( p == 4 ) { double area; if( r == 0 ) area = poly_area(p, Quad1_x, Quad1_y); else area = poly_area(p, Quad2_x, Quad2_y); return area; } for (int i = 0; i < 4; i++) // Quad one { quad[0][i].a.x = Quad1_x[i]; //Point A quad[0][i].a.y = Quad1_y[i]; quad[0][i].b.x = Quad1_x[(i + 1) % 4]; //Point B quad[0][i].b.y = Quad1_y[(i + 1) % 4]; quad[1][i].a.x = Quad2_x[i]; //Point A quad[1][i].a.y = Quad2_y[i]; quad[1][i].b.x = Quad2_x[(i + 1) % 4]; //Point B quad[1][i].b.y = Quad2_y[(i + 1) % 4]; } //----------------------------------------------------------- // Get the intersection line segments for the quads... //----------------------------------------------------------- for (int quadIndex = 0; quadIndex < 2; quadIndex++) { for (int i = 0; i < 4; i++) { struct FPOINT intersect[4]; double aArray[4]; int index = 0; //----------------------------------------------------------- // Check intersctions for every line in other quad... //----------------------------------------------------------- for (int n = 0; n < 4; n++) { if (!Parallel(quad[1 - quadIndex][n], quad[quadIndex][i])) { double a = LineTest(quad[1 - quadIndex][n], quad[quadIndex][i]); if (a >= 0.0 && a <= 1.0) { a = LineTest(quad[quadIndex][i], quad[1 - quadIndex][n]); intersect[index].x = GetX(quad[quadIndex][i], a); intersect[index].y = GetY(quad[quadIndex][i], a); aArray[index] = a; index++; } } } if (index == 0) continue; //----------------------------------------------------------- // Remove point doubles... //----------------------------------------------------------- for (int p = 0; p < index; p++) mask[p] = 1; for (int p = 0; p < index-1; p++) { if( mask[p] ) { for (int r = p+1; r < index; r++) if (mask[r] && PointEquals(intersect[r], intersect[p])) mask[r] = 0; } } cnt = 0; for (int p = 0; p < index; p++) if( mask[p] ) cnt++; index = cnt; if (index == 0) continue; //----------------------------------------------------------- // Sort the lists smallest a value to largest a value... //----------------------------------------------------------- struct FPOINT sortedIntersect[4]; sortedIntersect[0]; double sortedA[4]; sortedA[0]; for (int n = 0; n < index; n++) { int inum; for (inum = 0; inum < n; inum++) { if (aArray[n] < sortedA[inum]) { //----------------------------------------------------------- // shift list to fit segment in correct place... //----------------------------------------------------------- for (int p = n; p > inum; p--) { sortedA[p] = sortedA[p - 1]; sortedIntersect[p] = sortedIntersect[p - 1]; } sortedA[inum] = aArray[n]; sortedIntersect[inum] = intersect[n]; break; } } if (inum == n) { sortedA[inum] = aArray[n]; sortedIntersect[inum] = intersect[n]; } } //----------------------------------------------------------- // create line semgents and clip line segments to fit current quad and add to a list... //----------------------------------------------------------- for (int n = 0; n < index; n+=2) { if (sortedA[n] < 0.0) { if (sortedA[n + 1] > 1.0) { shapes[shapeIndex] = quad[quadIndex][i]; shapeIndex++; } else if (sortedA[n + 1] >= 0.0) { shapes[shapeIndex].a.x = GetX(quad[quadIndex][i], sortedA[n + 1]); shapes[shapeIndex].a.y = GetY(quad[quadIndex][i], sortedA[n + 1]); shapes[shapeIndex].b = quad[quadIndex][i].b; shapeIndex++; } } else if (sortedA[n] > 1.0) { if (sortedA[n + 1] <= 1.0 && sortedA[n + 1] >= 0.0) { shapes[shapeIndex].b.x = GetX(quad[quadIndex][i], sortedA[n + 1]); shapes[shapeIndex].b.y = GetY(quad[quadIndex][i], sortedA[n + 1]); shapes[shapeIndex].a = quad[quadIndex][i].a; shapeIndex++; } } else { if (sortedA[n + 1] > 1.0) { shapes[shapeIndex].b.x = GetX(quad[quadIndex][i], sortedA[n]); shapes[shapeIndex].b.y = GetY(quad[quadIndex][i], sortedA[n]); shapes[shapeIndex].a = quad[quadIndex][i].a; shapeIndex++; } else if (sortedA[n + 1] < 0.0) { shapes[shapeIndex].a.x = GetX(quad[quadIndex][i], sortedA[n]); shapes[shapeIndex].a.y = GetY(quad[quadIndex][i], sortedA[n]); shapes[shapeIndex].b = quad[quadIndex][i].b; shapeIndex++; } else { shapes[shapeIndex].a.x = GetX(quad[quadIndex][i], sortedA[n]); shapes[shapeIndex].a.y = GetY(quad[quadIndex][i], sortedA[n]); shapes[shapeIndex].b.x = GetX(quad[quadIndex][i], sortedA[n + 1]); shapes[shapeIndex].b.y = GetY(quad[quadIndex][i], sortedA[n + 1]); shapeIndex++; } } } } } //------------------------------------------------------ // if there are no intersections area is zero... //------------------------------------------------------ if (shapeIndex == 0) return 0.0; //------------------------------------------------------ // Remove lines with zero length //------------------------------------------------------ for (i = 0; i < shapeIndex; i++) mask[i] = 0; for (i = 0; i < shapeIndex; i++) if (PointEquals(shapes[i].a, shapes[i].b)) mask[i] = 1; int nshapes = shapeIndex; r = 0; for (p = 0; p < nshapes; p++) { if( mask[p] == 0 ) shapes[r++] = shapes[p]; else shapeIndex--; } //------------------------------------------------------ // Remove line doubles //------------------------------------------------------ for (i = 0; i < shapeIndex; i++) { for (n = 0; n < shapeIndex; n++) mask[n] = 0; for (n = 0; n < shapeIndex; n++) { if (LineEquals(shapes[i], shapes[n]) && (i != n)) { mask[n] = 1; } } int nshapes = shapeIndex; r = 0; for (p = 0; p < nshapes; p++) { if( mask[p] == 0 ) shapes[r++] = shapes[p]; else shapeIndex--; } } //------------------------------------------------------ // if there are no intersections area is zero... //------------------------------------------------------ if (shapeIndex == 0) return 0.0; //------------------------------------------------------ // find adjacent lines to points... //------------------------------------------------------ struct FPOINT points[16]; struct FPOINT lineadjacency[16]; int pointNumber = 0; for (int i = 0; i < shapeIndex; i++) { int alreadyDone = 0; for (int n = 0; n < pointNumber; n++) { if (PointEquals(points[n], shapes[i].a)) alreadyDone = 1; } if (!alreadyDone) { lineadjacency[pointNumber].x = i; points[pointNumber] = shapes[i].a; for (int n = 0; n < shapeIndex; n++) { if (n != i && (PointEquals(points[pointNumber], shapes[n].a) || PointEquals(points[pointNumber], shapes[n].b))) lineadjacency[pointNumber].y = n; } pointNumber++; } alreadyDone = 0; for (int n = 0; n < pointNumber; n++) { if (PointEquals(points[n], shapes[i].b)) alreadyDone = 1; } if (!alreadyDone) { lineadjacency[pointNumber].x = i; points[pointNumber] = shapes[i].b; for (int n = 0; n < shapeIndex; n++) { if (n != i && (PointEquals(points[pointNumber], shapes[n].a) || PointEquals(points[pointNumber], shapes[n].b))) lineadjacency[pointNumber].y = n; } pointNumber++; } } //------------------------------------------------------ // find connected points... //------------------------------------------------------ struct FPOINT adjacency[16]; int lineUsed[16]; for( int i = 0; i < 16; i++) lineUsed[i] = 0; for (int i = 0; i < pointNumber; i++) { if (!lineUsed[(int)lineadjacency[i].x]) { for (int n = 0; n < pointNumber; n++) { if (n != i) { if (lineadjacency[i].x == lineadjacency[n].x) { adjacency[i].x = n; adjacency[n].x = i; lineUsed[(int)lineadjacency[i].x] = 1; break; } else if (lineadjacency[i].x == lineadjacency[n].y) { adjacency[i].x = n; adjacency[n].y = i; lineUsed[(int)lineadjacency[i].x] = 1; break; } } } } if (!lineUsed[(int)lineadjacency[i].y]) { for (int n = 0; n < pointNumber; n++) { if (n != i) { if (lineadjacency[i].y == lineadjacency[n].x) { adjacency[i].y = n; adjacency[n].x = i; lineUsed[(int)lineadjacency[i].y] = 1; break; } else if (lineadjacency[i].y == lineadjacency[n].y) { adjacency[i].y = n; adjacency[n].y = i; lineUsed[(int)lineadjacency[i].y] = 1; break; } } } } } //------------------------------------------------------ // spit up complex figure into smaller triangles... //------------------------------------------------------ struct FPOINT triangles[16][3]; int triUsed[16]; for( int i = 0; i < 16; i++ ) triUsed[i] = 0; int currentTri = 0; int pointsRemain = pointNumber; while (pointsRemain > 2) { for (int i = 0; i < pointNumber; i++) { if (!triUsed[i]) { int n; for (n = 0; n < 2; n++) { struct LINESEGMENT test; test.a = points[i]; if (n == 0) test.b = points[(int)adjacency[i].x]; else test.b = points[(int)adjacency[i].y]; int less = 0; int greater = 0; for (int p = 0; p < shapeIndex; p++) { if (!Parallel(test, shapes[p])) { double a = LineTest(test, shapes[p]); if (a < 0.0) less++; if (a > 1.0) greater++; } } if (less % 2 == 1 || greater % 2 == 1) { // concave section, trainge not in shape so don't use it n = 3; continue; } } if (n == 3) continue; struct LINESEGMENT line; line.a = points[(int)adjacency[i].x]; line.b = points[(int)adjacency[i].y]; double a; for (n = 0; n < shapeIndex; n++) { if (!Parallel(line, shapes[n])) a = LineTest(line, shapes[n]); if (a > 0.0 && a < 1.0) n = shapeIndex + 1; } if (n == shapeIndex + 1) continue; // traingle intersects shape outline, don't use it //add triangle to list... triangles[currentTri][0] = points[i]; triangles[currentTri][1] = points[(int)adjacency[i].x]; triangles[currentTri][2] = points[(int)adjacency[i].y]; currentTri++; pointsRemain--; triUsed[i] = 1; // remove point from list and readjust adjacencies... if (adjacency[(int)adjacency[i].x].x == i) { if (adjacency[(int)adjacency[i].y].x == i) { adjacency[(int)adjacency[i].x].x = adjacency[i].y; adjacency[(int)adjacency[i].y].x = adjacency[i].x; } else { adjacency[(int)adjacency[i].x].x = adjacency[i].y; adjacency[(int)adjacency[i].y].y = adjacency[i].x; } } else { if (adjacency[(int)adjacency[i].y].x == i) { adjacency[(int)adjacency[i].x].y = adjacency[i].y; adjacency[(int)adjacency[i].y].x = adjacency[i].x; } else { adjacency[(int)adjacency[i].x].y = adjacency[i].y; adjacency[(int)adjacency[i].y].y = adjacency[i].x; } } break; } } } //------------------------------------------------------------------ // loop through all of the triangles created and add up the area... //------------------------------------------------------------------ double area = 0.0; //sum of the area of the triagles for (int i = 0; i < currentTri; i++) { double a = Length(triangles[i][0], triangles[i][1]); //sides lengths of the triagle double b = Length(triangles[i][1], triangles[i][2]); double c = Length(triangles[i][0], triangles[i][2]); double s = .5 * (a + b + c); // Heron's area of an SSS triagle formula if( s != 0. ) area += sqrt(s * (s - a) * (s - b) * (s - c)); // calculate area and add it to the sum of areas } return area; }