ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rcomb.c
Revision: 2.27
Committed: Sat Mar 1 01:37:49 2025 UTC (2 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.26: +9 -9 lines
Log Message:
refactor: Simplified some error messages

File Contents

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