| 1 | #ifndef lint | 
| 2 | static const char RCSid[] = "$Id: dircode.c,v 2.6 2003/06/21 14:48:53 greg Exp $"; | 
| 3 | #endif | 
| 4 | /* | 
| 5 | * Compute a 4-byte direction code (externals defined in rtmath.h). | 
| 6 | * | 
| 7 | * Mean accuracy is 0.0022 degrees, with a maximum error of 0.0058 degrees. | 
| 8 | */ | 
| 9 |  | 
| 10 | #include "rtmath.h" | 
| 11 |  | 
| 12 | #define DCSCALE         11585.2         /* (1<<13)*sqrt(2) */ | 
| 13 | #define FXNEG           01 | 
| 14 | #define FYNEG           02 | 
| 15 | #define FZNEG           04 | 
| 16 | #define F1X             010 | 
| 17 | #define F2Z             020 | 
| 18 | #define F1SFT           5 | 
| 19 | #define F2SFT           18 | 
| 20 | #define FMASK           0x1fff | 
| 21 |  | 
| 22 | int32 | 
| 23 | encodedir(dv)           /* encode a normalized direction vector */ | 
| 24 | FVECT   dv; | 
| 25 | { | 
| 26 | register int32  dc = 0; | 
| 27 | int     cd[3], cm; | 
| 28 | register int    i; | 
| 29 |  | 
| 30 | for (i = 0; i < 3; i++) | 
| 31 | if (dv[i] < 0.) { | 
| 32 | cd[i] = (int)(dv[i] * -DCSCALE); | 
| 33 | dc |= FXNEG<<i; | 
| 34 | } else | 
| 35 | cd[i] = (int)(dv[i] * DCSCALE); | 
| 36 | if (!(cd[0] | cd[1] | cd[2])) | 
| 37 | return(0);              /* zero normal */ | 
| 38 | if (cd[0] <= cd[1]) { | 
| 39 | dc |= F1X | cd[0] << F1SFT; | 
| 40 | cm = cd[1]; | 
| 41 | } else { | 
| 42 | dc |= cd[1] << F1SFT; | 
| 43 | cm = cd[0]; | 
| 44 | } | 
| 45 | if (cd[2] <= cm) | 
| 46 | dc |= F2Z | cd[2] << F2SFT; | 
| 47 | else | 
| 48 | dc |= cm << F2SFT; | 
| 49 | if (!dc)        /* don't generate 0 code normally */ | 
| 50 | dc = F1X; | 
| 51 | return(dc); | 
| 52 | } | 
| 53 |  | 
| 54 |  | 
| 55 | void | 
| 56 | decodedir(dv, dc)       /* decode a normalized direction vector */ | 
| 57 | register FVECT  dv;     /* returned */ | 
| 58 | register int32  dc; | 
| 59 | { | 
| 60 | double  d1, d2, der; | 
| 61 |  | 
| 62 | if (!dc) {              /* special code for zero normal */ | 
| 63 | dv[0] = dv[1] = dv[2] = 0.; | 
| 64 | return; | 
| 65 | } | 
| 66 | d1 = ((dc>>F1SFT & FMASK)+.5)*(1./DCSCALE); | 
| 67 | d2 = ((dc>>F2SFT & FMASK)+.5)*(1./DCSCALE); | 
| 68 | der = sqrt(1. - d1*d1 - d2*d2); | 
| 69 | if (dc & F1X) { | 
| 70 | dv[0] = d1; | 
| 71 | if (dc & F2Z) { dv[1] = der; dv[2] = d2; } | 
| 72 | else { dv[1] = d2; dv[2] = der; } | 
| 73 | } else { | 
| 74 | dv[1] = d1; | 
| 75 | if (dc & F2Z) { dv[0] = der; dv[2] = d2; } | 
| 76 | else { dv[0] = d2; dv[2] = der; } | 
| 77 | } | 
| 78 | if (dc & FXNEG) dv[0] = -dv[0]; | 
| 79 | if (dc & FYNEG) dv[1] = -dv[1]; | 
| 80 | if (dc & FZNEG) dv[2] = -dv[2]; | 
| 81 | } | 
| 82 |  | 
| 83 |  | 
| 84 | double | 
| 85 | dir2diff(dc1, dc2)              /* approx. radians^2 between directions */ | 
| 86 | int32   dc1, dc2; | 
| 87 | { | 
| 88 | FVECT   v1, v2; | 
| 89 |  | 
| 90 | decodedir(v1, dc1); | 
| 91 | decodedir(v2, dc2); | 
| 92 |  | 
| 93 | return(2. - 2.*DOT(v1,v2)); | 
| 94 | } | 
| 95 |  | 
| 96 |  | 
| 97 | double | 
| 98 | fdir2diff(dc1, v2)              /* approx. radians^2 between directions */ | 
| 99 | int32   dc1; | 
| 100 | register FVECT  v2; | 
| 101 | { | 
| 102 | FVECT   v1; | 
| 103 |  | 
| 104 | decodedir(v1, dc1); | 
| 105 |  | 
| 106 | return(2. - 2.*DOT(v1,v2)); | 
| 107 | } |