src/cv/bsdfpeaks.c

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