ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rcomb.c
Revision: 2.12
Committed: Tue May 21 17:39:17 2024 UTC (11 months, 1 week ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.11: +22 -8 lines
Log Message:
fix(rcomb): Switched to group-signaling to make multi-processing more reliable

File Contents

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