packages/services/gfx/mw/v2_0/src/mwin/winlib/graph3d.c

   1 /*
   2  * Copyright (c) 1999 Greg Haerr <greg@censoft.com>
   3  *
   4  * 3D Graphics Library for Micro-Windows
   5  */
   6 #define MWINCLUDECOLORS
   7 #include "windows.h"
   8 #include "device.h"
   9 #include "graph3d.h"
  10 #define USEBLIT         1       /* =1 to use memDC's*/
  11
  12 static int      nxpix;
  13 static int      nypix;
  14 static vec1     xscale;
  15 static vec1     yscale;
  16 static vec3     eye;
  17 static vec3     direct;
  18 static double   Q[5][5];
  19 static HDC      hdc;
  20 static HDC      hdcMem;
  21 static HBITMAP  hbmp, hbmpOrg;
  22
  23 /* setup eye, direction, calc observation matrix Q*/
  24 void
  25 look3(vec1 x, vec1 y, vec1 z)
  26 {
  27         eye.x = x;
  28         eye.y = y;
  29         eye.z = z;
  30         direct.x = -eye.x;
  31         direct.y = -eye.y;
  32         direct.z = -eye.z;
  33         findQ();
  34 }
  35
  36 void
  37 init3(HDC hDC, HWND memhwnd)
  38 {
  39         HBRUSH  hbr;
  40
  41         hdc = hDC;
  42         if(hdc) {
  43                 nxpix = hdc->hwnd->clirect.right - hdc->hwnd->clirect.left;
  44                 nypix = hdc->hwnd->clirect.bottom - hdc->hwnd->clirect.top;
  45                 xscale = (vec1)(nxpix-1) / nxpix * nxpix/2;
  46                 yscale = (vec1)(nypix-1) / nypix * nypix/2;
  47
  48                 if(memhwnd) {
  49                         hdcMem = CreateCompatibleDC(NULL);
  50                         if(hdcMem) {
  51                                 hbmp = CreateCompatibleBitmap(hdcMem,
  52                                         nxpix, nypix);
  53                                 hbmpOrg = SelectObject(hdcMem, hbmp);
  54                                 hdc = hdcMem;
  55                         }
  56                         hbr = (HBRUSH)GetClassLong(memhwnd, GCL_HBRBACKGROUND);
  57                         FillRect(hdc, NULL, hbr);
  58                 }
  59                 /* create pen for setcolor3() color override*/
  60                 SelectObject(hdc, CreatePen(PS_SOLID, 1, BLACK));
  61         }
  62 }
  63
  64 void
  65 paint3(HDC hDC)
  66 {
  67         if(hdcMem) {
  68                 BitBlt(hDC, 0, 0, nxpix, nypix, hdcMem, 0, 0, SRCCOPY);
  69                 DeleteObject(SelectObject(hdcMem, hbmpOrg));
  70                 DeleteDC(hdcMem);
  71         }
  72         hdcMem = NULL;
  73 }
  74
  75 int
  76 fx(vec1 x)
  77 {
  78         return (int)(x * xscale + nxpix*0.5 - 0.5);
  79 }
  80
  81 int
  82 fy(vec1 y)
  83 {
  84         return (int)(y * yscale + nypix*0.5 - 0.5);
  85 }
  86
  87 void
  88 moveto3(vec2 pt)
  89 {
  90         MoveToEx(hdc, fx(pt.x), fy(pt.y), NULL);
  91 }
  92
  93 void
  94 setcolor3(MWCOLORVAL c)
  95 {
  96         if(hdc)
  97                 hdc->pen->color = c;
  98 }
  99
 100 void
 101 lineto3(vec2 pt)
 102 {
 103         LineTo(hdc, fx(pt.x), fy(pt.y));
 104 }
 105
 106 void
 107 polyfill(int n, vec2 points[])
 108 {
 109         int     i;
 110         int     xoff, yoff;
 111         MWPOINT pv[MAXPOLY];
 112
 113         if(!hdc)
 114                 return;
 115
 116         /* calc window offset*/
 117         xoff = hdc->hwnd->clirect.left;
 118         yoff = hdc->hwnd->clirect.top;
 119
 120         /* only plot non-trivial polygons*/
 121         if(n > 2) {
 122                 for(i=0; i<n; ++i) {
 123                         pv[i].x = fx(points[i].x) + xoff;
 124                         pv[i].y = fy(points[i].y) + yoff;
 125                         /* fix: floating round error, y intercept difference
 126                          * with GdLine
 127                          */
 128                         /*pv[i].x = fx(points[i].x + xoff);*/
 129                         /*pv[i].y = fy(points[i].y + yoff);*/
 130                 }
 131                 GdSetForeground(GdFindColor(hdc->pen->color));
 132                 GdFillPoly(hdc->psd, n, pv);
 133         }
 134 }
 135
 136 void
 137 square(void)
 138 {
 139         vec2    pt0, pt1, pt2, pt3;
 140
 141         pt0.x = -1; pt0.y = 1;
 142         pt1.x = -1; pt1.y = -1;
 143         pt2.x = 1; pt2.y = -1;
 144         pt3.x = 1; pt3.y = 1;
 145         moveto3(pt0);
 146         lineto3(pt1);
 147         lineto3(pt2);
 148         lineto3(pt3);
 149         lineto3(pt0);
 150 }
 151
 152 void
 153 circle3(vec1 r)
 154 {
 155         vec1    theta = 0;
 156         vec1    thinc = 2*pi/100;
 157         int     i;
 158         vec2    pt;
 159
 160         pt.x = r;
 161         pt.y = 0.0;
 162         moveto3(pt);
 163
 164         for(i=0; i<100; ++i) {
 165                 theta = theta + thinc;
 166                 pt.x = r*cos(theta);
 167                 pt.y = r*sin(theta);
 168                 lineto3(pt);
 169         }
 170 }
 171
 172 void
 173 daisy(vec1 r,int points)
 174 {
 175         int     i, j;
 176         vec1    theta = 0;
 177         vec1    thinc;
 178         vec2    pt[100];
 179
 180         /* calculate n points on a circle*/
 181         thinc = 2*pi/points;
 182         for(i=0; i<points; ++i) {
 183                 pt[i].x = r*cos(theta);
 184                 pt[i].y = r*sin(theta);
 185                 theta += thinc;
 186         }
 187
 188         /* join point i to point j for all 0 <= i < j < n */
 189         for(i=0; i<points-1; ++i) {
 190                 for(j=i+1; j<points; ++j) {
 191                         moveto3(pt[i]);
 192                         lineto3(pt[j]);
 193                 }
 194         }
 195 }
 196
 197 void
 198 rose(vec1 r,int levels,int points)
 199 {
 200         int     i, j, m, n;
 201         vec1    r1, theta, thinc;
 202         vec2    inner[100];
 203         vec2    outer[100];
 204         vec2    triangle[3];
 205
 206         m = levels;
 207         n = points;
 208         thinc = 2*pi/n;
 209
 210         /* initial inner circle*/
 211         for(i=0; i<n; ++i) {
 212                 inner[i].x = 0.0;
 213                 inner[i].y = 0.0;
 214         }
 215
 216         /* loop thru m levels*/
 217         for(j=1; j<=m; ++j) {
 218                 theta = -j*pi/n;
 219                 r1 = r * (vec1)j/m;
 220
 221                 /* calc n points on outer circle*/
 222                 for(i=0; i<n; ++i) {
 223                         theta += thinc;
 224                         outer[i].x = r1*cos(theta);
 225                         outer[i].y = r1*sin(theta);
 226                 }
 227
 228                 /* construct/draw triangles with vertices on
 229                  * inner and outer circles
 230                  */
 231                 for(i=0; i<n; ++i) {
 232                         triangle[0] = outer[i];
 233                         triangle[1] = outer[(i+1) % n];
 234                         triangle[2] = inner[i];
 235
 236                         /* fill triangle in red*/
 237                         setcolor3(RED);
 238                         polyfill(3, triangle);
 239
 240 #if 1
 241                         /* outline triangle in white*/
 242                         setcolor3(WHITE);
 243                         moveto3(triangle[0]);
 244                         lineto3(triangle[1]);
 245                         lineto3(triangle[2]);
 246                         lineto3(triangle[0]);
 247 #endif
 248                 }
 249
 250                 /* copy points on outer circle to inner arrays*/
 251                 for(i=0; i<n; ++i)
 252                         inner[i] = outer[i];
 253         }
 254 }
 255
 256 /* draw a triangle with cordners v0, v1, v2*/
 257 void
 258 triangle(vec2 v0, vec2 v1, vec2 v2)
 259 {
 260         vec2    poly[3];
 261
 262         poly[0] = v0;
 263         poly[1] = v1;
 264         poly[2] = v2;
 265
 266         setcolor3(GREEN);
 267         polyfill(3, poly);
 268         setcolor3(BLACK);
 269         moveto3(poly[2]);
 270         lineto3(poly[0]);
 271         lineto3(poly[1]);
 272         lineto3(poly[2]);
 273 }
 274
 275 /* draw a quadrilateral with corners v0, v1, v2, v3*/
 276 void
 277 quadrilateral(vec2 v0, vec2 v1, vec2 v2, vec2 v3)
 278 {
 279         vec2    poly[4];
 280
 281         poly[0] = v0;
 282         poly[1] = v1;
 283         poly[2] = v2;
 284         poly[3] = v3;
 285         setcolor3(GREEN);
 286         polyfill(4, poly);
 287         setcolor3(BLACK);
 288         moveto3(poly[3]);
 289         lineto3(poly[0]);
 290         lineto3(poly[1]);
 291         lineto3(poly[2]);
 292         lineto3(poly[3]);
 293 }
 294
 295 /* find intersection of lines v0 to v1 and v2 to v3*/
 296 static int
 297 patch(vec2 v0, vec2 v1, vec2 v2, vec2 v3)
 298 {
 299         vec1    denom;
 300         vec1    mu;
 301         vec2    v4;
 302
 303         denom = (v1.x-v0.x)*(v3.y-v2.y) - (v1.y-v0.y)*(v3.x-v2.x);
 304         if(fabs(denom) > epsilon) {
 305                 mu = ((v2.x-v0.x)*(v3.y-v2.y) - (v2.y-v0.y)*(v3.x-v2.x))/denom;
 306
 307                 /* if intersection between lines v0 to v1 and v2 to v3,
 308                  * call it v4 and form triangles v0,v2,v4 and v1,v3,v4
 309                  */
 310                 if(mu >= 0 && mu <= 1) {
 311                         v4.x = (1-mu)*v0.x + mu*v1.x;
 312                         v4.y = (1-mu)*v0.y + mu*v1.y;
 313                         triangle(v0, v2, v4);
 314                         triangle(v1, v3, v4);
 315                         return 0;
 316                 }
 317         }
 318
 319         /* else find intersection of lines v0 to v2 and v1 to v3*/
 320         denom = (v2.x-v0.x)*(v3.y-v1.y) - (v2.y-v0.y)*(v3.x-v1.x);
 321         if(fabs(denom) > epsilon) {
 322                 mu = ((v1.x-v0.x)*(v3.y-v1.y) - (v1.y-v0.y)*(v3.x-v1.x))/denom;
 323
 324                 /* if intersection between v0 and v1, call it v4
 325                  * and form triangles v0,v1,v4 and v2,v3,v4
 326                  */
 327                 if(mu >= 0 && mu <= 1) {
 328                         v4.x = (1-mu)*v0.x + mu*v2.x;
 329                         v4.y = (1-mu)*v0.y + mu*v2.y;
 330                         triangle(v0, v1, v4);
 331                         triangle(v2, v3, v4);
 332                         return 0;
 333                 }
 334         }
 335
 336         /* there are no proper intersections so form quadrilateral v0,v1,v3,v2*/
 337         quadrilateral(v0, v1, v3, v2);
 338         return 1;
 339 }
 340
 341 /* plotted function*/
 342 static vec1
 343 plotfn(vec1 x, vec1 z)
 344 {
 345         vec1    t;
 346
 347         /* y = 4sin(sqrt(x*x+z*z))/sqrt(x*x+z*z) */
 348         t = sqrt(x*x + z*z);
 349         if(fabs(t) < epsilon)
 350                 return 4.0;
 351         return 4.0 * sin(t) / t;
 352 }
 353
 354 /* draw mathematical function plotfn*/
 355 void
 356 drawgrid(vec1 xmin, vec1 xmax, int nx, vec1 zmin, vec1 zmax, int nz)
 357 {
 358         int     i, j;
 359         vec1    xi, xstep, yij;
 360         vec1    zj, zstep;
 361         vec2    v[2][100];
 362         double  S[5][5];
 363
 364         /* scale it down*/
 365         scale3(1.0/(xmax-xmin)*2, 1.0/(xmax-xmin)*2, 1.0/(zmax-zmin), S);
 366         mult3(Q, S, Q);
 367
 368         /* grid from xmin to xmax in nx steps and zmin to xmax in nz steps*/
 369         xstep = (xmax-xmin)/nx;
 370         zstep = (zmax-zmin)/nz;
 371         xi = xmin;
 372         zj = zmin;
 373
 374         /* calc grid points on first fixed-z line, fine the y-height
 375          * and transfrorm the points (xi,yij,zj) into observed
 376          * position.  Observed first set stored in v[0,1..nx]
 377          */
 378         for(i=0; i<=nx; ++i) {
 379                 yij = plotfn(xi, zj);
 380                 v[0][i].x = Q[1][1]*xi + Q[1][2]*yij + Q[1][3]*zj;
 381                 v[0][i].y = Q[2][1]*xi + Q[2][2]*yij + Q[2][3]*zj;
 382                 xi += xstep;
 383         }
 384
 385         /* run thru consecutive fixed-z lines (the second set)*/
 386         for(j=0; j<nz; ++j) {
 387                 xi = xmin;
 388                 zj += zstep;
 389
 390                 /* calc grid points on this second set, find the
 391                  * y-height and transform the points (xi,yij,zj)
 392                  * into observed position.  Observed second set
 393                  * stored in v[1,0..nx]
 394                  */
 395                 for(i=0; i<=nx; ++i) {
 396                         yij = plotfn(xi, zj);
 397                         v[1][i].x = Q[1][1]*xi + Q[1][2]*yij + Q[1][3]*zj;
 398                         v[1][i].y = Q[2][1]*xi + Q[2][2]*yij + Q[2][3]*zj;
 399                         xi += xstep;
 400                 }
 401
 402                 /* run thru the nx patches formed by these two sets*/
 403                 for(i=0; i<nx; ++i)
 404                         patch(v[0][i], v[0][i+1], v[1][i], v[1][i+1]);
 405
 406                 /* copy second set into first set*/
 407                 for(i=0; i<=nx; ++i)
 408                         v[0][i] = v[1][i];
 409         }
 410 }
 411
 412 /* returns the angle whose tangent is y/x.
 413  * all anomalies such as x=0 are also checked
 414  */
 415 vec1
 416 angle(vec1 x, vec1 y)
 417 {
 418         if(fabs(x) < epsilon)
 419                 if(fabs(y) < epsilon)
 420                         return 0.0;
 421                 else
 422                         if(y > 0.0)
 423                                 return pi*0.5;
 424                         else return pi*1.5;
 425         else
 426                 if(x < 0.0)
 427                         return atan(y/x) + pi;
 428                 else return atan(y/x);
 429 }
 430
 431 /* calc 3d scaling matrix A giving scaling vector sx,sy,sz.
 432  * one unit on the x axis becomes sx units, one unit on y, sy,
 433  * and one unit on the z axis becomes sz units
 434  */
 435 void
 436 scale3(vec1 sx, vec1 sy, vec1 sz, double A[][5])
 437 {
 438         int     i, j;
 439
 440         for(i=1; i<5; ++i)
 441                 for(j=1; j<5; ++j)
 442                         A[i][j] = 0.0;
 443         A[1][1] = sx;
 444         A[2][2] = sy;
 445         A[3][3] = sz;
 446         A[4][4] = 1.0;
 447 }
 448
 449 /* calc 3d axes translation matrix A
 450  * origin translated by vectdor tx,ty,tz
 451  */
 452 void
 453 tran3(vec1 tx, vec1 ty, vec1 tz, double A[][5])
 454 {
 455         int     i, j;
 456
 457         for(i=1; i<5; ++i) {
 458                 for(j=1; j<5; ++j)
 459                         A[i][j] = 0.0;
 460                 A[i][i] = 1.0;
 461         }
 462         A[1][4] = -tx;
 463         A[2][4] = -ty;
 464         A[3][4] = -tz;
 465 }
 466
 467 /* calc 3d axes rotation matrix A.  The axes are
 468  * rotated anti-clockwise through an angle theta radians
 469  * about an axis specified by m: m=1 means x, m=2 y, m=3 z axis
 470  */
 471 void
 472 rot3(int m, vec1 theta, double A[][5])
 473 {
 474         int     i, j, m1, m2;
 475         vec1    c, s;
 476
 477         for(i=1; i<5; ++i)
 478                 for(j=1; j<5; ++j)
 479                         A[i][j] = 0.0;
 480         A[m][m] = 1.0;
 481         A[4][4] = 1.0;
 482         m1 = (m % 3) + 1;
 483         m2 = (m1 % 3) + 1;
 484         c = cos(theta);
 485         s = sin(theta);
 486         A[m1][m1] = c;
 487         A[m2][m2] = c;
 488         A[m1][m2] = s;
 489         A[m2][m1] = s;
 490 }
 491
 492 /* calc the matrix product C of two matrices A and B*/
 493 void
 494 mult3(double A[][5], double B[][5], double C[][5])
 495 {
 496         int     i, j, k;
 497         vec1    ab;
 498
 499         for(i=1; i<5; ++i)
 500                 for(j=1; j<5; ++j) {
 501                         ab = 0;
 502                         for(k=1; k<5; ++k)
 503                                 ab += A[i][k] * B[k][j];
 504                         C[i][j] = ab;
 505                 }
 506 }
 507
 508 /* calc observation matrix Q for given observer*/
 509 void
 510 findQ(void)
 511 {
 512         vec1    alpha, beta, gamma, v, w;
 513         double  E[5][5];
 514         double  F[5][5];
 515         double  G[5][5];
 516         double  H[5][5];
 517         double  U[5][5];
 518
 519         /* calc translation matrix F*/
 520         tran3(eye.x, eye.y, eye.z, F);
 521
 522         /* calc rotation matrix G*/
 523         alpha = angle(-direct.x, -direct.y);
 524         rot3(3, alpha, G);
 525
 526         /* calc rotation matrix H*/
 527         v = sqrt(direct.x*direct.x + direct.y*direct.y);
 528         beta = angle(-direct.z, v);
 529         rot3(2, beta, H);
 530
 531         /* calc rotation matrix U*/
 532         w = sqrt(v*v + direct.z*direct.z);
 533         gamma = angle(-direct.x*w, direct.y*direct.z);
 534         rot3(3, -gamma, U);
 535
 536         /* combine the transformations to find Q*/
 537         mult3(G, F, Q);
 538         mult3(H, Q, E);
 539         mult3(U, E, Q);
 540 }