Import Upstream version 3a1upstream/3a1

1 files changed, 536 insertions, 0 deletions

diff --git a/colorspc.scm b/colorspc.scm
new file mode 100644
index 0000000..3a88767
--- /dev/null
+++ b/colorspc.scm
@@ -0,0 +1,536 @@
+;;; "colorspc.scm" color-space conversions
+;Copyright 2001, 2002 Aubrey Jaffer
+;
+;Permission to copy this software, to modify it, to redistribute it,
+;to distribute modified versions, and to use it for any purpose is
+;granted, subject to the following restrictions and understandings.
+;
+;1.  Any copy made of this software must include this copyright notice
+;in full.
+;
+;2.  I have made no warranty or representation that the operation of
+;this software will be error-free, and I am under no obligation to
+;provide any services, by way of maintenance, update, or otherwise.
+;
+;3.  In conjunction with products arising from the use of this
+;material, there shall be no use of my name in any advertising,
+;promotional, or sales literature without prior written consent in
+;each case.
+
+(require 'logical)
+(require 'multiarg/and-)
+(require-if 'compiling 'sort)
+(require-if 'compiling 'ciexyz)
+;@
+(define (color:linear-transform matrix row)
+  (map (lambda (mrow) (apply + (map * mrow row)))
+       matrix))
+
+(define RGB709:into-matrix
+  '((  3.240479 -1.537150 -0.498535 )
+    ( -0.969256  1.875992  0.041556 )
+    (  0.055648 -0.204043  1.057311 )))
+
+;;; http://www.pima.net/standards/it10/PIMA7667/PIMA7667-2001.PDF gives
+;;; matrix identical to sRGB:from-matrix, but colors drift under
+;;; repeated conversions to and from CIEXYZ.  Instead use RGB709.
+
+(define RGB709:from-matrix
+  '((  0.412453  0.357580  0.180423 )
+    (  0.212671  0.715160  0.072169 )
+    (  0.019334  0.119193  0.950227 )))
+
+;; From http://www.cs.rit.edu/~ncs/color/t_convert.html
+;@
+(define (CIEXYZ->RGB709 XYZ)
+  (color:linear-transform RGB709:into-matrix XYZ))
+(define (RGB709->CIEXYZ rgb)
+  (color:linear-transform RGB709:from-matrix rgb))
+
+;;; From http://www.w3.org/Graphics/Color/sRGB.html
+
+(define sRGB-log
+  (lambda (sv)
+    (if (<= sv 0.00304)
+	(* 12.92 sv)
+	(+ -0.055 (* 1.055 (expt sv 10/24))))))
+(define sRGB-exp
+  (lambda (x)
+    (if (<= x 0.03928)
+	(/ x 12.92)
+	(expt (/ (+ 0.055 x) 1.055) 2.4))))
+
+;; Clipping as recommended by sRGB spec.
+;@
+(define (CIEXYZ->sRGB XYZ)
+  (map (lambda (sv)
+	 (inexact->exact (round (* 255 (sRGB-log (max 0 (min 1 sv)))))))
+       (color:linear-transform RGB709:into-matrix XYZ)))
+(define (sRGB->CIEXYZ sRGB)
+  (color:linear-transform
+   RGB709:from-matrix
+   (map sRGB-exp
+	(map (lambda (b8v) (/ b8v 255.0)) sRGB))))
+
+;;; sRGB values are sometimes written as 24-bit integers 0xRRGGBB
+;@
+(define (xRGB->sRGB xRGB)
+  (list (ash xRGB -16)
+	(logand (ash xRGB -8) 255)
+	(logand xRGB 255)))
+(define (sRGB->xRGB sRGB)
+  (apply + (map * sRGB '(#x10000 #x100 #x1))))
+;@
+(define (xRGB->CIEXYZ xRGB) (sRGB->CIEXYZ (xRGB->sRGB xRGB)))
+(define (CIEXYZ->xRGB xyz)   (sRGB->xRGB (CIEXYZ->sRGB xyz)))
+
+;;;  http://www.pima.net/standards/it10/PIMA7667/PIMA7667-2001.PDF
+;;;	    Photography ­ Electronic still picture imaging ­
+;;;		 Extended sRGB color encoding ­ e-sRGB
+
+(define e-sRGB-log
+  (lambda (sv)
+    (cond ((< sv -0.0031308)
+	   (- 0.055 (* 1.055 (expt (- sv) 10/24))))
+	  ((<= sv 0.0031308)
+	   (* 12.92 sv))
+	  (else (+ -0.055 (* 1.055 (expt sv 10/24)))))))
+(define e-sRGB-exp
+  (lambda (x)
+    (cond ((< x -0.04045)
+	   (- (expt (/ (- 0.055 x) 1.055) 2.4)))
+	  ((<= x 0.04045)
+	   (/ x 12.92))
+	  (else (expt (/ (+ 0.055 x) 1.055) 2.4)))))
+;@
+(define (CIEXYZ->e-sRGB n XYZ)
+  (define two^n-9 (ash 1 (- n 9)))
+  (define offset (* 3 (ash 1 (- n 3))))
+  (map (lambda (x)
+	 (+ (inexact->exact (round (* x 255 two^n-9))) offset))
+       (map e-sRGB-log
+	    (color:linear-transform
+	     RGB709:into-matrix
+	     XYZ))))
+;@
+(define (e-sRGB->CIEXYZ n rgb)
+  (define two^n-9 (ash 1 (- n 9)))
+  (define offset (* 3 (ash 1 (- n 3))))
+  (color:linear-transform
+   RGB709:from-matrix
+   (map e-sRGB-exp
+	(map (lambda (b8v) (/ (- b8v offset) 255.0 two^n-9))
+	     rgb))))
+;@
+(define (sRGB->e-sRGB n sRGB)
+  (define two^n-9 (ash 1 (- n 9)))
+  (define offset (* 3 (ash 1 (- n 3))))
+  (map (lambda (x) (+ offset (* two^n-9 x))) sRGB))
+;@
+(define (e-sRGB->sRGB n rgb)
+  (define two^n-9 (ash 1 (- n 9)))
+  (define offset (* 3 (ash 1 (- n 3))))
+  (map (lambda (x) (/ (- x offset) two^n-9)) rgb))
+;@
+(define (e-sRGB->e-sRGB n rgb m)
+  (define shft (- m n))
+  (cond ((zero? shft) rgb)
+	(else (map (lambda (x) (ash x shft)) rgb))))
+
+;;; From http://www.cs.rit.edu/~ncs/color/t_convert.html
+
+;;; CIE 1976 L*a*b* is based directly on CIE XYZ and is an attampt to
+;;; linearize the perceptibility of color differences. The non-linear
+;;; relations for L*, a*, and b* are intended to mimic the logarithmic
+;;; response of the eye. Coloring information is referred to the color
+;;; of the white point of the system, subscript n.
+
+;;;; L* is CIE lightness
+;;;     L* = 116 * (Y/Yn)^1/3 - 16    for Y/Yn > 0.008856
+;;;     L* = 903.3 * Y/Yn             otherwise
+
+(define (CIE:Y/Yn->L* Y/Yn)
+  (if (> Y/Yn 0.008856)
+      (+ -16 (* 116 (expt Y/Yn 1/3)))
+      (* 903.3 Y/Yn)))
+(define (CIE:L*->Y/Yn L*)
+  (cond ((<= L* (* 903.3 0.008856))
+	 (/ L* 903.3))
+	((<= L* 100.)
+	 (expt (/ (+ L* 16) 116) 3))
+	(else 1)))
+
+;;; a* = 500 * ( f(X/Xn) - f(Y/Yn) )
+;;; b* = 200 * ( f(Y/Yn) - f(Z/Zn) )
+;;;     where f(t) = t^1/3              for t > 0.008856
+;;;           f(t) = 7.787 * t + 16/116 otherwise
+
+(define (ab-log t)
+  (if (> t 0.008856)
+      (expt t 1/3)
+      (+ 16/116 (* t 7.787))))
+(define (ab-exp f)
+  (define f3 (expt f 3))
+  (if (> f3 0.008856)
+      f3
+      (/ (- f 16/116) 7.787)))
+;@
+(define (CIEXYZ->L*a*b* XYZ . white-point)
+  (apply (lambda (X/Xn Y/Yn Z/Zn)
+	   (list (CIE:Y/Yn->L* Y/Yn)
+		 (* 500 (- (ab-log X/Xn) (ab-log Y/Yn)))
+		 (* 200 (- (ab-log Y/Yn) (ab-log Z/Zn)))))
+	 (map / XYZ (if (null? white-point)
+			CIEXYZ:D65
+			(car white-point)))))
+
+;;; Here Xn, Yn and Zn are the tristimulus values of the reference white.
+;@
+(define (L*a*b*->CIEXYZ L*a*b* . white-point)
+  (apply (lambda (Xn Yn Zn)
+	   (apply (lambda (L* a* b*)
+		    (let* ((Y/Yn (CIE:L*->Y/Yn L*))
+			   (fY/Yn (ab-log Y/Yn)))
+		      (list (* Xn (ab-exp (+ fY/Yn (/ a* 500))))
+			    (* Yn Y/Yn)
+			    (* Zn (ab-exp (+ fY/Yn (/ b* -200)))))))
+		  L*a*b*))
+	 (if (null? white-point)
+	     CIEXYZ:D65
+	     (car white-point))))
+
+;;; XYZ to CIELUV
+
+;;; CIE 1976 L*u*u* (CIELUV) is based directly on CIE XYZ and is another
+;;; attampt to linearize the perceptibility of color differences.  L* is
+;;; CIE lightness as for L*a*b* above.  The non-linear relations for u*
+;;; and v* are:
+
+;;;     u* =  13 L* ( u' - un' )
+;;;     v* =  13 L* ( v' - vn' )
+
+;;; The quantities un' and vn' refer to the reference white or the light
+;;; source; for the 2° observer and illuminant C, un' = 0.2009, vn' =
+;;; 0.4610. Equations for u' and v' are given below:
+
+;;;     u' = 4 X / (X + 15 Y + 3 Z)
+;;;     v' = 9 Y / (X + 15 Y + 3 Z)
+
+(define (XYZ->uv XYZ)
+  (apply (lambda (X Y Z)
+	   (define denom (+ X (* 15 Y) (* 3 Z)))
+	   (if (zero? denom)
+	       '(4. 9.)
+	       (list (/ (* 4 X) denom)
+		     (/ (* 9 Y) denom))))
+	 XYZ))
+;@
+(define (CIEXYZ->L*u*v* XYZ . white-point)
+  (set! white-point (if (null? white-point)
+			CIEXYZ:D65
+			(car white-point)))
+  (let* ((Y/Yn (/ (cadr XYZ) (cadr white-point)))
+	 (L* (CIE:Y/Yn->L* Y/Yn)))
+    (cons L* (map (lambda (q) (* 13 L* q))
+		  (map - (XYZ->uv XYZ) (XYZ->uv white-point))))))
+
+;;; CIELUV to XYZ
+
+;;; The transformation from CIELUV to XYZ is performed as following:
+
+;;;     u' = u / ( 13 L* ) + un
+;;;     v' = v / ( 13 L* ) + vn
+;;;     X = 9 Y u' / 4 v'
+;;;     Z = ( 12 Y - 3 Y u' - 20 Y v' ) / 4 v'
+;@
+(define (L*u*v*->CIEXYZ L*u*v* . white-point)
+  (set! white-point (if (null? white-point)
+			CIEXYZ:D65
+			(car white-point)))
+  (apply (lambda (un vn)
+	   (apply (lambda (L* u* v*)
+		    (if (not (positive? L*))
+			'(0. 0. 0.)
+			(let* ((up (+ (/ u* 13 L*) un))
+			       (vp (+ (/ v* 13 L*) vn))
+			       (Y (* (CIE:L*->Y/Yn L*) (cadr white-point))))
+			  (list (/ (* 9 Y up) 4 vp)
+				Y
+				(/ (* Y (+ 12 (* -3 up) (* -20 vp))) 4 vp)))))
+		  L*u*v*))
+	 (XYZ->uv white-point)))
+
+;;; http://www.inforamp.net/~poynton/PDFs/coloureq.pdf
+
+(define pi (* 4 (atan 1)))
+(define pi/180 (/ pi 180))
+;@
+(define (L*a*b*->L*C*h lab)
+  (define h (/ (atan (caddr lab) (cadr lab)) pi/180))
+  (list (car lab)
+	(sqrt (apply + (map * (cdr lab) (cdr lab))))
+	(if (negative? h) (+ 360 h) h)))
+;@
+(define (L*C*h->L*a*b* lch)
+  (apply (lambda (L* C* h)
+	   (set! h (* h pi/180))
+	   (list L*
+		 (* C* (cos h))
+		 (* C* (sin h))))
+	 lch))
+;@
+(define (L*a*b*:DE* lab1 lab2)
+  (sqrt (apply + (map (lambda (x) (* x x)) (map - lab1 lab2)))))
+
+;;; http://www.colorpro.com/info/data/cie94.html
+
+(define (color:process-params parametric-factors)
+  (define ans
+    (case (length parametric-factors)
+      ((0) #f)
+      ((1) (if (list? parametric-factors)
+	       (apply color:process-params parametric-factors)
+	       (append parametric-factors '(1 1))))
+      ((2) (append parametric-factors '(1)))
+      ((3) parametric-factors)
+      (else (slib:error 'parametric-factors 'too-many parametric-factors))))
+  (and ans
+       (for-each (lambda (obj)
+		   (if (not (number? obj))
+		       (slib:error 'parametric-factors 'not 'number? obj)))
+		 ans))
+  ans)
+;@
+(define (L*C*h:DE*94 lch1 lch2 . parametric-factors)
+  (define C* (sqrt (* (cadr lch1) (cadr lch2)))) ;Geometric mean
+  (sqrt (apply + (map /
+		      (map (lambda (x) (* x x)) (map - lch1 lch2))
+		      (list 1		; S_l
+			    (+ 1 (* .045 C*)) ; S_c
+			    (+ 1 (* .015 C*))) ; S_h
+		      (or (color:process-params parametric-factors)
+			  '(1 1 1))))))
+
+;;; CMC-DE is designed only for small color-differences.  But try to do
+;;; something reasonable for large differences.  Use bisector (h*) of
+;;; the hue angles if separated by less than 90.o; otherwise, pick h of
+;;; the color with larger C*.
+;@
+(define (CMC-DE lch1 lch2 . parametric-factors)
+  (apply (lambda (L* C* h_)		;Geometric means
+	   (let ((ang1 (* pi/180 (caddr lch1)))
+		 (ang2 (* pi/180 (caddr lch2))))
+	     (cond ((>= 90 (abs (/ (atan (sin (- ang1 ang2))
+					 (cos (- ang1 ang2)))
+				   pi/180)))
+		    (set! h_ (/ (atan (+ (sin ang1) (sin ang2))
+				      (+ (cos ang1) (cos ang2)))
+				pi/180)))
+		   ((>= (cadr lch1) (cadr lch2)) (caddr lch1))
+		   (else (caddr lch2))))
+	   (let* ((C*^4 (expt C* 4))
+		  (f    (sqrt (/ C*^4 (+ C*^4 1900))))
+		  (T    (if (and (> h_ 164) (< h_ 345))
+			    (+ 0.56 (abs (* 0.2 (cos (* (+ h_ 168) pi/180)))))
+			    (+ 0.36 (abs (* 0.4 (cos (* (+ h_ 35) pi/180)))))))
+		  (S_l  (if (< L* 16)
+			    0.511
+			    (/ (* 0.040975 L*) (+ 1 (* 0.01765 L*)))))
+		  (S_c  (+ (/ (* 0.0638 C*) (+ 1 (* 0.0131 C*))) 0.638))
+		  (S_h  (* S_c (+ (* (+ -1 T) f) 1))))
+	     (sqrt (apply
+		    + (map /
+			   (map (lambda (x) (* x x)) (map - lch1 lch2))
+			   (list S_l S_c S_h)
+			   (or (color:process-params parametric-factors)
+			       '(2 1 1)))))))
+	 (map sqrt (map * lch1 lch2))))
+;@
+(define (XYZ:normalize-colors lst)
+  (define sum (apply max (map (lambda (XYZ) (apply + XYZ)) lst)))
+  (map (lambda (XYZ) (map (lambda (x) (/ x sum)) XYZ)) lst))
+;@
+(define (XYZ:normalize XYZ)
+  (car (XYZ:normalize-colors (list XYZ))))
+
+;;; Chromaticity
+;@
+(define (XYZ->chromaticity XYZ)
+  (define sum (apply + XYZ))
+  (list (/ (car XYZ) sum) (/ (cadr XYZ) sum)))
+;@
+(define (chromaticity->CIEXYZ x y)
+  (list x y (- 1 x y)))
+(define (chromaticity->whitepoint x y)
+  (list (/ x y) 1 (/ (- 1 x y) y)))
+;@
+(define (XYZ->xyY XYZ)
+  (define sum (apply + XYZ))
+  (if (zero? sum)
+      '(0 0 0)
+      (list (/ (car XYZ) sum) (/ (cadr XYZ) sum) (cadr XYZ))))
+;@
+(define (xyY->XYZ xyY)
+  (define x (car xyY))
+  (define y (cadr xyY))
+  (if (zero? y)
+      '(0 0 0)
+      (let ((Y/y (/ (caddr xyY) y)))
+	(list (* Y/y x) (caddr xyY) (* Y/y (- 1 x y))))))
+;@
+(define (xyY:normalize-colors lst . n)
+  (define (nthcdr n lst) (if (zero? n) lst (nthcdr (+ -1 n) (cdr lst))))
+  (define Ys (map caddr lst))
+  (set! n (if (null? n) 1 (car n)))
+  (let ((max-Y (if (positive? n)
+		   (* n (apply max Ys))
+		   (let ()
+		     (require 'sort)
+		     (apply max (nthcdr (- n) (sort Ys >=)))))))
+    (map (lambda (xyY)
+	   (let ((x (max 0 (car xyY)))
+		 (y (max 0 (cadr xyY))))
+	     (define sum (max 1 (+ x y)))
+	     (list (/ x sum)
+		   (/ y sum)
+		   (max 0 (min 1 (/ (caddr xyY) max-Y))))))
+	 lst)))
+
+;;;  http://www.aim-dtp.net/aim/technology/cie_xyz/cie_xyz.htm:
+;;;  Illuminant D65                           0.312713 0.329016
+;; (define CIEXYZ:D65 (chromaticity->whitepoint 0.312713 0.329016))
+;; (define CIEXYZ:D65 (chromaticity->whitepoint 0.3127 0.3290))
+;@
+(define CIEXYZ:D50 (chromaticity->whitepoint 0.3457 0.3585))
+
+;;; With its 16-bit resolution, e-sRGB-16 is extremely sensitive to
+;;; whitepoint.  Even the 6 digits of precision specified above is
+;;; insufficient to make (color->e-srgb 16 d65) ==> (57216 57216 57216)
+;@
+(define CIEXYZ:D65 (e-sRGB->CIEXYZ 16 '(57216 57216 57216)))
+
+;;; http://www.efg2.com/Lab/Graphics/Colors/Chromaticity.htm CIE 1931:
+;@
+(define CIEXYZ:A (chromaticity->whitepoint 0.44757 0.40745)) ; 2856.K
+(define CIEXYZ:B (chromaticity->whitepoint 0.34842 0.35161)) ; 4874.K
+(define CIEXYZ:C (chromaticity->whitepoint 0.31006 0.31616)) ; 6774.K
+(define CIEXYZ:E (chromaticity->whitepoint 1/3 1/3)) ; 5400.K
+
+;;; Converting spectra
+(define cie:x-bar #f)
+(define cie:y-bar #f)
+(define cie:z-bar #f)
+;@
+(define (load-ciexyz . path)
+  (let ((path (if (null? path)
+		  (in-vicinity (library-vicinity) "cie1931.xyz")
+		  (car path))))
+    (set! cie:x-bar (make-vector 80))
+    (set! cie:y-bar (make-vector 80))
+    (set! cie:z-bar (make-vector 80))
+    (call-with-input-file path
+      (lambda (iprt)
+	(do ((wlen 380 (+ 5 wlen))
+	     (idx 0 (+ 1 idx)))
+	    ((>= wlen 780))
+	  (let ((rlen (read iprt)))
+	    (if (not (eqv? wlen rlen))
+		(slib:error path 'expected wlen 'not rlen))
+	    (vector-set! cie:x-bar idx (read iprt))
+	    (vector-set! cie:y-bar idx (read iprt))
+	    (vector-set! cie:z-bar idx (read iprt))))))))
+;@
+(define (wavelength->XYZ wl)
+  (if (not cie:y-bar) (require 'ciexyz))
+  (set! wl (- (/ wl 5.e-9) 380/5))
+  (if (<= 0 wl (+ -1 400/5))
+      (let* ((wlf (inexact->exact (floor wl)))
+	     (res (- wl wlf)))
+	(define (interpolate vect idx res)
+	  (+ (* res (vector-ref vect idx))
+	     (* (- 1 res) (vector-ref vect (+ 1 idx)))))
+	(list (interpolate cie:x-bar wlf res)
+	      (interpolate cie:y-bar wlf res)
+	      (interpolate cie:z-bar wlf res)))
+      (slib:error 'wavelength->XYZ 'out-of-range wl)))
+(define (wavelength->CIEXYZ wl)
+  (XYZ:normalize (wavelength->XYZ wl)))
+(define (wavelength->chromaticity wl)
+  (XYZ->chromaticity (wavelength->XYZ wl)))
+;@
+(define (spectrum->XYZ . args)
+  (define x 0)
+  (define y 0)
+  (define z 0)
+  (if (not cie:y-bar) (require 'ciexyz))
+  (case (length args)
+    ((1)
+     (set! args (car args))
+     (do ((wvln 380.e-9 (+ 5.e-9 wvln))
+	  (idx 0 (+ 1 idx)))
+	 ((>= idx 80) (map (lambda (x) (/ x 80)) (list x y z)))
+       (let ((inten (args wvln)))
+	 (set! x (+ x (* (vector-ref cie:x-bar idx) inten)))
+	 (set! y (+ y (* (vector-ref cie:y-bar idx) inten)))
+	 (set! z (+ z (* (vector-ref cie:z-bar idx) inten))))))
+    ((3)
+     (let* ((vect (if (list? (car args)) (list->vector (car args)) (car args)))
+	    (vlen (vector-length vect))
+	    (x1 (cadr args))
+	    (x2 (caddr args))
+	    (xinc (/ (- x2 x1) (+ -1 vlen)))
+	    (x->j (lambda (x) (inexact->exact (round (/ (- x x1) xinc)))))
+	    (x->k (lambda (x) (inexact->exact (round (/ (- x 380.e-9) 5.e-9)))))
+	    (j->x (lambda (j) (+ x1 (* j xinc))))
+	    (k->x (lambda (k) (+ 380.e-9 (* k 5.e-9))))
+	    (xlo (max (min x1 x2) 380.e-9))
+	    (xhi (min (max x1 x2) 780.e-9))
+	    (jhi (x->j xhi))
+	    (khi (x->k xhi))
+	    (jinc (if (negative? xinc) -1 1)))
+       (if (<= (abs xinc) 5.e-9)
+	   (do ((wvln (j->x (x->j xlo)) (+ wvln (abs xinc)))
+		(jdx (x->j xlo) (+ jdx jinc)))
+	       ((>= jdx jhi)
+		(let ((nsmps (abs (- jhi (x->j xlo)))))
+		  (map (lambda (x) (/ x nsmps)) (list x y z))))
+	     (let ((ciedex (min 79 (x->k wvln)))
+		   (inten (vector-ref vect jdx)))
+	       (set! x (+ x (* (vector-ref cie:x-bar ciedex) inten)))
+	       (set! y (+ y (* (vector-ref cie:y-bar ciedex) inten)))
+	       (set! z (+ z (* (vector-ref cie:z-bar ciedex) inten)))))
+	   (do ((wvln (k->x (x->k xlo)) (+ wvln 5.e-9))
+		(kdx (x->k xlo) (+ kdx 1)))
+	       ((>= kdx khi)
+		(let ((nsmps (abs (- khi (x->k xlo)))))
+		  (map (lambda (x) (/ x nsmps)) (list x y z))))
+	     (let ((inten (vector-ref vect (x->j wvln))))
+	       (set! x (+ x (* (vector-ref cie:x-bar kdx) inten)))
+	       (set! y (+ y (* (vector-ref cie:y-bar kdx) inten)))
+	       (set! z (+ z (* (vector-ref cie:z-bar kdx) inten))))))))
+    (else (slib:error 'spectrum->XYZ 'wna args))))
+(define (spectrum->CIEXYZ . args)
+  (XYZ:normalize (apply spectrum->XYZ args)))
+(define (spectrum->chromaticity . args)
+  (XYZ->chromaticity (apply spectrum->XYZ args)))
+;@
+(define blackbody-spectrum
+  (let* ((c 2.998e8)
+	 (h 6.626e-34)
+	 (h*c (* h c))
+	 (k 1.381e-23)
+	 (pi*2*h*c*c (* 2 pi h*c c)))
+    (lambda (temp . span)
+      (define h*c/kT (/ h*c k temp))
+      (define pi*2*h*c*c*span (* pi*2*h*c*c (if (null? span) 1.e-9 (car span))))
+      (lambda (x)
+	(/ pi*2*h*c*c*span
+	   (expt x 5)
+	   (- (exp (/ h*c/kT x)) 1))))))
+;@
+(define (temperature->XYZ temp)
+  (spectrum->XYZ (blackbody-spectrum temp 5.e-9)))
+(define (temperature->CIEXYZ temp)
+  (XYZ:normalize (temperature->XYZ temp)))
+(define (temperature->chromaticity temp)
+  (XYZ->chromaticity (temperature->XYZ temp)))