ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rcomb.c
Revision: 2.18
Committed: Thu May 23 19:29:41 2024 UTC (10 days, 3 hours ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 2.17: +18 -6 lines
Error occurred while calculating annotation data.
Log Message:
fix(rcomb): Avoid inappropriate byte-swapping and improved final fclose() wait

File Contents

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