ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rcomb.c
Revision: 2.29
Committed: Fri Apr 4 18:06:48 2025 UTC (3 weeks, 6 days ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.28: +1 -5 lines
Log Message:
perf(rmtxop,rcomb): Switched default internal data type from double to float

File Contents

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