ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/Development/ray/src/util/rcomb.c
Revision: 2.35
Committed: Thu Oct 30 16:46:49 2025 UTC (5 days ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 2.34: +7 -2 lines
Log Message:
feat(rcomb): Added proper error message when input file cannot be opened

File Contents

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