ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rcomb.c
Revision: 2.16
Committed: Thu May 23 15:48:44 2024 UTC (11 months, 1 week ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.15: +17 -27 lines
Log Message:
fix(rcomb): Reverted to original synchronization method as signals proved unreliable for multi-processing

File Contents

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