src/cv/bsdfpeaks.c

#ifndef lint
static const char RCSid[] = "$Id: bsdfpeaks.c,v 2.1 2025/05/21 17:23:07 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) */
        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);

char    *progname;              /* needed by bsdfrep.c */

/* 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)
{
        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, 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[])
{
        RBFNODE *rbf;
        SDData  *sdp;
        FILE    *fp;
        int     ndirs;
        FWHM    *peaka;
        int     i;

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