Coons曲面

原理：
http://wenku.baidu.com/link?url=w2flm9and7b1cg_gPOHAi8N-TfFxeje0FVyHab4PsCpegJ9SapcDTvb22sRiBlBKz8tk96Ha_P9-XO0dO3sRe7gIRdPDtG2HV-P64xdpMyS

#define nPoints 30   //每段曲线的点数

//矩阵乘法
void mul(float *c, float *a, float *b, int row, int rc, int col)
{
    for(int i=0; i<row; i++)
    {
        for(int j=0; j<col; j++)
        {
            float t = 0;
            for(int k=0; k<rc; k++)
                t += a[i*rc+k] * b[k*col+j];
            c[i*col+j] = t;
        }
    }
}

void init(float *p, float *u, float *v, float *uv)
{
    for(int i=0; i<4; i++)
    {
        p[i] = rand() % 500;
        u[i] = (rand() % 2 ? -1 : 1) * (rand() % 10 + 1);
        v[i] = (rand() % 2 ? -1 : 1) * (rand() % 10 + 1);
        uv[i] = (rand() % 2 ? -1 : 1) * (rand() % 10 + 1);
    }
}

void setM(float *C, float *p, float *u, float *v, float *uv)
{
    float B[][4] = 
    {
        {p[0], p[1], v[0],  v[1]},
        {p[2], p[3], v[2],  v[3]},
        {u[0], u[1], uv[0], uv[1]},
        {u[2], u[3], uv[2], uv[3]}
    };

    for(int i=0; i<4; i++)
        for(int j=0; j<4; j++)
            C[i*4+j] = B[i][j];
}

//双三次Coons曲面
void Coons(CDC *pdc)
{
    CPen newPen(PS_SOLID, 1, RGB(0, 0, 255));
    pdc->SelectObject(&newPen);

    float p1[4], pu1[4], pv1[4], puv1[4];
    float p2[4], pu2[4], pv2[4], puv2[4];
    init(p1, pu1, pv1, puv1);
    init(p2, pu2, pv2, puv2);

    float delta = 1.0 / nPoints;
    float FG[1][4], FGt[4][1], Cx[4][4], Cy[4][4];
    float tmp1[nPoints], tmp2[nPoints];

    for(float u=0; u<nPoints; u++)
    {
        float t = u * delta;
        float t1 = t;
        float t2 = t1 * t;
        float t3 = t2 * t;

        FG[0][0] = 2*t3 - 3*t2 + 1;
        FG[0][1] = -2*t3 + 3*t2;
        FG[0][2] = t3 - 2*t2 + t1;
        FG[0][3] = t3 - t2;

        setM((float *)Cx, p1, pu1, pv1, puv1);
        setM((float *)Cy, p2, pu2, pv2, puv2);

        float tmpx[1][4], tmpy[1][4], x[nPoints], y[nPoints];
        mul((float *)tmpx, (float *)FG, (float *)Cx, 1, 4, 4);
        mul((float *)tmpy, (float *)FG, (float *)Cy, 1, 4, 4);

        for(int v=0; v<nPoints; v++)
        {
            t = v * delta;
            t1 = t;
            t2 = t1 * t;
            t3 = t2 * t;

            FGt[0][0] = 2*t3 - 3*t2 + 1;
            FGt[1][0] = -2*t3 + 3*t2;
            FGt[2][0] = t3 - 2*t2 + t1;
            FGt[3][0] = t3 - t2;

            float x1[1][1], y1[1][1];
            mul((float *)x1,  (float *)tmpx, (float *)FGt, 1, 4, 1);
            mul((float *)y1,  (float *)tmpy, (float *)FGt, 1, 4, 1);

            x[v] = x1[0][0];
            y[v] = y1[0][0];
        }

        pdc->MoveTo(x[0], y[0]);
        //连线
        for(int i=1; i<nPoints; i++)
        {
            pdc->LineTo(x[i], y[i]);
        }
        //连线
        for(i=1; i<nPoints; i++)
        {
            if(u!=0)
            {
                pdc->MoveTo(tmp1[i], tmp2[i]);  
                pdc->LineTo(x[i], y[i]);
            }
            tmp1[i] = x[i];
            tmp2[i] = y[i];
        }
    }
}