ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rcomb.c
Revision: 2.24
Committed: Mon Jun 10 18:20:06 2024 UTC (3 days, 13 hours ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 2.23: +6 -4 lines
Log Message:
fix(rcomb): pay better attention to -w option

File Contents

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