ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rcomb.c
Revision: 2.28
Committed: Fri Mar 28 00:06:36 2025 UTC (5 weeks, 1 day ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.27: +3 -1 lines
Log Message:
refactor: Made MAXCOMP redefinable and common between rmatrix tools

File Contents

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