src/cv/bsdfpeaks.c

#ifndef lint
static const char RCSid[] = "$Id: bsdfpeaks.c,v 2.5 2025/06/03 21:31:51 greg Exp $";
#endif
/*
 *  Compute minimum FWHM peak for each incident direction in SIR input.
 *  Report FWHM of corresponding peaks in XML representations if provided.
 */

#define _USE_MATH_DEFINES
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "bsdfrep.h"

typedef struct {
        float           peakv;          /* peak BSDF value */
        float           width;          /* smallest FWHM (deg) */
        const RBFNODE   *rbs;           /* incident system */
        int             ndx;            /* peak index for RBFVAL */
} FWHM;                 /* struct to hold peak value */

typedef double  eval_f(const FVECT vin, const FVECT vout, const void *p);

/* Comparison function to put larger peaks first */
int
cmpFWHM(const void *p0, const void *p1)
{
        float   diff = (*(const FWHM *)p0).peakv -
                        (*(const FWHM *)p1).peakv;

        if (diff > 0) return(-1);
        if (diff < 0) return(1);
        return(0);
}

/* BSDF evaluation function for RBF system */
double
rbf_eval(const FVECT vin, const FVECT vout, const void *p)
{
        /* XXX verify vin == p->invec ? */
        return(eval_rbfrep((const RBFNODE *)p, vout));
}

/* BSDF evaluation for XML input */
double
bsdf_eval(const FVECT vin, const FVECT vout, const void *p)
{
        SDValue sv;

        if (SDreportError(
                        SDevalBSDF(&sv, vin, vout, (const SDData *)p),
                        stderr))
                exit(1);

        return(sv.cieY);
}

/* Find full-width, half-maximum in radians around BSDF direction */
double
getFWHM(const FVECT vin, const FVECT vc, double rad0, eval_f *ev, const void *p)
{
        const double    peakv = (*ev)(vin, vc, p);
        double          rad1 = rad0;                    /* current radii */

        while (rad0 < M_PI/2.) {                        /* look for FWHM */
            FVECT       v0, vt;
            double      phi;
            v0[0] = 1; v0[1] = v0[2] = 0;
            geodesic(v0, vc, v0, rad0, GEOD_RAD);       /* use vc as pivot */
            for (phi = 0; phi < 2.*M_PI; phi += M_PI/18.) {
                spinvector(vt, v0, vc, phi);
                if ((*ev)(vin, vt, p) <= .5*peakv) {    /* found one side? */
                    FVECT       vt1;
                    while (rad1 < M_PI/2.) {            /* bracket peak */
                        geodesic(vt1, vt, vc, rad0+rad1, GEOD_RAD);
                        if ((*ev)(vin, vt1, p) <= .5*peakv)
                            return(rad0+rad1);          /* got both! */
                        rad1 *= 1.05;                   /* else bump rad1 */
                    }
                }
            }
            rad1 = rad0 *= 1.05;                        /* or expand search */
        }
        return(M_PI);                   /* failure return */
}

/* Get outgoing direction for the given FWHM record */
void
getOutDir(FVECT vo, FWHM *dp)
{
        const RBFVAL    *vp = dp->rbs->rbfa + dp->ndx;

        ovec_from_pos(vo, vp->gx, vp->gy);
}

/* Assign FWHM record for specified RBF system */
void
assignFWHM(FWHM *dp, const RBFNODE *rbf)
{
        FVECT   vo;
        int     j;
        double  rad;

        dp->rbs = rbf;
        dp->ndx = 0;            /* find peak outgoing */
        for (j = rbf->nrbf; --j; )
                if (rbf->rbfa[j].peak > rbf->rbfa[dp->ndx].peak)
                        dp->ndx = j;
                                /* record peak */
        getOutDir(vo, dp);
        dp->peakv = eval_rbfrep(rbf, vo);
                                /* get FWHM angle in degrees */
        dp->width = 180./M_PI * getFWHM(rbf->invec, vo,
                                        R2ANG(rbf->rbfa[dp->ndx].crad),
                                        rbf_eval, rbf);
}

/* Evaluate FWHM for each incident direction recorded in SIR */
int
main(int argc, char *argv[])
{
        const RBFNODE   *rbf;
        SDData          *sdp;
        FILE            *fp;
        int             ndirs;
        FWHM            *peaka;
        int             i;
                                                /* set global progname */
        fixargv0(argv[0]);
        if (argc < 2)
                goto userr;

        fp = fopen(argv[1], "rb");              /* load SIR input */
        if (fp == NULL) {
                fprintf(stderr, "%s: cannot open BSDF interpolant '%s'\n",
                                progname, argv[1]);
                return(1);
        }
        if (!load_bsdf_rep(fp))
                return(1);
        fclose(fp);
        for (i = 2; i < argc; i++)              /* check/load any XMLs */
                if (SDcacheFile(argv[i]) == NULL)
                        return(1);
        ndirs = 0;                              /* count input directions */
        for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
                if (rbf->nrbf <= 0) {
                        ndirs = 0;
                        break;
                }
                ++ndirs;
        }
        if (!ndirs) {
                fprintf(stderr, "%s: missing/bad RBFs in '%s'\n", progname, argv[1]);
                return(1);
        }
                                                /* print output header */
        printf("%d incident directions in '%s': %s -> %s\n", ndirs, argv[1],
                                input_orient>0 ? "Front" : "Back",
                                output_orient>0 ? "Front" : "Back");
        fputs("Incident (theta, phi)\tExiting (theta, phi)\tPeak\tFWHM", stdout);
        for (i = 2; i < argc; i++)
                printf("\t'%s'", argv[i]);
        fputc('\n', stdout);
                                                /* find SIR peaks */
        peaka = (FWHM *)malloc(sizeof(FWHM)*ndirs);
        if (peaka == NULL) return(1);
        for (i = 0, rbf = dsf_list; i < ndirs; i++, rbf = rbf->next)
                assignFWHM(&peaka[i], rbf);
                                                /* sort strong to weak */
        qsort(peaka, ndirs, sizeof(FWHM), cmpFWHM);

        for (i = 0; i < ndirs; i++) {           /* report FWHM for each incidence */
                FVECT   vout;
                int     j;
                getOutDir(vout, &peaka[i]);
                printf("%.0f %.0f\t%.0f %.0f", get_theta180(peaka[i].rbs->invec),
                                get_phi360(peaka[i].rbs->invec),
                                get_theta180(vout), get_phi360(vout));
                                                /* peak and FWHM from SIR */
                printf("\t%.2e\t%.1f", peaka[i].peakv, peaka[i].width);
                                                /* FWHM for each XML */
                for (j = 2; j < argc; j++) {
                        const SDData    *sd = SDcacheFile(argv[j]);
                        double  psa;
                        if (SDreportError(
                                        SDsizeBSDF(&psa, peaka[i].rbs->invec,
                                                NULL, SDqueryMin, sd),
                                        stderr))
                                return(1);

                        printf("\t%.1f", 180./M_PI * getFWHM(peaka[i].rbs->invec,
                                                vout, sqrt(psa/M_PI),
                                                bsdf_eval, sd));
                        SDfreeCache(sd);
                }
                fputc('\n', stdout);
        }
        /*                      we're exiting, anyway...
        SDfreeCache(NULL);
        clear_bsdf_rep();
        */
        return(0);
userr:
        fprintf(stderr, "Usage: %s bsdf.sir [bsdfrep1.xml ..]\n", progname);
        return(1);
}
Revision:	2.6
Committed:	Sat Jun 7 05:09:45 2025 UTC (10 hours, 36 minutes ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 2.5:	+1 -2 lines
Log Message:	refactor: Put some declarations into "paths.h" and included in "platform.h"
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.6	static const char RCSid[] = "$Id: bsdfpeaks.c,v 2.5 2025/06/03 21:31:51 greg Exp $";
3	greg	2.1	#endif
4			/*
5			* Compute minimum FWHM peak for each incident direction in SIR input.
6			* Report FWHM of corresponding peaks in XML representations if provided.
7			*/
8
9			#define _USE_MATH_DEFINES
10			#include <stdio.h>
11			#include <stdlib.h>
12			#include <math.h>
13			#include "bsdfrep.h"
14
15			typedef struct {
16	greg	2.3	float peakv; /* peak BSDF value */
17			float width; /* smallest FWHM (deg) */
18			const RBFNODE rbs; / incident system */
19			int ndx; /* peak index for RBFVAL */
20	greg	2.1	} FWHM; /* struct to hold peak value */
21
22	greg	2.2	typedef double eval_f(const FVECT vin, const FVECT vout, const void *p);
23	greg	2.1
24			/* Comparison function to put larger peaks first */
25			int
26			cmpFWHM(const void p0, const void p1)
27			{
28			float diff = ((const FWHM )p0).peakv -
29			((const FWHM )p1).peakv;
30
31			if (diff > 0) return(-1);
32			if (diff < 0) return(1);
33			return(0);
34			}
35
36			/* BSDF evaluation function for RBF system */
37			double
38	greg	2.2	rbf_eval(const FVECT vin, const FVECT vout, const void *p)
39	greg	2.1	{
40			/* XXX verify vin == p->invec ? */
41	greg	2.2	return(eval_rbfrep((const RBFNODE *)p, vout));
42	greg	2.1	}
43
44			/* BSDF evaluation for XML input */
45			double
46	greg	2.2	bsdf_eval(const FVECT vin, const FVECT vout, const void *p)
47	greg	2.1	{
48			SDValue sv;
49
50			if (SDreportError(
51	greg	2.2	SDevalBSDF(&sv, vin, vout, (const SDData *)p),
52	greg	2.1	stderr))
53			exit(1);
54
55			return(sv.cieY);
56			}
57
58			/* Find full-width, half-maximum in radians around BSDF direction */
59			double
60	greg	2.2	getFWHM(const FVECT vin, const FVECT vc, double rad0, eval_f ev, const void p)
61	greg	2.1	{
62			const double peakv = (*ev)(vin, vc, p);
63			double rad1 = rad0; /* current radii */
64
65			while (rad0 < M_PI/2.) { /* look for FWHM */
66			FVECT v0, vt;
67			double phi;
68			v0[0] = 1; v0[1] = v0[2] = 0;
69			geodesic(v0, vc, v0, rad0, GEOD_RAD); /* use vc as pivot */
70			for (phi = 0; phi < 2.*M_PI; phi += M_PI/18.) {
71			spinvector(vt, v0, vc, phi);
72			if ((ev)(vin, vt, p) <= .5peakv) { /* found one side? */
73			FVECT vt1;
74	greg	2.2	while (rad1 < M_PI/2.) { /* bracket peak */
75	greg	2.1	geodesic(vt1, vt, vc, rad0+rad1, GEOD_RAD);
76			if ((ev)(vin, vt1, p) <= .5peakv)
77			return(rad0+rad1); /* got both! */
78			rad1 = 1.05; / else bump rad1 */
79			}
80			}
81			}
82			rad1 = rad0 = 1.05; / or expand search */
83			}
84			return(M_PI); /* failure return */
85			}
86
87			/* Get outgoing direction for the given FWHM record */
88			void
89			getOutDir(FVECT vo, FWHM *dp)
90			{
91	greg	2.3	const RBFVAL *vp = dp->rbs->rbfa + dp->ndx;
92	greg	2.1
93			ovec_from_pos(vo, vp->gx, vp->gy);
94			}
95
96			/* Assign FWHM record for specified RBF system */
97			void
98	greg	2.3	assignFWHM(FWHM dp, const RBFNODE rbf)
99	greg	2.1	{
100			FVECT vo;
101			int j;
102			double rad;
103
104			dp->rbs = rbf;
105			dp->ndx = 0; /* find peak outgoing */
106			for (j = rbf->nrbf; --j; )
107			if (rbf->rbfa[j].peak > rbf->rbfa[dp->ndx].peak)
108			dp->ndx = j;
109			/* record peak */
110			getOutDir(vo, dp);
111			dp->peakv = eval_rbfrep(rbf, vo);
112			/* get FWHM angle in degrees */
113			dp->width = 180./M_PI * getFWHM(rbf->invec, vo,
114			R2ANG(rbf->rbfa[dp->ndx].crad),
115			rbf_eval, rbf);
116			}
117
118			/* Evaluate FWHM for each incident direction recorded in SIR */
119			int
120			main(int argc, char *argv[])
121			{
122	greg	2.3	const RBFNODE *rbf;
123			SDData *sdp;
124			FILE *fp;
125			int ndirs;
126			FWHM *peaka;
127			int i;
128	greg	2.5	/* set global progname */
129			fixargv0(argv[0]);
130	greg	2.1	if (argc < 2)
131			goto userr;
132
133			fp = fopen(argv[1], "rb"); /* load SIR input */
134			if (fp == NULL) {
135			fprintf(stderr, "%s: cannot open BSDF interpolant '%s'\n",
136			progname, argv[1]);
137			return(1);
138			}
139			if (!load_bsdf_rep(fp))
140			return(1);
141			fclose(fp);
142			for (i = 2; i < argc; i++) /* check/load any XMLs */
143			if (SDcacheFile(argv[i]) == NULL)
144			return(1);
145			ndirs = 0; /* count input directions */
146			for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
147			if (rbf->nrbf <= 0) {
148			ndirs = 0;
149			break;
150			}
151			++ndirs;
152			}
153			if (!ndirs) {
154			fprintf(stderr, "%s: missing/bad RBFs in '%s'\n", progname, argv[1]);
155			return(1);
156			}
157			/* print output header */
158			printf("%d incident directions in '%s': %s -> %s\n", ndirs, argv[1],
159			input_orient>0 ? "Front" : "Back",
160			output_orient>0 ? "Front" : "Back");
161			fputs("Incident (theta, phi)\tExiting (theta, phi)\tPeak\tFWHM", stdout);
162			for (i = 2; i < argc; i++)
163			printf("\t'%s'", argv[i]);
164			fputc('\n', stdout);
165			/* find SIR peaks */
166			peaka = (FWHM )malloc(sizeof(FWHM)ndirs);
167			if (peaka == NULL) return(1);
168			for (i = 0, rbf = dsf_list; i < ndirs; i++, rbf = rbf->next)
169			assignFWHM(&peaka[i], rbf);
170			/* sort strong to weak */
171			qsort(peaka, ndirs, sizeof(FWHM), cmpFWHM);
172
173			for (i = 0; i < ndirs; i++) { /* report FWHM for each incidence */
174			FVECT vout;
175			int j;
176			getOutDir(vout, &peaka[i]);
177			printf("%.0f %.0f\t%.0f %.0f", get_theta180(peaka[i].rbs->invec),
178			get_phi360(peaka[i].rbs->invec),
179			get_theta180(vout), get_phi360(vout));
180			/* peak and FWHM from SIR */
181			printf("\t%.2e\t%.1f", peaka[i].peakv, peaka[i].width);
182			/* FWHM for each XML */
183			for (j = 2; j < argc; j++) {
184			const SDData *sd = SDcacheFile(argv[j]);
185			double psa;
186			if (SDreportError(
187			SDsizeBSDF(&psa, peaka[i].rbs->invec,
188			NULL, SDqueryMin, sd),
189			stderr))
190			return(1);
191
192			printf("\t%.1f", 180./M_PI * getFWHM(peaka[i].rbs->invec,
193			vout, sqrt(psa/M_PI),
194	greg	2.2	bsdf_eval, sd));
195	greg	2.4	SDfreeCache(sd);
196	greg	2.1	}
197			fputc('\n', stdout);
198			}
199	greg	2.4	/* we're exiting, anyway...
200			SDfreeCache(NULL);
201			clear_bsdf_rep();
202			*/
203	greg	2.1	return(0);
204			userr:
205			fprintf(stderr, "Usage: %s bsdf.sir [bsdfrep1.xml ..]\n", progname);
206			return(1);
207			}