src/px/pcond2.c

#ifndef lint
static const char       RCSid[] = "$Id: pcond2.c,v 3.15 2005/09/08 19:08:52 greg Exp $";
#endif
/*
 * Input and output conditioning routines for pcond.
 *  Added white-balance adjustment 10/01 (GW).
 */

#include "pcond.h"
#include "warp3d.h"


RGBPRIMP        inprims = stdprims;     /* input primaries */
COLORMAT        inrgb2xyz;              /* convert input RGB to XYZ */
RGBPRIMP        outprims = stdprims;    /* output primaries */

double  (*lumf)() = rgblum;             /* input luminance function */
double  inpexp = 1.0;                   /* input exposure value */

char    *mbcalfile = NULL;              /* macbethcal mapping file */
char    *cwarpfile = NULL;              /* color space warping file */

static struct mbc {
        COLORMAT        cmat;
        float   xa[3][6], ya[3][6];
        COLOR   cmin, cmax;
}       mbcond;                         /* macbethcal conditioning struct */

static WARP3D   *cwarp;                 /* color warping structure */

static COLOR    *scanbuf;               /* scanline processing buffer */
static int      nread;                  /* number of scanlines processed */

static void sfscan(COLOR *sl, int len, double sf);
static void matscan(COLOR *sl, int len, COLORMAT mat);
static void mbscan(COLOR *sl, int len, struct mbc *mb);
static void cwscan(COLOR *sl, int len, WARP3D *wp);
static void getmbcalfile(char *fn, struct mbc *mb);


extern double
rgblum(         /* compute (scotopic) luminance of RGB color */
        COLOR   clr,
        int     scotopic
)
{
        if (scotopic)           /* approximate */
                return( WHTSEFFICACY * (colval(clr,RED)*.062 +
                        colval(clr,GRN)*.608 + colval(clr,BLU)*.330) );
        return( WHTEFFICACY * (colval(clr,RED)*inrgb2xyz[1][0] +
                        colval(clr,GRN)*inrgb2xyz[1][1] +
                        colval(clr,BLU)*inrgb2xyz[1][2]) );
}


extern double
cielum(         /* compute (scotopic) luminance of CIE color */
        COLOR   xyz,
        int     scotopic
)
{
        if (scotopic)           /* approximate */
                return(colval(xyz,CIEY) *
                        (1.33*(1. + (colval(xyz,CIEY)+colval(xyz,CIEZ))/
                        colval(xyz,CIEX)) - 1.68));
        return(colval(xyz,CIEY));
}


extern COLOR *
nextscan(void)                          /* read and condition next scanline */
{
        if (nread >= numscans(&inpres)) {
                if (cwarpfile != NULL)
                        free3dw(cwarp);
                free((void *)scanbuf);
                return(scanbuf = NULL);
        }
        if (what2do&DO_ACUITY)
                acuscan(scanbuf, nread);
        else if (freadscan(scanbuf, scanlen(&inpres), infp) < 0) {
                fprintf(stderr, "%s: %s: scanline read error\n",
                                progname, infn);
                exit(1);
        }
        if (what2do&DO_VEIL)                    /* add veiling */
                addveil(scanbuf, nread);
        if (what2do&DO_COLOR)                   /* scotopic color loss */
                scotscan(scanbuf, scanlen(&inpres));
        if (what2do&DO_LINEAR)                  /* map luminances */
                sfscan(scanbuf, scanlen(&inpres), scalef);
        else
                mapscan(scanbuf, scanlen(&inpres));
        if (mbcalfile != NULL)                  /* device color correction */
                mbscan(scanbuf, scanlen(&inpres), &mbcond);
        else if (cwarpfile != NULL)             /* device color space warp */
                cwscan(scanbuf, scanlen(&inpres), cwarp);
        else if ((lumf == cielum) | (inprims != outprims))
                matscan(scanbuf, scanlen(&inpres), mbcond.cmat);
        nread++;
        return(scanbuf);
}


extern COLOR *
firstscan(void)                         /* return first processed scanline */
{
        if (mbcalfile != NULL)          /* load macbethcal file */
                getmbcalfile(mbcalfile, &mbcond);
        else if (cwarpfile != NULL) {
                if ((cwarp = load3dw(cwarpfile, NULL)) == NULL)
                        syserror(cwarpfile);
        } else
                if (lumf == rgblum)
                        comprgb2rgbWBmat(mbcond.cmat, inprims, outprims);
                else
                        compxyz2rgbWBmat(mbcond.cmat, outprims);
        if (what2do&DO_ACUITY)
                initacuity();
        scanbuf = (COLOR *)malloc(scanlen(&inpres)*sizeof(COLOR));
        if (scanbuf == NULL)
                syserror("malloc");
        nread = 0;
        return(nextscan());
}


static void
sfscan(                 /* apply scalefactor to scanline */
        register COLOR  *sl,
        int     len,
        double  sf
)
{
        while (len--) {
                scalecolor(sl[0], sf);
                sl++;
        }
}


static double
greypoint(                      /* compute gamut mapping grey target */
        COLOR col
)
{
        COLOR   gryc;
        int     i;
                                /* improves saturated color rendering */
        copycolor(gryc, col);
        for (i = 3; i--; )
                if (gryc[i] > cwhite[i])
                        gryc[i] = cwhite[i];
                else if (gryc[i] < cblack[i])
                        gryc[i] = cblack[i];
        return((*lumf)(gryc,0));
}


static void
matscan(                        /* apply color matrix to scaline */
        register COLOR  *sl,
        int     len,
        COLORMAT        mat
)
{
        double  gryv;

        while (len--) {
                gryv = greypoint(sl[0]);
                colortrans(sl[0], mat, sl[0]);
                clipgamut(sl[0], gryv, CGAMUT, cblack, cwhite);
                sl++;
        }
}


static void
mbscan(                 /* apply macbethcal adj. to scaline */
        COLOR   *sl,
        int     len,
        register struct mbc     *mb
)
{
        double  d;
        register int    i, j;

        while (len--) {
                d = greypoint(sl[0]);
                colortrans(sl[0], mb->cmat, sl[0]);
                clipgamut(sl[0], d, CGAMUT, mb->cmin, mb->cmax);
                for (i = 0; i < 3; i++) {
                        d = colval(sl[0],i);
                        for (j = 0; j < 4 && mb->xa[i][j+1] <= d; j++)
                                ;
                        colval(sl[0],i) = ( (mb->xa[i][j+1] - d)*mb->ya[i][j] +
                                        (d - mb->xa[i][j])*mb->ya[i][j+1] ) /
                                        (mb->xa[i][j+1] - mb->xa[i][j]);
                }
                sl++;
        }
}


static void
cwscan(                 /* apply color space warp to scaline */
        COLOR   *sl,
        int     len,
        WARP3D  *wp
)
{
        int     rval;
        double  gryv;

        while (len--) {
                gryv = greypoint(sl[0]);
                rval = warp3d(sl[0], sl[0], wp);
                if (rval & W3ERROR)
                        syserror("warp3d");
                if (rval & W3BADMAP) {
                        fprintf(stderr, "%s: %s: bad color space map\n",
                                        progname, cwarpfile);
                        exit(1);
                }
                clipgamut(sl[0], gryv, CGAMUT, cblack, cwhite);
                sl++;
        }
}


static void
getmbcalfile(                   /* load macbethcal file */
        char    *fn,
        register struct mbc     *mb
)
{
        char    buf[128];
        FILE    *fp;
        int     inpflags = 0;
        register int    i;

        if ((fp = fopen(fn, "r")) == NULL)
                syserror(fn);
        while (fgets(buf, sizeof(buf), fp) != NULL) {
                if (!(inpflags & 01) &&
                                sscanf(buf,
                                "rxa(i) : select(i,%f,%f,%f,%f,%f,%f)",
                                &mb->xa[0][0], &mb->xa[0][1],
                                &mb->xa[0][2], &mb->xa[0][3],
                                &mb->xa[0][4], &mb->xa[0][5]
                                ) == 6)
                        inpflags |= 01;
                else if (!(inpflags & 02) &&
                                sscanf(buf,
                                "rya(i) : select(i,%f,%f,%f,%f,%f,%f)",
                                &mb->ya[0][0], &mb->ya[0][1],
                                &mb->ya[0][2], &mb->ya[0][3],
                                &mb->ya[0][4], &mb->ya[0][5]
                                ) == 6)
                        inpflags |= 02;
                else if (!(inpflags & 04) &&
                                sscanf(buf,
                                "gxa(i) : select(i,%f,%f,%f,%f,%f,%f)",
                                &mb->xa[1][0], &mb->xa[1][1],
                                &mb->xa[1][2], &mb->xa[1][3],
                                &mb->xa[1][4], &mb->xa[1][5]
                                ) == 6)
                        inpflags |= 04;
                else if (!(inpflags & 010) &&
                                sscanf(buf,
                                "gya(i) : select(i,%f,%f,%f,%f,%f,%f)",
                                &mb->ya[1][0], &mb->ya[1][1],
                                &mb->ya[1][2], &mb->ya[1][3],
                                &mb->ya[1][4], &mb->ya[1][5]
                                ) == 6)
                        inpflags |= 010;
                else if (!(inpflags & 020) &&
                                sscanf(buf,
                                "bxa(i) : select(i,%f,%f,%f,%f,%f,%f)",
                                &mb->xa[2][0], &mb->xa[2][1],
                                &mb->xa[2][2], &mb->xa[2][3],
                                &mb->xa[2][4], &mb->xa[2][5]
                                ) == 6)
                        inpflags |= 020;
                else if (!(inpflags & 040) &&
                                sscanf(buf,
                                "bya(i) : select(i,%f,%f,%f,%f,%f,%f)",
                                &mb->ya[2][0], &mb->ya[2][1],
                                &mb->ya[2][2], &mb->ya[2][3],
                                &mb->ya[2][4], &mb->ya[2][5]
                                ) == 6)
                        inpflags |= 040;
                else if (!(inpflags & 0100) &&
                                sscanf(buf,
                                "ro = %f*rn + %f*gn + %f*bn",
                                &mb->cmat[0][0], &mb->cmat[0][1],
                                &mb->cmat[0][2]) == 3)
                        inpflags |= 0100;
                else if (!(inpflags & 0200) &&
                                sscanf(buf,
                                "go = %f*rn + %f*gn + %f*bn",
                                &mb->cmat[1][0], &mb->cmat[1][1],
                                &mb->cmat[1][2]) == 3)
                        inpflags |= 0200;
                else if (!(inpflags & 0400) &&
                                sscanf(buf,
                                "bo = %f*rn + %f*gn + %f*bn",
                                &mb->cmat[2][0], &mb->cmat[2][1],
                                &mb->cmat[2][2]) == 3)
                        inpflags |= 0400;
        }
        if (inpflags != 0777) {
                fprintf(stderr,
                "%s: cannot grok macbethcal file \"%s\" (inpflags==0%o)\n",
                        progname, fn, inpflags);
                exit(1);
        }
        fclose(fp);
                                        /* compute gamut */
        for (i = 0; i < 3; i++) {
                colval(mb->cmin,i) = mb->xa[i][0] -
                                mb->ya[i][0] *
                                (mb->xa[i][1]-mb->xa[i][0]) /
                                (mb->ya[i][1]-mb->ya[i][0]);
                colval(mb->cmax,i) = mb->xa[i][4] +
                                (1.-mb->ya[i][4]) *
                                (mb->xa[i][5] - mb->xa[i][4]) /
                                (mb->ya[i][5] - mb->ya[i][4]);
        }
}
Revision:	3.16
Committed:	Thu Sep 8 21:16:44 2005 UTC (18 years, 7 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 3.15:	+2 -2 lines
Log Message:	Bug fix in last change
#	User	Rev	Content
1	greg	3.1	#ifndef lint
2	greg	3.16	static const char RCSid[] = "$Id: pcond2.c,v 3.15 2005/09/08 19:08:52 greg Exp $";
3	greg	3.1	#endif
4			/*
5			* Input and output conditioning routines for pcond.
6	greg	3.10	* Added white-balance adjustment 10/01 (GW).
7	greg	3.1	*/
8
9			#include "pcond.h"
10	greg	3.8	#include "warp3d.h"
11	greg	3.1
12
13			RGBPRIMP inprims = stdprims; /* input primaries */
14			COLORMAT inrgb2xyz; /* convert input RGB to XYZ */
15			RGBPRIMP outprims = stdprims; /* output primaries */
16
17			double (lumf)() = rgblum; / input luminance function */
18			double inpexp = 1.0; /* input exposure value */
19
20			char mbcalfile = NULL; / macbethcal mapping file */
21	greg	3.8	char cwarpfile = NULL; / color space warping file */
22	greg	3.1
23	greg	3.3	static struct mbc {
24	greg	3.7	COLORMAT cmat;
25	greg	3.3	float xa[3][6], ya[3][6];
26	greg	3.7	COLOR cmin, cmax;
27	greg	3.3	} mbcond; /* macbethcal conditioning struct */
28	greg	3.1
29	greg	3.8	static WARP3D cwarp; / color warping structure */
30
31	greg	3.1	static COLOR scanbuf; / scanline processing buffer */
32			static int nread; /* number of scanlines processed */
33
34	schorsch	3.12	static void sfscan(COLOR *sl, int len, double sf);
35			static void matscan(COLOR *sl, int len, COLORMAT mat);
36			static void mbscan(COLOR sl, int len, struct mbc mb);
37			static void cwscan(COLOR sl, int len, WARP3D wp);
38			static void getmbcalfile(char fn, struct mbc mb);
39	greg	3.1
40	schorsch	3.12
41			extern double
42			rgblum( /* compute (scotopic) luminance of RGB color */
43			COLOR clr,
44			int scotopic
45			)
46	greg	3.1	{
47			if (scotopic) /* approximate */
48			return( WHTSEFFICACY * (colval(clr,RED)*.062 +
49			colval(clr,GRN).608 + colval(clr,BLU).330) );
50			return( WHTEFFICACY * (colval(clr,RED)*inrgb2xyz[1][0] +
51			colval(clr,GRN)*inrgb2xyz[1][1] +
52			colval(clr,BLU)*inrgb2xyz[1][2]) );
53			}
54
55
56	schorsch	3.12	extern double
57			cielum( /* compute (scotopic) luminance of CIE color */
58			COLOR xyz,
59			int scotopic
60			)
61	greg	3.1	{
62			if (scotopic) /* approximate */
63			return(colval(xyz,CIEY) *
64			(1.33*(1. + (colval(xyz,CIEY)+colval(xyz,CIEZ))/
65			colval(xyz,CIEX)) - 1.68));
66			return(colval(xyz,CIEY));
67			}
68
69
70	schorsch	3.12	extern COLOR *
71			nextscan(void) /* read and condition next scanline */
72	greg	3.1	{
73			if (nread >= numscans(&inpres)) {
74	greg	3.8	if (cwarpfile != NULL)
75			free3dw(cwarp);
76	greg	3.10	free((void *)scanbuf);
77	greg	3.4	return(scanbuf = NULL);
78	greg	3.1	}
79	greg	3.2	if (what2do&DO_ACUITY)
80			acuscan(scanbuf, nread);
81			else if (freadscan(scanbuf, scanlen(&inpres), infp) < 0) {
82	greg	3.1	fprintf(stderr, "%s: %s: scanline read error\n",
83			progname, infn);
84			exit(1);
85			}
86			if (what2do&DO_VEIL) /* add veiling */
87	greg	3.2	addveil(scanbuf, nread);
88	greg	3.1	if (what2do&DO_COLOR) /* scotopic color loss */
89			scotscan(scanbuf, scanlen(&inpres));
90			if (what2do&DO_LINEAR) /* map luminances */
91			sfscan(scanbuf, scanlen(&inpres), scalef);
92			else
93			mapscan(scanbuf, scanlen(&inpres));
94			if (mbcalfile != NULL) /* device color correction */
95			mbscan(scanbuf, scanlen(&inpres), &mbcond);
96	greg	3.8	else if (cwarpfile != NULL) /* device color space warp */
97			cwscan(scanbuf, scanlen(&inpres), cwarp);
98	schorsch	3.11	else if ((lumf == cielum) \| (inprims != outprims))
99	greg	3.1	matscan(scanbuf, scanlen(&inpres), mbcond.cmat);
100	greg	3.2	nread++;
101	greg	3.1	return(scanbuf);
102			}
103
104
105	schorsch	3.12	extern COLOR *
106			firstscan(void) /* return first processed scanline */
107	greg	3.1	{
108			if (mbcalfile != NULL) /* load macbethcal file */
109			getmbcalfile(mbcalfile, &mbcond);
110	greg	3.8	else if (cwarpfile != NULL) {
111			if ((cwarp = load3dw(cwarpfile, NULL)) == NULL)
112			syserror(cwarpfile);
113			} else
114	greg	3.1	if (lumf == rgblum)
115	greg	3.10	comprgb2rgbWBmat(mbcond.cmat, inprims, outprims);
116	greg	3.1	else
117	greg	3.10	compxyz2rgbWBmat(mbcond.cmat, outprims);
118	greg	3.14	if (what2do&DO_ACUITY)
119			initacuity();
120	greg	3.1	scanbuf = (COLOR )malloc(scanlen(&inpres)sizeof(COLOR));
121			if (scanbuf == NULL)
122			syserror("malloc");
123			nread = 0;
124			return(nextscan());
125			}
126
127
128	schorsch	3.12	static void
129			sfscan( /* apply scalefactor to scanline */
130			register COLOR *sl,
131			int len,
132			double sf
133			)
134	greg	3.1	{
135			while (len--) {
136			scalecolor(sl[0], sf);
137			sl++;
138			}
139			}
140
141
142	greg	3.15	static double
143			greypoint( /* compute gamut mapping grey target */
144			COLOR col
145			)
146			{
147			COLOR gryc;
148			int i;
149			/* improves saturated color rendering */
150			copycolor(gryc, col);
151			for (i = 3; i--; )
152			if (gryc[i] > cwhite[i])
153			gryc[i] = cwhite[i];
154			else if (gryc[i] < cblack[i])
155			gryc[i] = cblack[i];
156	greg	3.16	return((*lumf)(gryc,0));
157	greg	3.15	}
158
159
160	schorsch	3.12	static void
161			matscan( /* apply color matrix to scaline */
162			register COLOR *sl,
163			int len,
164			COLORMAT mat
165			)
166	greg	3.1	{
167	greg	3.15	double gryv;
168
169	greg	3.1	while (len--) {
170	greg	3.15	gryv = greypoint(sl[0]);
171	greg	3.1	colortrans(sl[0], mat, sl[0]);
172	greg	3.15	clipgamut(sl[0], gryv, CGAMUT, cblack, cwhite);
173	greg	3.1	sl++;
174			}
175			}
176
177
178	schorsch	3.12	static void
179			mbscan( /* apply macbethcal adj. to scaline */
180			COLOR *sl,
181			int len,
182			register struct mbc *mb
183			)
184	greg	3.1	{
185			double d;
186			register int i, j;
187
188			while (len--) {
189	greg	3.15	d = greypoint(sl[0]);
190	greg	3.6	colortrans(sl[0], mb->cmat, sl[0]);
191	greg	3.15	clipgamut(sl[0], d, CGAMUT, mb->cmin, mb->cmax);
192	greg	3.1	for (i = 0; i < 3; i++) {
193			d = colval(sl[0],i);
194			for (j = 0; j < 4 && mb->xa[i][j+1] <= d; j++)
195			;
196			colval(sl[0],i) = ( (mb->xa[i][j+1] - d)*mb->ya[i][j] +
197			(d - mb->xa[i][j])*mb->ya[i][j+1] ) /
198			(mb->xa[i][j+1] - mb->xa[i][j]);
199			}
200	greg	3.8	sl++;
201			}
202			}
203
204
205	schorsch	3.12	static void
206			cwscan( /* apply color space warp to scaline */
207			COLOR *sl,
208			int len,
209			WARP3D *wp
210			)
211	greg	3.8	{
212			int rval;
213	greg	3.15	double gryv;
214	greg	3.8
215			while (len--) {
216	greg	3.15	gryv = greypoint(sl[0]);
217	greg	3.8	rval = warp3d(sl[0], sl[0], wp);
218			if (rval & W3ERROR)
219			syserror("warp3d");
220			if (rval & W3BADMAP) {
221			fprintf(stderr, "%s: %s: bad color space map\n",
222			progname, cwarpfile);
223			exit(1);
224			}
225	greg	3.15	clipgamut(sl[0], gryv, CGAMUT, cblack, cwhite);
226	greg	3.1	sl++;
227			}
228			}
229
230
231	schorsch	3.12	static void
232			getmbcalfile( /* load macbethcal file */
233			char *fn,
234			register struct mbc *mb
235			)
236	greg	3.1	{
237			char buf[128];
238			FILE *fp;
239			int inpflags = 0;
240	greg	3.7	register int i;
241	greg	3.1
242			if ((fp = fopen(fn, "r")) == NULL)
243			syserror(fn);
244			while (fgets(buf, sizeof(buf), fp) != NULL) {
245			if (!(inpflags & 01) &&
246			sscanf(buf,
247			"rxa(i) : select(i,%f,%f,%f,%f,%f,%f)",
248			&mb->xa[0][0], &mb->xa[0][1],
249			&mb->xa[0][2], &mb->xa[0][3],
250			&mb->xa[0][4], &mb->xa[0][5]
251			) == 6)
252			inpflags \|= 01;
253			else if (!(inpflags & 02) &&
254			sscanf(buf,
255			"rya(i) : select(i,%f,%f,%f,%f,%f,%f)",
256			&mb->ya[0][0], &mb->ya[0][1],
257			&mb->ya[0][2], &mb->ya[0][3],
258			&mb->ya[0][4], &mb->ya[0][5]
259			) == 6)
260			inpflags \|= 02;
261			else if (!(inpflags & 04) &&
262			sscanf(buf,
263			"gxa(i) : select(i,%f,%f,%f,%f,%f,%f)",
264			&mb->xa[1][0], &mb->xa[1][1],
265			&mb->xa[1][2], &mb->xa[1][3],
266			&mb->xa[1][4], &mb->xa[1][5]
267			) == 6)
268			inpflags \|= 04;
269			else if (!(inpflags & 010) &&
270			sscanf(buf,
271			"gya(i) : select(i,%f,%f,%f,%f,%f,%f)",
272			&mb->ya[1][0], &mb->ya[1][1],
273			&mb->ya[1][2], &mb->ya[1][3],
274			&mb->ya[1][4], &mb->ya[1][5]
275			) == 6)
276			inpflags \|= 010;
277			else if (!(inpflags & 020) &&
278			sscanf(buf,
279			"bxa(i) : select(i,%f,%f,%f,%f,%f,%f)",
280			&mb->xa[2][0], &mb->xa[2][1],
281			&mb->xa[2][2], &mb->xa[2][3],
282			&mb->xa[2][4], &mb->xa[2][5]
283			) == 6)
284			inpflags \|= 020;
285			else if (!(inpflags & 040) &&
286			sscanf(buf,
287			"bya(i) : select(i,%f,%f,%f,%f,%f,%f)",
288			&mb->ya[2][0], &mb->ya[2][1],
289			&mb->ya[2][2], &mb->ya[2][3],
290			&mb->ya[2][4], &mb->ya[2][5]
291			) == 6)
292			inpflags \|= 040;
293			else if (!(inpflags & 0100) &&
294			sscanf(buf,
295	gwlarson	3.9	"ro = %frn + %fgn + %f*bn",
296	greg	3.1	&mb->cmat[0][0], &mb->cmat[0][1],
297			&mb->cmat[0][2]) == 3)
298			inpflags \|= 0100;
299			else if (!(inpflags & 0200) &&
300			sscanf(buf,
301	gwlarson	3.9	"go = %frn + %fgn + %f*bn",
302	greg	3.1	&mb->cmat[1][0], &mb->cmat[1][1],
303			&mb->cmat[1][2]) == 3)
304			inpflags \|= 0200;
305			else if (!(inpflags & 0400) &&
306			sscanf(buf,
307	gwlarson	3.9	"bo = %frn + %fgn + %f*bn",
308	greg	3.1	&mb->cmat[2][0], &mb->cmat[2][1],
309			&mb->cmat[2][2]) == 3)
310			inpflags \|= 0400;
311			}
312			if (inpflags != 0777) {
313			fprintf(stderr,
314			"%s: cannot grok macbethcal file \"%s\" (inpflags==0%o)\n",
315			progname, fn, inpflags);
316			exit(1);
317			}
318			fclose(fp);
319	greg	3.7	/* compute gamut */
320			for (i = 0; i < 3; i++) {
321			colval(mb->cmin,i) = mb->xa[i][0] -
322			mb->ya[i][0] *
323			(mb->xa[i][1]-mb->xa[i][0]) /
324			(mb->ya[i][1]-mb->ya[i][0]);
325			colval(mb->cmax,i) = mb->xa[i][4] +
326			(1.-mb->ya[i][4]) *
327			(mb->xa[i][5] - mb->xa[i][4]) /
328			(mb->ya[i][5] - mb->ya[i][4]);
329			}
330	greg	3.1	}