程序 构造平面三角形类，并使该类具有判断任意两三角形是否相交

<pre t="code" l="cpp">class Vector2D { //2D向量类

public:

double x, y

Vector2D(double x_ = 0, double y_ = 0): x(x_), y(y_) {}

}

bool isZero(const Vector2D vec, double tolerance = 0.05) const {

return std::abs(vec.x) <tolerance std::abs(vec.y) <tolerance

}

Vector2D tripleProduct(const Vector2D v1, const Vector2D v2, const Vector2D v3) { //三角乘积

double p = v1.x * v2.y - v1.y * v2.x

return Vector2D(-p * v3.y, p * v3.x)

}

double dot(const Vector2D v1, cosnt Vector2D v2) { //点乘

return v1.x * v2.x + v1.y * v2.y

}

Vector2D negate(const Vector2D vec) { //向量取反

return Vector2D(-vec.x, -vec.y)

}

Vector2D sub(const Vector2D v1, const Vector2D v2) { //向量相减

return Vector2D(v1.x - v2.x, v1.y - v2.y)

}

Vector2D getFarthestPoint(const Vector2D dir, const Vector2D poly[], int len) { //求一个凸多边形poly在某个方向dir上最远的一点

Vector2D bestPoint = poly[0]

double bestProj = dot(bestPoint, dir)

for (int i = 1i <len++i) {

Vector2D curr = poly[i]

double proj = dot(curr, dir)

if (proj >bestProj) {

bestPoint = curr

bestProj = proj

}

}

return bestPoint

}

Vector2D getSupportPoint(const Vector2D dir, const Vector2D poly1[], int len1, const Vector2D poly2[], int len2) { //GJK算法的一部分, 求两个凸多边形的闵可夫斯基差在某个方向dir上最远的一点

Vector2D v1 = getFarthestPoint(dir, poly1, len1),

v2 = getFarthestPoint(negate(dir), poly2, len2)

return sub(v1, v2)

}

bool isIntersect(const Vector2D poly1[], int len1, const Vector2D poly2[], int len2) {

Vector2D simplexA, simplexB, simplexC, dir = Vector2D(-1, -1)

simplexA = getSupportPoint(dir, poly1, len1, poly2, len2)

if (dot(simplexA, dir) <= 0) {

return false

}

dir = negate(simplexA)

simplexB = getSupportPoint(dir, poly1, len1, poly2, len2)

if (dot(simplexB, dir) <= 0) {

return false

}

Vector2D ab = sub(simplexB, simplexA)

dir = tripleProduct(ab, negate(simplexA), ab)

for (int i = 25i--) { //25是最大迭代次数, 一般迭代3次以内就能完成

if (isZero(dir)) {

return true

}

simplexC = getSupportPoint(dir, poly1, len1, poly2, len2)

if (dot(simplexC, dir) <= 0) {

return false

}

Vector2D ba = sub(simplexA, simplexB),

ac = sub(simplexC, simplexA),

bc = sub(simplexC, simplexB),

acPerp = tripleProduct(ac, negate(ba), ac),

bcPerp = tripleProduct(bc, ba, bc)

if (dot(acPerp, simplexA) >0) {

simplexB = simplexC

dir = negate(acPerp)

} else if (dot(bcPerp, simplexB) >0) {

simplexA = simplexC

dir = negate(bcPerp)

} else {

return true

}

}

return false

}

这个是我把自己的代码里的一段截出来改了下变量名和参数什么的, 不确定有没哪里漏了改的, 如果运行有错误可以自己修一下

isIntersect可以判断任意两个凸多边形是否相交, 用的GJK算法三角形和矩形存在一个Vector2D数组里, 用isIntersect(poly1, poly1len, poly2, poly2len)就能判断

改一下getFarthestPoint还能变成任意两个凸图形是否相交

用c语言编写单纯形法的程序怎么写

#include<stdio.h>

#include<math.h>

intm //记录约束条件方程组的个数

intn //记录未知量的个数

floatM=1000000.0

floatA[100][100] //用于记录方程组的数目和系数

floatC[100] //用于存储目标函数中各个变量的系数

floatb[100] //用于存储常约束条件中的常数

floatCB[100]//用于存储基变量的系数

floatseta[100] //存放出基与入基的变化情况

floatcn[100] //存储检验数矩阵

floatx[100]

intnum[100] //用于存放出基与进基变量的情况

floatZ=0 //记录目标函数值

voidshuru()

voidprint()

intmincz()

intfind_line(int a)

voidexchange(int a,int b)

intmain()

{

int i,j=0

int p,q,temp //q：换入，p：换出

shuru()

printf("\n--------------------------------------------------------------------------\n")

printf(" \tCB\tXB\tb\t")

for(i=0i<ni++)

printf(" X(%d)\t",i+1)

for(i=0i<ni++)

x[i]=0

printf("\n")

while(1) {

q=mincz()

if(q==-1) {

print()

printf("\n所得解已经是最优解!\n")

printf("\n最优解为:\n")

for(j=0j<mj++) {

temp=num[j]-1

x[temp]=b[j]

}

for(i=0i<ni++) {

printf("x%d=%.2f",i+1,x[i])

Z=Z+x[i]*C[i]

}

printf("Z=%.2f",Z)

break

}

print()

p=find_line(q)

printf("\np=%d,q=%d",p,q)

if(q==-1) break

exchange(p,q)

}

return 0

}

intmincz()

{

int i,k=0

int flag=0//检验数标记

float min=0

for(i=0i<ni++)

if(cn[i]>=0)

flag=1

else {

flag=0

break

}

if(flag==1)

return -1

//进行到此处，说明存在<0的检验数

//找到最小的检验数，作为换入变量

for(i=0i<ni++) {

if(min>cn[i]) {

min=cn[i]

k=i

}

}

return k

}

intfind_line(int a)

{

int i,k,j

int flag=0

float min

k=a

for(i=0i<mi++)

if(A[i][k]<=0)

flag=1

else {

flag=0

break

}

if(flag==1) {

printf("\n该线性规划无最优解!\n")

return -1

}

for(i=0i<mi++) {

if(A[i][k]>0)

seta[i]=b[i]/A[i][k]

else seta[i]=M

}

min=M

for(i=0i<mi++) {

if(min>=seta[i]) {

min=seta[i]

j=i

}

}

num[j]=k+1

CB[j]=C[k]

return j

}

voidexchange(int p,int q)

{

int i,j,c,l

float temp1,temp2,temp3

c=p//行号，换出

l=q//列号，换入

temp1=A[c][l] //A[c][l]主元

b[c]=b[c]/temp1

for(j=0j<nj++)

A[c][j]=A[c][j]/temp1 //主元化为1

for(i=0i<mi++) {

if(i!=c)

if(A[i][l]!=0) {

temp2=A[i][l]

b[i]=b[i]-b[c]*temp2

//主元所在列，其余元素化为0

for(j=0j<nj++)

A[i][j]=A[i][j]-A[c][j]*temp2

}

}

temp3=cn[l]

for(i=0i<ni++)

cn[i]=cn[i]-A[c][i]*temp3

}

voidprint()

{

int i,j=0

printf("\n--------------------------------------------------------------------------\n")

for(i=0i<mi++) {

printf("%8.2f\tX(%d) %8.2f",CB[i],num[i],b[i])

for(j=0j<nj++)

printf("%8.2f ",A[i][j])

printf("\n")

}

printf("\n--------------------------------------------------------------------------\n")

printf("\t\t\t")

for(i=0i<ni++)

printf(" %8.2f",cn[i])

printf("\n--------------------------------------------------------------------------\n")

}

voidshuru()

{

int i,j//循环变量

int k

printf("请输入线性规划问题的约束条件个数M：")

scanf("%d",&m)

printf("请输入线性规划问题的决策变量个数N：")

scanf("%d",&n)

printf("\n请输入方程组的系数矩阵A(%d行%d列):\n",m,n)

for(i=0i<mi++)

for(j=0j<nj++)

scanf("%f",&A[i][j])

printf("\n请输入初始基变量的数字代码矩阵:\n")

for(i=0i<mi++)

scanf("%d",&num[i])

printf("\n请输入方程组右边的值矩阵b:\n")

for(i=0i<mi++)

scanf("%f",&b[i])

printf("\n请输入目标函数各个变量的系数所构成的系数阵C:\n")

for(i=0i<ni++)

scanf("%f",&C[i])

for(i=0i<ni++)

cn[i]=-C[i]

for(i=0i<mi++) {

k=num[i]-1

CB[i]=C[k]

}

}

class Vector2D { //2D向量类

public:

    double x, y

    Vector2D(double x_ = 0, double y_ = 0): x(x_), y(y_) {}

}

bool isZero(const Vector2D &vec, double tolerance = 0.05) const {

            return std::abs(vec.x) < tolerance && std::abs(vec.y) < tolerance

}

Vector2D tripleProduct(const Vector2D &v1, const Vector2D &v2, const Vector2D &v3) { //三角乘积

    double p = v1.x * v2.y - v1.y * v2.x

    return Vector2D(-p * v3.y, p * v3.x)

}

double dot(const Vector2D &v1, cosnt Vector2D &v2) { //点乘

    return v1.x * v2.x + v1.y * v2.y

}

Vector2D negate(const Vector2D &vec) { //向量取反

    return Vector2D(-vec.x, -vec.y)

}

Vector2D sub(const Vector2D &v1, const Vector2D &v2) { //向量相减

    return Vector2D(v1.x - v2.x, v1.y - v2.y)

}

Vector2D getFarthestPoint(const Vector2D &dir, const Vector2D poly[], int len) { //求一个凸多边形poly在某个方向dir上最远的一点

    Vector2D bestPoint = poly[0]

    double bestProj = dot(bestPoint, dir)

    for (int i = 1 i < len ++i) {

        Vector2D curr = poly[i]

        double proj = dot(curr, dir)

        if (proj > bestProj) {

            bestPoint = curr

            bestProj = proj

        }

    }

    return bestPoint

}

Vector2D getSupportPoint(const Vector2D &dir, const Vector2D poly1[], int len1, const Vector2D poly2[], int len2) { //GJK算法的一部分, 求两个凸多边形的闵可夫斯基差在某个方向dir上最远的一点

    Vector2D v1 = getFarthestPoint(dir, poly1, len1),

             v2 = getFarthestPoint(negate(dir), poly2, len2)

    return sub(v1, v2)

}

bool isIntersect(const Vector2D poly1[], int len1, const Vector2D poly2[], int len2) {

    Vector2D simplexA, simplexB, simplexC, dir = Vector2D(-1, -1)

    simplexA = getSupportPoint(dir, poly1, len1, poly2, len2)

    if (dot(simplexA, dir) <= 0) {

        return false

    }

    dir = negate(simplexA)

    simplexB = getSupportPoint(dir, poly1, len1, poly2, len2)

    if (dot(simplexB, dir) <= 0) {

        return false

    }

    Vector2D ab = sub(simplexB, simplexA)

    dir = tripleProduct(ab, negate(simplexA), ab)

    for (int i = 25 i--) { //25是最大迭代次数, 一般迭代3次以内就能完成

        if (isZero(dir)) {

            return true

        }

        simplexC = getSupportPoint(dir, poly1, len1, poly2, len2)

        if (dot(simplexC, dir) <= 0) {

            return false

        }

        Vector2D ba = sub(simplexA, simplexB),

                ac = sub(simplexC, simplexA),

                bc = sub(simplexC, simplexB),

                acPerp = tripleProduct(ac, negate(ba), ac),

                bcPerp = tripleProduct(bc, ba, bc)

        if (dot(acPerp, simplexA) > 0) {

            simplexB = simplexC

            dir = negate(acPerp)

        } else if (dot(bcPerp, simplexB) > 0) {

            simplexA = simplexC

            dir = negate(bcPerp)

        } else {

            return true

        }

    }

    return false

}

这个是我把自己的代码里的一段截出来改了下变量名和参数什么的, 不确定有没哪里漏了改的, 如果运行有错误可以自己修一下

isIntersect可以判断任意两个凸多边形是否相交, 用的GJK算法

三角形和矩形存在一个Vector2D数组里, 用isIntersect(poly1, poly1len, poly2, poly2len)就能判断

改一下getFarthestPoint还能变成任意两个凸图形是否相交

---