ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/rcomb.c
Revision: 2.28
Committed: Fri Mar 28 00:06:36 2025 UTC (5 weeks ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.27: +3 -1 lines
Log Message:
refactor: Made MAXCOMP redefinable and common between rmatrix tools

File Contents

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