Return to 3d.c CVS log

Up to [Research Unix] / researchv10no / cmd / view2d / term

Annotation of researchv10no/cmd/view2d/term/3d.c, revision 1.1

1.1     ! root        1: /*% 3cc -c %
        !             2:  * Fixed-point 3d coordinate transform routines for machines with
        !             3:  * 16 bit shorts and 32 bit longs.
        !             4:  */
        !             5: #include       <CC/jerq.h>
        !             6: #include       <CC/3d.h>
        !             7: #include       <CC/font.h>
        !             8: static Hcoordl cur;
        !             9: static Point cen;
        !            10: Point Ap, Bp;
        !            11: static int wid;
        !            12: /*
        !            13:  * angles are scaled so that PI=>32768
        !            14:  * Since unsigned shorts are 16 bits, 2*PI overflows and yields
        !            15:  * 0, which seems appropriate.
        !            16:  */
        !            17: angle cvtangle[361]={
        !            18: 0,     182,    364,    546,    728,    910,    1092,   1274,
        !            19: 1456,  1638,   1820,   2002,   2184,   2366,   2548,   2730,
        !            20: 2912,  3094,   3276,   3458,   3640,   3822,   4004,   4187,
        !            21: 4369,  4551,   4733,   4915,   5097,   5279,   5461,   5643,
        !            22: 5825,  6007,   6189,   6371,   6553,   6735,   6917,   7099,
        !            23: 7281,  7463,   7645,   7827,   8009,   8192,   8374,   8556,
        !            24: 8738,  8920,   9102,   9284,   9466,   9648,   9830,   10012,
        !            25: 10194, 10376,  10558,  10740,  10922,  11104,  11286,  11468,
        !            26: 11650, 11832,  12014,  12196,  12379,  12561,  12743,  12925,
        !            27: 13107, 13289,  13471,  13653,  13835,  14017,  14199,  14381,
        !            28: 14563, 14745,  14927,  15109,  15291,  15473,  15655,  15837,
        !            29: 16019, 16201,  16384,  16566,  16748,  16930,  17112,  17294,
        !            30: 17476, 17658,  17840,  18022,  18204,  18386,  18568,  18750,
        !            31: 18932, 19114,  19296,  19478,  19660,  19842,  20024,  20206,
        !            32: 20388, 20571,  20753,  20935,  21117,  21299,  21481,  21663,
        !            33: 21845, 22027,  22209,  22391,  22573,  22755,  22937,  23119,
        !            34: 23301, 23483,  23665,  23847,  24029,  24211,  24393,  24576,
        !            35: 24758, 24940,  25122,  25304,  25486,  25668,  25850,  26032,
        !            36: 26214, 26396,  26578,  26760,  26942,  27124,  27306,  27488,
        !            37: 27670, 27852,  28034,  28216,  28398,  28580,  28763,  28945,
        !            38: 29127, 29309,  29491,  29673,  29855,  30037,  30219,  30401,
        !            39: 30583, 30765,  30947,  31129,  31311,  31493,  31675,  31857,
        !            40: 32039, 32221,  32403,  32585,  32768,  32950,  33132,  33314,
        !            41: 33496, 33678,  33860,  34042,  34224,  34406,  34588,  34770,
        !            42: 34952, 35134,  35316,  35498,  35680,  35862,  36044,  36226,
        !            43: 36408, 36590,  36772,  36955,  37137,  37319,  37501,  37683,
        !            44: 37865, 38047,  38229,  38411,  38593,  38775,  38957,  39139,
        !            45: 39321, 39503,  39685,  39867,  40049,  40231,  40413,  40595,
        !            46: 40777, 40960,  41142,  41324,  41506,  41688,  41870,  42052,
        !            47: 42234, 42416,  42598,  42780,  42962,  43144,  43326,  43508,
        !            48: 43690, 43872,  44054,  44236,  44418,  44600,  44782,  44964,
        !            49: 45147, 45329,  45511,  45693,  45875,  46057,  46239,  46421,
        !            50: 46603, 46785,  46967,  47149,  47331,  47513,  47695,  47877,
        !            51: 48059, 48241,  48423,  48605,  48787,  48969,  49152,  49334,
        !            52: 49516, 49698,  49880,  50062,  50244,  50426,  50608,  50790,
        !            53: 50972, 51154,  51336,  51518,  51700,  51882,  52064,  52246,
        !            54: 52428, 52610,  52792,  52974,  53156,  53339,  53521,  53703,
        !            55: 53885, 54067,  54249,  54431,  54613,  54795,  54977,  55159,
        !            56: 55341, 55523,  55705,  55887,  56069,  56251,  56433,  56615,
        !            57: 56797, 56979,  57161,  57344,  57526,  57708,  57890,  58072,
        !            58: 58254, 58436,  58618,  58800,  58982,  59164,  59346,  59528,
        !            59: 59710, 59892,  60074,  60256,  60438,  60620,  60802,  60984,
        !            60: 61166, 61348,  61531,  61713,  61895,  62077,  62259,  62441,
        !            61: 62623, 62805,  62987,  63169,  63351,  63533,  63715,  63897,
        !            62: 64079, 64261,  64443,  64625,  64807,  64989,  65171,  65353,
        !            63: 65536, 
        !            64: };
        !            65: static fract sintab[128+2]={
        !            66: 0,     201,    402,    603,    803,    1004,   1205,   1405,
        !            67: 1605,  1805,   2005,   2204,   2404,   2602,   2801,   2998,
        !            68: 3196,  3393,   3589,   3785,   3980,   4175,   4369,   4563,
        !            69: 4756,  4948,   5139,   5329,   5519,   5708,   5896,   6083,
        !            70: 6269,  6455,   6639,   6822,   7005,   7186,   7366,   7545,
        !            71: 7723,  7900,   8075,   8249,   8423,   8594,   8765,   8934,
        !            72: 9102,  9268,   9434,   9597,   9759,   9920,   10079,  10237,
        !            73: 10393, 10548,  10701,  10853,  11002,  11150,  11297,  11442,
        !            74: 11585, 11726,  11866,  12003,  12139,  12273,  12406,  12536,
        !            75: 12665, 12791,  12916,  13038,  13159,  13278,  13395,  13510,
        !            76: 13622, 13733,  13842,  13948,  14053,  14155,  14255,  14353,
        !            77: 14449, 14543,  14634,  14723,  14810,  14895,  14978,  15058,
        !            78: 15136, 15212,  15286,  15357,  15426,  15492,  15557,  15618,
        !            79: 15678, 15735,  15790,  15842,  15892,  15940,  15985,  16028,
        !            80: 16069, 16107,  16142,  16175,  16206,  16234,  16260,  16284,
        !            81: 16305, 16323,  16339,  16353,  16364,  16372,  16379,  16382,
        !            82: 16384, 16382,  
        !            83: };
        !            84: fract isin(x)
        !            85: register angle x;
        !            86: {
        !            87:        register fract a, b;
        !            88:        register i, f=0;
        !            89:        if(x>=PI){
        !            90:                x-=PI;
        !            91:                f=1;
        !            92:        }
        !            93:        if(x>PI/2)
        !            94:                x=PI-x;
        !            95:        i=x>>7;
        !            96:        a=sintab[i];
        !            97:        b=sintab[i+1]-a;
        !            98:        i=x&0177;
        !            99:        a+=(b*i+64)/128;
        !           100:        return(f?-a:a);
        !           101: }
        !           102: fract icos(x)
        !           103: register angle x;
        !           104: {
        !           105:        return(isin(x+PI/2));
        !           106: }
        !           107: long isqrt(a)
        !           108: register long a;
        !           109: {
        !           110:        register long x, y;
        !           111:        if(a<=0)return(0);
        !           112:        for(y=a,x=1;y!=0;y>>=2,x<<=1);
        !           113:        while((y=(a/x+x)>>1)<x)x=y;
        !           114:        return(x);
        !           115: }
        !           116: #define        NSTACK  32
        !           117: Hcoord hcoord(x, y, z, w)
        !           118: fract x, y, z, w;
        !           119: {
        !           120:        Hcoord r;
        !           121:        r.x=x;
        !           122:        r.y=y;
        !           123:        r.z=z;
        !           124:        r.w=w;
        !           125:        return(r);
        !           126: }
        !           127: matrix stack3[NSTACK]={
        !           128:        ONE,0,0,0,
        !           129:        0,ONE,0,0,
        !           130:        0,0,ONE,0,
        !           131:        0,0,0,ONE,
        !           132: };
        !           133: matrix *sp3=stack3;            /* matrix stack pointer */
        !           134: /*
        !           135:  * Initialize the matrix stack.
        !           136:  * v is the viewport into which to map the image.  (d/s) is the distance to the
        !           137:  * screen, whose half-width is 1. (equivalently, d is the distance to a screen of
        !           138:  * half-size s.)
        !           139:  */
        !           140: Bitmap *viewport;
        !           141: init3(v, d, s)
        !           142: Bitmap *v;
        !           143: {
        !           144:        sp3=stack3;
        !           145:        ident(*sp3);
        !           146:        (*sp3)[0][0]=(*sp3)[2][2]=d;
        !           147:        (*sp3)[1][1]=s;
        !           148:        (*sp3)[3][3]=s;
        !           149:        viewport=v;
        !           150:        cen=div(add(v->rect.origin, v->rect.corner), 2);
        !           151:        wid=min(v->rect.corner.x-v->rect.origin.x,
        !           152:                v->rect.corner.y-v->rect.origin.y)/2;
        !           153: }
        !           154: /*
        !           155:  * Push a copy of the top of the matrix stack onto the stack
        !           156:  */
        !           157: push3(){
        !           158:        register i, j;
        !           159:        ;if(sp3==stack3+NSTACK-1)
        !           160:                return(-1);
        !           161:        for(i=0;i!=4;i++) for(j=0;j!=4;j++)
        !           162:                sp3[1][i][j]=sp3[0][i][j];
        !           163:        sp3++;
        !           164:        return(0);
        !           165: }
        !           166: /*
        !           167:  * pop the top matrix off the stack
        !           168:  */
        !           169: pop3(){
        !           170:        if(sp3==stack3)
        !           171:                return(-1);
        !           172:        --sp3;
        !           173:        return(0);
        !           174: }
        !           175: /*
        !           176:  * Multiply the top matrix on the stack by rotation of angle theta about
        !           177:  * the given axis.
        !           178:  */
        !           179: rot3(theta, axis)
        !           180: angle theta;
        !           181: {
        !           182:        rotsc3(isin(theta), icos(theta), axis);
        !           183: }
        !           184: /*
        !           185:  * rotate, given the sin and cosine of the rotation angle.
        !           186:  */
        !           187: rotsc3(s, c, axis)
        !           188: short s, c;
        !           189: {
        !           190:        matrix m;
        !           191:        register i=(axis+1)%3, j=(axis+2)%3;
        !           192:        ident(m);
        !           193:        m[i][i] = c;
        !           194:        m[i][j] = -s;
        !           195:        m[j][i] = s;
        !           196:        m[j][j] = c;
        !           197:        xform3(m);
        !           198: }
        !           199: /*
        !           200:  * Scale by a factor of [xyz]/w about axis [xyz]
        !           201:  */
        !           202: scale3(p)
        !           203: Hcoord p;
        !           204: {
        !           205:        matrix m;
        !           206:        ident(m);
        !           207:        m[0][0]=p.x;
        !           208:        m[1][1]=p.y;
        !           209:        m[2][2]=p.z;
        !           210:        m[3][3]=p.w;
        !           211:        xform3(m);
        !           212: }
        !           213: /*
        !           214:  * translate by (dx/dw, dy/dw, dz/dw)
        !           215:  */
        !           216: tran3(p)
        !           217: Hcoord p;
        !           218: {
        !           219:        matrix m;
        !           220:        ident(m);
        !           221:        m[0][0]=p.w;
        !           222:        m[1][1]=p.w;
        !           223:        m[2][2]=p.w;
        !           224:        m[3][3]=p.w;
        !           225:        m[0][3]=p.x;
        !           226:        m[1][3]=p.y;
        !           227:        m[2][3]=p.z;
        !           228:        xform3(m);
        !           229: }
        !           230: /*
        !           231:  * Store an identity matrix in m
        !           232:  */
        !           233: ident(m)
        !           234: matrix m;
        !           235: {
        !           236:        register fract *s = &m[0][0];
        !           237:        *s++=ONE;*s++=0;*s++=0;*s++=0;
        !           238:        *s++=0;*s++=ONE;*s++=0;*s++=0;
        !           239:        *s++=0;*s++=0;*s++=ONE;*s++=0;
        !           240:        *s++=0;*s++=0;*s++=0;*s++=ONE;
        !           241: }
        !           242: /*
        !           243:  * Multiply the top of the matrix stack by m
        !           244:  */
        !           245: xform3(m)
        !           246: matrix m;
        !           247: {
        !           248:        register i, j, k;
        !           249:        long sum, max=0;
        !           250:        long tmp[4][4];
        !           251:        char buf[128];
        !           252: 
        !           253:        for(i=0;i!=4;i++) for(j=0;j!=4;j++){
        !           254:                sum=0;
        !           255:                for(k=0;k!=4;k++)
        !           256:                        sum+=(long)(*sp3)[i][k]*m[k][j];
        !           257:                tmp[i][j]=sum;
        !           258:                if(sum<0)
        !           259:                        sum = -sum;
        !           260:                if(sum>max)
        !           261:                        max=sum;
        !           262:        }
        !           263:        if((k = max/(ONE-1)) == 0) k = 1;
        !           264:        for(i=0;i!=4;i++) for(j=0;j!=4;j++)
        !           265:                (*sp3)[i][j]=tmp[i][j]/k;
        !           266: }
        !           267: long dot(a, b)
        !           268: Hcoord a, b;
        !           269: {
        !           270:        return((long)a.x*b.x+(long)a.y*b.y+(long)a.z*b.z);
        !           271: }
        !           272: Hcoord unitize(a)
        !           273: Hcoord a;
        !           274: {
        !           275:        short l=isqrt(dot(a, a));
        !           276:        Hcoord v;
        !           277:        if(l==0){
        !           278:                v=a;
        !           279:                v.w=0;
        !           280:        }
        !           281:        else{
        !           282:                if(a.w<0) l = -l;
        !           283:                v.x=muldiv(a.x, ONE, l);
        !           284:                v.y=muldiv(a.y, ONE, l);
        !           285:                v.z=muldiv(a.z, ONE, l);
        !           286:                v.w=ONE;
        !           287:        }
        !           288:        return(v);
        !           289: }
        !           290: Hcoord cross(a, b)
        !           291: Hcoord a, b;
        !           292: {
        !           293:        Hcoord r;
        !           294:        r.x=((long)a.y*b.z-(long)a.z*b.y)/ONE;
        !           295:        r.y=((long)a.z*b.x-(long)a.x*b.z)/ONE;
        !           296:        r.z=((long)a.x*b.y-(long)a.y*b.x)/ONE;
        !           297:        r.w=ONE;
        !           298:        return(r);
        !           299: }
        !           300: /*
        !           301:  * multiply the top of the maxrix stack by a view-pointing transformation
        !           302:  * with the eyepoint at e, looking in direction (lx, ly, lz),
        !           303:  * with up direction (ux, uy, uz)
        !           304:  */
        !           305: look3(e, l, u)
        !           306: Hcoord e, l, u;
        !           307: {
        !           308:        matrix m;
        !           309:        fract len;
        !           310:        Hcoord r;
        !           311: 
        !           312:        /*printf("look3(e=(%d,%d,%d,%d), l=(%d,%d,%d,%d), u=(%d,%d,%d,%d))\n",
        !           313:                e.x,e.y,e.z,e.w,l.x,l.y,l.z,l.w,u.x,u.y,u.z,u.w);*/
        !           314:        l=unitize(l);
        !           315:        if(l.w==0){
        !           316:                l.y=ONE;
        !           317:                l.w=ONE;
        !           318:        }
        !           319:        /* Adjust u to be perpendicular to l */
        !           320:        u=unitize(u);
        !           321:        len=dot(u, l)/ONE;
        !           322:        u.x-=muldiv(len, l.x, ONE);
        !           323:        u.y-=muldiv(len, l.y, ONE);
        !           324:        u.z-=muldiv(len, l.z, ONE);
        !           325:        u=unitize(u);
        !           326:        if(u.w==0){
        !           327:                u.z=ONE;
        !           328:                u.w=ONE;
        !           329:        }
        !           330:        r=cross(l, u);
        !           331:        /*printf("r=(%d,%d,%d,%d), l=(%d,%d,%d,%d), u=(%d,%d,%d,%d))\n",
        !           332:                r.x,r.y,r.z,r.w,l.x,l.y,l.z,l.w,u.x,u.y,u.z,u.w);*/
        !           333:        /* make the matrix to transform from (rlu) space to (xyz) space */
        !           334:        ident(m);
        !           335:        m[0][0]=r.x; m[0][1]=r.y; m[0][2]=r.z;
        !           336:        m[1][0]=l.x; m[1][1]=l.y; m[1][2]=l.z;
        !           337:        m[2][0]=u.x; m[2][1]=u.y; m[2][2]=u.z;
        !           338:        xform3(m);
        !           339:        tran3(hcoord(e.x, e.y, e.z, -e.w));
        !           340: }
        !           341: Hcoordl xformp(p)
        !           342: Hcoord p;
        !           343: {
        !           344:        Hcoordl r;
        !           345:        r.x=(long)(*sp3)[0][0]*p.x+(long)(*sp3)[0][1]*p.y
        !           346:                +(long)(*sp3)[0][2]*p.z+(long)(*sp3)[0][3]*p.w;
        !           347:        r.y=(long)(*sp3)[1][0]*p.x+(long)(*sp3)[1][1]*p.y
        !           348:                +(long)(*sp3)[1][2]*p.z+(long)(*sp3)[1][3]*p.w;
        !           349:        r.z=(long)(*sp3)[2][0]*p.x+(long)(*sp3)[2][1]*p.y
        !           350:                +(long)(*sp3)[2][2]*p.z+(long)(*sp3)[2][3]*p.w;
        !           351:        r.w=(long)(*sp3)[3][0]*p.x+(long)(*sp3)[3][1]*p.y
        !           352:                +(long)(*sp3)[3][2]*p.z+(long)(*sp3)[3][3]*p.w;
        !           353:        return(r);
        !           354: }
        !           355: move3(p)
        !           356: Hcoord p;
        !           357: {
        !           358:        cur=xformp(p);
        !           359:        /*printf("move3(%d,%d,%d,%d)\n", p.x, p.y, p.z, p.w);*/
        !           360: }
        !           361: line3(p)
        !           362: Hcoord p;
        !           363: {
        !           364:        Hcoordl new;
        !           365:        /*printf("line3(%d,%d,%d,%d)\n", p.x, p.y, p.z, p.w);*/
        !           366:        new=xformp(p);
        !           367:        clipanddraw(cur, new);
        !           368:        cur=new;
        !           369: }
        !           370: /*
        !           371:  * clip points to the viewing volume, then do perspective division
        !           372:  * and mapping to screen coordinates and draw lines.  We only clip against the
        !           373:  * eye plane, trusting segment to clip in x and z.  There probably should
        !           374:  * be some overflow precaution taken in the multiplication by n/d.
        !           375:  */
        !           376: int clip3scale = 0;
        !           377: 
        !           378: static clipanddraw(a, b)
        !           379: Hcoordl a, b;
        !           380: {
        !           381:        Point A, B;
        !           382:        Hcoordl t;
        !           383:        int clip=0;
        !           384:        long n, d;
        !           385:        long m;
        !           386:        long as, bs;
        !           387: #define                OV      (32767)         /* 2**(15.5) */
        !           388:        /*
        !           389:         * we clip against the plane p=(0,1,0,-1)  [this corresponds
        !           390:         * to the euclidean plane y=1.  This is picked rather than y=0
        !           391:         * to avoid a divide check in the perspective divide]
        !           392:         * First we check to see if one endpoint is on each side of the plane.
        !           393:         * and if necessary switch them so that a is behind and b is in front.
        !           394:         * then we replace a by a+(b-a)*n/d, where p.(a+(b-a)*n/d)=0
        !           395:         * i.e. a.y+(b.y-a.y)*n/d-a.w-(b.w-a.w)*n/d=0
        !           396:         *      (b.y-a.y-b.w+a.w)*n/d=a.w-a.y
        !           397:         *      n/d=(a.w-a.y)/(b.y-a.y-b.w+a.w)
        !           398:         */
        !           399:        /*printf("clipanddraw(a=(%d,%d,%d,%d), b=(%d,%d,%d,%d)\n",
        !           400:                a.x,a.y,a.z,a.w,b.x,b.y,b.z,b.w);*/
        !           401:        if(a.w>0?a.y<a.w:a.y>a.w)       /* eqv. a.y/a.w<1, but avoids divide */
        !           402:                clip++;
        !           403:        if(b.w>0?b.y<b.w:b.y>b.w){
        !           404:                if(clip)                /* both points behind y=1 plane */
        !           405:                {
        !           406:                        /*printf("clipped out of existence\n");*/
        !           407:                        return;
        !           408:                }
        !           409:                t=a;
        !           410:                a=b;
        !           411:                b=t;
        !           412:                clip++;
        !           413:        }
        !           414:        /* check for overflow */
        !           415:        m = max(max(abs(a.x), abs(a.y)), max(abs(a.z), abs(b.x)));
        !           416:        m = max(m, max(max(abs(b.y), abs(b.z)), max(abs(a.w), abs(b.w))));
        !           417:        if(m > OV)
        !           418:        {
        !           419:                m = (m+OV-1)/OV;
        !           420:                a.x /= m;
        !           421:                a.y /= m;
        !           422:                a.z /= m;
        !           423:                a.w /= m;
        !           424:                b.x /= m;
        !           425:                b.y /= m;
        !           426:                b.z /= m;
        !           427:                b.w /= m;
        !           428:        }
        !           429:        if(clip){
        !           430:                n=a.w-a.y;
        !           431:                d=b.y-a.y-b.w+a.w;
        !           432:                if(d==0)                        /* never happens? */
        !           433:                        return;
        !           434:                a.x+=(b.x-a.x)*n/d;
        !           435:                a.y+=(b.y-a.y)*n/d;
        !           436:                a.z+=(b.z-a.z)*n/d;
        !           437:                a.w+=(b.w-a.w)*n/d;     /* don't use w below */
        !           438:                /*printf("clipped to a=(%d,%d,%d,%d), b=(%d,%d,%d,%d)\n",
        !           439:                        a.x,a.y,a.z,a.w,b.x,b.y,b.z,b.w);*/
        !           440:        }
        !           441:        if(clip3scale)
        !           442:        {
        !           443:                as = a.w*clip3scale;
        !           444:                bs = b.w*clip3scale;
        !           445:        }
        !           446:        else
        !           447:        {
        !           448:                as = a.y;
        !           449:                bs = b.y;
        !           450:        }
        !           451:        A=add(Pt((short)muldiv(a.x, wid, as), (short)muldiv(-a.z, wid, as)), cen);
        !           452:        B=add(Pt((short)muldiv(b.x, wid, bs), (short)muldiv(-b.z, wid, bs)), cen);
        !           453: #define        safe(P) if(P.x<-16000)P.x= -16000;if(P.x>16000)P.x=16000;if(P.y<-16000)P.y= -16000;if(P.y>16000)P.y=16000
        !           454:        safe(A);
        !           455:        safe(B);
        !           456:        Ap = A;
        !           457:        Bp = B;
        !           458:        segment(viewport, A, B, F_OR);
        !           459:        /*printf("---->drew points (%d,%d) - (%d,%d)\n", A.x, A.y, B.x, B.y);*/
        !           460: }

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.