ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rcomb.c
Revision: 2.31
Committed: Fri Apr 18 23:59:03 2025 UTC (13 days, 14 hours ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 2.30: +2 -2 lines
Log Message:
refactor(rmtxop,rcomb,pvsum): Removed BSDF library dependency where unneeded

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.31 static const char RCSid[] = "$Id: rcomb.c,v 2.30 2025/04/04 18:18:06 greg Exp $";
3 greg 2.1 #endif
4     /*
5     * General component matrix combiner, operating on a row at a time.
6 greg 2.18 *
7     * Multi-processing mode under Unix creates children that each work
8 greg 2.19 * on one input row at a time, fed by the original process. Final conversion
9 greg 2.20 * and output to stdout is sorted by last child while its siblings send it
10 greg 2.19 * their record calculations.
11 greg 2.1 */
12    
13     #include <math.h>
14     #include "platform.h"
15 greg 2.10 #include "rtprocess.h"
16 greg 2.1 #include "rtio.h"
17     #include "rmatrix.h"
18     #include "calcomp.h"
19    
20     #ifndef M_PI
21     #define M_PI 3.14159265358979323846
22     #endif
23    
24     /* Unary matrix operation(s) */
25     typedef struct {
26     double cmat[MAXCOMP*MAXCOMP]; /* component transformation */
27     double sca[MAXCOMP]; /* scalar coefficients */
28     const char *csym; /* symbolic coefficients */
29     short clen; /* number of coefficients */
30     short nsf; /* number of scalars */
31     } RUNARYOP;
32    
33     /* Input matrix */
34     typedef struct {
35     const char *inspec; /* input specification */
36     RUNARYOP preop; /* transform operation */
37     RMATRIX imx; /* input matrix header info */
38     RMATRIX *rmp; /* active single-row matrix */
39     FILE *infp; /* open input stream */
40     } ROPMAT;
41    
42     ROPMAT *mop = NULL; /* allocated input array */
43     int nall = 0; /* number allocated */
44     int nmats = 0; /* number of actual inputs */
45    
46     RMATRIX *mcat = NULL; /* final concatenation */
47     int mcat_last = 0; /* goes after trailing ops? */
48    
49 greg 2.2 int in_nrows; /* number of input rows (or 0) */
50     #define in_ncols (mop[0].rmp->ncols) /* number of input columns */
51 greg 2.1 #define in_ncomp (mop[0].rmp->ncomp) /* input #components */
52    
53     extern int nowarn; /* turn off warnings? */
54    
55     int cur_row; /* current input/output row */
56     int cur_col; /* current input/output column */
57     int cur_chan; /* if we're looping channels */
58    
59 greg 2.10 SUBPROC *cproc = NULL; /* child process array */
60     int nchildren = 0; /* # of child processes */
61     int inchild = -1; /* our child ID (-1: parent) */
62    
63 greg 2.22 extern int checksymbolic(ROPMAT *rop);
64 greg 2.1
65 greg 2.22 int
66 greg 2.1 split_input(ROPMAT *rop)
67     {
68     if (rop->rmp == &rop->imx && !(rop->rmp = rmx_copy(&rop->imx))) {
69     fputs("Out of memory in split_input()\n", stderr);
70     return(0);
71     }
72     rmx_reset(rop->rmp);
73     return(1);
74     }
75    
76     /* Check/set transform based on a reference input file */
77 greg 2.22 int
78 greg 2.1 checkreffile(ROPMAT *rop)
79     {
80     static const char *curRF = NULL;
81     static RMATRIX refm;
82     const int nc = rop->imx.ncomp;
83     int i;
84    
85     if (!curRF || strcmp(rop->preop.csym, curRF)) {
86     FILE *fp = fopen(rop->preop.csym, "rb");
87     if (!rmx_load_header(&refm, fp)) {
88     fprintf(stderr, "%s: cannot read info header\n",
89     rop->preop.csym);
90     curRF = NULL;
91     if (fp) fclose(fp);
92     return(0);
93     }
94     fclose(fp);
95     curRF = rop->preop.csym;
96     }
97     if (refm.ncomp == 3) {
98     rop->preop.csym = (refm.dtype == DTxyze) ? "XYZ" : "RGB";
99     return(checksymbolic(rop));
100     }
101     if (refm.ncomp == 2) {
102     fprintf(stderr, "%s: cannot convert to 2 components\n",
103     curRF);
104     return(0);
105     }
106     if (refm.ncomp == 1) {
107     rop->preop.csym = "Y"; /* XXX big assumption */
108     return(checksymbolic(rop));
109     }
110     if (refm.ncomp == nc &&
111     !memcmp(refm.wlpart, rop->imx.wlpart, sizeof(refm.wlpart)))
112     return(1); /* nothing to do */
113    
114     if ((nc <= 3) | (nc > MAXCSAMP) | (refm.ncomp > MAXCSAMP)) {
115     fprintf(stderr, "%s: cannot resample from %d to %d components\n",
116     curRF, nc, refm.ncomp);
117     return(0);
118     }
119     if (!split_input(rop)) /* get our own struct */
120     return(0);
121     rop->preop.clen = refm.ncomp * nc; /* compute spec to ref */
122    
123     for (i = 0; i < nc; i++) {
124     SCOLOR scstim, scresp;
125     int j;
126     memset(scstim, 0, sizeof(COLORV)*nc);
127     scstim[i] = 1.f;
128     convertscolor(scresp, refm.ncomp, refm.wlpart[0], refm.wlpart[3],
129     scstim, nc, rop->imx.wlpart[0], rop->imx.wlpart[3]);
130     for (j = refm.ncomp; j-- > 0; )
131     rop->preop.cmat[j*nc + i] = scresp[j];
132     }
133     /* remember new spectral params */
134     memcpy(rop->rmp->wlpart, refm.wlpart, sizeof(rop->rmp->wlpart));
135     rop->rmp->ncomp = refm.ncomp;
136     return(1);
137     }
138    
139     /* Compute conversion row from spectrum to one channel of RGB */
140 greg 2.22 void
141 greg 2.1 rgbrow(ROPMAT *rop, int r, int p)
142     {
143     const int nc = rop->imx.ncomp;
144     const float * wlp = rop->imx.wlpart;
145     int i;
146    
147     for (i = nc; i--; ) {
148     int nmEnd = wlp[0] + (wlp[3] - wlp[0])*i/nc;
149     int nmStart = wlp[0] + (wlp[3] - wlp[0])*(i+1)/nc;
150     COLOR crgb;
151     spec_rgb(crgb, nmStart, nmEnd);
152     rop->preop.cmat[r*nc+i] = crgb[p];
153     }
154     }
155    
156     /* Compute conversion row from spectrum to one channel of XYZ */
157 greg 2.22 void
158 greg 2.1 xyzrow(ROPMAT *rop, int r, int p)
159     {
160     const int nc = rop->imx.ncomp;
161     const float * wlp = rop->imx.wlpart;
162     int i;
163    
164     for (i = nc; i--; ) {
165     int nmEnd = wlp[0] + (wlp[3] - wlp[0])*i/nc;
166     int nmStart = wlp[0] + (wlp[3] - wlp[0])*(i+1)/nc;
167     COLOR cxyz;
168     spec_cie(cxyz, nmStart, nmEnd);
169     rop->preop.cmat[r*nc+i] = cxyz[p];
170     }
171     }
172    
173     /* Use the spectral sensitivity function to compute matrix coefficients */
174 greg 2.22 void
175 greg 2.25 sensrow(ROPMAT *rop, int r, double (*sf)(const SCOLOR sc, int ncs, const float wlpt[4]))
176 greg 2.1 {
177     const int nc = rop->imx.ncomp;
178     int i;
179    
180     for (i = nc; i--; ) {
181     SCOLOR sclr;
182     memset(sclr, 0, sizeof(COLORV)*nc);
183     sclr[i] = 1.f;
184     rop->preop.cmat[r*nc+i] = (*sf)(sclr, nc, rop->imx.wlpart);
185     }
186     }
187    
188     /* Check/set symbolic transform */
189 greg 2.22 int
190 greg 2.1 checksymbolic(ROPMAT *rop)
191     {
192     const int nc = rop->imx.ncomp;
193     const int dt = rop->imx.dtype;
194 greg 2.3 double cf = 1;
195 greg 2.1 int i, j;
196     /* check suffix => reference file */
197     if (strchr(rop->preop.csym, '.') > rop->preop.csym)
198     return(checkreffile(rop));
199    
200     if (nc < 3) {
201     fprintf(stderr, "%s: -c '%s' requires at least 3 components\n",
202     rop->inspec, rop->preop.csym);
203     return(0);
204     }
205     rop->preop.clen = strlen(rop->preop.csym) * nc;
206     if (rop->preop.clen > MAXCOMP*MAXCOMP) {
207     fprintf(stderr, "%s: -c '%s' results in too many components\n",
208     rop->inspec, rop->preop.csym);
209     return(0);
210     }
211     for (j = 0; rop->preop.csym[j]; j++) {
212     int comp = 0;
213     switch (rop->preop.csym[j]) {
214     case 'B':
215 greg 2.3 case 'b':
216 greg 2.1 ++comp;
217     /* fall through */
218     case 'G':
219 greg 2.3 case 'g':
220 greg 2.1 ++comp;
221     /* fall through */
222     case 'R':
223 greg 2.3 case 'r':
224     if (rop->preop.csym[j] <= 'Z')
225     cf = 1./WHTEFFICACY;
226 greg 2.1 if (dt == DTxyze) {
227     for (i = 3; i--; )
228 greg 2.3 rop->preop.cmat[j*nc+i] = cf*xyz2rgbmat[comp][i];
229 greg 2.1 } else if (nc == 3)
230     rop->preop.cmat[j*nc+comp] = 1.;
231     else
232     rgbrow(rop, j, comp);
233     break;
234     case 'Z':
235 greg 2.3 case 'z':
236 greg 2.1 ++comp;
237     /* fall through */
238     case 'Y':
239 greg 2.3 case 'y':
240 greg 2.1 ++comp;
241     /* fall through */
242     case 'X':
243 greg 2.3 case 'x':
244     if ((rop->preop.csym[j] <= 'Z') & (dt != DTxyze))
245     cf = WHTEFFICACY;
246 greg 2.1 if (dt == DTxyze) {
247     rop->preop.cmat[j*nc+comp] = 1.;
248     } else if (nc == 3) {
249     for (i = 3; i--; )
250     rop->preop.cmat[j*nc+i] =
251     rgb2xyzmat[comp][i];
252     } else if (comp == CIEY)
253     sensrow(rop, j, scolor2photopic);
254     else
255     xyzrow(rop, j, comp);
256    
257 greg 2.3 for (i = nc*(cf != 1); i--; )
258     rop->preop.cmat[j*nc+i] *= cf;
259 greg 2.1 break;
260     case 'S': /* scotopic (il)luminance */
261 greg 2.3 cf = WHTSCOTOPIC;
262     /* fall through */
263     case 's':
264 greg 2.1 sensrow(rop, j, scolor2scotopic);
265 greg 2.3 for (i = nc*(cf != 1); i--; )
266     rop->preop.cmat[j*nc+i] *= cf;
267 greg 2.1 break;
268     case 'M': /* melanopic (il)luminance */
269 greg 2.3 cf = WHTMELANOPIC;
270     /* fall through */
271     case 'm':
272 greg 2.1 sensrow(rop, j, scolor2melanopic);
273 greg 2.3 for (i = nc*(cf != 1); i--; )
274     rop->preop.cmat[j*nc+i] *= cf;
275 greg 2.1 break;
276     case 'A': /* average component */
277 greg 2.3 case 'a':
278 greg 2.1 for (i = nc; i--; )
279     rop->preop.cmat[j*nc+i] = 1./(double)nc;
280     break;
281     default:
282     fprintf(stderr, "%s: -c '%c' unsupported\n",
283     rop->inspec, rop->preop.csym[j]);
284     return(0);
285     }
286     }
287     if (!split_input(rop)) /* get our own struct */
288     return(0);
289     memcpy(rop->rmp->wlpart, WLPART, sizeof(rop->rmp->wlpart));
290     rop->rmp->ncomp = rop->preop.clen / nc;
291     /* decide on output type */
292 greg 2.4 if (!strcasecmp(rop->preop.csym, "XYZ")) {
293 greg 2.1 if (dt <= DTspec)
294     rop->rmp->dtype = DTxyze;
295 greg 2.4 } else if (!strcasecmp(rop->preop.csym, "RGB")) {
296 greg 2.1 if (dt <= DTspec)
297     rop->rmp->dtype = DTrgbe;
298     } else if (rop->rmp->dtype == DTspec)
299     rop->rmp->dtype = DTfloat;
300     return(1);
301     }
302    
303 greg 2.22 int
304 greg 2.1 get_component_xfm(ROPMAT *rop)
305     {
306     int i, j;
307    
308     if (rop->rmp != &rop->imx) { /* reset destination matrix */
309     rmx_free(rop->rmp);
310     rop->rmp = &rop->imx;
311     }
312     if (rop->preop.csym && /* symbolic transform? */
313     !checksymbolic(rop))
314     return(0);
315     /* undo exposure? */
316     if (fabs(1. - bright(rop->rmp->cexp)) > .025) {
317     if (rop->rmp->ncomp == 1)
318     rop->rmp->cexp[RED] = rop->rmp->cexp[GRN] =
319     rop->rmp->cexp[BLU] = bright(rop->rmp->cexp);
320     if (rop->preop.nsf <= 0) {
321     rop->preop.nsf = i = rop->rmp->ncomp;
322     while (i--)
323     rop->preop.sca[i] = 1.;
324     }
325     if (rop->preop.nsf == 1) {
326     if (rop->rmp->ncomp == 3) {
327     rop->preop.sca[2] = rop->preop.sca[1] =
328     rop->preop.sca[0];
329     rop->preop.nsf = 3;
330     } else
331     rop->preop.sca[0] /= bright(rop->rmp->cexp);
332     }
333     if (rop->preop.nsf == 3) {
334     opcolor(rop->preop.sca, /=, rop->rmp->cexp);
335     } else if (rop->preop.nsf > 3) { /* punt */
336     double mult = 1./bright(rop->rmp->cexp);
337     for (i = rop->preop.nsf; i--; )
338     rop->preop.sca[i] *= mult;
339     }
340     setcolor(rop->rmp->cexp, 1., 1., 1.);
341     }
342     if (rop->preop.clen > 0) { /* use component transform? */
343     if (rop->preop.clen % rop->imx.ncomp) {
344     fprintf(stderr, "%s: -c must have N x %d coefficients\n",
345     rop->inspec, rop->imx.ncomp);
346     return(0);
347     }
348     if (rop->preop.nsf > 0) { /* scale transform, instead */
349     if (rop->preop.nsf == 1) {
350     for (i = rop->preop.clen; i--; )
351     rop->preop.cmat[i] *= rop->preop.sca[0];
352     } else if (rop->preop.nsf*rop->imx.ncomp != rop->preop.clen) {
353     fprintf(stderr, "%s: -s must have one or %d factors\n",
354     rop->inspec,
355     rop->preop.clen/rop->imx.ncomp);
356     return(0);
357     } else {
358     for (i = rop->preop.nsf; i--; )
359     for (j = rop->imx.ncomp; j--; )
360     rop->preop.cmat[i*rop->imx.ncomp+j]
361     *= rop->preop.sca[i];
362     }
363     }
364     rop->preop.nsf = 0; /* now folded in */
365     if (!split_input(rop)) /* get our own struct */
366     return(0);
367     rop->rmp->ncomp = rop->preop.clen / rop->imx.ncomp;
368     if ((rop->rmp->ncomp > 3) & (rop->rmp->dtype <= DTspec)) {
369     rop->rmp->dtype = DTfloat; /* probably not actual spectrum */
370     memcpy(rop->rmp->wlpart, WLPART, sizeof(rop->rmp->wlpart));
371     }
372     } else if (rop->preop.nsf > 0) { /* else use scalar(s)? */
373     if (rop->preop.nsf == 1) {
374     for (i = rop->rmp->ncomp; --i; )
375     rop->preop.sca[i] = rop->preop.sca[0];
376     rop->preop.nsf = rop->rmp->ncomp;
377     } else if (rop->preop.nsf != rop->rmp->ncomp) {
378     fprintf(stderr, "%s: -s must have one or %d factors\n",
379     rop->inspec, rop->rmp->ncomp);
380     return(0);
381     }
382     }
383     return(1);
384     }
385    
386 greg 2.22 int
387 greg 2.1 apply_op(RMATRIX *dst, const RMATRIX *src, const RUNARYOP *ro)
388     {
389     if (ro->clen > 0) {
390     RMATRIX *res = rmx_transform(src, dst->ncomp, ro->cmat);
391     if (!res) {
392     fputs("Error in call to rmx_transform()\n", stderr);
393     return(0);
394     }
395     if (!rmx_transfer_data(dst, res, 0))
396     return(0);
397     rmx_free(res);
398     } else if (dst != src)
399 greg 2.9 memcpy(dst->mtx, src->mtx, rmx_array_size(dst));
400 greg 2.1 if (ro->nsf == dst->ncomp)
401     rmx_scale(dst, ro->sca);
402     return(1);
403     }
404    
405 greg 2.22 int
406 greg 2.1 open_input(ROPMAT *rop)
407     {
408     int outtype;
409    
410     if (!rop || !rop->inspec || !rop->inspec[0])
411     return(0);
412     if (rop->inspec == stdin_name)
413     rop->infp = stdin;
414     else if (rop->inspec[0] == '!')
415     rop->infp = popen(rop->inspec+1, "r");
416     else
417     rop->infp = fopen(rop->inspec, "rb");
418    
419     if (!rmx_load_header(&rop->imx, rop->infp)) {
420     fprintf(stderr, "Bad header from: %s\n", rop->inspec);
421     return(0);
422     }
423     return(get_component_xfm(rop));
424     }
425    
426     /* Return nominal wavelength associated with input component (return nm) */
427 greg 2.22 double
428 greg 2.1 l_wavelength(char *nam)
429     {
430     double comp = argument(1);
431    
432     if ((comp < -.5) | (comp >= in_ncomp+.5)) {
433     errno = EDOM;
434     return(.0);
435     }
436     if (comp < .5) /* asking for #components? */
437     return(in_ncomp);
438    
439     if (in_ncomp == 3) { /* special case for RGB */
440     const int w0 = (int)(comp - .5);
441     return(mop[0].rmp->wlpart[w0] +
442     (comp-.5)*(mop[0].rmp->wlpart[w0+1] -
443     mop[0].rmp->wlpart[w0]));
444     }
445     return(mop[0].rmp->wlpart[0] + /* general case, even div. */
446     (comp-.5)/(double)in_ncomp *
447     (mop[0].rmp->wlpart[3] - mop[0].rmp->wlpart[0]));
448     }
449    
450     /* Return ith input with optional channel selector */
451 greg 2.22 double
452 greg 2.1 l_chanin(char *nam)
453     {
454     double inp = argument(1);
455     int mi, chan;
456    
457     if ((mi = (int)(inp-.5)) < 0 || mi >= nmats) {
458     errno = EDOM;
459     return(.0);
460     }
461     if (inp < .5) /* asking for #inputs? */
462     return(nmats);
463    
464     if (nargum() >= 2) {
465     double cval = argument(2);
466     if (cval < .5 || (chan = (int)(cval-.5)) >= in_ncomp) {
467     errno = EDOM;
468     return(.0);
469     }
470     } else
471     chan = cur_chan;
472    
473     return(mop[mi].rmp->mtx[cur_col*in_ncomp + chan]);
474     }
475    
476 greg 2.22 int
477 greg 2.1 initialize(RMATRIX *imp)
478     {
479     int i;
480     /* XXX struct is zeroed coming in */
481     setcolor(imp->cexp, 1.f, 1.f, 1.f);
482     for (i = 0; i < nmats; i++) { /* open each input */
483     int restype;
484     if (!open_input(&mop[i]))
485     return(0);
486     restype = mop[i].rmp->dtype;
487     if (!imp->dtype || (restype = rmx_newtype(restype, imp->dtype)) > 0)
488     imp->dtype = restype;
489 greg 2.24 else if (!nowarn)
490 greg 2.1 fprintf(stderr, "%s: warning - data type mismatch\n",
491     mop[i].inspec);
492     if (!i) {
493     imp->ncols = mop[0].rmp->ncols;
494     imp->ncomp = mop[0].rmp->ncomp;
495     memcpy(imp->wlpart, mop[0].rmp->wlpart, sizeof(imp->wlpart));
496 greg 2.2 } else if ((mop[i].rmp->ncols != imp->ncols) |
497     (mop[i].rmp->ncomp != imp->ncomp) |
498     ((in_nrows > 0) & (mop[i].rmp->nrows > 0) &
499     (mop[i].rmp->nrows != in_nrows))) {
500 greg 2.1 fprintf(stderr, "%s: mismatch in size or #components\n",
501     mop[i].inspec);
502     return(0);
503     } /* XXX should check wlpart? */
504 greg 2.2 if (in_nrows <= 0)
505     in_nrows = imp->nrows = mop[i].rmp->nrows;
506 greg 2.1 } /* set up .cal environment */
507     esupport |= E_VARIABLE|E_FUNCTION|E_RCONST;
508     esupport &= ~(E_OUTCHAN|E_INCHAN);
509     varset("PI", ':', M_PI);
510     varset("nfiles", ':', nmats);
511     varset("nrows", ':', in_nrows);
512     varset("ncols", ':', in_ncols);
513     varset("ncomp", ':', in_ncomp);
514     varset("R", ':', 1.);
515     varset("G", ':', 2.);
516     varset("B", ':', 3.);
517     funset("wl", 1, ':', l_wavelength);
518     funset("ci", 1, '=', l_chanin);
519     scompile("ri(i)=ci(i,R);gi(i)=ci(i,G);bi(i)=ci(i,B)", NULL, 0);
520     return(1);
521     }
522    
523 greg 2.22 void
524 greg 2.1 output_headinfo(FILE *fp)
525     {
526     int i;
527    
528     for (i = 0; i < nmats; i++) {
529     const char *cp = mop[i].imx.info;
530     fputs(mop[i].inspec, fp);
531     fputs(":\n", fp);
532     if (!cp) continue;
533     while (*cp) {
534     if (*cp == '\n') {
535     cp++; /* avoid inadvertant terminus */
536     continue;
537     }
538     fputc('\t', fp); /* indent this input's info */
539     do
540     putc(*cp, fp);
541     while (*cp++ != '\n');
542     }
543     }
544     }
545    
546 greg 2.22 int
547 greg 2.10 spawned_children(int np)
548 greg 2.1 {
549 greg 2.10 int i, rv;
550    
551     #if defined(_WIN32) || defined(_WIN64)
552     if (np > 1) {
553 greg 2.24 if (!nowarn)
554     fputs("Warning: only one process under Windows\n", stderr);
555 greg 2.10 np = 1;
556     } else
557     #endif
558 greg 2.18 if ((in_nrows > 0) & (np*4 > in_nrows))
559     np = in_nrows/4;
560 greg 2.10 /* we'll be doing a row at a time */
561 greg 2.1 for (i = 0; i < nmats; i++) {
562 greg 2.10 mop[i].imx.nrows = 1;
563 greg 2.1 if (!rmx_prepare(&mop[i].imx))
564     goto memerror;
565     if (mop[i].rmp != &mop[i].imx) {
566     mop[i].rmp->nrows = 1;
567     if (!rmx_prepare(mop[i].rmp))
568     goto memerror;
569     }
570     }
571 greg 2.20 /* prep output row buffer(s) */
572 greg 2.23 if (mop[nmats].preop.clen > 0) {
573 greg 2.10 if (!split_input(&mop[nmats])) /* need separate buffer */
574 greg 2.23 goto memerror;
575     mop[nmats].rmp->ncomp = mop[nmats].preop.clen /
576     mop[nmats].imx.ncomp;
577 greg 2.1 }
578 greg 2.10 mop[nmats].imx.nrows = 1;
579     if (!rmx_prepare(&mop[nmats].imx))
580     goto memerror;
581 greg 2.20 if (mop[nmats].rmp != &mop[nmats].imx) {
582     mop[nmats].rmp->nrows = 1;
583     if (!rmx_prepare(mop[nmats].rmp))
584     goto memerror;
585     }
586 greg 2.18 if (np <= 1) { /* single process return */
587 greg 2.10 #ifdef getc_unlocked
588     for (i = 0; i < nmats; i++)
589     flockfile(mop[i].infp);
590     flockfile(stdout);
591     #endif
592     return(0);
593     }
594     fflush(stdout); /* flush header & spawn children */
595 greg 2.20 nchildren = np + 1; /* extra child to sequence output */
596 greg 2.19 cproc = (SUBPROC *)malloc(sizeof(SUBPROC)*nchildren);
597 greg 2.10 if (!cproc)
598     goto memerror;
599 greg 2.20 for (i = nchildren; i--; ) cproc[i] = sp_inactive;
600     cproc[nchildren-1].flags |= PF_FILT_OUT;
601 greg 2.23 /* start each child from parent */
602     for (i = 0; i < nchildren; i++)
603     if ((rv = open_process(&cproc[i], NULL)) <= 0)
604     break; /* child breaks here */
605 greg 2.10 if (rv < 0) {
606 greg 2.23 perror("fork"); /* WTH? */
607 greg 2.17 close_processes(cproc, i);
608 greg 2.10 exit(1);
609     }
610 greg 2.23 if (i != nchildren-1) { /* last child is sole reader */
611     int j = i;
612     while (j-- > 0) {
613     close(cproc[j].r);
614     cproc[j].r = -1;
615 greg 2.21 }
616 greg 2.23 }
617     if (rv > 0)
618 greg 2.20 return(1); /* parent return value */
619 greg 2.23
620     inchild = i; /* else set our child index */
621 greg 2.20 while (i-- > 0) /* only parent writes siblings */
622 greg 2.11 close(cproc[i].w);
623 greg 2.20
624 greg 2.23 i = nmats; /* close matrix streams (carefully) */
625     while (i-- > 0) {
626     if (mop[i].infp != stdin) {
627     close(fileno(mop[i].infp)); /* avoid lseek() */
628     fclose(mop[i].infp); /* ! pclose() */
629     }
630     mop[i].infp = NULL;
631     }
632     fpurge(stdin); /* discard previously buffered input */
633    
634 greg 2.22 if (inchild == nchildren-1)
635     return(-1); /* output process return value */
636 greg 2.20
637 greg 2.23 i = nmats; /* get matrix rows from parent */
638 greg 2.20 while (i-- > 0) {
639     mop[i].infp = stdin;
640     mop[i].imx.dtype = DTrmx_native;
641     mop[i].imx.pflags &= ~RMF_SWAPIN;
642     }
643 greg 2.10 #ifdef getc_unlocked
644     flockfile(stdin);
645     #endif
646 greg 2.20 mop[nmats].rmp->dtype = DTrmx_native;
647     return(0); /* worker child return value */
648 greg 2.10 memerror:
649     fputs("Out of memory in spawned_children()\n", stderr);
650     exit(1);
651     }
652    
653 greg 2.22 int
654 greg 2.19 parent_loop(void)
655 greg 2.10 {
656     int i;
657    
658 greg 2.20 rmx_reset(&mop[nmats].imx); /* not touching output side */
659     if (mop[nmats].rmp != &mop[nmats].imx) {
660     rmx_free(mop[nmats].rmp);
661     mop[nmats].rmp = &mop[nmats].imx;
662     }
663 greg 2.10 #ifdef getc_unlocked
664 greg 2.20 for (i = 0; i < nmats; i++) /* we handle matrix inputs */
665 greg 2.10 flockfile(mop[i].infp);
666     #endif
667 greg 2.19 /* load & send rows to kids */
668 greg 2.10 for (cur_row = 0; (in_nrows <= 0) | (cur_row < in_nrows); cur_row++) {
669 greg 2.20 int wfd = cproc[cur_row % (nchildren-1)].w;
670 greg 2.11 for (i = 0; i < nmats; i++)
671 greg 2.10 if (!rmx_load_row(mop[i].imx.mtx, &mop[i].imx, mop[i].infp)) {
672     if (cur_row > in_nrows) /* unknown #input rows? */
673     break;
674 greg 2.27 fprintf(stderr, "%s: load error at row %d\n",
675 greg 2.10 mop[i].inspec, cur_row);
676     return(0);
677     }
678 greg 2.11 if (i < nmats)
679     break;
680     for (i = 0; i < nmats; i++)
681 greg 2.19 if (writebuf(wfd, mop[i].imx.mtx, rmx_array_size(&mop[i].imx))
682 greg 2.23 != rmx_array_size(&mop[i].imx)) {
683 greg 2.27 fprintf(stderr, "%s: write error at row %d\n",
684 greg 2.23 mop[i].inspec, cur_row);
685 greg 2.19 return(0);
686 greg 2.23 }
687 greg 2.10 }
688 greg 2.23 i = close_processes(cproc, nchildren); /* collect family */
689 greg 2.20 free(cproc); cproc = NULL; nchildren = 0;
690 greg 2.10 if (i < 0) {
691 greg 2.24 if (!nowarn)
692 greg 2.27 fputs("Warning: lost child process\n", stderr);
693 greg 2.10 return(1);
694     }
695     if (i > 0) {
696     fprintf(stderr, "Child exited with status %d\n", i);
697     return(0);
698     }
699 greg 2.20 return(1); /* return success! */
700 greg 2.10 }
701    
702 greg 2.22 int
703 greg 2.19 combine_input(void)
704 greg 2.10 {
705 greg 2.20 const int row0 = (inchild >= 0)*inchild;
706     const int rstep = nchildren ? nchildren-1 : 1;
707 greg 2.10 ROPMAT *res = &mop[nmats];
708     int set_r, set_c;
709     RMATRIX *tmp = NULL;
710     int co_set;
711     int i;
712    
713 greg 2.23 if (mcat_last && !(tmp = rmx_alloc(1, res->imx.ncols, res->rmp->ncomp))) {
714     fputs("Out of buffer space in combine_input()\n", stderr);
715     return(0);
716     }
717 greg 2.1 /* figure out what the user set */
718     co_set = fundefined("co");
719     if (!co_set)
720     co_set = -vardefined("co");
721     if (!co_set & (in_ncomp == 3) && vardefined("ro") &&
722     vardefined("go") && vardefined("bo")) {
723     scompile("co(p)=select(p,ro,go,bo)", NULL, 0);
724     co_set = 1;
725     }
726     if (co_set) { /* set if user wants, didn't set */
727     set_r = varlookup("r") != NULL && !vardefined("r");
728     set_c = varlookup("c") != NULL && !vardefined("c");
729     } else /* save a little time */
730     set_r = set_c = 0;
731 greg 2.16 /* read/process row-by-row */
732 greg 2.10 for (cur_row = row0; (in_nrows <= 0) | (cur_row < in_nrows); cur_row += rstep) {
733 greg 2.1 RMATRIX *mres = NULL;
734 greg 2.14 for (i = 0; i < nmats; i++)
735 greg 2.1 if (!rmx_load_row(mop[i].imx.mtx, &mop[i].imx, mop[i].infp)) {
736 greg 2.5 if (cur_row > in_nrows) /* unknown #input rows? */
737 greg 2.10 break;
738 greg 2.27 fprintf(stderr, "%s: load error at row %d\n",
739 greg 2.1 mop[i].inspec, cur_row);
740     return(0);
741     }
742 greg 2.10 if (i < nmats)
743     break;
744     for (i = 0; i < nmats; i++)
745 greg 2.1 if (!apply_op(mop[i].rmp, &mop[i].imx, &mop[i].preop))
746     return(0);
747     if (set_r) varset("r", '=', cur_row);
748     for (cur_col = 0; cur_col < in_ncols; cur_col++) {
749     if (set_c) varset("c", '=', cur_col);
750     for (cur_chan = 0; cur_chan < in_ncomp; cur_chan++) {
751     const int ndx = cur_col*in_ncomp + cur_chan;
752     eclock++;
753     if (!co_set) { /* just summing elements? */
754     res->imx.mtx[ndx] = 0;
755     for (i = nmats; i--; )
756     res->imx.mtx[ndx] += mop[i].rmp->mtx[ndx];
757     } else if (co_set > 0) {
758     double dchan = cur_chan+1;
759     res->imx.mtx[ndx] = funvalue("co", 1, &dchan);
760     } else
761     res->imx.mtx[ndx] = varvalue("co");
762     }
763     } /* final conversions */
764     if (!mcat) {
765     if (!apply_op(res->rmp, &res->imx, &res->preop))
766     return(0);
767     } else if (mcat_last) {
768     if (!apply_op(tmp, &res->imx, &res->preop))
769     return(0);
770     mres = rmx_multiply(tmp, mcat);
771     if (!mres)
772     goto multerror;
773     if (!rmx_transfer_data(res->rmp, mres, 0))
774     return(0);
775     } else /* mcat && !mcat_last */ {
776     mres = rmx_multiply(&res->imx, mcat);
777     if (!mres)
778     goto multerror;
779     if (!apply_op(res->rmp, mres, &res->preop))
780     return(0);
781     }
782     rmx_free(mres); mres = NULL;
783     if (!rmx_write_data(res->rmp->mtx, res->rmp->ncomp,
784 greg 2.23 res->rmp->ncols, res->rmp->dtype, stdout) ||
785     (inchild >= 0 && fflush(stdout) == EOF)) {
786 greg 2.27 fprintf(stderr, "Conversion/write error at row %d\n",
787 greg 2.23 cur_row);
788 greg 2.16 return(0);
789 greg 2.23 }
790 greg 2.1 }
791 greg 2.20 return(inchild >= 0 || fflush(stdout) != EOF);
792 greg 2.1 multerror:
793 greg 2.27 fputs("Unexpected matrix multiply error\n", stderr);
794 greg 2.1 return(0);
795     }
796    
797 greg 2.22 int
798     output_loop(void)
799     {
800     const size_t row_size = rmx_array_size(mop[nmats].rmp);
801 greg 2.23 int cur_child = 0;
802 greg 2.22 int i = nmats;
803    
804 greg 2.23 while (i-- > 0) { /* free input buffers */
805 greg 2.22 rmx_reset(&mop[i].imx);
806     if (mop[i].rmp != &mop[i].imx) {
807     rmx_free(mop[i].rmp);
808     mop[i].rmp = &mop[i].imx;
809     }
810     }
811 greg 2.23 if (mop[nmats].rmp != &mop[nmats].imx) /* output is split? */
812     rmx_reset(&mop[nmats].imx);
813 greg 2.22 #ifdef getc_unlocked
814     flockfile(stdout); /* we own this, now */
815     #endif
816     for ( ; ; ) { /* loop until no more */
817     ssize_t rv;
818     rv = readbuf(cproc[cur_child].r, mop[nmats].rmp->mtx, row_size);
819     if (!rv) /* out of rows? */
820     break;
821     if (rv != row_size) {
822 greg 2.27 fputs("Read error\n", stderr);
823 greg 2.22 return(0);
824     } /* do final conversion */
825     if (!rmx_write_data(mop[nmats].rmp->mtx, mop[nmats].rmp->ncomp,
826     mop[nmats].rmp->ncols, mop[nmats].rmp->dtype, stdout)) {
827 greg 2.27 fputs("Conversion/write error\n", stderr);
828 greg 2.22 return(0);
829     }
830     cur_child++;
831 greg 2.23 cur_child *= (cur_child < inchild); /* loop over workers */
832 greg 2.22 }
833     return(fflush(stdout) != EOF);
834     }
835    
836     int
837 greg 2.1 get_factors(double da[], int n, char *av[])
838     {
839     int ac;
840    
841     for (ac = 0; ac < n && isflt(av[ac]); ac++)
842     da[ac] = atof(av[ac]);
843     return(ac);
844     }
845    
846 greg 2.22 void
847 greg 2.1 resize_inparr(int n2alloc)
848     {
849     int i;
850    
851 greg 2.6 if (n2alloc == nall)
852     return;
853 greg 2.26 for (i = nall; i-- > n2alloc; ) {
854 greg 2.1 rmx_reset(&mop[i].imx);
855     if (mop[i].rmp != &mop[i].imx)
856     rmx_free(mop[i].rmp);
857     }
858     mop = (ROPMAT *)realloc(mop, n2alloc*sizeof(ROPMAT));
859     if (mop == NULL) {
860     fputs("Out of memory in resize_inparr()\n", stderr);
861     exit(1);
862     }
863 greg 2.6 if (n2alloc > nall)
864     memset(mop+nall, 0, (n2alloc-nall)*sizeof(ROPMAT));
865 greg 2.1 nall = n2alloc;
866     }
867    
868     /* Load one or more matrices and operate on them, sending results to stdout */
869     int
870     main(int argc, char *argv[])
871     {
872    
873     int outfmt = DTfromHeader;
874     const char *defCsym = NULL;
875     int echoheader = 1;
876     int stdin_used = 0;
877 greg 2.10 int nproc = 1;
878 greg 2.1 const char *mcat_spec = NULL;
879     int n2comp = 0;
880     uby8 comp_ndx[128];
881     int i;
882     /* get starting input array */
883     mop = (ROPMAT *)calloc(nall=2, sizeof(ROPMAT));
884     /* get options and arguments */
885     for (i = 1; i < argc; i++)
886     if (argv[i][0] != '-' || !argv[i][1]) {
887     if (argv[i][0] == '-') {
888     if (stdin_used++) goto stdin_error;
889     mop[nmats].inspec = stdin_name;
890     } else
891     mop[nmats].inspec = argv[i];
892     if (!mop[nmats].preop.csym)
893     mop[nmats].preop.csym = defCsym;
894     if (++nmats >= nall)
895     resize_inparr(nmats + (nmats>>2) + 2);
896     } else {
897     int n = argc-1 - i;
898     switch (argv[i][1]) { /* get option */
899     case 'w':
900     nowarn = !nowarn;
901     break;
902     case 'h':
903     echoheader = !echoheader;
904     break;
905 greg 2.10 case 'n':
906     nproc = atoi(argv[++i]);
907     if (nproc <= 0)
908     goto userr;
909     break;
910 greg 2.1 case 'e':
911     if (!n) goto userr;
912     comp_ndx[n2comp++] = i++;
913     break;
914     case 'f':
915     switch (argv[i][2]) {
916     case '\0':
917     if (!n) goto userr;
918     comp_ndx[n2comp++] = i++;
919     break;
920     case 'd':
921     outfmt = DTdouble;
922     break;
923     case 'f':
924     outfmt = DTfloat;
925     break;
926     case 'a':
927     outfmt = DTascii;
928     break;
929     case 'c':
930     outfmt = DTrgbe;
931     break;
932     default:
933     goto userr;
934     }
935     break;
936     case 's':
937     if (n > MAXCOMP) n = MAXCOMP;
938     i += mop[nmats].preop.nsf =
939     get_factors(mop[nmats].preop.sca,
940     n, argv+i+1);
941     if (mop[nmats].preop.nsf <= 0) {
942     fprintf(stderr, "%s: -s missing arguments\n",
943     argv[0]);
944     goto userr;
945     }
946     break;
947     case 'C':
948 greg 2.23 mcat_last = 0;
949 greg 2.1 if (!n || isflt(argv[i+1]))
950     goto userr;
951     defCsym = mop[nmats].preop.csym = argv[++i];
952     mop[nmats].preop.clen = 0;
953     break;
954     case 'c':
955 greg 2.23 mcat_last = 0;
956 greg 2.1 if (n && !isflt(argv[i+1])) {
957     mop[nmats].preop.csym = argv[++i];
958     mop[nmats].preop.clen = 0;
959     break;
960     }
961     if (n > MAXCOMP*MAXCOMP) n = MAXCOMP*MAXCOMP;
962     i += mop[nmats].preop.clen =
963     get_factors(mop[nmats].preop.cmat,
964     n, argv+i+1);
965     if (mop[nmats].preop.clen <= 0) {
966     fprintf(stderr, "%s: -c missing arguments\n",
967     argv[0]);
968     goto userr;
969     }
970     mop[nmats].preop.csym = NULL;
971     break;
972     case 'm':
973 greg 2.23 mcat_last = 1;
974 greg 2.1 if (!n) goto userr;
975     if (argv[++i][0] == '-' && !argv[i][1]) {
976     if (stdin_used++) goto stdin_error;
977     mcat_spec = stdin_name;
978     } else
979     mcat_spec = argv[i];
980     break;
981     default:
982     fprintf(stderr, "%s: unknown option '%s'\n",
983     argv[0], argv[i]);
984     goto userr;
985     }
986     }
987     if (!nmats) {
988     fprintf(stderr, "%s: need at least one input matrix\n", argv[0]);
989     goto userr;
990     }
991     resize_inparr(nmats+1); /* extra matrix at end for result */
992     mop[nmats].inspec = "trailing_ops";
993     /* load final concatenation matrix */
994 greg 2.31 if (mcat_spec && !(mcat = rmx_load(mcat_spec))) {
995 greg 2.1 fprintf(stderr, "%s: error loading concatenation matrix: %s\n",
996     argv[0], mcat_spec);
997     return(1);
998     }
999     /* get/check inputs, set constants */
1000     if (!initialize(&mop[nmats].imx))
1001     return(1);
1002    
1003     for (i = 0; i < n2comp; i++) /* user .cal files and expressions */
1004     if (argv[comp_ndx[i]][1] == 'f') {
1005     char *fpath = getpath(argv[comp_ndx[i]+1],
1006     getrlibpath(), 0);
1007     if (fpath == NULL) {
1008     fprintf(stderr, "%s: cannot find file '%s'\n",
1009     argv[0], argv[comp_ndx[i]+1]);
1010     return(1);
1011     }
1012     fcompile(fpath);
1013     } else /* (argv[comp_ndx[i]][1] == 'e') */
1014     scompile(argv[comp_ndx[i]+1], NULL, 0);
1015    
1016     /* get trailing color transform */
1017     if (!get_component_xfm(&mop[nmats]))
1018     return(1);
1019     /* adjust output dimensions and #components */
1020     if (mcat) {
1021     if (mop[nmats].imx.ncols != mcat->nrows) {
1022     fprintf(stderr,
1023     "%s: number of input columns does not match number of rows in '%s'\n",
1024     argv[0], mcat_spec);
1025     return(1);
1026     }
1027     if (mcat->ncomp != (mcat_last ? mop[nmats].rmp->ncomp : mop[nmats].imx.ncomp)) {
1028     fprintf(stderr,
1029     "%s: number of components does not match those in '%s'\n",
1030     argv[0], mcat_spec);
1031     return(1);
1032     }
1033     if (!split_input(&mop[nmats]))
1034     return(1);
1035     mop[nmats].rmp->ncols = mcat->ncols;
1036     }
1037 greg 2.30 #if DTrmx_native==DTfloat
1038     if (outfmt == DTdouble)
1039     fprintf(stderr, "%s: warning - writing float result as double\n", argv[0]);
1040     #endif
1041 greg 2.1 newheader("RADIANCE", stdout); /* write output header */
1042     if (echoheader)
1043     output_headinfo(stdout);
1044     printargs(argc, argv, stdout);
1045     fputnow(stdout);
1046     mop[nmats].rmp->dtype = rmx_write_header(mop[nmats].rmp, outfmt, stdout);
1047     if (!mop[nmats].rmp->dtype) {
1048     fprintf(stderr, "%s: unsupported output format\n", argv[0]);
1049     return(1);
1050     }
1051 greg 2.7 doptimize(1); /* optimize definitions */
1052 greg 2.22 i = spawned_children(nproc); /* create multiple processes if requested */
1053     if (i > 0) /* running in parent process? */
1054 greg 2.10 return(parent_loop() ? 0 : 1);
1055 greg 2.22 if (i < 0) /* running in output process? */
1056     return(output_loop() ? 0 : 1);
1057     /* else we are a worker process */
1058 greg 2.10 return(combine_input() ? 0 : 1);
1059 greg 2.1 stdin_error:
1060     fprintf(stderr, "%s: %s used for more than one input\n",
1061     argv[0], stdin_name);
1062     return(1);
1063     userr:
1064     fprintf(stderr,
1065 greg 2.10 "Usage: %s [-h][-f{adfc}][-n nproc][-e expr][-f file][-s sf .. | -c ce ..] m1 .. -m mcat > mres\n",
1066 greg 2.1 argv[0]);
1067     return(1);
1068     }