src/common/dircode.c

#ifndef lint
static const char RCSid[] = "$Id: dircode.c,v 2.12 2020/03/06 19:05:19 greg Exp $";
#endif
/*
 * Compute a 4-byte direction code (externals defined in rtmath.h).
 *
 * Mean accuracy is 0.0022 degrees, with a maximum error of 0.0058 degrees.
 */

#include "rtmath.h"

#define DCSCALE         11584.5         /* ((1<<13)-.5)*sqrt(2) */
#define FXNEG           01
#define FYNEG           02
#define FZNEG           04
#define F1X             010
#define F2Z             020
#define F1SFT           5
#define F2SFT           18
#define FMASK           0x1fff

int32
encodedir(FVECT dv)             /* encode a normalized direction vector */
{
        int32   dc = 0;
        int     cd[3], cm;
        int     i;

        for (i = 0; i < 3; i++)
                if (dv[i] < 0.) {
                        cd[i] = (int)(dv[i] * -DCSCALE + .5);
                        dc |= FXNEG<<i;
                } else
                        cd[i] = (int)(dv[i] * DCSCALE + .5);
        if (!(cd[0] | cd[1] | cd[2]))
                return(0);              /* zero normal */
        if (cd[0] <= cd[1]) {
                dc |= F1X | cd[0] << F1SFT;
                cm = cd[1];
        } else {
                dc |= cd[1] << F1SFT;
                cm = cd[0];
        }
        if (cd[2] <= cm)
                dc |= F2Z | cd[2] << F2SFT;
        else
                dc |= cm << F2SFT;
        if (!dc)        /* don't generate 0 code normally */
                dc = F1X;
        return(dc);
}

#if 0           /* original version for reference */

void
decodedir(FVECT dv, int32 dc)   /* decode a normalized direction vector */
{
        double  d1, d2, der;

        if (!dc) {              /* special code for zero normal */
                dv[0] = dv[1] = dv[2] = 0.;
                return;
        }
        d1 = (dc>>F1SFT & FMASK)*(1./DCSCALE);
        d2 = (dc>>F2SFT & FMASK)*(1./DCSCALE);
        der = sqrt(1. - d1*d1 - d2*d2);
        if (dc & F1X) {
                dv[0] = d1;
                if (dc & F2Z) { dv[1] = der; dv[2] = d2; }
                else { dv[1] = d2; dv[2] = der; }
        } else {
                dv[1] = d1;
                if (dc & F2Z) { dv[0] = der; dv[2] = d2; }
                else { dv[0] = d2; dv[2] = der; }
        }
        if (dc & FXNEG) dv[0] = -dv[0];
        if (dc & FYNEG) dv[1] = -dv[1];
        if (dc & FZNEG) dv[2] = -dv[2];
}

#else   

void
decodedir(FVECT dv, int32 dc)   /* decode a normalized direction vector */
{
        static const short      itab[4][3] = {
                                        {1,0,2},{0,1,2},{1,2,0},{0,2,1}
                                };
        static const RREAL      neg[2] = {1., -1.};
        const int               ndx = ((dc & F2Z) != 0)<<1 | ((dc & F1X) != 0);
        double                  d1, d2, der;

        if (!dc) {              /* special code for zero normal */
                dv[0] = dv[1] = dv[2] = 0.;
                return;
        }
        d1 = (dc>>F1SFT & FMASK)*(1./DCSCALE);
        d2 = (dc>>F2SFT & FMASK)*(1./DCSCALE);
        der = sqrt(1. - d1*d1 - d2*d2);
        dv[itab[ndx][0]] = d1;
        dv[itab[ndx][1]] = d2;
        dv[itab[ndx][2]] = der;
        dv[0] *= neg[(dc&FXNEG)!=0];
        dv[1] *= neg[(dc&FYNEG)!=0];
        dv[2] *= neg[(dc&FZNEG)!=0];
}

#endif

double
dir2diff(int32 dc1, int32 dc2)  /* approx. radians^2 between directions */
{
        FVECT   v1, v2;

        if (dc1 == dc2)
                return 0.;

        decodedir(v1, dc1);
        decodedir(v2, dc2);

        return(2. - 2.*DOT(v1,v2));
}

double
fdir2diff(int32 dc1, FVECT v2)  /* approx. radians^2 between directions */
{
        FVECT   v1;

        decodedir(v1, dc1);

        return(2. - 2.*DOT(v1,v2));
}
Revision:	2.13
Committed:	Tue Apr 7 18:20:52 2020 UTC (5 years, 6 months ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad5R4, rad5R3, rad6R0, HEAD
Changes since 2.12:	+2 -2 lines
Log Message:	Corrected constant that was causing bad encodings for diagonals
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.13	static const char RCSid[] = "$Id: dircode.c,v 2.12 2020/03/06 19:05:19 greg Exp $";
3	greg	2.1	#endif
4			/*
5	greg	2.7	* Compute a 4-byte direction code (externals defined in rtmath.h).
6	greg	2.3	*
7			* Mean accuracy is 0.0022 degrees, with a maximum error of 0.0058 degrees.
8	greg	2.1	*/
9
10	greg	2.7	#include "rtmath.h"
11	greg	2.1
12	greg	2.13	#define DCSCALE 11584.5 /* ((1<<13)-.5)sqrt(2) /
13	greg	2.1	#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	greg	2.5	int32
23	greg	2.8	encodedir(FVECT dv) /* encode a normalized direction vector */
24	greg	2.1	{
25	greg	2.8	int32 dc = 0;
26	greg	2.1	int cd[3], cm;
27	greg	2.8	int i;
28	greg	2.1
29			for (i = 0; i < 3; i++)
30			if (dv[i] < 0.) {
31	greg	2.11	cd[i] = (int)(dv[i] * -DCSCALE + .5);
32	greg	2.1	dc \|= FXNEG<<i;
33			} else
34	greg	2.11	cd[i] = (int)(dv[i] * DCSCALE + .5);
35	greg	2.6	if (!(cd[0] \| cd[1] \| cd[2]))
36			return(0); /* zero normal */
37	greg	2.1	if (cd[0] <= cd[1]) {
38			dc \|= F1X \| cd[0] << F1SFT;
39			cm = cd[1];
40			} else {
41			dc \|= cd[1] << F1SFT;
42			cm = cd[0];
43			}
44			if (cd[2] <= cm)
45			dc \|= F2Z \| cd[2] << F2SFT;
46			else
47			dc \|= cm << F2SFT;
48	greg	2.6	if (!dc) /* don't generate 0 code normally */
49	greg	2.2	dc = F1X;
50	greg	2.1	return(dc);
51			}
52
53	greg	2.9	#if 0 /* original version for reference */
54	greg	2.1
55			void
56	greg	2.8	decodedir(FVECT dv, int32 dc) /* decode a normalized direction vector */
57	greg	2.1	{
58			double d1, d2, der;
59
60	greg	2.6	if (!dc) { /* special code for zero normal */
61			dv[0] = dv[1] = dv[2] = 0.;
62			return;
63			}
64	greg	2.12	d1 = (dc>>F1SFT & FMASK)*(1./DCSCALE);
65			d2 = (dc>>F2SFT & FMASK)*(1./DCSCALE);
66	greg	2.1	der = sqrt(1. - d1d1 - d2d2);
67			if (dc & F1X) {
68			dv[0] = d1;
69			if (dc & F2Z) { dv[1] = der; dv[2] = d2; }
70			else { dv[1] = d2; dv[2] = der; }
71			} else {
72			dv[1] = d1;
73			if (dc & F2Z) { dv[0] = der; dv[2] = d2; }
74			else { dv[0] = d2; dv[2] = der; }
75			}
76			if (dc & FXNEG) dv[0] = -dv[0];
77			if (dc & FYNEG) dv[1] = -dv[1];
78			if (dc & FZNEG) dv[2] = -dv[2];
79			}
80
81	greg	2.9	#else
82
83			void
84			decodedir(FVECT dv, int32 dc) /* decode a normalized direction vector */
85			{
86			static const short itab[4][3] = {
87			{1,0,2},{0,1,2},{1,2,0},{0,2,1}
88			};
89			static const RREAL neg[2] = {1., -1.};
90			const int ndx = ((dc & F2Z) != 0)<<1 \| ((dc & F1X) != 0);
91			double d1, d2, der;
92
93			if (!dc) { /* special code for zero normal */
94			dv[0] = dv[1] = dv[2] = 0.;
95			return;
96			}
97	greg	2.11	d1 = (dc>>F1SFT & FMASK)*(1./DCSCALE);
98			d2 = (dc>>F2SFT & FMASK)*(1./DCSCALE);
99	greg	2.9	der = sqrt(1. - d1d1 - d2d2);
100			dv[itab[ndx][0]] = d1;
101			dv[itab[ndx][1]] = d2;
102			dv[itab[ndx][2]] = der;
103			dv[0] *= neg[(dc&FXNEG)!=0];
104			dv[1] *= neg[(dc&FYNEG)!=0];
105			dv[2] *= neg[(dc&FZNEG)!=0];
106			}
107
108			#endif
109	greg	2.1
110			double
111	greg	2.8	dir2diff(int32 dc1, int32 dc2) /* approx. radians^2 between directions */
112	greg	2.1	{
113			FVECT v1, v2;
114
115	greg	2.10	if (dc1 == dc2)
116			return 0.;
117
118	greg	2.1	decodedir(v1, dc1);
119			decodedir(v2, dc2);
120
121			return(2. - 2.*DOT(v1,v2));
122			}
123
124			double
125	greg	2.8	fdir2diff(int32 dc1, FVECT v2) /* approx. radians^2 between directions */
126	greg	2.1	{
127			FVECT v1;
128
129			decodedir(v1, dc1);
130
131			return(2. - 2.*DOT(v1,v2));
132			}