README 660 401 24 4154 4650403730 4763 ------------------------------------------------------------------------------ "Illumination and Color in Computer Generated Imagery" Source for code given in the Appendices (c)Roy Hall - 1988,1989 ------------------------------------------------------------------------------ This code was developed to test and demonstrate concepts discussed in "Illumination and Color for Computer Generated Imagery" (1989 from Springer-Verlag). While the code was used in the generation of the diagrams and images presented in the text, it may still contain bugs. The user is cautioned to independently verify that the code is suitable to his/her needs. This code is distributed on an "as is" basis with all faults and without any implied or expressed warranties or support from either the author, Roy Hall, or from Springer-Verlag Publishers. I am making this code available to the user community. You are free to use it as you wish. In return I request two things. First, that this README file always stays with any copy of the source, more specifically, that proper credit is given to the original author. Second, that any corrections, improvements, etc. are passed back to the author for inclusion in subsequent distribution and editions of the book. Please note that the code was written for clarity of presentation, not for performance. You will find that optimizations can be easily made. Please do not report these as I will not incorporate them into future releases. Note also that this code relects current work I am involved in and may differ slightly from that presented in the text. There have been some additions to the distribution (16 July 90): Corrections: list of corrections for the book and the code Slectral/: a directory containing a few spectral curves to help get you started Roy Hall Program of Computer Graphics 120 Rand Hall Cornell University Ithaca, NY 14853 roy@wisdom.graphics.cornell.edu Contents: README Corrections geo.h geo.c F.h F.c D.h D.c G.h G.c clr.h clr.c clr_clip.c clr_sample.c illum_mod.h illum_mod.c Spectral/ ome additions to the distribution (16 July 90): Corrections: list of corrections for the book and the code Slectral/: a directory containing a few spectral curves to help get you started Roy Hall Program of Computer Graphics 120 Rand Hall Cornell University Ithaca, NY 14853 roy@wisdom.graphics.cornell.edu Contents: README Corrections geo.h geo.c F.h F.c Corrections 660 401 24 11565 4656154650 6356 ----------------------------------------------------------------- Corrections since the first edition of the book ----------------------------------------------------------------- Confusion about emissivity -- this usually refers to energy, I mean intensity. p.40 - last line in 'Emissive Illumination', 'emissivity' should read 'emissive intensity' p.40 - 2nd line in 'A Generalized Illumination Expression', 'emissivity' should read 'emissive intensity' p.40 - footnote 18 should read: Emissive intensity refers to the intensity leaving the surface in a given direction at any wavelength. However, in computer graphics we are typically concerned only with the visible spectrum. p.95 - Eq (4.22) - function of lambda is missing for the ambient illumination in the ambient term p.99 - 2nd line in 'Radiosity for Diffuse Environments', 'emissivity' should read 'emissive intensity' p.154 - 'emissivity' should read 'emissive intensity', revised some dependencies - emissive intensity - add lambda,theta, phi intensity - add lambda coherent reflectance - add sigma coherent trans - add sigma attanuation - add lambda fix spelling on coherent reflection p.275 - Add to the Sancer (1969) reference the journal name: IEEE Transactions on Antennas and Propagation. p.280 - 'emissivity' should read 'emissive intensity' ---- p.37 - equation 2.46 -- denomenator should be cos(theta sub i)cos(theta sub r) p.46 - end of second paragraph - 'radial position' should read 'circumferential position' p.116 - 1st para, Rogers ref should be Rogers 1985 p.120 - 2nd para, chapter reference 3.3, Color Spaces ... should be Section 3.5, Color Spaces ... p.125 - sect 5.1.3, 2nd par, 1st line: add 'and' after 'primaries of the monitor' p.125 - sect 5.1.3, 3rd par (after bullets), 1st line: 'results negative values' should be 'results in negative values' p.126 - 6th para, 5th line: 'maintains hue and saturation and modifies the chromaticity' should read 'maintains hue and saturation and modifies the intensity' p.127 - 3rd par, 3rd line: 'image bands of represent' should read 'image bands of color represent' p.133 - There is an error in the RGB to YIQ matrix. it should be: 0.299 0.587 0.114 0.596 -0.275 -0.321 0.212 -0.523 0.311 - footnote 15 -- missing '(' before the R in Y(R-Y)(B-Y) p.249 - 5th para: 'plane that passes through the block point' to 'plane that passes through the black point' problems in the code: p.174 - change: * geo_getv - generate a vector to: * geo_line - generate a line p.209 - 2nd line "ignored if a NULL" instead of "ignored is a NULL" p.210 - remove 'surface has' in the first line of Notes. p.227 - change 'CLR_xyz_lab' to 'CLR_xyz_to_lab' change 'CLR_xyz_luv' to 'CLR_xyz_to_luv' p.228 - change rgb_to_yiq matrix (lines 86-89 of clr.c) to: static double rgb_yiq_mat[3][3] = { {0.299, 0.587, 0.114}, {0.596, -0.275, -0.321}, {0.212, -0.523, 0.311}}; p.288 - change the data for 480nm entry in CIEXYZ (line 121 of clr.c) to: {480, 0.0956, 0.1390, 0.8130}, {485, 0.0580, 0.1693, 0.6162}, p.230 - note after XYZ_done label: 'Y os a sampled' should be 'Y of a sampled' p.232 - a routine 'CLR_blackbody' has been inserted after 'CLR_read_mtl' and before 'CLR_add_spect' p.243 - The samp_to_ACC matrix is incorrect, should be: static double samp_to_ACC[3][4] = {{0.00000, 0.18892, 0.67493, 0.19253}, {0.00000, 0.31824, 0.00000, -0.46008}, {0.54640, 0.00000, 0.00000, 0.00000}}; p.244 - first line in the CLR_init_samples routine should be: CLR_exit_samples (); immediately before the: if(method == CLR_SAMPLE_MEYER) { p.248 - in notes (couple lines down): 'resulting a a step' should read 'resulting in a step' Added a CLR_xyz_to_spec routine that uses a method described by Glassner (How to Derive a Spectrum from an RGB Triplet, IEEE CG&A, July 1989). problems in figures: p.59 - annotate the start and end of the functions p.81 - figure 4.11 the image is reversed (left to right) - when corrected, the left image should be labelled 'Gouraud shading' and the right image should be labelled 'Phong shading' p.132 - figure 5.12 should be a black and white print p.164 - figure II.5 row1, column 3 Cka and Ckd should be 0.35 p.171 - figure III.2 - the reversed normal caption is wrong -- should be: N'' = -N sqrt(1 - |S'|^2) get rid of the divide by N.V end of the functions p.81 - figure 4.11 the image is reversed (left to right) - when corrected, the left image should be labellgeo.h 660 401 24 2200 4623574116 5023 /* **************************************************************** * geo.h * **************************************************************** * include file for the geometric utilities */ #ifndef GEO_H #define GEO_H #ifndef TRUE #define TRUE 1 #endif TRUE #ifndef FALSE #define FALSE 0 #endif FALSE /* common geometric constructs */ typedef struct {double i, j, k;} DIR_VECT; typedef struct {double x, y, z;} POINT3; typedef struct {double x, y, z, h;} POINT4; typedef struct {POINT3 start; DIR_VECT dir; } LINE; /* geometric manipulation routines */ double geo_dot(); /* vector dot product */ DIR_VECT geo_cross(); /* vector cross product */ double geo_norm(); /* vector normalize */ double geo_line(); /* vector between two points */ DIR_VECT *geo_rfl(); /* reflected vector */ DIR_VECT *geo_rfr(); /* refracted vector */ DIR_VECT *geo_H(); /* H vector */ DIR_VECT *geo_Ht(); /* H' vector */ #endif GEO_H /* ************************************************************* */ geo_dot(); /* vector dot product */ DIR_VECT geo_cross(); /* vector cross product */ double geo_norm(); /* vector normalize */ double geo_line(); /* vector between two points */ DIR_VECT *geo_rfl(); /* reflected vector */ DIR_VECT *geo_rfr(); /* refracted vector */ DIR_VECT *geo_H(); /* H vector */ DIgeo.c 660 401 24 16406 4623574116 5053 /* **************************************************************** * geo.c * **************************************************************** * MODULE PURPOSE: * This module contains routines that preform common * geometric manipulations for image synthesis. * * MODULE CONTENTS: * geo_dot - vector dot product * geo_cross - vector cross product * geo_norm - vector normalization * geo_line - generate a line * geo_rfl - compute the reflected vector * geo_rfr - compute the refracted vector * geo_H - compute H vector * geo_Ht - compute the H' vector */ #include #include #include "geo.h" /* **************************************************************** * double geo_dot (v0, v1) * DIR_VECT *v0, *v1 (in) normalized direction vectors * * Returns the dot product of two direction vectors, ref * Eq.(\*(rK). Note that the dot product is the cosine * between the direction vectors because the direction * vectors are normalized. */ double geo_dot (v0, v1) DIR_VECT *v0, *v1; { return (v0->i * v1->i) + (v0->j * v1->j) + (v0->k * v1->k); } /* **************************************************************** * DIR_VECT geo_cross (v0, v1) * DIR_VECT *v0, *v1 (in) normalized direction vectors * * Returns the cross product of two direction vectors, refer * to Eq.(\*(rL). */ DIR_VECT geo_cross (v0, v1) DIR_VECT *v0, *v1; { DIR_VECT v2; v2.i = (v0->j * v1->k) - (v0->k * v1->j); v2.j = (v0->k * v1->i) - (v0->i * v1->k); v2.k = (v0->i * v1->j) - (v0->j * v1->i); return v2; } /* **************************************************************** * double geo_norm (*v0) * DIR_VECT *v0 (mod) vector to be normalized * * Returns the length of the vector before normalization. * Returns zero if the vector could not be normalized. * Note that the input vector is modified. */ double geo_norm (v0) DIR_VECT *v0; { double len; if ((len = geo_dot(v0, v0)) <= 0.0) return FALSE; len = sqrt((double)len); v0->i /= len; v0->j /= len; v0->k /= len; return TRUE; } /* **************************************************************** * double geo_line (*p0, *p1, *v) * POINT3 *p0, *p1 (in) from, to points * LINE *v (out) generated line * * Returns the distance from p0 to p1. Returns 0.0 * on error (p0=p1). The line is expressed in parametric * form with a start point (p0) and a normalized direction * vector. The vector direction of the line is normalized. */ double geo_line (p0, p1, v) POINT3 *p0, *p1; LINE *v; { v->start = *p0; v->dir.i = p1->x - p0->x; v->dir.j = p1->y - p0->y; v->dir.k = p1->z - p0->z; return geo_norm(&v->dir); } /* **************************************************************** * DIR_VECT *geo_rfl (V, N) * DIR_VECT *V (in) incident vector * DIR_VECT *N (in) surface normal * * Returns the reflected direction vector. The reflected * direction is computed using the method given by Whitted * (1980), refer to Eq.(\*(rM). */ DIR_VECT *geo_rfl (V, N) DIR_VECT *V; DIR_VECT *N; { double N_dot_V; static DIR_VECT rfl; N_dot_V = geo_dot (N,V); rfl.i = (2.0 * N_dot_V * N->i) - V->i; rfl.j = (2.0 * N_dot_V * N->j) - V->j; rfl.k = (2.0 * N_dot_V * N->k) - V->k; return &rfl; } /* **************************************************************** * DIR_VECT *geo_rfr (V, N, ni, nt) * DIR_VECT *V (in) incident vector * DIR_VECT *N (in) surface normal * double ni (in) index of refraction for the * material on the front of the * interface (same side as N) * double nt (in) index of refraction for the * material on the back of the * interface (opposite size as N) * * Returns the refracted vector, if there complete internal * refracted (no refracted vector) then a NULL vector is * NULL is returned. The vector is computed using the * method given by Hall (1983) with enhancements as * described in Appendix III. */ DIR_VECT *geo_rfr (V, N, ni, nt) DIR_VECT *V; DIR_VECT *N; double ni; double nt; { static DIR_VECT T; /* the refracted vector */ DIR_VECT sin_T; /* sin vect of the refracted vect */ DIR_VECT cos_V; /* cos vect of the incident vect */ double len_sin_T; /* length of sin T squared */ double n_mult; /* ni over nt */ double N_dot_V; double N_dot_T; if ((N_dot_V = geo_dot(N,V)) > 0.0) n_mult = ni / nt; else n_mult = nt / ni; cos_V.i = N_dot_V * N->i; cos_V.j = N_dot_V * N->j; cos_V.k = N_dot_V * N->k; sin_T.i = (n_mult) * (cos_V.i - V->i); sin_T.j = (n_mult) * (cos_V.j - V->j); sin_T.k = (n_mult) * (cos_V.k - V->k); if ((len_sin_T = geo_dot(&sin_T, &sin_T)) >= 1.0) return NULL; /* internal reflection */ N_dot_T = sqrt(1.0 - len_sin_T); if (N_dot_V < 0.0) N_dot_T = -N_dot_T; T.i = sin_T.i - (N->i * N_dot_T); T.j = sin_T.j - (N->j * N_dot_T); T.k = sin_T.k - (N->k * N_dot_T); return &T; } /* **************************************************************** * DIR_VECT *geo_H (L, V) * DIR_VECT *L (in) incident vector * DIR_VECT *V (in) reflected vector * * Returns H, NULL on error (if L+H = 0). Refer * to Eq.(\*(rO). */ DIR_VECT *geo_H (L, V) DIR_VECT *L; DIR_VECT *V; { static DIR_VECT H; H.i = L->i + V->i; H.j = L->j + V->j; H.k = L->k + V->k; if (!geo_norm(&H)) return NULL; return &H; } /* **************************************************************** * DIR_VECT *geo_Ht(L, T, ni, nt) * DIR_VECT *L (in) incident vector * DIR_VECT *T (in) transmitted hector * double ni (in) incident index * double nt (in) transmitted index * * Returns H' oriented to the same side of the surface * as L computed using the method suggested by * Hall (1983). Returns NULL on error (if the angle * between V and L is less than the critical angle). */ DIR_VECT *geo_Ht (L, T, ni, nt) DIR_VECT *L; DIR_VECT *T; double ni; double nt; { float L_dot_T; float divisor; static DIR_VECT Ht; L_dot_T = -(geo_dot(L,T)); /* check for special cases */ if (ni == nt) { if (L_dot_T == 1.0) return L; else return NULL; } if (ni < nt) { if (L_dot_T < ni/nt) return NULL; divisor = (nt / ni) - 1.0; Ht.i = -(((L->i + T->i) / divisor) + T->i); Ht.j = -(((L->j + T->j) / divisor) + T->j); Ht.k = -(((L->k + T->k) / divisor) + T->k); } else { if (L_dot_T < nt/ni) return NULL; divisor = (ni / nt) - 1.0; Ht.i = ((L->i + T->i) / divisor) + L->i; Ht.j = ((L->j + T->j) / divisor) + L->j; Ht.k = ((L->k + T->k) / divisor) + L->k; } (void)geo_norm(&Ht); return &Ht; } /* ************************************************************* */ ni/nt) return NULL; divisor = (nt / ni) - 1.0; Ht.i = -(((L->i + T->i) / divisor) + T->i); Ht.j = -(((L->j + T->j) / divisor) + T->j); Ht.k = -(((L->k + T->k) / divisor) + T->k); } else { if (L_dot_T < nt/ni) return NULL; divisor F.h 660 401 24 1032 4623574120 4433 /* **************************************************************** * F.h * **************************************************************** * include file for the microfacet distribution functions */ #ifndef F_H #define F_H double F_d_pl(); double F_d_pp(); double F_d_R(); double F_c_pl(); double F_c_pp(); double F_c_R(); double F_approx_n(); double F_approx_k(); double *F_approx_Fr(); double *F_approx_Fr_Ft(); #endif F_H /* ************************************************************* */ ***************************************** * F.h * **************************************************************** * include file for the microfacet distribution functions */ #ifndef F_H #define F_H double F_d_pl(); double F_d_pp(); double F_d_R(); double F_c_pl(); double F_c_pp(); double F_c_R(); double F_approx_n(); double F_approx_k(); double *F_approx_Fr(); double *F_approx_Fr_Ft(); #endif F_H /* ***************************************F.c 660 401 24 23574 4623574116 4472 /* **************************************************************** * F.c * **************************************************************** * MODULE PURPOSE: * This module contains routines to for fresnel calculations. * * MODULE CONTENTS: * F_d_pl - reflectance for a dielectric, * parallel polarized, * F_d_pp - reflectance for a dielectric, * perpendicular polarized * F_d_R - reflectance for a dielectric * F_c_pl - reflectance for a conductor, * parallel polarized, * F_c_pp - reflectance for a conductor, * perpendicular polarized * F_c_R - reflectance for a conductor * F_approx_n - approximate the average index of * refraction for a material * F_approx_k - approximate the average absoption * index for a material * F_approx_Fr - approximate conductor reflectance * F_approx_Fr_Ft - approximate dielectric reflectance * and transmittance * * * ASSUMPTIONS: * The following are defined in math.h * HUGE largest floating pont value * M_PI_2 pi/2 * M_PI_4 pi/4 */ #include #include #include "geo.h" #include "F.h" /* **************************************************************** * double F_d_pl (N, L, T, ni, nt) * DIR_VECT *N, *L, *T (in) - N, L, T vectors * double ni (in) - incident index of refraction * double nt (in) - transmitted index of refraction * * Returns the reflectance of a dielectric for light with * polarization parallel to the plane of incidence (plane of the * N and L vector), refer to Eq.(\*(eJ). N and L are assumed to * be to the same side of the surface. */ double F_d_pl (N, L, T, ni, nt) DIR_VECT *N, *L, *T; double ni, nt; { double amplitude, N_dot_L, N_dot_T; N_dot_L = geo_dot(N, L); N_dot_T = geo_dot(N, T); amplitude = ((nt * N_dot_L) + (ni * N_dot_T)) / ((nt * N_dot_L) - (ni * N_dot_T)); return amplitude * amplitude; } /* **************************************************************** * double F_d_pp (N, L, T, ni, nt) * DIR_VECT *N, *L, *T (in) - N, L, T vectors * double ni (in) - incident index of refraction * double nt (in) - transmitted index of refraction * * Returns the reflectance of a dielectric for light with * polarization perpendicular to the plane of incidence * (plane of the N and L vector), refer to Eq.(\*(eK). N and L * are assumed to be to the same side of the surface. */ double F_d_pp (N, L, T, ni, nt) DIR_VECT *N, *L, *T; double ni, nt; { double amplitude, N_dot_L, N_dot_T; N_dot_L = geo_dot(N, L); N_dot_T = geo_dot(N, T); amplitude = ((ni * N_dot_L) + (nt * N_dot_T)) / ((ni * N_dot_L) - (nt * N_dot_T)); return amplitude * amplitude; } /* **************************************************************** * double F_d_R (N, L, T, ni, nt) * DIR_VECT *N, *L, *T (in) - N, L, T vectors * double ni (in) incident index of refraction * double nt (in) transmitted index of refraction * * Returns the average reflectance for a dielectric. Refer to * Eq.(\*(eL). N and L are assumed to be to the same side of the * surface. */ double F_d_R (N, L, T, ni, nt) DIR_VECT *N, *L, *T; double ni, nt; { return (F_d_pl (N, L, T, ni, nt) + F_d_pp (N, L, T, ni, nt)) / 2.0; } /* **************************************************************** * double F_c_pl (N, L, nt, k) * DIR_VECT *N, *L (in) - N, L vectors * double nt (in) - index of refraction * double k (in) - absorption coefficient * * Returns the reflectance of a conductor for light with * polarization parallel to the plane of incidence (plane of the * N and L vector), refer to Eq.(\*(eO). N and L are assumed to * be to the same side of the surface. */ double F_c_pl (N, L, nt, k) DIR_VECT *N, *L; double nt, k; { double tmp_f, N_dot_L; N_dot_L = geo_dot(N, L); tmp_f = ((nt * nt) + (k * k)) * N_dot_L * N_dot_L; return (tmp_f - (2.0 * nt * N_dot_L) + 1) / (tmp_f + (2.0 * nt * N_dot_L) + 1); } /* **************************************************************** * double F_c_pp (N, L, nt, k) * DIR_VECT *N, *L (in) - N, L vectors * double nt (in) - index of refraction * double k (in) - absorption coefficient * * Returns the reflectance of a conductor for light with * polarization perpendicular to the plane of incidence (plane of * the N and L vector), refer to Eq.(\*(eP). N and L are assumed * to be to the same side of the surface. */ double F_c_pp (N, L, nt, k) DIR_VECT *N, *L; double nt, k; { double tmp_f, N_dot_L; N_dot_L = geo_dot(N, L); tmp_f = (nt * nt) + (k * k); return (tmp_f - (2.0 * nt * N_dot_L) + (N_dot_L * N_dot_L)) / (tmp_f + (2.0 * nt * N_dot_L) + (N_dot_L * N_dot_L)); } /* **************************************************************** * double F_c_R (N, L, nt, k) * DIR_VECT *N, *L (in) - N, L vectors * double nt (in) - index of refraction * double k (in) - absorption coefficient * * Returns the average reflectance for a conductor. Refer to * Eq.(\*(eL). N and L are assumed to be to the same side of * the surface. */ double F_c_R (N, L, nt, k) DIR_VECT *N, *L; double nt, k; { return (F_c_pl (N, L, nt, k) + F_c_pp (N, L, nt, k))/2.0; } /* **************************************************************** * double F_approx_n (Ro, n_samp, mtl) * double *Ro (out) - average reflectance * int n_samp (in) - number of material samples * double *mtl (in) - material spectral curve * * Returns the approximate n, and loads Ro. This gives the * correct n for a single mtl reflectance if the material is a * dielectric, and the correct n assuming k=0 for a conductor, * refer to Eq.(\*(eN). */ double F_approx_n (Ro, n_samp, mtl) double *Ro; int n_samp; double *mtl; { int ct; *Ro = 0.0; for (ct=n_samp; --ct>=0; ) *Ro += *mtl++; *Ro /= (double)n_samp; return (1.0 + sqrt(*Ro)) / (1.0 - sqrt(*Ro)); } /* **************************************************************** * double F_approx_k (Ro, n_samp, mtl) * double *Ro (out) - average reflectance * int n_samp (in) - number of material samples * double *mtl (in) - material spectral curve * * Returns the approximate k, and loads Ro. This assumes n=1.0 * and is applicable to conductors only. Refer to Eq.(\*(eQ). */ double F_approx_k (Ro, n_samp, mtl) double *Ro; int n_samp; double *mtl; { int ct; *Ro = 0.0; for (ct=n_samp; --ct>=0; ) *Ro += *mtl++; *Ro /= (double)n_samp; return 2.0 * sqrt(*Ro / (1.0 - *Ro)); } /* **************************************************************** * double *F_approx_Fr (N, L, Ro, n, k, n_samp, mtl, Fr) * DIR_VECT *N, *L (in) - N, L vectors * double Ro (in) - average reflectance * double n (in) - measured or approximate index * of refraction * double k (in) - measured or approximate * absorption coefficient * int n_samp (in) number of material samples * double *mtl (in) material spectral curve * double *Fr (out) reflectance * * Loads the vector Fr with the approximate reflectance for each * color sample point for a conductor, refer to Eq.(\*(eR). * Returns the pointer to Fr. The measured values for n and k * should be used if they are available. Otherwise, assume k=0 * and approximate n using F_approx_n(), or assume n=1 and * approximate k using F_approx_k(). */ double *F_approx_Fr (N, L, Ro, n, k, n_samp, mtl, Fr) DIR_VECT *N, *L; double Ro, n, k; int n_samp; double *mtl, *Fr; { double R_theta; /* ave R for N and L*/ double factor, *R_out; R_theta = F_c_R (N, L, n, k); factor = (R_theta - Ro) / (1.0 - Ro); for (R_out=Fr ;--n_samp>=0; Fr++, mtl++) if ((*Fr = *mtl + ((1.0 - *mtl) * factor)) < 0.0) *Fr = 0.0; return R_out; } /* **************************************************************** * double *F_approx_Fr_Ft (N, L, T, Ro, ni, nt, * n_samp, *mtl_p, Fr, Ft) * DIR_VECT *N, *L, *T (in) - N, L vectors * double Ro (in) average reflectance * double ni (in) ave n for incident material * double nt (in) ave n for transmitted material * int n_samp (in) number of material samples * double *mtl_p (in) pseudo material curve * double *Fr (out) reflectance * double *Tr (out) transmittance * * Loads the vector Fr with the approximate reflectance and Ft * with the approximate transmittance for each sample point for * a conductor, refer to Eq.(\*(eR). Returns the pointer to Fr. * The measured value for n should be used if it is available, * otherwise use F_approx_n() to approximate the index of * refraction. */ double *F_approx_Fr_Ft (N, L, T, Ro, ni, nt, n_samp, mtl_p, Fr, Ft) DIR_VECT *N, *L, *T; double Ro, ni, nt; int n_samp; double *mtl_p, *Fr, *Ft; { double R_theta; /* ave R for N and L*/ double factor, *R_out; if (T == NULL) /* internal reflection */ for (R_out=Fr; --n_samp>=0; *Fr++ = 1.0, *Ft++ = 0.0) ; else { R_theta = F_d_R (N, L, T, ni, nt); factor = (R_theta - Ro) / (1.0 - Ro); for (R_out=Fr; --n_samp>=0; mtl_p++) { if ((*Fr = *mtl_p + ((1.0 - *mtl_p) * factor)) < 0.0) *Fr = 0.0; *Ft++ = 1.0 - *Fr++; } } return R_out; } /* ************************************************************* */ { double R_theta; /* ave R for N and L*/ double factor, *R_out; if (T == NULL) /* internal reflection *D.h 660 401 24 1304 4623574117 4441 /* **************************************************************** * D.h * **************************************************************** * include file for the microfacet distribution functions */ #ifndef D_H #define D_H /* Microfacet distribution routines */ double D_phong_init(); double D_phong(); double D_blinn_init(); double D_blinn(); double D_gaussian_init(); double D_gaussian(); double D_reitz_init(); double D_reitz(); double D_cook_init(); double D_cook(); double D_g(); double D_tau(); double D_sigma(); double D_m(); double D_coherent(); double D_incoherent(); double D_Vxz(); #endif D_H /* ************************************************************* */ crofacet distribution functions */ #ifndef D_H #define D_H /* Microfacet distribution routines */ double D_phong_init(); double D_phong(); double D_blinn_init(); double D_blinn(); double D_gaussian_init(); double D_gaussian(); double D_reitz_init(); double D_reitz(); double D_cook_init(); double D_cook(); doubleD.c 660 401 24 33764 4623574117 4473 /* **************************************************************** * D.c * **************************************************************** * MODULE PURPOSE: * This module contains microfacet distribution routines and * various support routines. One of support functions includes * computing the coefficients corresponding to beta for each * function (as described in Section 4.3.1-Improved Specular * Reflection). Beta is the angle between H and N where the * function is equal to half the value at N = H. Beta is in * radians. * * MODULE CONTENTS: * D_phong_init - compute Ns given beta * D_phong - Phong cosine power function * D_blinn_init - compute Ns given beta * D_blinn - Blinn cosine power function * D_gaussian_init - compute C1 given beta * D_gaussian - Gaussian distribution * D_reitz_init - compute c2 given beta * D_reitz - Trowbridge-Reitz * D_cook_init - compute m given beta * D_cook - Cook model * D_g - compute apparent roughness * D_tau - compute correlation distance * D_sigma - compute rms roughness * D_m - compute slope * D_coherent - coherent roughness attenuation * (per Beckmann) * D_incoherent - incoherent microfacet attenuation * (per Beckmann and Bahar) * D_Vxy - Vxy used in D_incoherent() * * ASSUMPTIONS: * > The following are defined in math.h * HUGE largest floating pont value * M_PI pi * M_PI_2 pi/2 * M_PI_4 pi/4 * M_SQRT2 root of 2 * * > The direction vectors N, V, and L are assumed to be oriented * to the same side of the surface. */ #include #include #include "geo.h" #include "D.h" /* **************************************************************** * double D_phong_init (beta) * double beta (in) - angle between N and * H where function = .5 * Returns Ns given beta for the Phong (1975) specular function, * per footnote to Eq.(\*(rH) */ double D_phong_init (beta) double beta; { if (beta <= 0.0) return HUGE; if (beta >= M_PI_4) return 0.0; return -(log(2.0) / log(cos(2.0 * beta))); } /* **************************************************************** * double D_phong (N, L, V, Ns) * DIR_VECT *N, *L, *V (in) - N, L, V vectors * double Ns (in) - specular exponent * * Returns the value of the Phong (1975) specular function given * surface normal, incident light direction, view direction, and * the specular exponent. Refer to Eqs.(\*(rU) and (\*(rI). */ double D_phong (N, L, V, Ns) DIR_VECT *N, *L, *V; double Ns; { double Rv_dot_L; if ((Rv_dot_L = geo_dot(geo_rfl(V,N),L)) < 0.0) return 0.0; return pow (Rv_dot_L, Ns); } /* **************************************************************** * double D_blinn_init (beta) * double beta (in) - angle between N and * H where function = .5 * Returns Ns given beta for the cosine power distribution * function presented by Blinn (1977). Refer to Eq.(\*(rHa) */ double D_blinn_init (beta) double beta; { if (beta <= 0.0) return HUGE; if (beta >= M_PI_2) return 0.0; return -(log(2.0) / log(cos(beta))); } /* **************************************************************** * double D_blinn (N, L, V, Ns) * DIR_VECT *N, *L, *V (in) - N, L, V vectors * double Ns (in) - specular exponent * * Returns the cosine power distribution function presented by * Blinn (1977) given surface normal, incident light direction, * view direction, and the specular exponent. Refer to Eq.(\*(e5) */ double D_blinn (N, L, V, Ns) DIR_VECT *N, *L, *V; double Ns; { double N_dot_H; if ((N_dot_H = geo_dot(N,geo_H(V,L))) < 0.0) return 0.0; return pow (N_dot_H, Ns); } /* **************************************************************** * double D_gaussian_init (beta) * double beta (in) - angle between N and * H where function = .5 * Returns C1 given beta for the Gaussian distribution presented * by Blinn (1977). Refer to Eq.(\*(rHb) */ double D_gaussian_init (beta) double beta; { if (beta <= 0.0) return HUGE; return sqrt(log(2.0)) / beta; } /* **************************************************************** * double D_gaussian (N, L, V, C1) * DIR_VECT *N, *L, *V (in) - N, L, V vectors * double C1 (in) - shape coefficient * * Returns the Gaussian distribution function presented by * Blinn (1977) given surface normal, incident light direction, * view direction, and the specular exponent. Refer to Eq.(\*(eA) */ double D_gaussian (N, L, V, C1) DIR_VECT *N, *L, *V; double C1; { double tmp; double N_dot_H; if ((N_dot_H = geo_dot(N,geo_H(V,L))) < 0.0) return 0.0; tmp = acos(N_dot_H) * C1; return exp (-(tmp * tmp)); } /* **************************************************************** * double D_reitz_init (beta) * double beta (in) - angle between N and * H where function = .5 * Returns C2 squared given beta for the Trowbridge-Reitz * distribution function presented by Blinn (1977). Refer to * Eq.(\*(rHc). C2 is squared for computational efficiency later. */ double D_reitz_init (beta) double beta; { double cos_beta; if (beta <= 0.0) return 0.0; cos_beta = cos(beta); return ((cos_beta * cos_beta) - 1.0) / ((cos_beta * cos_beta) - sqrt(2.0)); } /* **************************************************************** * double D_reitz (N, L, V, C2_2) * DIR_VECT *N, *L, *V (in) - N, L, V vectors * double C2_2 (in) - C2 squared * * Returns the Trowbridge-Reitz distribution function presented * by Blinn (1977) given surface normal, incident light * direction, view direction, and the specular exponent. Refer * to Eq.(\*(eB) */ double D_reitz (N, L, V, C2_2) DIR_VECT *N, *L, *V; double C2_2; { double tmp; double N_dot_H; if ((N_dot_H = geo_dot(N,geo_H(V,L))) < 0.0) return 0.0; tmp = C2_2 / ((N_dot_H * N_dot_H * (C2_2 - 1.0)) + 1.0); return tmp * tmp; } /* **************************************************************** * double D_cook_init (beta) * double beta (in) - angle between N and * H where function = .5 * Returns m squared corresponding to beta, refer to Eq.(\*(rT). * m is squared for computational efficiency in later use. */ double D_cook_init (beta) double beta; { double tan_beta; if (beta <= 0.0) return 0.0; if (beta >= M_PI_2) return HUGE; tan_beta = tan(beta); return -(tan_beta * tan_beta) / log(pow(cos(beta),4.0)/2.0); } /* **************************************************************** * double D_cook (N, L, V, m_2) * DIR_VECT *N, *L, *V (in) - N, L, V vectors * double m_2 (in) - m squared * * Returns microfacet distribution probability (Cook 1982) derived * from (Beckmann 1963). Refer to Eq.(\*(rX). */ double D_cook (N, L, V, m_2) DIR_VECT *N, *L, *V; double m_2; { double tmp; double N_dot_H; if ((N_dot_H = geo_dot(N,geo_H(V,L))) < 0.0) return 0.0; tmp = -(1.0 - (N_dot_H * N_dot_H)) / (N_dot_H * N_dot_H * m_2); return exp(tmp)/(4.0 * M_PI * m_2 * pow(N_dot_H,4.0)); } /* **************************************************************** * double D_g (N, L, V, sigma, lambda) * DIR_VECT *N, *L, *V (in) - N, L, V vectors * double sigma (in) RMS roughness, nm * double lambda (in) wavelength, nm * * Returns the apparent roughness, g, as given by Beckmann (1963). * Refer to Eq.(\*(rY). */ double D_g (N, L, V, sigma, lambda) DIR_VECT *N, *L, *V; double sigma, lambda; { double tmp; tmp = geo_dot(N,L) + geo_dot(N,V); return (4.0 * M_PI * M_PI * sigma * sigma * tmp * tmp) / (lambda * lambda); } /* **************************************************************** * double D_tau (sigma, m) * double sigma (in) rms roughness (nm) * double m (in) rms slope * * Returns tau (correlation distance) in nm. The relationship of * roughness, slope, and correlation distance given by * Beckmann (1963) is assumed, refer to Eq.(\*(rD). */ double D_tau (sigma, m) double sigma, m; { return (2.0 * sigma) / m; } /* **************************************************************** * double D_sigma (m, tau) * double m (in) rms slope * double tau (in) correlation distance (nm) * * Returns sigma (rms roughness). The relationship of roughness, * slope, and correlation distance given by Beckmann (1963) is * assumed, refer to Eq.(\*(rD). */ double D_sigma (m, tau) double m, tau; { return (m * tau)/ 2.0; } /* **************************************************************** * double D_m (sigma, tau) * double sigma (in) rms roughness (nm) * double tau (in) correlation distance (nm) * * Returns m (rms slope). The relationship of roughness, slope, * and correlation distance given by Beckmann (1963) is assumed, * refer to Eq.(\*(rD). */ double D_m (sigma, tau) double sigma, tau; { return (2.0 * sigma) / tau; } /* **************************************************************** * double D_coherent (g) * double g (in) apparent roughness * * Returns the coherent roughness attenuation given * by Beckmann (1963), refer to Eq.(\*(rZ). */ double D_coherent (g) double g; { return exp(-g); } /* **************************************************************** * double D_incoherent (N, L, V, m, g, lambda, tau) * DIR_VECT *N, *L, *V (in) N, L, V vectors * double m (in) m (slope) * double g (in) apparent roughness * double lambda (in) wavelength (nm) * double tau (in) correlation distance (nm) * * Returns the microfacet distribution function given by * Beckmann (1963) as interpreted by Koestner (1986), refer to * Eq.(\*(rC). */ double D_incoherent (N, L, V, m, g, lambda, tau) DIR_VECT *N, *L, *V; double m, g, lambda, tau; { DIR_VECT *H; double D1, D2; double denom; double N_dot_H, N_dot_H_2; if (g <= 0.0) return 0.0; if ((H = geo_H(V,L)) == NULL) return 0.0; if ((N_dot_H = geo_dot(N,H)) <= 0.0) return 0.0; N_dot_H_2 = N_dot_H * N_dot_H; denom = 4.0 * M_PI * m * m * N_dot_H_2 * N_dot_H_2; /* if g < 8.0, evaluate the series expansion, Eq.(\*(rCb) * terminating when a term falls below 0.00001. For small g * the evaluation terminates quickly. For large g the series * converges slowly. If g is near 8.0 the first terms of the * series will be less that the termination criteria. This * requires an additional termination criteria that the values * of the terms are decreasing at the time the termination * criteria is reached. */ if (g < 8.0) { double inc, ct, ct_factorial, g_pow; double last_inc, Vxz_t_over_4; Vxz_t_over_4 = (tau * tau * D_Vxz(N,L,V,lambda)) / 4.0; D1 = 0.0; ct = 1.0; ct_factorial = 1.0; g_pow = g * g; inc = 0.0; do { last_inc = inc; inc = (g_pow / (ct * ct_factorial)) * exp(-(g + Vxz_t_over_4/ct)); D1 += inc; ct += 1.0; ct_factorial *= ct; g_pow *= g; } while (((inc/denom) > 0.00001) || (inc > last_inc)); D1 /= denom; /* if g < 5.0, only the series expansion is used. If * 5.0 < g < 8.0, then we are in a transition range * for interpolation between the series expansion and * the convergence expression. */ if (g < 5.0) return D1; } /* if g > 5.0, evaluate the convergence expression, Eq.(\*(rCc) * (this routine would have returned earlier if G < 5.0). */ { double tmp; tmp = (N_dot_H_2 - 1.0) / (N_dot_H_2 * m * m); D2 = exp(tmp) / denom; /* if g > 8.0, only the convergence expression is used, * otherwise we are in a transition zone between the * convergent expression and series expansion. */ if (g > 8.0) return D2; } /* linear interpolation between D1 and D2 for * 5.0 < g < 8.0 */ return (((8.0 - g) * D1) + ((g - 5.0) * D2)) / 3.0; } /* **************************************************************** * double D_Vxz (N, L, V, lambda) * DIR_VECT *N, *L, *V (in) N, L, V vectors * double lambda (in) wavelength (nm) * * Returns the value of Vxz per Eq.(\*(rCd). */ double D_Vxz (N, L, V, lambda) DIR_VECT *N, *L, *V; double lambda; { double N_dot_L, N_dot_V, cos_phi; DIR_VECT sin_i, sin_r; double len_2_i, len_2_r; /* compute the sine squared of the incident angle */ N_dot_L = geo_dot(N, L); len_2_i = 1.0 - (N_dot_L * N_dot_L); /* compute the sine squared of the reflected angle */ N_dot_V = geo_dot(N, V); len_2_r = 1.0 - (N_dot_V * N_dot_V); /* if the sine of either the incident or the * reflected angle is zero, then the middle * term is zero. */ if ((len_2_i > 0.0) && (len_2_r > 0.0)) { /* the two sine vectors are of length equal to the sine of * the angles. The dot product is the cosine times the * lengths of the vectors (sines of the angles). */ sin_i.i = L->i - (N_dot_L * N->i); sin_i.j = L->j - (N_dot_L * N->j); sin_i.k = L->k - (N_dot_L * N->k); sin_r.i = V->i - (N_dot_V * N->i); sin_r.j = V->j - (N_dot_V * N->j); sin_r.k = V->k - (N_dot_V * N->k); cos_phi = geo_dot (&sin_i, &sin_r); return (4.0 * M_PI * M_PI * (len_2_i + (2.0 * cos_phi) + len_2_r)) / (lambda * lambda); } else { return (4.0 * M_PI * M_PI * (len_2_i + len_2_r)) / (lambda * lambda); } } /* ************************************************************* */ _i.i = L->i G.h 660 401 24 635 4623574122 4426 /* **************************************************************** * G.h * **************************************************************** * include file for the geometric attenuation functions */ #ifndef G_H #define G_H /* Microfacet distribution routines */ double G_torrance(); double G_sancer(); #endif G_H /* ************************************************************* */ * lambda); } } /* ************************************************************* */ _i.i = L->i G.c 660 401 24 5376 4623574117 4454 /* **************************************************************** * G.c * **************************************************************** * MODULE PURPOSE: * This module contains geometric attenuation * functions * * MODULE CONTENTS: * D_torrance - function by Torrance and Sparrow * D_Sancer - function by Sancer */ #include #include "geo.h" #define ROOT_PI 1.7724538509055159 /* **************************************************************** * double G_torrance (N, L, V) * DIR_VECT *N, *L, *V (in) - N, L, V vectors * double dummy (in) - just to make the * calls identical * * Returns geometric attenuation (Torrance and Sparrow * 1967), refer to Eq.(\*(eS). The direction vectors * N, V, and L are assumed to be oriented to the same * side of the surface. */ double G_torrance (N, L, V, dummy) DIR_VECT *N, *L, *V; double dummy; { double N_dot_H, V_dot_H; double N_dot_V, N_dot_L; double g1, g2; DIR_VECT *H; H = geo_H(V,L); N_dot_H = geo_dot (N,H); V_dot_H = geo_dot (V,H); N_dot_V = geo_dot (N,V); N_dot_L = geo_dot (N,L); if ((g1 = (2.0 * N_dot_H * N_dot_V) / V_dot_H) > 1.0) g1 = 1.0; if ((g2 = (2.0 * N_dot_H * N_dot_L) / V_dot_H) < g1) return g2; else return g1; } /* **************************************************************** * double G_sancer (N, L, V, m_2) * DIR_VECT *N, *L, *V (in) - N, L, V vectors * double m_2 (in) - m squared * * Returns geometric attenuation (Sancer 1969), refer to Eq.(\*(rT). * The direction vectors N, V, and L are assumed to be oriented * to the same side of the surface. * * NOTE: * > m can be related to beta using D_cook_init() * > This implementation assumes slope relates roughness and * correlation distance as described by Beckmann. This * explains the missing factor of 2 with m_2 when compared * to Eq.(\*(rT). */ double G_sancer (N, L, V, m_2) DIR_VECT *N, *L, *V; double m_2; { double N_dot_L, N_dot_V; double c1, c2, root_c1, root_c2, ci, cr; if (m_2 <= 0.0) return 1.0; if ((N_dot_L = geo_dot(N,L)) <= 0.0) return 0.0; c1 = (N_dot_L * N_dot_L) / (m_2 * (1.0 - (N_dot_L * N_dot_L))); root_c1 = sqrt(c1); ci = (exp(-c1) / (2.0 * ROOT_PI * root_c1)) - (erfc(root_c1) / 2.0); if ((N_dot_V = geo_dot(N,V)) <= 0.0) return 0.0; c2 = (N_dot_V * N_dot_V) / (m_2 * (1.0 - (N_dot_V * N_dot_V))); root_c2 = sqrt(c2); cr = (exp(-c2) / (2.0 * ROOT_PI * root_c2)) - (erfc(root_c2) / 2.0); return 1.0/(1.0 + ci + cr); } /* ************************************************************* */ )) <= 0.0) return 0.0; c1 = (N_dot_L * N_dot_L) / (m_2 * (1.0 - (N_dot_L * N_dot_L))); root_c1 = sqrt(c1); ci = (exp(-c1) / (2.0 * ROOT_PI * root_c1)) - (erfc(root_c1) / 2.0); if ((N_dot_V = geo_dot(N,V)) <= 0.0) return 0.0; c2 = (N_dotclr.h 660 401 24 3310 4630566317 5036 /* **************************************************************** * clr.h * **************************************************************** * include file for the geometric utilities */ #ifndef CLR_H #define CLR_H #ifndef TRUE #define TRUE 1 #endif TRUE #ifndef FALSE #define FALSE 0 #endif FALSE #define CLR_SAMPLE_MEYER 0 #define CLR_SAMPLE_HALL 1 #define CLR_SAMPLE_LIN_RAMP 2 /* common geometric constructs */ typedef struct {double r, g, b;} CLR_RGB; typedef struct {double x, y, z;} CLR_XYZ; typedef struct {double l, a, b;} CLR_LAB; typedef struct {double l, u, v;} CLR_LUV; /* color routine declarations */ int CLR_init(); int CLR_read_mtl(); int CLR_add_spect(); int CLR_mult_spect(); int CLR_scale_spect(); double CLR_area_spect(); CLR_XYZ CLR_spect_to_xyz(); int CLR_xyz_to_spect(); CLR_RGB CLR_spect_to_rgb(); int CLR_get_rgb(); int CLR_get_min_wl(); int CLR_get_max_wl(); int CLR_get_xyz_rgb(); int CLR_get_rgb_xyz(); int CLR_get_yiq_rgb(); int CLR_get_rgb_yiq(); int CLR_rgb_to_aux_rgb(); CLR_LAB CLR_xyz_to_lab(); CLR_LUV CLR_xyz_to_luv(); int CLR_t_concat(); int CLR_t_inverse(); CLR_RGB CLR_clamp_rgb(); CLR_RGB CLR_scale_rgb(); CLR_RGB CLR_clip_rgb(); int CLR_init_samples(); int CLR_num_samples(); int CLR_spect_to_sample(); int CLR_get_sample_rgb(); int CLR_get_sample_xyz(); int CLR_reconstruct(); int CLR_exit_samples(); #endif CLR_H /* ************************************************************* */ LAB CLR_xyz_to_lab(); CLR_LUV CLR_xyz_to_luv(); int CLR_t_concat(); int CLR_t_inverse(); CLR_RGB CLR_clamp_rgb(); CLR_RGB CLR_scale_rgb(); CLR_RGB CLR_clip_rgb(); int CLR_init_samples(); int CLR_num_samples(); int CLR_spect_to_sample(); int Cclr.c 660 401 24 77367 4650366572 5103 /* **************************************************************** * clr.c * **************************************************************** * MODULE PURPOSE: * This module contains routines for color * transformations, color space manipulations, * and spectral sampling. * * MODULE CONTENTS: * CLR_init - initialized the color routines * CLR_read_mtl - read a material file * CLR_blackbody - returns the spectral curve for * a blackbody at the specified * temperature. * CLR_add_spect - add two spectral curves * CLR_mult_spect - multiply two spectral curves * CLR_scale_spect - scale a spectral curve * CLR_area_spect - get the area under a curve * CLR_spect_to_xyz - sample a curve to xyz * CLR_xyz_to_spect - reconstruct a spectral curve * CLR_spect_to_rgb - sample a curve to rgb * CLR_get_CIEXYZ - get the CIEXYZ matching functions * CLR_get_rgb - returns the color space primaries * CLR_get_min_wl - returns the max wavelength * CLR_get_max_wl - returns the min wavelength * CLR_get_xyz_rgb - returns the xyz to rgb matrix * CLR_get_rgb_xyz - returns the rgb to xyz matrix * CLR_get_rgb_yiq - returns the rgb to yiq matrix * CLR_get_yiq_rgb - returns the yiq to rgb matrix * CLR_rgb_to_aux_rgb - returns the matrix from the * current color space to an * auxiliary color space * CLR_xyz_to_lab - transforms from xyz to Lab * CLR_xyz_to_luv - transforms from xyz to Luv * CLR_t_concat - concatenate color transforms * CLR_t_inverse - return an inverse transform * CLR_exit - finish with the color routines * * NOTES: * > Spectral curve manipulations use the lower and upper * bounds established at init. Curves are at 1nm * increments and are arrays of doubles, (max - min + 1) * elements long. * > Matrix scaling is so that an RGB of 1,1,1 transforms * to an XYZ with a Y value of 1.0 * > The XYZ and RGB sampling from spectral curves are * scaled so that that a mirror reflector (1.0 for all * wavelengths) has a Y value of 1.0. * > Wavelength bounds of 380nm to 780nm (the visible * spectrum) are appropriate for most applications. */ #include #include #include "clr.h" #define RED 0 #define GREEN 1 #define BLUE 2 #define WHITE 3 #define LOAD_MAT(a,b) \ { int ii,jj; \ for (ii=0; ii<=2; ii++) \ for (jj=0; jj<=2; jj++) a[ii][jj] = b[ii][jj]; \ } static int init = FALSE; static int min_wl = 380; static int max_wl = 780; static double xyz_rgb_mat[3][3]; static double rgb_xyz_mat[3][3]; static double *X_tristim = NULL; static double *Y_tristim = NULL; static double *Z_tristim = NULL; static double *work_curve = NULL; static double *reconstruct[3] = NULL; static double XYZ_to_recon[3][3]; static double xyz_scale = 1.0; static CLR_XYZ rgb_primary[4]; static int clr__cspace_to_xyz (); static CLR_XYZ white; static CLR_LUV ref_luv; /* This is the RGB to YIQ matrix and YIQ to RGB matrix as * given in Sect 5.2 NTSC and RGB Video */ static double rgb_yiq_mat[3][3] = { {0.299, 0.587, 0.114}, {0.596, -0.275, -0.321}, {0.212, -0.523, 0.311}}; static double yiq_rgb_mat[3][3] = { {1.0000, 0.9557, 0.6199}, {1.0000, -0.2716, -0.6469}, {1.0000, -1.1082, 1.7051}}; /* This is the NTSC primaries with D6500 white point for use as * the default initialization as given in sect 5.1.1 Color * Correction for Display. */ static CLR_XYZ rgb_NTSC[4] = { {0.670, 0.330, 0.000}, /* red */ {0.210, 0.710, 0.080}, /* green */ {0.140, 0.080, 0.780}, /* blue */ {0.313, 0.329, 0.358}}; /* white */ /* This is the data for the CIEXYZ curves take from Judd and * Wyszecki (1975), table 2.6, these are for the 1931 standard * observer with a 2-degree visual field. */ static double CIEXYZ[81][4] = { {380, 0.0014, 0.0000, 0.0065}, {385, 0.0022, 0.0001, 0.0105}, {390, 0.0042, 0.0001, 0.0201}, {395, 0.0076, 0.0002, 0.0362}, {400, 0.0143, 0.0004, 0.0679}, {405, 0.0232, 0.0006, 0.1102}, {410, 0.0435, 0.0012, 0.2074}, {415, 0.0776, 0.0022, 0.3713}, {420, 0.1344, 0.0040, 0.6456}, {425, 0.2148, 0.0073, 1.0391}, {430, 0.2839, 0.0116, 1.3856}, {435, 0.3285, 0.0168, 1.6230}, {440, 0.3483, 0.0230, 1.7471}, {445, 0.3481, 0.0298, 1.7826}, {450, 0.3362, 0.0380, 1.7721}, {455, 0.3187, 0.0480, 1.7441}, {460, 0.2908, 0.0600, 1.6692}, {465, 0.2511, 0.0739, 1.5281}, {470, 0.1954, 0.0910, 1.2876}, {475, 0.1421, 0.1126, 1.0419}, {480, 0.0956, 0.1390, 0.8130}, {485, 0.0580, 0.1693, 0.6162}, {490, 0.0320, 0.2080, 0.4652}, {495, 0.0147, 0.2586, 0.3533}, {500, 0.0049, 0.3230, 0.2720}, {505, 0.0024, 0.4073, 0.2123}, {510, 0.0093, 0.5030, 0.1582}, {515, 0.0291, 0.6082, 0.1117}, {520, 0.0633, 0.7100, 0.0782}, {525, 0.1096, 0.7932, 0.0573}, {530, 0.1655, 0.8620, 0.0422}, {535, 0.2257, 0.9149, 0.0298}, {540, 0.2904, 0.9540, 0.0203}, {545, 0.3597, 0.9803, 0.0134}, {550, 0.4334, 0.9950, 0.0087}, {555, 0.5121, 1.0000, 0.0057}, {560, 0.5945, 0.9950, 0.0039}, {565, 0.6784, 0.9786, 0.0027}, {570, 0.7621, 0.9520, 0.0021}, {575, 0.8425, 0.9154, 0.0018}, {580, 0.9163, 0.8700, 0.0017}, {585, 0.9786, 0.8163, 0.0014}, {590, 1.0263, 0.7570, 0.0011}, {595, 1.0567, 0.6949, 0.0010}, {600, 1.0622, 0.6310, 0.0008}, {605, 1.0456, 0.5668, 0.0006}, {610, 1.0026, 0.5030, 0.0003}, {615, 0.9384, 0.4412, 0.0002}, {620, 0.8544, 0.3810, 0.0002}, {625, 0.7514, 0.3210, 0.0001}, {630, 0.6424, 0.2650, 0.0000}, {635, 0.5419, 0.2170, 0.0000}, {640, 0.4479, 0.1750, 0.0000}, {645, 0.3608, 0.1382, 0.0000}, {650, 0.2835, 0.1070, 0.0000}, {655, 0.2187, 0.0816, 0.0000}, {660, 0.1649, 0.0610, 0.0000}, {665, 0.1212, 0.0446, 0.0000}, {670, 0.0874, 0.0320, 0.0000}, {675, 0.0636, 0.0232, 0.0000}, {680, 0.0468, 0.0170, 0.0000}, {685, 0.0329, 0.0119, 0.0000}, {690, 0.0227, 0.0082, 0.0000}, {695, 0.0158, 0.0057, 0.0000}, {700, 0.0114, 0.0041, 0.0000}, {705, 0.0081, 0.0029, 0.0000}, {710, 0.0058, 0.0021, 0.0000}, {715, 0.0041, 0.0015, 0.0000}, {720, 0.0029, 0.0010, 0.0000}, {725, 0.0020, 0.0007, 0.0000}, {730, 0.0014, 0.0005, 0.0000}, {735, 0.0010, 0.0004, 0.0000}, {740, 0.0007, 0.0002, 0.0000}, {745, 0.0005, 0.0002, 0.0000}, {750, 0.0003, 0.0001, 0.0000}, {755, 0.0002, 0.0001, 0.0000}, {760, 0.0002, 0.0001, 0.0000}, {765, 0.0001, 0.0000, 0.0000}, {770, 0.0001, 0.0000, 0.0000}, {775, 0.0001, 0.0000, 0.0000}, {780, 0.0000, 0.0000, 0.0000} }; /* **************************************************************** * CLR_init (min, max, rgb) * int min, max (in) wavelength bounds (nm) * CLR_XYZ rgb[4] (in) the primaries. If NULL, * then the NTSC primaries with * a D6500 white point are used. * * Initializes the color routines for a spectral range from * 'min'nm to 'max'nm and an RGB color space given by 'rgb' * Returns TRUE if successful and FALSE on error. */ int CLR_init (min, max, rgb) int min; int max; CLR_XYZ rgb[4]; { int color, ct; double len; char *malloc(); if (init) return FALSE; init = TRUE; min_wl = min; max_wl = max; /* load primaries and build transformations, use the * defaults is rgb==NULL */ for (color=RED; color<=WHITE; color++) { if (rgb == NULL) rgb_primary[color] = rgb_NTSC[color]; else rgb_primary[color] = rgb[color]; rgb_primary[color].z = 1.0 - rgb_primary[color].x - rgb_primary[color].y; } if (!clr__cspace_to_xyz(rgb_primary, rgb_xyz_mat)) goto error; if (!CLR_t_inverse(rgb_xyz_mat, xyz_rgb_mat)) goto error; /* build reference info for L*a*b*, L*u*v* transforms */ white.x = (100.0 * rgb_primary[WHITE].x) / rgb_primary[WHITE].y; white.y = 100.0; white.z = (100.0 * rgb_primary[WHITE].z) / rgb_primary[WHITE].y; ref_luv.u = (4.0 * white.x) / (white.x + (15.0 * white.y) + (3.0 * white.z)); ref_luv.v = (9.0 * white.y) / (white.x + (15.0 * white.y) + (3.0 * white.z)); /* build the XYZ sampling curves - allocate space for the * curves and interpolate the CIEXYZ table into continuous * curves at 1nm increments. */ { int ii, nm; double *x, *y, *z; double x_cur, y_cur, z_cur; double x_inc, y_inc, z_inc; if ((work_curve = (double *)malloc((unsigned) (sizeof(double) * (max_wl-min_wl+1)))) == NULL) goto error; if ((x = X_tristim = (double *)malloc((unsigned) (sizeof(double) * (max_wl-min_wl+1)))) == NULL) goto error; if ((y = Y_tristim = (double *)malloc((unsigned) (sizeof(double) * (max_wl-min_wl+1)))) == NULL) goto error; if ((z = Z_tristim = (double *)malloc((unsigned) (sizeof(double) * (max_wl-min_wl+1)))) == NULL) goto error; for (nm=min_wl; nm<380; nm++) *x++ = *y++ = *z++ = 0.0; for (ct=0; ct<80; ct++) { x_cur = CIEXYZ[ct][1]; y_cur = CIEXYZ[ct][2]; z_cur = CIEXYZ[ct][3]; x_inc = (CIEXYZ[ct+1][1] - x_cur) / 5.0; y_inc = (CIEXYZ[ct+1][2] - y_cur) / 5.0; z_inc = (CIEXYZ[ct+1][3] - z_cur) / 5.0; for (ii=0; ii<5; ii++, nm++) { if (nm > max_wl) goto XYZ_done; if (nm >= min_wl) { *x++ = x_cur; *y++ = y_cur; *z++ = z_cur; } x_cur += x_inc; y_cur += y_inc; z_cur += z_inc; } } if (nm <= max_wl) { *x++ = x_cur; *y++ = y_cur; *z++ = z_cur; nm++; } for ( ;nm<=max_wl; nm++) *x++ = *y++ = *z++ = 0.0; } XYZ_done: /* determine the scaling factor to be used in sampling spectral * curves to bring Y of a sampled identity curve to 1. */ xyz_scale = 1.0 / CLR_area_spect (Y_tristim); /* allocate and build the reconstruction curves and the * transformation from XYZ to the reconstruction weights */ { double reconst_to_XYZ[3][3]; CLR_XYZ xyz; if ((reconstruct[0] = (double *)malloc((unsigned) (sizeof(double) * (3 * (max_wl-min_wl+1))))) == NULL) goto error; reconstruct[1] = reconstruct[0] + (max_wl-min_wl+1); reconstruct[2] = reconstruct[1] + (max_wl-min_wl+1); for (ct=min; ct<=max; ct++) { *(reconstruct[0] + ct - min) = 1.0; *(reconstruct[1] + ct - min) = sin((double)((2.0 * M_PI / 400.0) * ((double)ct - 380.0))); *(reconstruct[2] + ct - min) = cos((double)((2.0 * M_PI / 400.0) * ((double)ct - 380.0))); } for (ct=0; ct<=2; ct++) { xyz = CLR_spect_to_xyz (reconstruct[ct]); reconst_to_XYZ[ct][0] = xyz.x; reconst_to_XYZ[ct][1] = xyz.y; reconst_to_XYZ[ct][2] = xyz.z; } CLR_t_inverse (reconst_to_XYZ, XYZ_to_recon); } return TRUE; error: CLR_exit(); return FALSE; } /* **************************************************************** * CLR_read_mtl (file, conduct, n, k, curve) * char *file (in) material file name * int *conduct (out) conductor - set to FALSE if not * specified in the material file * double *n, *k (out) n and k, - set to 0.0 if not * specified in the material file * double *curve (out) curve array * * Returns TRUE if the material was read, FALSE if the material * could not be found or an error was detected in the material * file. NULL pointers can be given as arguments for 'conduct', * 'n', 'k', and/or 'curve' if the properties are not of interest. * * The material file must be in the format given in Sect III.7.1 * Material Data. */ int CLR_read_mtl (file, conduct, n, k, curve) char *file; int *conduct; double *n, *k; double *curve; { FILE *fp, *fopen(); char in_str[200]; int flag = 1; int cur_wl, last_wl, ind_wl; double cur_val, last_val, inc_val; if (!init) return FALSE; if ((fp=fopen(file,"r")) == NULL) return FALSE; if (n != NULL) *n = 0.0; if (k != NULL) *k = 0.0; if (conduct != NULL) *conduct = FALSE; if (curve == NULL) flag = -1; ind_wl = min_wl; while (fgets(in_str,200,fp)) { switch (in_str[0]) { case '#': /* comment */ break; case 'k': /* absorption coefficient */ if (k != NULL) (void)sscanf (in_str, "k %lf", k); break; case 'n': /* index of refraction */ if (n != NULL) (void)sscanf (in_str, "n %lf", n); break; case 'c': /* conductor */ if (conduct != NULL) *conduct = TRUE; break; case 'd': /* dielectric */ if (conduct != NULL) *conduct = FALSE; break; default: /* spectral data */ if ((flag != -1) && (sscanf(in_str, "%d %lf", &cur_wl, &cur_val) == 2)) { if (flag == 1) { flag = 0; while (ind_wl <= cur_wl) { *curve++ = cur_val; if (++ind_wl > max_wl) { flag == -1; break; } } } else { if (cur_wl < last_wl) { (void)fclose(fp); return FALSE; } inc_val = (cur_val - last_val) / (cur_wl - last_wl); while (last_wl < ind_wl) { last_val += inc_val; last_wl++; } while (ind_wl <= cur_wl) { *curve++ = last_val; last_val += inc_val; if (++ind_wl > max_wl) { flag = -1; break; } } } last_wl = cur_wl; last_val = cur_val; } } } while (ind_wl <= max_wl) { *curve++ = last_val; ind_wl++; } (void)fclose(fp); return TRUE; } /* **************************************************************** * CLR_blackbody (temp, curve) * double temp (in) temperature (degrees K) * double *curve (out) curve array * * Fills the 'curve' array with the illumination values for an * ideal blackbody radiator at the specified temperature. * Following the convention used for normalizing other illuminant * curves, the blackbody curve is normalized so that the value * at 560nm is 1.0. The formulation is taken from Judd and * Wysecki (1975) p.106. NOTE: c1 and c2 have been adjusted * to account for the wavelength to be in nm. */ int CLR_blackbody (temp, curve) double temp; double *curve; { double norm; double c1 = 3.74150E29; double c2 = 1.4388E7; double wl; int ct; if (temp < 0.0) return FALSE; norm = (pow((double)560.0, (double)5.0) * (exp (c2 / (560.0 * temp)) - 1.0)) / c1; for (ct=max_wl-min_wl+1, wl=min_wl; --ct>=0; wl += 1.0) { *curve++ = (norm * c1) / (pow(wl, (double)5.0) * (exp (c2 / (wl * temp)) - 1.0)); } return TRUE; } /* **************************************************************** * CLR_add_spect (c1, c2, c3) * double *c1, *c2 (in) - input curves * double *c3 (out) - output curve * * Add spectral curve 'c1' to spectral curve 'c2' to get spectral * curve 'c3'. Returns TRUE if successful, FALSE if the CLR_ * routines are not initialized. */ CLR_add_spect (c1,c2,c3) double *c1, *c2, *c3; { int ct; if (!init) return FALSE; for (ct=max_wl-min_wl+1; --ct>=0; ) *c3++ = *c1++ + *c2++; return TRUE; } /* **************************************************************** * CLR_mult_spect (c1, c2, c3) * double *c1, *c2 (in) - input curves * double *c3 (out) - output curve * * Multiply spectral curve 'c1' by spectral curve 'c2' to get * spectral curve 'c3'. Returns TRUE if successful, FALSE if * the CLR_ routines are not initialized. */ CLR_mult_spect (c1,c2,c3) double *c1, *c2, *c3; { int ct; if (!init) return FALSE; for (ct=max_wl-min_wl+1; --ct>=0; ) *c3++ = *c1++ * *c2++; return TRUE; } /* **************************************************************** * CLR_scale_spect (c1, scale, c3) * double *c1 (in) - input curve * double scale (in) - scale * double *c3 (out) - output curve * * Multiply spectral curve 'c1' by 'scale' to get spectral curve * 'c3'. Returns TRUE if successful, FALSE if the CLR_ routines * are not initialized. */ CLR_scale_spect (c1,scale,c3) double *c1, scale, *c3; { int ct; if (!init) return FALSE; for (ct=max_wl-min_wl+1; --ct>=0; ) *c3++ = *c1++ * scale; return TRUE; } /* **************************************************************** * double CLR_area_spect (c1) * double *c1 (in) - input curve * * Returns the are under the spectral curve 'c1'. Returns 0 if * the CLR_ routines are not initialized. */ double CLR_area_spect (c1) double *c1; { int ct; double area = 0.0; if (!init) return FALSE; for (ct=max_wl-min_wl+1; --ct>=0; ) area += *c1++; return area; } /* **************************************************************** * CLR_XYZ CLR_spect_to_xyz(spectral) * double *spectral (in) - the spectral curve * * Returns the sample values if successful, and a sample value * of (0,0,0) if the CLR_ routines are not initialized. * * Multiplies the spectral curve by each of the sampling curves * then integrates the resulting curves. The XYZ values are then * normalized such that an identity material has Y=1. */ CLR_XYZ CLR_spect_to_xyz (spectral) double *spectral; { CLR_XYZ xyz; if (!init) {xyz.x = xyz.y = xyz.z = 0.0; return xyz;} (void)CLR_mult_spect(spectral, X_tristim, work_curve); xyz.x = xyz_scale * CLR_area_spect(work_curve); (void)CLR_mult_spect(spectral, Y_tristim, work_curve); xyz.y = xyz_scale * CLR_area_spect(work_curve); (void)CLR_mult_spect(spectral, Z_tristim, work_curve); xyz.z = xyz_scale * CLR_area_spect(work_curve); return xyz; } /* **************************************************************** * CLR_xyz_to_spect(xyz,spectral) * CLR_XYZ xyz; (in) - the XYZ value of the color * double *spectral (out) - the reconstructed spectral curve * * Reconstructs a spectral curve using the fitting functions * described in Glassner (1989). */ CLR_xyz_to_spect(xyz, spectral) CLR_XYZ xyz; double *spectral; { double recon[3]; double *r0, *r1, *r2; int ct; recon[0] = (xyz.x * XYZ_to_recon[0][0]) + (xyz.y * XYZ_to_recon[1][0]) + (xyz.z * XYZ_to_recon[2][0]); recon[1] = (xyz.x * XYZ_to_recon[0][1]) + (xyz.y * XYZ_to_recon[1][1]) + (xyz.z * XYZ_to_recon[2][1]); recon[2] = (xyz.x * XYZ_to_recon[0][2]) + (xyz.y * XYZ_to_recon[1][2]) + (xyz.z * XYZ_to_recon[2][2]); r0 = reconstruct[0]; r1 = reconstruct[1]; r2 = reconstruct[2]; for (ct=max_wl-min_wl+1; --ct>=0; ) *spectral++ = (recon[0] * *r0++) + (recon[1] * *r1++) + (recon[2] * *r2++); return TRUE; } /* **************************************************************** * CLR_RGB CLR_spect_to_rgb (spectral) * double *spectral (in) - the spectral curve * * Returns the sample values if successful, and a sample value of * (0,0,0) if the CLR_ routines are not initialized. * * The curve is first sampled to XYZ then transformed to RGB */ CLR_RGB CLR_spect_to_rgb (spectral) double *spectral; { CLR_XYZ xyz; CLR_RGB rgb; if (!init) {rgb.r = rgb.g = rgb.b = 0.0; return rgb;} xyz = CLR_spect_to_xyz (spectral); rgb.r = (xyz_rgb_mat[0][0] * xyz.x) + (xyz_rgb_mat[0][1] * xyz.y) + (xyz_rgb_mat[0][2] * xyz.z); rgb.g = (xyz_rgb_mat[1][0] * xyz.x) + (xyz_rgb_mat[1][1] * xyz.y) + (xyz_rgb_mat[1][2] * xyz.z); rgb.b = (xyz_rgb_mat[2][0] * xyz.x) + (xyz_rgb_mat[2][1] * xyz.y) + (xyz_rgb_mat[2][2] * xyz.z); return rgb; } /* **************************************************************** * CLR_get_CIEXYZ (X,Y,Z) * double *X,*Y,*Z (mod) arrays to hold the curves * * Copies the XYZ sampling curves into the user supplied buffers. * Returns TRUE is successful and FALSE if the CLR_ routines are * not initialized. */ CLR_get_CIEXYZ(X,Y,Z) double *X,*Y,*Z; { int ct; double *x,*y,*z; if (!init) return FALSE; x = X_tristim; y = Y_tristim; z = Z_tristim; for (ct=max_wl-min_wl+1; --ct>=0; ) { *X++ = *x++; *Y++ = *y++; *Z++ = *z++; } return TRUE; } /* **************************************************************** * CLR_get_rgb (rgb) * CLR_XYZ rgb[4] (mod) the primaries * * Fills 'rgb' with the primaries the CLR_ routines are initialized * to. Returns TRUE if successful and FALSE if the CLR_ routines * are not initialized. */ CLR_get_rgb (rgb) CLR_XYZ *rgb; { int ct; if (!init) return FALSE; for (ct=0; ct<=3; ct++) rgb[ct] = rgb_primary[ct]; return TRUE; } /* **************************************************************** * CLR_get_min_wl () * * Returns the minimum wavelength bound for which the CLR_ routines * are initialized, returns -1 if the CLR_ routines are not * initialized. */ CLR_get_min_wl() { if (!init) return -1; return min_wl; } /* **************************************************************** * CLR_get_max_wl () * * Returns the minimum wavelength bound for which the CLR_ routines * are initialized, returns -1 if the CLR_ routines are not * initialized. */ CLR_get_max_wl() { if (!init) return -1; return max_wl; } /* **************************************************************** * CLR_get_xyz_rgb (mat) * double mat[3][3] (mod) - matrix to be loaded * * Copies the CIEXYZ to RGB transformation into the user supplied * matrix. Returns TRUE if successful and FALSE if the CLR_ * routines are not initialized. */ int CLR_get_xyz_rgb(mat) double mat[][3]; { if (!init) return FALSE; LOAD_MAT(mat, xyz_rgb_mat); return TRUE; } /* **************************************************************** * CLR_get_rgb_xyz (mat) * double mat[3][3] (mod) - matrix to be loaded * * Copies the RGB to CIEXYZ transformation into the user supplied * matrix. Returns TRUE if successful and FALSE if the CLR_ * routines are not initialized. */ int CLR_get_rgb_xyz(mat) double mat[][3]; { if (!init) return FALSE; LOAD_MAT(mat, rgb_xyz_mat); return TRUE; } /* **************************************************************** * CLR_get_yiq_rgb (mat) * double mat[3][3] (mod) - matrix to be loaded * * Copies the YIQ to RGB transformation into the user supplied * matrix. Returns TRUE if successful and FALSE if the CLR_ * routines are not initialized. */ int CLR_get_yiq_rgb(mat) double mat[][3]; { if (!init) return FALSE; LOAD_MAT(mat, yiq_rgb_mat); return TRUE; } /* **************************************************************** * CLR_get_rgb_yiq (mat) * double mat[3][3] (mod) - matrix to be loaded * * Copies the RGB to YIQ transformation into the user supplied * matrix. Returns TRUE if successful and FALSE if the CLR_ * routines are not initialized. */ int CLR_get_rgb_yiq(mat) double mat[][3]; { if (!init) return FALSE; LOAD_MAT(mat, rgb_yiq_mat); return TRUE; } /* **************************************************************** * CLR_rgb_to_aux_rgb (to, from, rgb_aux) * double to[3][3] (mod) - matrix to aux rgb * double from[3][3] (mod) - matrix from aux rgb * CLR_XYZ rgb_aux[4] (in) - the color space definition, * 3 primaries and white * * Creates the transformations between RGB color space the CLR_ * routines are initialized for and another RGB color space as * specified in 'aux_rgb'. The 'to' and 'from' matrices are * filled with the resulting transformations. TRUE is returned * if successful, FALSE is returned if the CLR_ routines are not * initialized or if there is a singularity detected when the * transformation is being generated. * * The technique used is described in Sect 3.2, Colorimetry and * the RGB Monitor. */ int CLR_rgb_to_aux_rgb(to,from,rgb_aux) double to[][3], from[][3]; CLR_XYZ rgb_aux[]; { CLR_XYZ rgb_tmp[4]; double rgb_xyz_aux[3][3], xyz_rgb_aux[3][3]; int color; if (!init) return FALSE; /* normalize the chromaticities of the auxiliary primaries * and white point. */ for (color=RED; color<=WHITE; color++) { rgb_tmp[color] = rgb_aux[color]; rgb_tmp[color].z = 1.0 - rgb_aux[color].x - rgb_aux[color].y; } /* get the transform between XYZ and the auxiliary RGB */ if (!clr__cspace_to_xyz(rgb_tmp, rgb_xyz_aux)) return FALSE; if (!CLR_t_inverse(rgb_xyz_aux,xyz_rgb_aux)) return FALSE; /* concatenate with the transforms for the RGB color space * that the CLR_ routine set is initialized to */ (void)CLR_t_concat(xyz_rgb_aux, rgb_xyz_mat, to); (void)CLR_t_concat(xyz_rgb_mat, rgb_xyz_aux, from); return TRUE; } /* **************************************************************** * CLR_LAB CLR_xyz_to_lab (xyz) * CLR_XYZ xyz (in) - xyz color (.01<=Y<=1) * * Returns L*a*b* given XYZ. Returns (0,0,0) if the CLR_ routines * are not initialized. This transformation is described in * Section 3.3, Alternate Color Representations, Eq.(\*(sC), * and is taken from Judd and Wyszecki (1975) pg.320. The * transformation is scaled for (.01 < Y < 1) for consistency * within the CLR_ routine set. */ CLR_LAB CLR_xyz_to_lab(xyz) CLR_XYZ xyz; { CLR_LAB lab; if (!init) {lab.l = lab.a = lab.b = 0.0; return lab;} xyz.x *= 100.0; xyz.y *= 100.0; xyz.z *= 100.0; lab.l = 25.0 * pow(((100.0*xyz.y)/white.y),.33333) - 16.0; lab.a = 500.0 * (pow(xyz.x/white.x,.33333) - pow(xyz.y/white.y,.33333)); lab.b = 200.0 * (pow(xyz.y/white.y,.33333) - pow(xyz.z/white.z,.33333)); return lab; } /* **************************************************************** * CLR_LUV CLR_xyz_to_luv (xyz) * CLR_XYZ xyz (in) - xyz color (.01<=Y<=1) * * Returns L*u*v* given XYZ. Returns (0,0,0) if the CLR_ routines * are not initialized. This transformation is described in * Section 3.3, Alternate Color Representations, Eq.(\*(sD), * and is taken from Judd and Wyszecki (1975) pg.328. The * transformation is scaled for (.01 < Y < 1) for consistency * within the CLR_ routine set. */ CLR_LUV CLR_xyz_to_luv(xyz) CLR_XYZ xyz; { CLR_LUV luv; double u, v; if (!init) {luv.l = luv.u = luv.v = 0.0; return luv;} xyz.x *= 100.0; xyz.y *= 100.0; xyz.z *= 100.0; luv.l = 25.0 * pow(((100.0*xyz.y)/white.y),.33333) - 16.0; u = (4.0 * xyz.x) / (xyz.x + (15.0 * xyz.y) + (3.0 * xyz.z)); v = (9.0 * xyz.y) / (xyz.x + (15.0 * xyz.y) + (3.0 * xyz.z)); luv.u = 13.0 * luv.l * (u - ref_luv.u); luv.v = 13.0 * luv.l * (v - ref_luv.v); return luv; } /* **************************************************************** * CLR_t_concat (m1, m2, m3) * double m1[3][3] (in) - matrix * double m2[3][3] (in) - matrix to concat * double m3[3][3] (in) - concatenated matrix * * Concatenate 'm1' to 'm2' resulting in 'm3'. In use, suppose you * have an XYZ to RGB and an RGB to YIQ matrix. Concatenate the * RGB to YIQ matrix to the XYZ to RGB matrix to get an XYZ to * YIQ matrix. Returns TRUE. */ int CLR_t_concat(m1,m2,m3) double m1[][3], m2[][3], m3[][3]; { double t1[3][3], t2[3][3]; LOAD_MAT(t1,m1); LOAD_MAT(t2,m2); { int ii,jj; for (ii=0; ii<=2; ii++) for (jj=0; jj<=2; jj++) m3[ii][jj] = (t1[ii][0] * t2[0][jj]) + (t1[ii][1] * t2[1][jj]) + (t1[ii][2] * t2[2][jj]); } return TRUE; } /* **************************************************************** * CLR_t_inverse (mat, inv_mat) * double mat[3][3] (in) - matrix to be inverted * double inv_mat[3][3] (mod) - inverted matrix * * Inverts 'mat' using Gaussian elimination. Returns TRUE if * successful and FALSE if there is a singularity. */ int CLR_t_inverse (mat, inv_mat) double mat[][3]; double inv_mat[][3]; { int ii, jj, kk; double tmp_mat[3][3], tmp_d; for (ii=0; ii<=2; ii++) for (jj=0; jj<=2; jj++) { tmp_mat[ii][jj] = mat[ii][jj]; inv_mat[ii][jj] = (ii==jj ? 1.0 : 0.0); } for (ii=0; ii<=2; ii++) { for (jj=ii+1, kk=ii; jj<=2; jj++) if (fabs(tmp_mat[jj][ii]) > fabs(tmp_mat[kk][ii])) kk = jj; /* check for singularity */ if (tmp_mat[kk][ii] == 0.0) return FALSE; /* pivot - switch rows kk and ii */ if (kk != ii) for (jj=0; jj<=2; jj++) { tmp_d = tmp_mat[ii][jj]; tmp_mat[ii][jj] = tmp_mat[kk][jj]; tmp_mat[kk][jj] = tmp_d; tmp_d = inv_mat[ii][jj]; inv_mat[ii][jj] = inv_mat[kk][jj]; inv_mat[kk][jj] = tmp_d; } /* normalize the row - make the diagonal 1 */ for (tmp_d = 1.0 / tmp_mat[ii][ii], jj=0; jj<=2; jj++) { tmp_mat[ii][jj] *= tmp_d; inv_mat[ii][jj] *= tmp_d; } /* zero the non-diagonal terms in this column */ for (jj=0; jj<=2; jj++) if (jj != ii) for (tmp_d = -tmp_mat[jj][ii], kk=0; kk<=2; kk++) { tmp_mat[jj][kk] += tmp_mat[ii][kk] * tmp_d; inv_mat[jj][kk] += inv_mat[ii][kk] * tmp_d; } } return TRUE; } /* **************************************************************** * clr__cspace_to_xyz (cspace, t_mat) * CLR_XYZ cspace[4] (in) - the color space definition, * 3 primaries and white * double t_mat[3][3] (mod) - the color transformation * * Builds the transformation from a set of primaries to the CIEXYZ * color space. This is the basis for the generation of the color * transformations in the CLR_ routine set. The method used is * that detailed in Sect 3.2 Colorimetry and the RGB monitor. * Returns FALSE if there is a singularity. */ static int clr__cspace_to_xyz (cspace, t_mat) CLR_XYZ cspace[]; double t_mat[][3]; { int ii, jj, kk, tmp_i, ind[3]; double mult, white[3], scale[3]; /* normalize the white point to Y=1 */ if (cspace[WHITE].y <= 0.0) return FALSE; white[0] = cspace[WHITE].x / cspace[WHITE].y; white[1] = 1.0; white[2] = cspace[WHITE].z / cspace[WHITE].y; for (ii=0; ii<=2; ii++) { t_mat[0][ii] = cspace[ii].x; t_mat[1][ii] = cspace[ii].y; t_mat[2][ii] = cspace[ii].z; ind[ii] = ii; } /* gaussian elimination with partial pivoting */ for (ii=0; ii<2; ii++) { for (jj=ii+1; jj<=2; jj++) if (fabs(t_mat[ind[jj]][ii]) > fabs(t_mat[ind[ii]][ii])) { tmp_i=ind[jj]; ind[jj]=ind[ii]; ind[ii]=tmp_i; } if (t_mat[ind[ii]][ii] == 0.0) return FALSE; for (jj=ii+1; jj<=2; jj++) { mult = t_mat[ind[jj]][ii] / t_mat[ind[ii]][ii]; for (kk=ii+1; kk<=2; kk++) t_mat[ind[jj]][kk] -= t_mat[ind[ii]][kk] * mult; white[ind[jj]] -= white[ind[ii]] * mult; } } if (t_mat[ind[2]][2] == 0.0) return FALSE; /* back substitution to solve for scale */ scale[ind[2]] = white[ind[2]] / t_mat[ind[2]][2]; scale[ind[1]] = (white[ind[1]] - (t_mat[ind[1]][2] * scale[ind[2]])) / t_mat[ind[1]][1]; scale[ind[0]] = (white[ind[0]] - (t_mat[ind[0]][1] * scale[ind[1]]) - (t_mat[ind[0]][2] * scale[ind[2]])) / t_mat[ind[0]][0]; /* build matrix */ for (ii=0; ii<=2; ii++) { t_mat[0][ii] = cspace[ii].x * scale[ii]; t_mat[1][ii] = cspace[ii].y * scale[ii]; t_mat[2][ii] = cspace[ii].z * scale[ii]; } return TRUE; } /* **************************************************************** * CLR_exit() * * Completes use of the CLR_ routine set and frees any * allocated space */ CLR_exit() { if (!init) return TRUE; (void)CLR_exit_samples(); if (X_tristim != NULL) free((char *)X_tristim); if (Y_tristim != NULL) free((char *)Y_tristim); if (Z_tristim != NULL) free((char *)Z_tristim); if (work_curve != NULL) free((char *)work_curve); X_tristim = Y_tristim = Z_tristim = work_curve = NULL; init = FALSE; return TRUE; } /* ************************************************************* */ use of the CLR_ routine set and frees any * allocated space */ CLR_exit() { if (!init) return TRUE; (void)CLR_exit_samples(); if (X_tristim != NULL) free((char *)X_tristim); if (Y_tristim != NULL) free((char *)Y_tristim); if (Z_tristim != NUclr_clip.c 660 401 24 12716 4623574150 6066 /* **************************************************************** * clr_clip.c * **************************************************************** * MODULE PURPOSE: * This module contains routines for bringing colors outside * the displayable gamut into the displayable gamut. * * MODULE CONTENTS: * CLR_clamp_rgb - clamps rgb values to 0.0-1.0 range * CLR_scale_rgb - scales rgb values to 0.0-1.0 range * CLR_clip_rgb - clips rgb values to 0.0-1.0 range */ #include #include #include "clr.h" /* define the max and min scaling values for equal intensity * clipping. These are defined to be slightly inside of the * 0 to 1 range to avoid numerical problems that occur when * colors are on or very close to the clipping boundary. */ #define MAX_CLIP 0.9999 /* 1 - MAX_CLIP > 1 res step */ #define MIN_CLIP 0.0001 /* less than 1 res step */ /* **************************************************************** * CLR_clamp_rgb (rgb) * CLR_RGB rgb (in) - the input rgb * * Returns clamped color. Values greater than 1 are clamped to 1, * values less than 0 are clamped to 0. */ CLR_RGB CLR_clamp_rgb(in_rgb) CLR_RGB in_rgb; { CLR_RGB rgb; rgb = in_rgb; if (rgb.r < 0.0) rgb.r = 0.0; else if (rgb.r > 1.0) rgb.r = 1.0; if (rgb.g < 0.0) rgb.g = 0.0; else if (rgb.g > 1.0) rgb.g = 1.0; if (rgb.b < 0.0) rgb.b = 0.0; else if (rgb.b > 1.0) rgb.b = 1.0; return rgb; } /* **************************************************************** * CLR_scale_rgb (rgb) * CLR_RGB rgb (in) - the input rgb * * Returns scaled color. Values less than 0 are clamped to zero. * If any values are greater than 1, all color components are * scaled so the offending value is equal to 1. */ CLR_RGB CLR_scale_rgb(in_rgb) CLR_RGB in_rgb; { double scale=1.0; double tmp_scale; CLR_RGB rgb; rgb = in_rgb; if (rgb.r < 0.0) rgb.r = 0.0; else if ((rgb.r > 1.0) && ((tmp_scale = (1.0 / rgb.r)) < scale)) scale = tmp_scale; if (rgb.g < 0.0) rgb.g = 0.0; else if ((rgb.g > 1.0) && ((tmp_scale = (1.0 / rgb.g)) < scale)) scale = tmp_scale; if (rgb.b < 0.0) rgb.b = 0.0; else if ((rgb.b > 1.0) && ((tmp_scale = (1.0 / rgb.b)) < scale)) scale = tmp_scale; if (scale < 1.0) { rgb.r *= scale; rgb.g *= scale; rgb.b *= scale; } return rgb; } /* **************************************************************** * CLR_RGB CLR_clip_rgb(in_rgb) * CLR_RGB in_rgb (in) - the input rgb * * Returns the clipped color - clipping is by desaturating the * color by shifting it towards the neutral axis in a plane * perpendicular to the neutral axis. The neutral axis is taken * as the vector from the monitor black to the monitor white. */ CLR_RGB CLR_clip_rgb(in_rgb) CLR_RGB in_rgb; { CLR_RGB diff, out_rgb; double axis_rgb, diff_mult, tmp_mult; /* check to see if clipping is required */ if ((in_rgb.r < 0.0) || (in_rgb.r > 1.0) || (in_rgb.g < 0.0) || (in_rgb.g > 1.0) || (in_rgb.b < 0.0) || (in_rgb.b > 1.0)) { /* clipping is required, determine the distance from * the origin to the equal intensity plane containing * the color. The distance is normalized the origin * at color (0,0,0) and a distance of 1 at (1,1,1) */ axis_rgb = (in_rgb.r + in_rgb.g + in_rgb.b) * .333333; /* check for the intensity plane of the color being * outside the displayable range -- if it is, set * color to either black or white. */ if (axis_rgb <= MIN_CLIP) /* this is not a visible color -- it should not * have been computed */ out_rgb.r = out_rgb.g = out_rgb.b = 0.0; else if (axis_rgb >= MAX_CLIP) /* This is way beyond white in intensity, set it * to the white point. */ out_rgb.r = out_rgb.g = out_rgb.b = 1.0; else { /* the intensity plane is within the displayable * range. Compute the vector from the neutral * axis to the color on it's intensity plane. * The intersection of the neutral axis and the * intensity plane is at r=g=b=axis_rgb. */ diff.r = in_rgb.r - axis_rgb; diff.g = in_rgb.g - axis_rgb; diff.b = in_rgb.b - axis_rgb; /* determine the relative length of the vector * to the edge of the displayable color gamut. */ diff_mult = 1.0; if (in_rgb.r > 1.0) { if ((tmp_mult = (MAX_CLIP - axis_rgb) / diff.r) < diff_mult) diff_mult = tmp_mult; } else if (in_rgb.r < 0.0) { if ((tmp_mult = (MIN_CLIP - axis_rgb) / diff.r) < diff_mult) diff_mult = tmp_mult; } if (in_rgb.g > 1.0) { if ((tmp_mult = (MAX_CLIP - axis_rgb) / diff.g) < diff_mult) diff_mult = tmp_mult; } else if (in_rgb.g < 0.0) { if ((tmp_mult = (MIN_CLIP - axis_rgb) / diff.g) < diff_mult) diff_mult = tmp_mult; } if (in_rgb.b > 1.0) { if ((tmp_mult = (MAX_CLIP - axis_rgb) / diff.b) < diff_mult) diff_mult = tmp_mult; } else if (in_rgb.b < 0.0) { if ((tmp_mult = (MIN_CLIP - axis_rgb) / diff.b) < diff_mult) diff_mult = tmp_mult; } /* determine the location of the color at the * edge of the displayable color gamut. */ out_rgb.r = axis_rgb + (diff.r * diff_mult); out_rgb.g = axis_rgb + (diff.g * diff_mult); out_rgb.b = axis_rgb + (diff.b * diff_mult); } return out_rgb; } else return in_rgb; } /* ************************************************************* */ ; } else if (in_rgb.b < 0.0) { if ((tmp_mulclr_sample.c 660 401 24 34202 4650402570 6406 /* **************************************************************** * clr_sample.c * **************************************************************** * MODULE PURPOSE: * This module contains the routines for spectral sampling. * The operations include defining the sampling space, * sampling spectral curves, and transformations from spectral * sample space to RGB or XYZ. * * MODULE CONTENTS: * CLR_init_samples - initialize spectral sampling * CLR_num_samples - returns the number of samples * CLR_spect_to_sample - sample a spectral curve * CLR_get_sample_rgb - get the sample to RGB matrix * CLR_get_sample_xyz - get the sample to XYZ matrix * CLR_reconstruct - reconstruct a spectral curve * CLR_exit_samples - finish with spectral sampling * * NOTES: * > The CLR_ routines must be initialized (CLR_init()) before * using any of the sampling routines. * * > When the CLR_SAMPLE_HALL method is selected, sampling * uses abutting, non-overlapping, box sample and * reconstruction functions are described in Section * 3.5.2 Spectral Sampling Approaches, and in Hall (1983). * This provides continuous sampling between the low bound * of the first sample and the high bound of the last * sample. For applications requiring discrete isolated * samples, user modification of these routines * is required. * * > When the CLR_SAMPLE_MEYER method is used, 4 samples are * used as described in Meyer (1988). * * > When using CLR_SAMPLE_LIN_RAMP method is used, sampling * uses box convolution filters and reconstruction uses * overlapping linear ramped triangles. */ #include #include "clr.h" static int sample_type = -1; static int samples = 0; static double *sample_func = NULL; static double *reconst_func = NULL; static double XYZ_to_ACC[3][3] = {{-0.0177, 1.0090, 0.0073}, {-1.5370, 1.0821, 0.3209}, { 0.1946, -0.2045, 0.5264}}; static double ACC_to_XYZ[3][3]; static double samp_to_ACC[3][4] = {{0.00000, 0.18892, 0.67493, 0.19253}, {0.00000, 0.31824, 0.00000, -0.46008}, {0.54640, 0.00000, 0.00000, 0.00000}}; /* These are the bounds for the Hall method optimized for the * Macbeth Colorchecker chart illuminated by D5500, D6500, D7500, * and standard illuminant C */ static int Hall_bounds[10][11] = { {544, 627, 000, 000, 000, 000, 000, 000, 000, 000, 000}, {407, 535, 652, 000, 000, 000, 000, 000, 000, 000, 000}, {414, 494, 580, 658, 000, 000, 000, 000, 000, 000, 000}, {414, 494, 562, 595, 656, 000, 000, 000, 000, 000, 000}, {418, 479, 511, 563, 595, 654, 000, 000, 000, 000, 000}, {418, 479, 510, 555, 581, 604, 655, 000, 000, 000, 000}, {419, 472, 494, 517, 556, 580, 604, 655, 000, 000, 000}, {419, 473, 495, 517, 552, 573, 591, 613, 659, 000, 000}, {419, 468, 486, 504, 522, 552, 571, 588, 609, 656, 000}, {415, 427, 473, 494, 515, 545, 564, 580, 595, 615, 659} }; extern char *malloc(); /* **************************************************************** * CLR_init_samples (method, num_samples, sample_bounds) * int method (in) - sampling method: * CLR_SAMPLE_MEYER Meyer (1988) * CLR_SAMPLE_HALL Hall (1983) * int num_samples (in) - number of sample functions * int sample_bounds[] (in) - boundaries of the sampling * functions. There must be * num_samples+1 bounds arranged in * ascending order in this array. If * a NULL pointer is passed in and * num_samples <= 10, then a default * sampling is used. * * For the CLR_SAMPLE_HALL method the bound wavelength is included * in the sampling function. For example, using 3 samples with * bounds at (411, 491, 571, 651), the actual samples are 411-490, * 491-570, and 571-650. * * The CLR_SAMPLE_MEYER method uses a prescribed sampling with 4 * samples. The num_samples and sample_bounds arguments are * ignored. * * Returns TRUE if successful, FALSE if sample bounds are not valid * or sampling is previously initialized to some other value. */ CLR_init_samples (method, num_samples, sample_bounds) int method; int num_samples; int *sample_bounds; { int ct; CLR_exit_samples (); if (method == CLR_SAMPLE_MEYER) { samples = 4; sample_type = method; CLR_t_inverse (XYZ_to_ACC, ACC_to_XYZ); } else if (method == CLR_SAMPLE_HALL) { /* build a set of sampling and reconstruction curves */ int min_wl, max_wl, cur_wl; double *cur_s, *cur_r, fill; if (num_samples <= 0) return FALSE; if (sample_bounds == NULL) sample_bounds = Hall_bounds[num_samples - 1]; if ((min_wl = CLR_get_min_wl()) < 0) return FALSE; max_wl = CLR_get_max_wl(); if ((sample_func = (double *)malloc((unsigned)(sizeof(double) * num_samples * (max_wl - min_wl + 1)))) == NULL) goto error; if ((reconst_func = (double *)malloc((unsigned)(sizeof(double) * num_samples * (max_wl - min_wl + 1)))) == NULL) goto error; cur_s = sample_func; cur_r = reconst_func; for (ct=0; ct 456) || (max_wl < 457)) *sample++ = 0.0; else *sample++ = spectral[456-min_wl] + (0.4 * (spectral[457-min_wl] - spectral[456-min_wl])); if ((min_wl > 490) || (max_wl < 491)) *sample++ = 0.0; else *sample++ = spectral[490-min_wl] + (0.9 * (spectral[491-min_wl] - spectral[490-min_wl])); if ((min_wl > 557) || (max_wl < 558)) *sample++ = 0.0; else *sample++ = spectral[557-min_wl] + (0.7 * (spectral[558-min_wl] - spectral[557-min_wl])); if ((min_wl > 631) || (max_wl < 632)) *sample++ = 0.0; else *sample++ = spectral[631-min_wl] + (0.4 * (spectral[632-min_wl] - spectral[631-min_wl])); } else if ((sample_type == CLR_SAMPLE_HALL) || (sample_type == CLR_SAMPLE_LIN_RAMP)) { double *cur_s, *spec; cur_s = sample_func; for (cur_samp=0; cur_samp=0; ) *sample += *spec++ * *cur_s++; } } else return FALSE; return TRUE; } /* **************************************************************** * CLR_get_sample_rgb (matrix) * double *matrix (mod) - matrix to be filled. * * Returns the matrix for conversion from the sampled space to the * RGB the CLR_ routines have been initialized for. The matrix * is a 3 x num_samples matrix. */ CLR_get_sample_rgb(matrix) double *matrix; { double *xyz = NULL; double xyz_rgb[3][3]; int ct; /* get the XYZ matrix, then transform it into RGB */ if ((xyz = (double *)malloc((unsigned)(3 * sizeof(double) * samples))) == NULL) goto error; if (!CLR_get_sample_xyz(xyz)) goto error; if (!CLR_get_xyz_rgb(xyz_rgb)) goto error; for (ct=0; ct=0; ) *spec++ = 0.0; cur_r = reconst_func; for (cur_samp=0; cur_samp=0; ) *spec++ += *sample * *cur_r++; } } else return FALSE; return TRUE; } /* **************************************************************** * CLR_exit_samples() * * Complete use of spectral sampling, free any allocated space. */ CLR_exit_samples () { sample_type = -1; samples = 0; if (sample_func != NULL) { free((char *)sample_func); sample_func = NULL; } if (reconst_func != NULL) { free((char *)reconst_func); reconst_func = NULL; } return TRUE; } ; --ct>=0; ) *spec++ += *sample * *cur_r++; } } else return FALSE; return TRUE; } /* **************************************************************** * CLR_exit_samples() * * Complete use of spectral sampling, free any allocated space. */ CLR_exit_samples () { sample_type = -1; samples = 0; if (sample_func != NULL) { free((char *)sample_func); sampillum_mod.h 660 401 24 7575 4623574153 6257 /* **************************************************************** * illum_mod.h * **************************************************************** * include file for the illumination models */ #ifndef ILLUM_H #define ILLUM_H #include "geo.h" /* The material structure maintains the description of a material * interface. An interface between a conductor or dielectric and * air is characterized by loading the properties of the material * and an index of refraction of 1 for the outside material. An * interface between two materials is characterized by generating * a pseudo-material as described in appendix III.5.1, Approxima- * tion Methodology. * * In filling the material structure, the reflected direction is * the 'outside' of the material. That is, for an interface * between air and glass, for example, the the reflected direction * or 'outside' is air (Ar_spectral = NULL, nr=1.0), and the * transmitted direction is glass (nt=1.5, etc.) * * Individual spectral components of the material are characterized * by the sampled spectral values and a multiplier to scale these * values (as discussed in appendix II.1.5, Selecting Ka, Ka, Ks, * Ft and Fr) */ typedef struct { double Ka_scale; /* ambient multiplier */ double *Ka_spectral; /* sampled ambient spectral curve */ double Kd_scale; /* diffuse multiplier */ double *Kd_spectral; /* sampled diffuse spectral curve */ double Ks_scale; /* specular multiplier */ double *Ks_spectral; /* sampled specular spectral curve */ double Kt_scale; /* transmitted multiplier */ double *Kt_spectral; /* sampled specular transmitted */ /* curve. The Kt_ properties are */ /* used for the Whitted model only */ double *Ar_spectral; /* sampled filter spectral curve. */ double *At_spectral; /* sampled filter spectral curve */ /* These are used in the Hall model */ /* only. Filter attenuation is */ /* ignored if a NULL pointer is */ /* specified. */ double beta; /* roughness indicator */ double (*D)(); /* microfacet distribution */ double D_coeff; /* the coefficient for the */ /* microfacet distribution function */ /* computed from by microfacet */ /* distribution init function */ double (*G)(); /* geometric attenuation function */ double G_coeff; /* the coefficient for the geometric */ /* attenuation function (m_2 for */ /* the Sancer function, unused for */ /* the Torrance-Sparrow function) */ double Ro; /* average reflectance */ double nr; /* index of refraction (incident) */ double nt; /* index of refraction (transmitted) */ double k; /* absorption coefficient */ int conductor; /* flags the specular material as a */ /* conductor */ int transparent; /* flags whether the material is */ /* transparent --- note, composite */ /* materials have a dielectric */ /* specular component but may not */ /* be transparent */ int r_is_air; /* flags that the 'outside' or */ /* reflect side of the interface is */ /* air */ } ILLUM_MATL; typedef struct { double I_scale; /* illumination multiplier */ double *I_spectral; /* sampled source spectral curve */ POINT3 center; /* center of the light source */ } ILLUM_LGT; int IM_init(); double *IM_bouknight(); double *IM_phong(); double *IM_blinn(); double *IM_whitted(); double *IM_hall(); int IM_exit(); #endif CLR_H /* ************************************************************* */ is */ /* air */ } ILLUM_MATL; typedef struct { double I_scale; /* illumination muillum_mod.c 660 401 24 55746 4623574123 6272 /* **************************************************************** * illum_mod.c * **************************************************************** * MODULE PURPOSE: * This module contains the source code for the illumination * models discussed for in Section 4.2, Incremental Shading, * Empirical Models, and in Section 4.3, Ray Tracing, * Transitional Models. * * MODULE CONTENTS: * IM_init - initialize the illumination models * IM_bouknight - Bouknight (1970) illumination model * IM_phong - Phong (1975) illumination model * IM_blinn - Blinn (1976) illumination model * IM_whitted - Whitted (1980) illumination model * IM_hall - Hall (1983) illumination model * IM_exit - finish with the illumination models * * NOTES: * > The illumination model is called once a * point on a surface has been identified and the list * of illuminating light sources has been generated. * * > There exists a routine that the illumination models can * call to get a color from any direction. Specifically * this is used for inquiring about the reflected or * transmitted directions in the ray tracing models. * This routine is passed the view vector for which the * color is required. * * > a common calling convention is used for ease in * interchanging the illumination model routines. Each * routine is passed the location and orientation of the * surface, a view vector, an interface material, a * description of the lighting, and an array to receive * the computed color. * * The orientation of the surface is given by the surface * normal which is ALWAYS directed to the 'outside' or * reflected side of the material. * * The view vector is specified by start position and * direction vector. During visibility computations the * view vector is typically directed towards the surface. * The direction cosines MUST be negated prior to calling * the illumination model for consistency with the vector * conventions used. * * See 'illum_mod.h' for material details. * * The light vector is a list of pointers to ILLUM_LGT * structures terminated by a NULL pointer. The first * entry is taken as the ambient illumination. Only * light that is incident from the air side of a material * can provide illumination. * * > These models assume that the material structure is * correctly loaded and that the surface is facing the * viewer (N.V > 0) for illumination models that do not * consider transparency. */ #include #include #include "illum_mod.h" #include "F.h" static int num_samples = 0; static int (*get_color)() = NULL; static double *Fr_buffer = NULL; static double *Ft_buffer = NULL; /* **************************************************************** * int IM_init (num_clr_samples, get_clr_routine) * int num_clr_samples (in) - number of color samples * int (*get_clr_routine)() (in) - routine to call to get * the color from some direction * * Initializes the illumination model routine set, returns TRUE * if successful and FALSE upon failure. */ IM_init (num_clr_samples, get_clr_routine) int num_clr_samples, (*get_clr_routine)(); { char *malloc(); if (((num_samples = num_clr_samples) <= 0) || ((Fr_buffer = (double *)malloc((unsigned)(num_samples * sizeof(double)))) == NULL) || ((Ft_buffer = (double *)malloc((unsigned)(num_samples * sizeof(double)))) == NULL)) { (void)IM_exit(); return FALSE; } get_color = get_clr_routine; return TRUE; } /* **************************************************************** * double *IM_bouknight (surface, V, matl, lgts, color) * LINE *surface (in) - surface for color computation * LINE *V (in) - view vector * ILLUM_MATL *matl (in) - material properties * ILLUM_LGT **lgts (in) - illuminating sources * double *color (mod) - array to receive the color * * Evaluates the color using the Bouknight (1970) illumination * model as described in Eq.(\*(sE). Returns 'color' upon * success and NULL upon failure. */ double *IM_bouknight (surface, V, matl, lgts, color) LINE *surface, *V; ILLUM_MATL *matl; ILLUM_LGT **lgts; double *color; { int ct; double N_dot_L; LINE L; /* load the ambient illumination */ for (ct=0; ctKa_scale * matl->Ka_spectral[ct] * (*lgts)->I_scale * (*lgts)->I_spectral[ct]; } lgts++; /* load the diffuse component of the illumination. Loop * through the lights and compute (N.L). If it is positive * then the surface is illuminated. */ while (*lgts != NULL) { if ((geo_line(&(surface->start),&((*lgts)->center),&L) > 0) && ((N_dot_L=geo_dot(&(surface->dir),&(L.dir))) > 0)) { /* The surface is illuminated by this light. Loop through * the color samples and sum the diffuse contribution for * this light to the the color. */ for (ct=0; ctKd_scale * matl->Kd_spectral[ct] * N_dot_L * (*lgts)->I_scale * (*lgts)->I_spectral[ct]; } } lgts++; } return color; } /* **************************************************************** * double *IM_phong (surface, V, matl, lgts, color) * LINE *surface (in) - surface for color computation * LINE *V (in) - view vector * ILLUM_MATL *matl (in) - material properties * ILLUM_LGT **lgts (in) - illuminating sources * double *color (mod) - array to receive the color * * Evaluates the color using the Phong (1975) illumination * model as described in Eq.(\*(rU). Returns 'color' upon * success and NULL upon failure. * * The actual Phong model results when the microfacet distribution * D_phong() is used, and matl->Ks_spectral is the identity * material. * * Using the microfacet distribution D_blinn() gives Blinn's * interpretation of the Phong model per Eq.(\*(rI). */ double *IM_phong (surface, V, matl, lgts, color) LINE *surface, *V; ILLUM_MATL *matl; ILLUM_LGT **lgts; double *color; { int ct; double N_dot_L, D; LINE L; /* load the ambient illumination */ for (ct=0; ctKa_scale * matl->Ka_spectral[ct] * (*lgts)->I_scale * (*lgts)->I_spectral[ct]; } lgts++; /* load the diffuse and specular illumination components. * Loop through the lights and compute (N.L). If it is * positive then the surface is illuminated. */ while (*lgts != NULL) { if ((geo_line(&(surface->start),&((*lgts)->center),&L) > 0) && ((N_dot_L=geo_dot(&(surface->dir),&(L.dir))) > 0)) { /* The surface is illuminated. Compute the microfacet * distribution function. */ D = (*(matl->D))(&(surface->dir), &(L.dir), &(V->dir), matl->D_coeff); /* Loop through the color samples and sum the diffuse * and specular contribution for this light to the the * color. */ for (ct=0; ctKd_scale * matl->Kd_spectral[ct]) + (D * matl->Ks_scale * matl->Ks_spectral[ct])) * (*lgts)->I_scale * (*lgts)->I_spectral[ct]; } } lgts++; } return color; } /* **************************************************************** * double *IM_blinn (surface, V, matl, lgts, color) * LINE *surface (in) - surface for color computation * LINE *V (in) - view vector * ILLUM_MATL *matl (in) - material properties * ILLUM_LGT **lgts (in) - illuminating sources * double *color (mod) - array to receive the color * * Evaluates the color using the Blinn (1977) illumination * model as described in Eq.(\*(rV). Returns 'color' upon * success and NULL upon failure. The microfacet distribution * functions D_blinn(), D_gaussian(), and D_reitz are the three * functions presented by Blinn. The geometric attenuation * function G_torrance() is the attenuation function used by Blinn. * If matl->G is NULL then the geometric attenuation is omitted. * The Fresnel reflectance is approximated using the Cook (1983) * technique. */ double *IM_blinn (surface, V, matl, lgts, color) LINE *surface, *V; ILLUM_MATL *matl; ILLUM_LGT **lgts; double *color; { int ct; double N_dot_L, N_dot_V, D, G, *Fr; LINE L; DIR_VECT *T, *H; /* load the ambient illumination */ for (ct=0; ctKa_scale * matl->Ka_spectral[ct] * (*lgts)->I_scale * (*lgts)->I_spectral[ct]; } lgts++; /* load the diffuse and specular illumination components. * Loop through the lights and compute (N.L). If it is * positive then the surface is illuminated. */ while (*lgts != NULL) { if ((geo_line(&(surface->start),&((*lgts)->center),&L) > 0) && ((N_dot_L=geo_dot(&(surface->dir),&(L.dir))) > 0)) { /* The surface is illuminated. Compute the microfacet * distribution, geometric attenuation, Fresnel * reflectance and (N.V) for the specular function. */ D = (*(matl->D))(&(surface->dir), &(L.dir), &(V->dir), matl->D_coeff); if (matl->G == NULL) G = 1.0; else G = (*(matl->G))(&(surface->dir), &(L.dir), &(V->dir), matl->G_coeff); H = geo_H(&(L.dir), &(V->dir)); if (matl->conductor) { Fr = F_approx_Fr(H, &(L.dir), matl->Ro, matl->nt, matl->k, num_samples, matl->Ks_spectral, Fr_buffer); } else { T = geo_rfr(&(V->dir),H, matl->nr, matl->nt); Fr = F_approx_Fr_Ft(H, &(L.dir), T, matl->Ro, matl->nr, matl->nt, num_samples, matl->Ks_spectral, Fr_buffer, Ft_buffer); } /* Loop through the color samples and sum the diffuse * and specular contribution for this light to the the * color. Note the threshold on N_dot_V to prevent * divide by zero at grazing. */ if ((N_dot_V=geo_dot(&(surface->dir),&(V->dir))) > 0.0001){ for (ct=0; ctKd_scale * matl->Kd_spectral[ct]) + ((D * G * matl->Ks_scale * Fr[ct]) / N_dot_V)) * (*lgts)->I_scale * (*lgts)->I_spectral[ct]; } } } lgts++; } return color; } /* **************************************************************** * double *IM_whitted (surface, V, matl, lgts, color) * LINE *surface (in) - surface for color computation * LINE *V (in) - view vector * ILLUM_MATL *matl (in) - material properties * ILLUM_LGT *lgts (in) - illuminating sources * double *color (mod) - array to receive the color * * Evaluates the color using the Whitted (1980) illumination * model as described in Eq.(\*(e6). Returns 'color' upon * success and NULL upon failure. * * The actual Whitted model results when the microfacet * distribution D_blinn() is used, and when both matl->Ks_spectral * and matl->Kt_spectral are the identity material. * * The matl->Kt_scale and matl->Kt_spectral are required for this * illumination model only. */ double *IM_whitted (surface, V, matl, lgts, color) LINE *surface, *V; ILLUM_MATL *matl; ILLUM_LGT **lgts; double *color; { int inside = FALSE; int ct; double N_dot_L, D; LINE L; /* figure out whether we are on the reflected or transmitted * side of the material interface (outside or inside). If * we are inside a material, then there is no illumination * from lights and such - skip directly to reflection and * refraction contribution. */ if ((geo_dot(&(surface->dir),&(V->dir)) < 0) || (!matl->r_is_air)) { for (ct=0; ctr_is_air) goto rfl_rfr; /* load the ambient illumination */ for (ct=0; ctKa_scale * matl->Ka_spectral[ct] * (*lgts)->I_scale * (*lgts)->I_spectral[ct]; } lgts++; /* load the diffuse and specular illumination components. * Loop through the lights and compute (N.L). If it is * positive then the surface is illuminated. */ while (*lgts != NULL) { if ((geo_line(&(surface->start),&((*lgts)->center),&L) > 0) && ((N_dot_L=geo_dot(&(surface->dir),&(L.dir))) > 0)) { /* The surface is illuminated. Compute the microfacet * distribution function. */ D = (*(matl->D))(&(surface->dir), &(L.dir), &(V->dir), matl->D_coeff); /* Loop through the color samples and sum the diffuse * and specular contribution for this light to the the * color. */ for (ct=0; ctKd_scale * matl->Kd_spectral[ct]) + (D * matl->Ks_scale * matl->Ks_spectral[ct])) * (*lgts)->I_scale * (*lgts)->I_spectral[ct]; } } lgts++; } /* compute the contribution from the reflection and * refraction directions. Get a buffer to hold the * computed colors, then process the reflected direction * and the refracted direction. */ rfl_rfr: if (get_color != NULL) { char *malloc(); double *Ir_or_It; LINE V_new; DIR_VECT *T; if ((Ir_or_It = (double *)malloc((unsigned)(num_samples * sizeof(double)))) == NULL) return color; /* get the reflected vector then ask for the color from * the reflected direction. If there is a color, then * sum it with the current color */ V_new.start = surface->start; V_new.dir = *(geo_rfl(&(V->dir),&(surface->dir))); if ((*get_color)(&V_new, Ir_or_It) != NULL) for (ct=0; ctKs_scale * matl->Ks_spectral[ct]; /* if the material is transparent then get the refracted * vector and ask for the color from the refracted * direction. If there is a color, then sum it with * the current color */ if (matl->transparent) { V_new.start = surface->start; if (inside) T = geo_rfr(&(V->dir),&(surface->dir), matl->nt, matl->nr); else T = geo_rfr(&(V->dir),&(surface->dir), matl->nr, matl->nt); if (T != NULL) { V_new.dir = *T; if ((*get_color)(&V_new, Ir_or_It) != NULL) for (ct=0; ctKt_scale * matl->Kt_spectral[ct]; } } (void)free((char *)Ir_or_It); } return color; } /* **************************************************************** * double *IM_hall (surface, V, matl, lgts, color) * LINE *surface (in) - surface for color computation * LINE *V (in) - view vector * ILLUM_MATL *matl (in) - material properties * ILLUM_LGT *lgts (in) - illuminating sources * double *color (mod) - array to receive the color * * Evaluates the color using the Hall (1983) illumination * model as described in Eq.(\*(sF). Returns 'color' upon * success and NULL upon failure. * * The actual Hall model results when the microfacet * distribution D_blinn() is used, and matl->D = NULL. * * The transmittance is computed from the reflectance, so * matl->Kt_scale and matl->Kt_spectral are not used in the model. */ double *IM_hall (surface, V, matl, lgts, color) LINE *surface, *V; ILLUM_MATL *matl; ILLUM_LGT **lgts; double *color; { int inside = FALSE; int ct; double N_dot_L, D, G, *Fr; LINE L; DIR_VECT *T, *H, *Ht, pseudo_V; /* figure out whether we are on the reflected or transmitted * side of the material interface (outside or inside). If * we are inside a material, then there may be illumination * from outside. */ if (geo_dot(&(surface->dir),&(V->dir)) < 0) { for (ct=0; ctr_is_air) goto rfl_rfr; /* load the ambient illumination if we are not inside * inside the material */ if (!inside) { for (ct=0; ctKa_scale * matl->Ka_spectral[ct] * (*lgts)->I_scale * (*lgts)->I_spectral[ct]; } } lgts++; /* load the diffuse and specular illumination components. * Loop through the lights and compute (N.L). If it is * positive then the surface is illuminated. If it is * negative, then there may be transmitted illumination. */ while (*lgts != NULL) { if ((geo_line(&(surface->start),&((*lgts)->center),&L) > 0) && ((N_dot_L=geo_dot(&(surface->dir),&(L.dir))) > 0) && !inside) { /* The surface is illuminated. Compute the microfacet * distribution, geometric attenuation, Fresnel * reflectance and (N.V) for the specular function. */ D = (*(matl->D))(&(surface->dir), &(L.dir), &(V->dir), matl->D_coeff); if (matl->G == NULL) G = 1.0; else G = (*(matl->G))(&(surface->dir), &(L.dir), &(V->dir), matl->G_coeff); H = geo_H(&(L.dir), &(V->dir)); if (matl->conductor) { Fr = F_approx_Fr(H, &(L.dir), matl->Ro, matl->nt, matl->k, num_samples, matl->Ks_spectral, Fr_buffer); } else { T = geo_rfr(&(L.dir), H, matl->nr, matl->nt); Fr = F_approx_Fr_Ft(H, &(L.dir), T, matl->Ro, matl->nr, matl->nt, num_samples, matl->Ks_spectral, Fr_buffer, Ft_buffer); } /* Loop through the color samples and sum the diffuse * and specular contribution for this light to the the * color. */ for (ct=0; ctKd_scale * matl->Kd_spectral[ct]) + (D * G * matl->Ks_scale * Fr[ct])) * (*lgts)->I_scale * (*lgts)->I_spectral[ct]; } } else if ((N_dot_L > 0) && inside) { /* We are inside and the light is outside. Compute * the transmitted contribution from the light */ if ((Ht = geo_Ht(&(L.dir),&(V->dir),matl->nr, matl->nt)) != NULL) { /* The microfacet distribution functions could * only be equated when cast in terms of the primary * vectors L, V, and N. A pseudo_V vector is required * so that any of the distribution functions can be * applied. Ht is the vector bisector between of the * angle between L and pseudo_V, thus pseudo_V can be * computed by reflecting L about Ht. Refer to the * text for details. */ pseudo_V = *(geo_rfl(&(L.dir), Ht)); D = (*(matl->D))(&(surface->dir), &(L.dir), &pseudo_V, matl->D_coeff); Fr = F_approx_Fr_Ft(Ht, &(L.dir), &(V->dir), matl->Ro, matl->nr, matl->nt, num_samples, matl->Ks_spectral, Fr_buffer, Ft_buffer); if (matl->G == NULL) G = 1.0; else { /* To include the geometric attenuation, the view * vector direction must be reversed so that it * is to the same side of the surface as the * normal, see text for details. */ pseudo_V.i = -V->dir.i; pseudo_V.j = -V->dir.j; pseudo_V.k = -V->dir.k; G = (*(matl->G))(&(surface->dir), &(L.dir), &pseudo_V, matl->G_coeff); } for (ct=0; ctKs_scale * Ft_buffer[ct]) * (*lgts)->I_scale * (*lgts)->I_spectral[ct]; } } } lgts++; } /* compute the contribution from the reflection and * refraction directions. Get a buffer to hold the * computed colors, then process the reflected direction * and the refracted direction. */ rfl_rfr: if (get_color != NULL) { char *malloc(); double *Ir_or_It; LINE V_new; DIR_VECT pseudo_N; if ((Ir_or_It = (double *)malloc((unsigned)(num_samples * sizeof(double)))) == NULL) return color; /* Determine the Fresnel reflectance and transmittance. * If we are inside the material, then a pseudo normal * is required that faces to the same side of the * interface as the view vector. */ if (matl->conductor) { Fr = F_approx_Fr(&(surface->dir), &(V->dir), matl->Ro, matl->nt, matl->k, num_samples, matl->Ks_spectral, Fr_buffer); } else if (inside) { T = geo_rfr(&(V->dir),&(surface->dir), matl->nt, matl->nr); pseudo_N.i = -surface->dir.i; pseudo_N.j = -surface->dir.j; pseudo_N.k = -surface->dir.k; Fr = F_approx_Fr_Ft(&pseudo_N, &(V->dir), T, matl->Ro, matl->nt, matl->nr, num_samples, matl->Ks_spectral, Fr_buffer, Ft_buffer); } else { T = geo_rfr(&(V->dir),&(surface->dir), matl->nr, matl->nt); Fr = F_approx_Fr_Ft(&(surface->dir), &(V->dir), T, matl->Ro, matl->nr, matl->nt, num_samples, matl->Ks_spectral, Fr_buffer, Ft_buffer); } /* get the reflected vector then ask for the color from * the reflected direction. If there is a color, then * sum it with the current color */ V_new.start = surface->start; V_new.dir = *(geo_rfl(&(V->dir),&(surface->dir))); if ((*get_color)(&V_new, Ir_or_It) != NULL) { for (ct=0; ctKs_scale * Fr[ct]; } /* if the material is transparent then get the refracted * vector and ask for the color from the refracted * direction. If there is a color, then sum it with * the current color */ if (matl->transparent && (T != NULL)) { V_new.start = surface->start; V_new.dir = *T; if ((*get_color)(&V_new, Ir_or_It) != NULL) for (ct=0; ctKt_scale * Ft_buffer[ct]; } (void)free((char *)Ir_or_It); } /* If we are inside a material that has a filter attenuation * then apply the attenuation to the color. */ if ( ((!inside) && ((Fr = matl->Ar_spectral) != NULL)) || ((inside) && ((Fr = matl->At_spectral) != NULL)) ) { double dist; dist = sqrt ( ((surface->start.x - V->start.x) * (surface->start.x - V->start.x)) + ((surface->start.y - V->start.y) * (surface->start.y - V->start.y)) + ((surface->start.z - V->start.z) * (surface->start.z - V->start.z)) ); for (ct=0; ct