ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rcomb.c
Revision: 2.14
Committed: Wed May 22 15:38:04 2024 UTC (3 weeks, 4 days ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.13: +4 -5 lines
Log Message:
perf(rcomb): Minor optimization of signal blocking during child input

File Contents

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