home *** CD-ROM | disk | FTP | other *** search
/ Club Amiga de Montreal - CAM / CAM_CD_1.iso / files / 038.lha / BCubic.c < prev    next >
Encoding:
C/C++ Source or Header  |  1987-05-16  |  13.8 KB  |  613 lines
  1. #! /bin/sh
  2. # This is a shell archive, meaning:
  3. # 1. Remove everything above the #! /bin/sh line.
  4. # 2. Save the resulting text in a file.
  5. # 3. Execute the file with /bin/sh (not csh) to create:
  6. #    bcubic.c
  7. #    bcubic.uue
  8. # This archive created: Wed May 13 10:11:10 1987
  9. export PATH; PATH=/bin:/usr/bin:$PATH
  10. echo shar: "extracting 'bcubic.c'" '(13600 characters)'
  11. if test -f 'bcubic.c'
  12. then
  13.     echo shar: "will not over-write existing file 'bcubic.c'"
  14. else
  15. cat << \!Funky!Stuff! > 'bcubic.c'
  16.  
  17. /*
  18.  *  BCUBIC.C        BEZIER CUBIC (3D) CURVED SURFACES
  19.  *
  20.  *  (C)Copyright 1987 by Matthew Dillon, All Rights Reserved.
  21.  *  Permission to redistribute for non-profit only.
  22.  *  Permission to embellish the source however you wish and redistribute.
  23.  *
  24.  *  Compiled with Aztec.  You MUST use the +L option (32 bit ints) or you
  25.  *  are sunk.  You might have to specify more expression space for cc:
  26.  *  (-E1000 should do it).   Addtionaly, the source expects all Amiga
  27.  *  symbols to be precompiled, and MY.LIB to exist.  A good programmer
  28.  *  would not have a problem removing the MY.LIB dependancies.
  29.  *
  30.  *  Usage:
  31.  *    Prop. Gadget changes granularity.  Use the SELECT button to move
  32.  *    control points on the X and Y axis.  Use the MENU button to move
  33.  *    control points on the Z axis.
  34.  *
  35.  *    SELECT: right-left is X axis movement
  36.  *        up-down    is Y axis movement
  37.  *    MENU:    up-down    is Z axis movement
  38.  *
  39.  *    The 3D plane is a linear projection with the Z axis straight up,
  40.  *    the X axis going to the right, and the Y axis going diagonal
  41.  *    (lower-left to upper-right).
  42.  *
  43.  *  NOTE:
  44.  *    This program has been incredibly optimized.  Do not attempt to
  45.  *    gleam the matrix theory from the source unless you are familar
  46.  *    with Bezier Cubic equations.
  47.  */
  48.  
  49. #include <typedefs.h>
  50. #include <xmisc.h>
  51.  
  52. #define ONE 512
  53. #define SSF 9
  54. #define SSFD    (SSF*2)
  55. #define MINGRAN 10            /* must be at least 10  */
  56. #define CSIZE    (ONE/MINGRAN+2)
  57.  
  58. typedef unsigned long ulong;
  59. typedef unsigned short uword;
  60. typedef unsigned char  ubyte;
  61.  
  62. typedef struct PropInfo XPI;
  63. typedef struct Image    IM;
  64.  
  65. extern IMESS *GetMsg();
  66. extern char *malloc();
  67. extern GADGET Gadgets[];
  68.  
  69. #define MYGADGETS   (WINDOWSIZING|WINDOWDRAG|WINDOWDEPTH|WINDOWCLOSE)
  70.  
  71. NW Nw = {
  72.     64, 64, 400, 120,
  73.     0, 1,
  74.     NEWSIZE|MOUSEBUTTONS|MOUSEMOVE|CLOSEWINDOW|GADGETDOWN|GADGETUP,
  75.     MYGADGETS|REPORTMOUSE|ACTIVATE|NOCAREREFRESH|RMBTRAP,
  76.     0, 0, (UBYTE *)"Bezier Cubic (3D curved surfaces), By Matthew Dillon", NULL, NULL,
  77.     32, 64, -1, -1, WBENCHSCREEN
  78. };
  79.  
  80. WIN *Win;
  81. RP  *Rp;
  82. short Ux, Uy, Mx, My;
  83. short SStep = 512/16;
  84. short TStep = 512/16;
  85.  
  86. main(ac, av)
  87. char *av[];
  88. {
  89.     register IMESS *mess;
  90.     char notdone = 1;
  91.     short mx, my, basex, basey, mm, mrmb = 0;
  92.     short gg, gy;
  93.     short pt = -1;
  94.     XPI *po;
  95.  
  96.     disablebreak();
  97.     exiterr(!init_cubic(), "");
  98.     init_gadgets(&Nw, &po);
  99.     exiterr(!openlibs(INTUITION_LIB|GRAPHICS_LIB), "unable to open libs");
  100.     exiterr(!(Win = OpenWindow(&Nw)), "unable to open window");
  101.     Rp = Win->RPort;
  102.     SetAPen(Rp, 3);
  103.     SetDrMd(Rp, JAM2);
  104.  
  105.     uparams();
  106.     precalculate();
  107.     plotcurve();
  108.     plotcontrolpts();
  109.     while (notdone) {
  110.     WaitPort(Win->UserPort);
  111.     mm = 0;
  112.     gg = 0;
  113.     while (mess = GetMsg(Win->UserPort)) {
  114.         switch(mess->Class) {
  115.         case CLOSEWINDOW:
  116.         notdone = 0;
  117.         break;
  118.         case NEWSIZE:
  119.         uparams();
  120.         precalculate();
  121.         clearwindow();
  122.         plotcurve();
  123.         plotcontrolpts();
  124.         break;
  125.         case GADGETUP:
  126.         case GADGETDOWN:
  127.         if (mess->IAddress != (APTR)&Gadgets[0]) {
  128.             do_help();
  129.             break;
  130.         }
  131.         gy = po->VertPot / 256;
  132.         gg = mess->Class;
  133.         break;
  134.         case MOUSEBUTTONS:
  135.         switch(mess->Code) {
  136.         case SELECTDOWN:
  137.             pt = findpoint(mess->MouseX, mess->MouseY);
  138.             basex = mess->MouseX;
  139.             basey = mess->MouseY;
  140.             break;
  141.         case SELECTUP:
  142.             pt = -1;
  143.             break;
  144.         case MENUDOWN:
  145.             pt = findpoint(mess->MouseX, mess->MouseY);
  146.             basex = mess->MouseX;
  147.             basey = mess->MouseY;
  148.             mrmb = 1;
  149.             break;
  150.         case MENUUP:
  151.             mrmb = 0;
  152.             pt = -1;
  153.             break;
  154.         }
  155.         break;
  156.         case MOUSEMOVE:
  157.         if (gg == GADGETDOWN) {
  158.             gy = po->VertPot / 256;
  159.             break;
  160.         }
  161.         if (pt >= 0) {
  162.             mm = 1;
  163.             mx = mess->MouseX;
  164.             my = mess->MouseY;
  165.         }
  166.         break;
  167.         default:
  168.         break;
  169.         }
  170.         ReplyMsg(mess);
  171.     }
  172.     if (gg) {
  173.         char buf[32];
  174.         mm = 0;
  175.         if (gg == GADGETUP)
  176.         gg = 0;
  177.         if (gy + 10 >= 0 && gy + 10 != SStep) {
  178.         char buf[32];
  179.         SStep = gy + 10;
  180.         TStep = SStep;
  181.         Move(Rp, 32, 16);
  182.         sprintf (buf, "Step: %-3ld", SStep);
  183.         Text(Rp, buf, strlen(buf));
  184.         precalculate();
  185.         clearwindow();
  186.         plotcurve();
  187.         plotcontrolpts();
  188.         }
  189.     }
  190.     if (mm && pt >= 0) {
  191.         if (mrmb) {
  192.         recalculate_one(pt, 0, 0, basey - my);
  193.         basey = my;
  194.         } else {
  195.         recalculate_one(pt, mx - basex, basey - my, 0);
  196.         basex = mx;
  197.         basey = my;
  198.         }
  199.         clearwindow();
  200.         plotcurve();
  201.         plotcontrolpts();
  202.     }
  203.     }
  204.     exiterr(1, NULL);
  205. }
  206.  
  207. exiterr(n, str)
  208. char *str;
  209. {
  210.     if (n) {
  211.     if (str)
  212.         puts(str);
  213.     if (Win)
  214.         CloseWindow(Win);
  215.     closelibs(-1);
  216.     exit(1);
  217.     }
  218. }
  219.  
  220. uparams()
  221. {
  222.     Ux = Win->BorderLeft;
  223.     Uy = Win->BorderTop;
  224.     Mx = Win->Width - Win->BorderRight;
  225.     My = Win->Height- Win->BorderBottom;
  226. }
  227.  
  228. long    Pmatrix[3][4][4] = {     0,50,100,150,          /* X    */
  229.                  0,50,100,150,
  230.                  0,50,100,150,
  231.                  0,50,100,150,
  232.  
  233.                  0, 0, 0, 0,        /* Y    */
  234.                 16,16,16,16,
  235.                 32,32,32,32,
  236.                 48,48,48,48,
  237.  
  238.                  0,16,16, 0,        /* Z    */
  239.                 16,32,32,16,
  240.                 16,32,32,16,
  241.                  0,16,16, 0
  242.                 };
  243.  
  244. long    Mmatrix[4][4] = {  -1,    3, -3,    1,
  245.                 3, -6,  3,    0,
  246.                -3,    3,  0,    0,
  247.                 1,    0,  0,    0
  248.             };
  249.  
  250. long    Mmatrix_trans[4][4];
  251.  
  252. long    Amatrix[CSIZE][4];  /* # entries used depends on step rate  */
  253. long    Bmatrix[CSIZE][4];  /* # entries used depends on step rate  */
  254. long    *CXarray;        /* each contains CSIZE*CSIZE*4 bytes    */
  255. long    *CYarray;
  256. long    *CZarray;
  257. short    Amax, Bmax;
  258.  
  259.  
  260. init_cubic()
  261. {
  262.     register short i, j;
  263.  
  264.     CXarray = (long *)malloc(CSIZE*CSIZE*4);
  265.     CYarray = (long *)malloc(CSIZE*CSIZE*4);
  266.     CZarray = (long *)malloc(CSIZE*CSIZE*4);
  267.     if (!CXarray || !CYarray || !CZarray) {
  268.     printf("unable to allocate %ld bytes\n", CSIZE*CSIZE*4*3);
  269.     return(0);
  270.     }
  271.     /* Transpose Mmatrix    */
  272.     for (i = 0; i < 4; ++i) {
  273.     for (j = 0; j < 4; ++j)
  274.         Mmatrix_trans[i][j] = Mmatrix[j][i];
  275.     }
  276.     return(1);
  277. }
  278.  
  279.  
  280. /*
  281.  *  return index of point closest to the given 2D window pixel coordinate.
  282.  *
  283.  *  To do this, each of the 16 control points must be transformed into their
  284.  *  plot coordinates and the range to the specified window pixel taken.
  285.  *  the point with the smallest range wins.
  286.  *  (NOTE: Since we are using the ranges for comparison only, there is
  287.  *   no need to take the square root)
  288.  */
  289.  
  290. findpoint(wx, wy)
  291. {
  292.     register short i, pt;
  293.     register ulong minrange = -1;
  294.     ulong range;
  295.     short x, y;
  296.  
  297.     pt = 0;
  298.     for (i = 0; i < 16; ++i) {
  299.     translate(Pmatrix[0][0][i], Pmatrix[1][0][i], Pmatrix[2][0][i], &x, &y);
  300.     x -= wx;
  301.     y -= wy;
  302.     range = x*x + y*y;
  303.     if (range < minrange) {
  304.         pt = i;
  305.         minrange = range;
  306.     }
  307.     }
  308.     return(pt);
  309. }
  310.  
  311. /*
  312.  *  This need be called only when the step size changes... it calculates
  313.  *  the A and B matrices from the following equation:
  314.  *
  315.  *    Overall equation:   C = SMPM.T. (Where X. means transpose of X)
  316.  *    S = { s^3 s^2 s^1 1 }    s ranging 0 to 1
  317.  *    T = { t^3 t^2 t^1 1 }    t ranging 0 to 1
  318.  *
  319.  *    Thus, a given step size will require calculation of a certain
  320.  *    number of S and T matrices.  We might as well do the multiplication
  321.  *    with M as well and stick the results in A and B for each step
  322.  *
  323.  *    A = SM
  324.  *    B = M.T.
  325.  *
  326.  *    Note: As far as mmult_l() goes, any matrix with either the number
  327.  *    of rows = 1 or number of columns = 1 can be transposed
  328.  *    inherently (that is, it takes no work to transpose it).
  329.  */
  330.  
  331. precalculate()
  332. {
  333.     long Smatrix[4], Tmatrix[4];
  334.  
  335.     register short s, t, i;
  336.  
  337.     Smatrix[3] = Tmatrix[3] = ONE;
  338.     for (s = i = 0; s <= ONE; s += SStep) {
  339.     Smatrix[2] = s;
  340.     Smatrix[1] = (s * s) >> SSF;
  341.     Smatrix[0] = (Smatrix[1] * s) >> SSF;
  342.     mmult_l(Smatrix,Mmatrix,Amatrix[i],1,4,4);
  343.     ++i;
  344.     if (s != ONE && s + SStep > ONE)
  345.         s = ONE - SStep;
  346.     }
  347.     Amax = i;
  348.     for (t = i = 0; t <= ONE; t += TStep) {
  349.     Tmatrix[2] = t;
  350.     Tmatrix[1] = (t * t) >> SSF;
  351.     Tmatrix[0] = (Tmatrix[1] * t) >> SSF;
  352.     mmult_l(Mmatrix_trans,Tmatrix,Bmatrix[i],4,4,1);
  353.     ++i;
  354.     if (t != ONE && t + TStep > ONE)
  355.         t = ONE - TStep;
  356.     }
  357.     Bmax = i;
  358.     recalculate_all();
  359. }
  360.  
  361. /*
  362.  *  Recalculate all global points from the semi-compiled matrices
  363.  *  A and B, and the P matrix (the 16 control points).
  364.  *
  365.  *  C = APB (matrix multiplication), for each dimension.  C is a 1x1
  366.  *  matrix, which means that it is a simple number.
  367.  */
  368.  
  369. recalculate_all()
  370. {
  371.     register short t, i, s, idx;
  372.     register long *A, *B;
  373.     long C[3];        /* 3D result Coordinates    */
  374.  
  375.     for (s = 0; s < Amax; ++s) {
  376.     A = Amatrix[s];
  377.     idx = s * Bmax;
  378.     for (t = 0; t < Bmax; ++t, ++idx) {
  379.         B = Bmatrix[t];
  380.         for (i = 0; i < 3; ++i) {
  381.         long T[4];
  382.         mmult_l(A,Pmatrix[i],T,1,4,4);
  383.         C[i] = T[0] * B[0] +
  384.                T[1] * B[1] +
  385.                T[2] * B[2] +
  386.                T[3] * B[3];
  387.         }
  388.         CXarray[idx] = C[0];
  389.         CYarray[idx] = C[1];
  390.         CZarray[idx] = C[2];
  391.     }
  392.     }
  393. }
  394.  
  395. /*
  396.  *  At this point, we have calculated all the global points and stuck
  397.  *  them in a *huge* table.  If we now want to change one of the control
  398.  *  points, we only have to recalculate one of the matrix elements.
  399.  *  (But still have to loop through the entire grid).
  400.  *
  401.  *  Note that there is a possible optimization here, since this scheme
  402.  *  allows us to determine which line segments will change when we do the
  403.  *  actual plotting.
  404.  */
  405.  
  406. recalculate_one(pt, dx, dy, dz)
  407. {
  408.     register short t, i, idx, col;
  409.     register long  temp;
  410.     long     arow;
  411.     short    s, row;
  412.  
  413.     col = pt & 3;
  414.     row = pt >> 2;
  415.  
  416.     Pmatrix[0][0][pt] += dx;    /* not required until a recalc_all()    */
  417.     Pmatrix[1][0][pt] += dy;
  418.     Pmatrix[2][0][pt] += dz;
  419.  
  420.     for (s = 0; s < Amax; ++s) {
  421.     arow = Amatrix[s][row];
  422.     idx = s * Bmax;
  423.     for (t = 0; t < Bmax; ++t, ++idx) {
  424.         temp = Bmatrix[t][col] * arow;
  425.         CXarray[idx] += temp * dx;
  426.         CYarray[idx] += temp * dy;
  427.         CZarray[idx] += temp * dz;
  428.     }
  429.     }
  430. }
  431.  
  432.  
  433. plotcurve()
  434. {
  435.     register short i, s, t, idx;
  436.     uword x, y;
  437.     static short Corr[CSIZE][2];
  438.  
  439.     for (s = 0; s < Amax; ++s) {
  440.     idx = s * Bmax;
  441.     translate(CXarray[idx]>>SSFD, CYarray[idx]>>SSFD, CZarray[idx]>>SSFD, &x, &y);
  442.     Move(Rp, x, y);
  443.     for (i = t = 0; t < Bmax; ++t, ++idx, ++i)
  444.         translate(CXarray[idx]>>SSFD, CYarray[idx]>>SSFD, CZarray[idx]>>SSFD, &Corr[i][0], &Corr[i][1]);
  445.     PolyDraw(Rp, i, Corr);
  446.     }
  447.     for (t = 0; t < Bmax; ++t) {
  448.     idx = t;
  449.     translate(CXarray[idx]>>SSFD, CYarray[idx]>>SSFD, CZarray[idx]>>SSFD, &x, &y);
  450.     Move(Rp, x, y);
  451.     for (s = i = 0; s < Amax; ++s, idx += Bmax, ++i)
  452.         translate(CXarray[idx]>>SSFD, CYarray[idx]>>SSFD, CZarray[idx]>>SSFD, &Corr[i][0], &Corr[i][1]);
  453.     PolyDraw(Rp, i, Corr);
  454.     }
  455. }
  456.  
  457.  
  458. plotcontrolpts()
  459. {
  460.     register short i;
  461.     short x, y;
  462.  
  463.     SetAPen(Rp, 1);
  464.     for (i = 0; i < 16; ++i) {
  465.     translate(Pmatrix[0][0][i], Pmatrix[1][0][i], Pmatrix[2][0][i], &x, &y);
  466.     Move(Rp, x - 2, y + 0);
  467.     Draw(Rp, x + 2, y + 0);
  468.     Move(Rp, x + 0, y - 2);
  469.     Draw(Rp, x + 0, y + 2);
  470.     }
  471.     SetAPen(Rp, 3);
  472. }
  473.  
  474. translate(x, y, z, wx, wy)
  475. uword *wx, *wy;
  476. {
  477.     *wx = Ux + 50 + x + y/2;
  478.     *wy = My - 50 - (z + y/2);
  479. }
  480.  
  481. /*
  482.  *  MATRIX ROUTINES
  483.  *
  484.  *  mmult_l(a,b,d,n1,n2,n3)
  485.  *
  486.  *  a    n1 x n2
  487.  *  b    n2 x n3
  488.  *  d    n1 x n3
  489.  *
  490.  *  NOTE: This isn't the most efficient way to handle matrix multiplication.
  491.  */
  492.  
  493. mmult_l(a,b,d,n1,n2,n3)
  494. long *a, *b, *d;
  495. {
  496.     register short i, j, k;
  497.     register long *an2 = a;
  498.     register long *bn3;
  499.     register long *dn3 = d;
  500.     register long sum;
  501.  
  502.     for (i = 0; i < n1; ++i, dn3 += n3, an2 += n2) {
  503.     for (j = 0; j < n3; ++j) {
  504.         sum = 0;
  505.         for (k = 0, bn3 = b; k < n2; ++k, bn3 += n3)
  506.         sum += *(an2+k) * *(bn3+j);
  507.         *(dn3+j) = sum;
  508.     }
  509.     }
  510. }
  511.  
  512.  
  513. /*
  514.  *  MISC
  515.  */
  516.  
  517. clearwindow()
  518. {
  519.     SetAPen(Rp, 0);
  520.     RectFill(Rp, Ux, Uy, Mx - 1, My - 1);
  521.     SetAPen(Rp, 3);
  522. }
  523.  
  524. /*
  525.  *  GADGET ROUTINES!    ------------------------------------------------
  526.  */
  527.  
  528. #define NG(nn)    &Gadgets[nn+1]
  529. #define G_YGLOB 1
  530. #define G_XGLOB 2
  531.  
  532. XPI Props[] = {
  533.     { AUTOKNOB|FREEVERT , 0, 0, 0x1FFF, 0x1FFF }
  534. };
  535.  
  536. short pairs[] = { 0,0,10,0,10,10,0,10 };
  537.  
  538. struct Border Border[] = {
  539.     { 0,0,1,1,JAM2,4,pairs,NULL }
  540. };
  541.  
  542. ITEXT Itext[] = {
  543.     { 0,1,JAM2,1,2,NULL,(ubyte *)"?",NULL }
  544. };
  545.  
  546. IM Images[] = {
  547.     { 0,0,2,1,1, NULL, 1, 0, NULL },
  548. };
  549.  
  550. GADGET Gadgets[] = {
  551.     {
  552.     &Gadgets[1], -15, 22, 15, -19, GADGIMAGE|GADGHCOMP|GRELRIGHT|GRELHEIGHT,
  553.     GADGIMMEDIATE|RIGHTBORDER|RELVERIFY,PROPGADGET,
  554.     (APTR)&Images[0],NULL,NULL,0,(APTR)&Props[0], G_YGLOB, 0
  555.     },
  556.     {
  557.     NULL,         -15, 11, 15, 11,  GRELRIGHT,
  558.     RELVERIFY|RIGHTBORDER, BOOLGADGET,
  559.     (APTR)&Border[0],NULL,&Itext[0],0,0,0,0
  560.     }
  561. };
  562.  
  563. GADGET *Gc;
  564. long GUx, GUy;
  565.  
  566. init_gadgets(nw, ppo)
  567. NW *nw;
  568. XPI **ppo;
  569. {
  570.     nw->FirstGadget = &Gadgets[0];
  571.     *ppo = &Props[0];
  572. }
  573.  
  574.  
  575. do_help()
  576. {
  577.     long fh;
  578.     static char *help[] = {
  579. "",
  580. "(C)Copyright 1987 by Matthew Dillon, All Rights Reserved",
  581. "Freely distributable for non-profit only",
  582. "",
  583. "Bezier Cubic Surfaces.  Allows you to fool around with the Bezier Cubic",
  584. "equation for 3D curved surface generation.  The Bezier form uses 16",
  585. "control points to define the surface.",
  586. "",
  587. "To move a control point along the Z axis, position the mouse over the",
  588. "control point in question and move it UP or DOWN with the MENU BUTTON",
  589. "held down.  To move a control point along the X or Y axis, position the",
  590. "mouse over the control point in question and move it LEFT/RIGHT or",
  591. "UP/DOWN with the SELECT BUTTON held down.  The Y axis is on the diagonal,",
  592. "in a lower-left to upper-right going direction.",
  593. "",
  594. "The prop. Gadget is used to define the granualarity of the surface.",
  595. "",
  596. "(Return to continue)",
  597.     NULL
  598.     };
  599.  
  600.     fh = Open("con:0/0/640/200/HELP", 1006);
  601.     if (fh) {
  602.     register short i;
  603.     char c;
  604.     for (i = 0; help[i]; ++i) {
  605.         Write(fh, help[i], strlen(help[i]));
  606.         Write(fh, "\n", 1);
  607.     }
  608.     Read(fh, &c, 1);
  609.     Close(fh);
  610.     }
  611. }
  612.  
  613.