src/common/depthcodec.c

#ifndef lint
static const char RCSid[] = "$Id: depthcodec.c,v 2.12 2022/03/03 16:09:31 greg Exp $";
#endif
/*
 * Routines to encode/decoded 16-bit depths
 */

#include "copyright.h"

#include <stdlib.h>
#include <math.h>
#include <ctype.h>
#include "rtio.h"
#include "fvect.h"
#include "depthcodec.h"


#ifndef depth2code
/* Encode depth as 16-bit signed integer */
int
depth2code(double d, double dref)
{
        if (d <= .0)
                return -32768;

        if (d > dref)
                return (int)(32768.001 - 32768.*dref/d) - 1;

        return (int)(32767.*d/dref - 32768.999);
}
#endif


#ifndef code2depth
/* Decode depth from 16-bit signed integer */
double
code2depth(int c, double dref)
{
        if (c <= -32768)
                return .0;

        if (c >= 32767)
                return FHUGE;

        if (c >= -1)
                return dref*32768./(32766.5 - c);

        return dref*(32768.5 + c)*(1./32767.);
}
#endif


/* Set codec defaults */
void
set_dc_defaults(DEPTHCODEC *dcp)
{
        memset(dcp, 0, sizeof(DEPTHCODEC));
        dcp->finp = stdin;
        dcp->inpname = "<stdin>";
        dcp->format = 'a';
        dcp->refdepth = 1.;
        dcp->depth_unit[0] = '1';
        dcp->vw = stdview;
        dcp->res.rt = PIXSTANDARD;
        if (!progname) progname = "depth_codec";
}


/* process header line */
static int
headline(char *s, void *p)
{
        DEPTHCODEC      *dcp = (DEPTHCODEC *)p;
        int             rv;

        if (formatval(dcp->inpfmt, s))  /* don't pass format */
                return 0;

        if ((rv = isbigendian(s)) >= 0) {
                dcp->swapped = (nativebigendian() != rv);
                return 0;
        }
        if (!strncmp(s, "NCOMP=", 6)) {
                if (atoi(s+6) != 1) {
                        if (dcp->hdrflags & HF_STDERR) {
                                fputs(dcp->inpname, stderr);
                                fputs(": NCOMP must equal 1\n", stderr);
                        }
                        return -1;
                }
                return 0;
        }
        if (!strncmp(s, "NROWS=", 6)) {
                dcp->res.yr = atoi(s+6);
                return 0;
        }
        if (!strncmp(s, "NCOLS=", 6)) {
                dcp->res.xr = atoi(s+6);
                return 0;
        }
                                                /* check for reference depth */
        if (!strncmp(s, DEPTHSTR, LDEPTHSTR)) {
                char    *cp;
                strlcpy(dcp->depth_unit, s+LDEPTHSTR, sizeof(dcp->depth_unit));
                cp = dcp->depth_unit;
                while (*cp) cp++;
                while (cp > dcp->depth_unit && isspace(cp[-1])) cp--;
                *cp = '\0';
                dcp->refdepth = atof(dcp->depth_unit);
                if (dcp->refdepth <= .0) {
                        if (dcp->hdrflags & HF_STDERR) {
                                fputs(dcp->inpname, stderr);
                                fputs(": bad reference depth in input header\n",
                                                stderr);
                        }
                        return -1;
                }
                if (dcp->hdrflags & HF_ENCODE)
                        return 0;       /* will add this later */
        } else if (!strncmp(s, "SAMP360=", 8))
                dcp->gotview--;
        else if (isview(s))             /* get view params */
                dcp->gotview += (sscanview(&dcp->vw, s) > 0);
        if (dcp->hdrflags & HF_HEADOUT)
                fputs(s, stdout);       /* copy to standard output */
        return 1;
}


/* Load/copy header */
int
process_dc_header(DEPTHCODEC *dcp, int ac, char *av[])
{
        if (dcp->hdrflags & HF_HEADIN &&
                        getheader(dcp->finp, headline, dcp) < 0) {
                if (dcp->hdrflags & HF_STDERR) {
                        fputs(dcp->inpname, stderr);
                        fputs(": bad header\n", stderr);
                }
                return 0;
        }
        dcp->gotview *= (dcp->gotview > 0);
                                        /* get resolution string? */
        if (dcp->hdrflags & HF_RESIN &&
                        (dcp->res.xr <= 0) | (dcp->res.yr <= 0) &&
                        !fgetsresolu(&dcp->res, dcp->finp)) {
                if (dcp->hdrflags & HF_STDERR) {
                        fputs(dcp->inpname, stderr);
                        fputs(": bad resolution string\n", stderr);
                }
                return 0;
        }
        if (dcp->hdrflags & HF_HEADOUT) {       /* finish header */
                if (!(dcp->hdrflags & HF_HEADIN))
                        newheader("RADIANCE", stdout);
                if (ac > 0)
                        printargs(ac, av, stdout);
                if (dcp->hdrflags & HF_ENCODE) {
                        fputs(DEPTHSTR, stdout);
                        fputs(dcp->depth_unit, stdout);
                        fputc('\n', stdout);
                        fputformat(DEPTH16FMT, stdout);
                } else {
                        fputs("NCOMP=1\n", stdout);
                        if ((dcp->hdrflags & (HF_RESIN|HF_RESOUT)) == HF_RESIN)
                                printf("NCOLS=%d\nNROWS=%d\n",
                                                scanlen(&dcp->res),
                                                numscans(&dcp->res));
                        switch (dcp->format) {
                        case 'a':
                                fputformat("ascii", stdout);
                                break;
                        case 'f':
                                fputendian(stdout);
                                fputformat("float", stdout);
                                break;
                        case 'd':
                                fputendian(stdout);
                                fputformat("double", stdout);
                                break;
                        }
                }
                fputc('\n', stdout);
        }
        if (dcp->hdrflags & HF_RESOUT)  /* put resolution string? */
                fputsresolu(&dcp->res, stdout);

        dcp->dstart = dcp->curpos = ftell(dcp->finp);
        return 1;
}


/* Check that we have what we need to decode depths */
int
check_decode_depths(DEPTHCODEC *dcp)
{
        if (dcp->hdrflags & HF_ENCODE) {
                if (dcp->hdrflags & HF_STDERR) {
                        fputs(progname, stderr);
                        fputs(": wrong header mode for decode\n", stderr);
                }
                return 0;
        }
        if (dcp->inpfmt[0] && strcmp(dcp->inpfmt, DEPTH16FMT)) {
                if (dcp->hdrflags & HF_STDERR) {
                        fputs(dcp->inpname, stderr);
                        fputs(": unexpected input format: ", stderr);
                        fputs(dcp->inpfmt, stderr);
                        fputc('\n', stderr);
                }
                return 0;
        }
        return 1;
}


/* Check that we have what we need to decode world positions */
int
check_decode_worldpos(DEPTHCODEC *dcp)
{
        char    *err;

        if (!check_decode_depths(dcp))
                return 0;
        if ((dcp->res.xr <= 0) | (dcp->res.yr <= 0)) {
                if (dcp->hdrflags & HF_STDERR) {
                        fputs(progname, stderr);
                        fputs(": missing map resolution\n", stderr);
                }
                return 0;
        }
        if (!dcp->gotview) {
                if (dcp->hdrflags & HF_STDERR) {
                        fputs(dcp->inpname, stderr);
                        fputs(": missing view\n", stderr);
                }
                return 0;
        }
        if ((err = setview(&dcp->vw)) != NULL) {
                if (dcp->hdrflags & HF_STDERR) {
                        fputs(dcp->inpname, stderr);
                        fputs(": input view error: ", stderr);
                        fputs(err, stderr);
                        fputc('\n', stderr);
                }
                return 0;
        }
        return 1;       
}


/* Decode next depth pixel */
double
decode_depth_next(DEPTHCODEC *dcp)
{
        int     c;

        if (dcp->use_last) {
                dcp->use_last = 0;
                return code2depth(dcp->last_dc, dcp->refdepth);
        }
        c = getint(2, dcp->finp);

        if (c == EOF && feof(dcp->finp))
                return -1.;

        dcp->last_dc = c;
        dcp->curpos += 2;

        return code2depth(c, dcp->refdepth);
}


/* Compute world position from depth */
int
compute_worldpos(FVECT wpos, DEPTHCODEC *dcp, int x, int y, double d)
{
        RREAL   loc[2];
        FVECT   rdir;

        if (d >= FHUGE*.99)
                goto badval;

        pix2loc(loc, &dcp->res, x, y);

        if (viewray(wpos, rdir, &dcp->vw, loc[0], loc[1]) >= -FTINY) {
                VSUM(wpos, wpos, rdir, d);
                return 1;
        }
badval:
        VCOPY(wpos, dcp->vw.vp);
        return 0;
}


/* Decode the next world position */
int
decode_worldpos_next(FVECT wpos, DEPTHCODEC *dcp)
{
        const long      n = (dcp->curpos - dcp->dstart)>>1;
        int             x, y;
        double          d;

        d = decode_depth_next(dcp);
        if (d < -FTINY)
                return -1;

        x = scanlen(&dcp->res);
        y = n / x;
        x = n - (long)y*x;

        return compute_worldpos(wpos, dcp, x, y, d);
}


/* Seek to the indicated pixel position */
int
seek_dc_pix(DEPTHCODEC *dcp, int x, int y)
{
        long    seekpos;

        if ((dcp->res.xr <= 0) | (dcp->res.yr <= 0)) {
                if (dcp->hdrflags & HF_STDERR) {
                        fputs(progname, stderr);
                        fputs(": need map resolution to seek\n", stderr);
                }
                return -1;
        }
        if ((x < 0) | (y < 0) ||
                        (x >= scanlen(&dcp->res)) | (y >= numscans(&dcp->res))) {
                if (dcp->hdrflags & HF_STDERR) {
                        fputs(dcp->inpname, stderr);
                        fputs(": warning - pixel index off map\n", stderr);
                }
                return 0;
        }
        seekpos = dcp->dstart + 2*((long)y*scanlen(&dcp->res) + x);

        if (seekpos == dcp->curpos-2) {
                dcp->use_last++;        /* avoids seek/read */
                return 1;
        }
        if (seekpos != dcp->curpos &&
                        fseek(dcp->finp, seekpos, SEEK_SET) == EOF) {
                if (dcp->hdrflags & HF_STDERR) {
                        fputs(dcp->inpname, stderr);
                        fputs(": seek error\n", stderr);
                }
                return -1;
        }
        dcp->curpos = seekpos;
        dcp->use_last = 0;
        return 1;
}


/* Read and decode depth for the given pixel */
double
decode_depth_pix(DEPTHCODEC *dcp, int x, int y)
{
        int     rval = seek_dc_pix(dcp, x, y);

        if (rval < 0)
                return -1.;
        if (!rval)
                return .0;

        return decode_depth_next(dcp);
}


/* Read and decode the world position at the given pixel */
int
get_worldpos_pix(FVECT wpos, DEPTHCODEC *dcp, int x, int y)
{
        double  d = decode_depth_pix(dcp, x, y);

        if (d < -FTINY)
                return -1;
        
        return compute_worldpos(wpos, dcp, x, y, d);
}
Revision:	2.13
Committed:	Thu Sep 22 21:45:28 2022 UTC (20 months, 1 week ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad5R4, HEAD
Changes since 2.12:	+23 -3 lines
Log Message:	feat(rcode_depth,rcode_norm): Allow NROWS,NCOLS in header rather than res string
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.13	static const char RCSid[] = "$Id: depthcodec.c,v 2.12 2022/03/03 16:09:31 greg Exp $";
3	greg	2.1	#endif
4			/*
5			* Routines to encode/decoded 16-bit depths
6			*/
7
8			#include "copyright.h"
9
10			#include <stdlib.h>
11			#include <math.h>
12			#include <ctype.h>
13			#include "rtio.h"
14			#include "fvect.h"
15			#include "depthcodec.h"
16
17
18	greg	2.3	#ifndef depth2code
19	greg	2.1	/* Encode depth as 16-bit signed integer */
20			int
21			depth2code(double d, double dref)
22			{
23			if (d <= .0)
24			return -32768;
25
26			if (d > dref)
27	greg	2.3	return (int)(32768.001 - 32768.*dref/d) - 1;
28	greg	2.1
29	greg	2.8	return (int)(32767.*d/dref - 32768.999);
30	greg	2.1	}
31			#endif
32
33
34	greg	2.3	#ifndef code2depth
35	greg	2.1	/* Decode depth from 16-bit signed integer */
36			double
37			code2depth(int c, double dref)
38			{
39			if (c <= -32768)
40			return .0;
41
42			if (c >= 32767)
43			return FHUGE;
44	greg	2.5
45	greg	2.7	if (c >= -1)
46	greg	2.5	return dref*32768./(32766.5 - c);
47
48	greg	2.7	return dref(32768.5 + c)(1./32767.);
49	greg	2.1	}
50	greg	2.3	#endif
51	greg	2.1
52
53			/* Set codec defaults */
54			void
55			set_dc_defaults(DEPTHCODEC *dcp)
56			{
57			memset(dcp, 0, sizeof(DEPTHCODEC));
58			dcp->finp = stdin;
59			dcp->inpname = "<stdin>";
60			dcp->format = 'a';
61			dcp->refdepth = 1.;
62			dcp->depth_unit[0] = '1';
63			dcp->vw = stdview;
64			dcp->res.rt = PIXSTANDARD;
65			if (!progname) progname = "depth_codec";
66			}
67
68
69			/* process header line */
70			static int
71			headline(char s, void p)
72			{
73			DEPTHCODEC dcp = (DEPTHCODEC )p;
74	greg	2.4	int rv;
75	greg	2.1
76			if (formatval(dcp->inpfmt, s)) /* don't pass format */
77			return 0;
78	greg	2.4
79			if ((rv = isbigendian(s)) >= 0) {
80			dcp->swapped = (nativebigendian() != rv);
81			return 0;
82			}
83	greg	2.13	if (!strncmp(s, "NCOMP=", 6)) {
84			if (atoi(s+6) != 1) {
85			if (dcp->hdrflags & HF_STDERR) {
86			fputs(dcp->inpname, stderr);
87			fputs(": NCOMP must equal 1\n", stderr);
88			}
89			return -1;
90			}
91			return 0;
92			}
93			if (!strncmp(s, "NROWS=", 6)) {
94			dcp->res.yr = atoi(s+6);
95			return 0;
96			}
97			if (!strncmp(s, "NCOLS=", 6)) {
98			dcp->res.xr = atoi(s+6);
99			return 0;
100			}
101			/* check for reference depth */
102	greg	2.1	if (!strncmp(s, DEPTHSTR, LDEPTHSTR)) {
103			char *cp;
104			strlcpy(dcp->depth_unit, s+LDEPTHSTR, sizeof(dcp->depth_unit));
105			cp = dcp->depth_unit;
106			while (*cp) cp++;
107			while (cp > dcp->depth_unit && isspace(cp[-1])) cp--;
108			*cp = '\0';
109			dcp->refdepth = atof(dcp->depth_unit);
110			if (dcp->refdepth <= .0) {
111	greg	2.2	if (dcp->hdrflags & HF_STDERR) {
112			fputs(dcp->inpname, stderr);
113			fputs(": bad reference depth in input header\n",
114			stderr);
115			}
116	greg	2.1	return -1;
117			}
118	greg	2.9	if (dcp->hdrflags & HF_ENCODE)
119			return 0; /* will add this later */
120	greg	2.6	} else if (!strncmp(s, "SAMP360=", 8))
121			dcp->gotview--;
122			else if (isview(s)) /* get view params */
123	greg	2.1	dcp->gotview += (sscanview(&dcp->vw, s) > 0);
124			if (dcp->hdrflags & HF_HEADOUT)
125			fputs(s, stdout); /* copy to standard output */
126			return 1;
127			}
128
129
130			/* Load/copy header */
131			int
132			process_dc_header(DEPTHCODEC dcp, int ac, char av[])
133			{
134			if (dcp->hdrflags & HF_HEADIN &&
135			getheader(dcp->finp, headline, dcp) < 0) {
136	greg	2.2	if (dcp->hdrflags & HF_STDERR) {
137			fputs(dcp->inpname, stderr);
138			fputs(": bad header\n", stderr);
139			}
140			return 0;
141	greg	2.1	}
142	greg	2.6	dcp->gotview *= (dcp->gotview > 0);
143	greg	2.12	/* get resolution string? */
144	greg	2.13	if (dcp->hdrflags & HF_RESIN &&
145			(dcp->res.xr <= 0) \| (dcp->res.yr <= 0) &&
146			!fgetsresolu(&dcp->res, dcp->finp)) {
147	greg	2.12	if (dcp->hdrflags & HF_STDERR) {
148			fputs(dcp->inpname, stderr);
149			fputs(": bad resolution string\n", stderr);
150			}
151			return 0;
152			}
153	greg	2.1	if (dcp->hdrflags & HF_HEADOUT) { /* finish header */
154			if (!(dcp->hdrflags & HF_HEADIN))
155			newheader("RADIANCE", stdout);
156			if (ac > 0)
157			printargs(ac, av, stdout);
158			if (dcp->hdrflags & HF_ENCODE) {
159			fputs(DEPTHSTR, stdout);
160			fputs(dcp->depth_unit, stdout);
161			fputc('\n', stdout);
162			fputformat(DEPTH16FMT, stdout);
163	greg	2.12	} else {
164			fputs("NCOMP=1\n", stdout);
165			if ((dcp->hdrflags & (HF_RESIN\|HF_RESOUT)) == HF_RESIN)
166			printf("NCOLS=%d\nNROWS=%d\n",
167			scanlen(&dcp->res),
168			numscans(&dcp->res));
169	greg	2.1	switch (dcp->format) {
170			case 'a':
171			fputformat("ascii", stdout);
172			break;
173			case 'f':
174	greg	2.4	fputendian(stdout);
175	greg	2.1	fputformat("float", stdout);
176			break;
177			case 'd':
178	greg	2.4	fputendian(stdout);
179	greg	2.1	fputformat("double", stdout);
180			break;
181			}
182	greg	2.12	}
183	greg	2.1	fputc('\n', stdout);
184			}
185	greg	2.12	if (dcp->hdrflags & HF_RESOUT) /* put resolution string? */
186	greg	2.1	fputsresolu(&dcp->res, stdout);
187
188			dcp->dstart = dcp->curpos = ftell(dcp->finp);
189			return 1;
190			}
191
192
193			/* Check that we have what we need to decode depths */
194			int
195			check_decode_depths(DEPTHCODEC *dcp)
196			{
197			if (dcp->hdrflags & HF_ENCODE) {
198	greg	2.2	if (dcp->hdrflags & HF_STDERR) {
199			fputs(progname, stderr);
200			fputs(": wrong header mode for decode\n", stderr);
201			}
202	greg	2.1	return 0;
203			}
204			if (dcp->inpfmt[0] && strcmp(dcp->inpfmt, DEPTH16FMT)) {
205	greg	2.2	if (dcp->hdrflags & HF_STDERR) {
206			fputs(dcp->inpname, stderr);
207			fputs(": unexpected input format: ", stderr);
208			fputs(dcp->inpfmt, stderr);
209			fputc('\n', stderr);
210			}
211	greg	2.1	return 0;
212			}
213			return 1;
214			}
215
216
217			/* Check that we have what we need to decode world positions */
218			int
219			check_decode_worldpos(DEPTHCODEC *dcp)
220			{
221			char *err;
222
223			if (!check_decode_depths(dcp))
224			return 0;
225			if ((dcp->res.xr <= 0) \| (dcp->res.yr <= 0)) {
226	greg	2.2	if (dcp->hdrflags & HF_STDERR) {
227			fputs(progname, stderr);
228			fputs(": missing map resolution\n", stderr);
229			}
230	greg	2.1	return 0;
231			}
232			if (!dcp->gotview) {
233	greg	2.2	if (dcp->hdrflags & HF_STDERR) {
234			fputs(dcp->inpname, stderr);
235			fputs(": missing view\n", stderr);
236			}
237	greg	2.1	return 0;
238			}
239			if ((err = setview(&dcp->vw)) != NULL) {
240	greg	2.2	if (dcp->hdrflags & HF_STDERR) {
241			fputs(dcp->inpname, stderr);
242			fputs(": input view error: ", stderr);
243			fputs(err, stderr);
244			fputc('\n', stderr);
245			}
246	greg	2.1	return 0;
247			}
248			return 1;
249			}
250
251
252			/* Decode next depth pixel */
253			double
254			decode_depth_next(DEPTHCODEC *dcp)
255			{
256	greg	2.5	int c;
257
258			if (dcp->use_last) {
259			dcp->use_last = 0;
260			return code2depth(dcp->last_dc, dcp->refdepth);
261			}
262			c = getint(2, dcp->finp);
263	greg	2.1
264			if (c == EOF && feof(dcp->finp))
265			return -1.;
266
267	greg	2.5	dcp->last_dc = c;
268	greg	2.1	dcp->curpos += 2;
269
270			return code2depth(c, dcp->refdepth);
271			}
272
273
274			/* Compute world position from depth */
275			int
276			compute_worldpos(FVECT wpos, DEPTHCODEC *dcp, int x, int y, double d)
277			{
278			RREAL loc[2];
279			FVECT rdir;
280
281	greg	2.11	if (d >= FHUGE*.99)
282			goto badval;
283
284	greg	2.1	pix2loc(loc, &dcp->res, x, y);
285
286	greg	2.11	if (viewray(wpos, rdir, &dcp->vw, loc[0], loc[1]) >= -FTINY) {
287			VSUM(wpos, wpos, rdir, d);
288			return 1;
289	greg	2.1	}
290	greg	2.11	badval:
291			VCOPY(wpos, dcp->vw.vp);
292			return 0;
293	greg	2.1	}
294
295
296			/* Decode the next world position */
297			int
298			decode_worldpos_next(FVECT wpos, DEPTHCODEC *dcp)
299			{
300			const long n = (dcp->curpos - dcp->dstart)>>1;
301			int x, y;
302			double d;
303
304			d = decode_depth_next(dcp);
305			if (d < -FTINY)
306			return -1;
307
308			x = scanlen(&dcp->res);
309			y = n / x;
310			x = n - (long)y*x;
311
312			return compute_worldpos(wpos, dcp, x, y, d);
313			}
314
315
316			/* Seek to the indicated pixel position */
317			int
318			seek_dc_pix(DEPTHCODEC *dcp, int x, int y)
319			{
320			long seekpos;
321
322			if ((dcp->res.xr <= 0) \| (dcp->res.yr <= 0)) {
323	greg	2.2	if (dcp->hdrflags & HF_STDERR) {
324			fputs(progname, stderr);
325			fputs(": need map resolution to seek\n", stderr);
326			}
327	greg	2.1	return -1;
328			}
329			if ((x < 0) \| (y < 0) \|\|
330			(x >= scanlen(&dcp->res)) \| (y >= numscans(&dcp->res))) {
331	greg	2.2	if (dcp->hdrflags & HF_STDERR) {
332			fputs(dcp->inpname, stderr);
333			fputs(": warning - pixel index off map\n", stderr);
334			}
335	greg	2.1	return 0;
336			}
337			seekpos = dcp->dstart + 2((long)yscanlen(&dcp->res) + x);
338
339	greg	2.5	if (seekpos == dcp->curpos-2) {
340			dcp->use_last++; /* avoids seek/read */
341			return 1;
342			}
343	greg	2.1	if (seekpos != dcp->curpos &&
344			fseek(dcp->finp, seekpos, SEEK_SET) == EOF) {
345	greg	2.2	if (dcp->hdrflags & HF_STDERR) {
346			fputs(dcp->inpname, stderr);
347			fputs(": seek error\n", stderr);
348			}
349	greg	2.1	return -1;
350			}
351			dcp->curpos = seekpos;
352	greg	2.10	dcp->use_last = 0;
353	greg	2.1	return 1;
354			}
355
356
357			/* Read and decode depth for the given pixel */
358			double
359			decode_depth_pix(DEPTHCODEC *dcp, int x, int y)
360			{
361			int rval = seek_dc_pix(dcp, x, y);
362
363			if (rval < 0)
364			return -1.;
365			if (!rval)
366			return .0;
367
368			return decode_depth_next(dcp);
369			}
370
371
372			/* Read and decode the world position at the given pixel */
373			int
374			get_worldpos_pix(FVECT wpos, DEPTHCODEC *dcp, int x, int y)
375			{
376			double d = decode_depth_pix(dcp, x, y);
377
378			if (d < -FTINY)
379			return -1;
380
381			return compute_worldpos(wpos, dcp, x, y, d);
382			}