Return to 3d.c CVS log

Up to [Research Unix] / researchv9 / jerq / src / lib / jj

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

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