color-model.rkt
1 #lang scheme/unit 2 (require "sig.rkt") 3 4 (import) 5 (export framework:color-model^) 6 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 7 ;;; ;;; 8 ;;; matrix ops ;;; 9 ;;; ;;; 10 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 11 12 ;; matrix inversion using cramer's rule 13 14 ;; submatrix : (list-of (list-of num)) int int -> (list-of (list-of num)) 15 ;; submatrix "crosses out" row i and column j from the matrix, returning a new one 16 17 (define (submatrix source i j) 18 (let row-loop ([row 0]) 19 (cond 20 [(eq? row (length source)) null] 21 [(eq? row i) (row-loop (+ row 1))] 22 [else 23 (cons 24 (let col-loop ([col 0]) 25 (cond 26 [(eq? col (length (car source))) null] 27 [(eq? col j) (col-loop (+ col 1))] 28 [else 29 (cons (list-ref (list-ref source row) col) 30 (col-loop (+ col 1)))])) 31 (row-loop (+ row 1)))]))) 32 33 ;;(equal? (submatrix test-matrix 1 2) 34 ;; '((1 2 6) (7 8 4))) 35 36 ;; det : (list-of (list-of num)) -> num 37 38 (define (det matrix) 39 (if (null? matrix) 40 1 41 (let loop ([row 0] [sign 1]) 42 (if (= row (length matrix)) 43 0 44 (+ (* sign 45 (list-ref (list-ref matrix row) 0) 46 (det (submatrix matrix row 0))) 47 (loop (+ row 1) (- sign))))))) 48 49 ;;(define square-test-matrix '((3 20 3) (37 0 8) (2 1 4))) 50 51 ;;(= (det square-test-matrix) -2553) 52 53 ;; invert : (list-of (list-of num)) -> (list-of (list-of num)) 54 55 (define (matrix-invert matrix) 56 (let-values ([(width height) (matrix-dimension matrix)]) 57 (when (not (= width height)) 58 (error 'invert "matrix is not square: ~s" matrix)) 59 (let ([delta-inv (/ 1 (det matrix))]) 60 (let row-loop ([row 0] [sign 1]) 61 (if (= row (length matrix)) 62 null 63 (cons 64 (let col-loop ([col 0] [sign sign]) 65 (if (= col (length (car matrix))) 66 null 67 (cons (* delta-inv 68 sign 69 (det (submatrix matrix col row))) 70 (col-loop (+ col 1) (- sign))))) 71 (row-loop (+ row 1) (- sign)))))))) 72 73 ;;(equal? (matrix-invert square-test-matrix) 74 ;; '((8/2553 77/2553 -160/2553) (44/851 -2/851 -29/851) (-1/69 -1/69 20/69))) 75 76 ;; matrix-dimension : (list-of (list-of num)) -> (values num num) 77 ;; takes a matrix, returns width and height 78 79 (define (matrix-dimension matrix) 80 (when (not (pair? matrix)) 81 (error 'matrix-dimension "matrix argument is not a list: ~s" matrix)) 82 (let ([height (length matrix)]) 83 (when (= height 0) 84 (error 'matrix-dimension "matrix argument is empty: ~s" matrix)) 85 (when (not (pair? (car matrix))) 86 (error 'matrix-dimension "matrix row is not a list: ~s" (car matrix))) 87 (let ([width (length (car matrix))]) 88 (when (= width 0) 89 (error 'matrix-dimension "matrix argument has width 0: ~s" matrix)) 90 (let loop ([rows matrix]) 91 (if (null? rows) 92 (values width height) 93 (begin 94 (when (not (pair? (car rows))) 95 (error 'matrix-dimension "row is not a list: ~s" (car rows))) 96 (when (not (= width (length (car rows)))) 97 (error 'matrix-dimension "rows have different widths: ~s and ~s" width (length (car rows)))) 98 (loop (cdr rows)))))))) 99 100 ;; transpose : (list-of (list-of num)) -> (list-of (list-of num)) 101 (define (transpose vector) (apply map list vector)) 102 103 104 ;; test code 105 ;;(equal? (transpose '((3 2 1) (9 8 7))) '((3 9) (2 8) (1 7))) 106 107 ;; inner-product : (list-of num) (list-of num) -> num 108 (define (inner-product a b) 109 (foldl + 0 (map * a b))) 110 111 ;; test code 112 ;; (= (inner-product '(4 1 3) '(0 3 4)) 15) 113 114 ;; matrix-multiply: (list-of (list-of num)) (list-of (list-of num)) -> (list-of (list-of num)) 115 ;; multiplies the two matrices. 116 (define (matrix-multiply a b) 117 (let-values ([(width-a height-a) (matrix-dimension a)] 118 [(width-b height-b) (matrix-dimension b)]) 119 (when (not (= width-a height-b)) 120 (error 'matrix-multiply "matrix dimensions do not match for multiplication")) 121 (let ([b-t (transpose b)]) 122 (map (λ (row) 123 (map (λ (col) 124 (inner-product row col)) 125 b-t)) 126 a)))) 127 128 ;; test code 129 ;; (equal? (matrix-multiply '((1 2 3 4) (9 8 3 2)) '((0) (2) (0) (3))) 130 ;; '((16) (22))) 131 132 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 133 ;;; ;;; 134 ;;; color model ;;; 135 ;;; ;;; 136 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 137 138 ;; ITU reccommendation phosphors: 139 140 ;; red green blue 141 ;;x 0.64 0.29 0.15 142 ;;y 0.33 0.60 0.06 143 ;; 144 ;; white point: 145 ;; c : x-w = 0.313, y-w = 0.329, big-y-w = 100.0 146 147 (define x-r 0.64) 148 (define y-r 0.33) 149 (define x-g 0.29) 150 (define y-g 0.60) 151 (define x-b 0.15) 152 (define y-b 0.06) 153 154 (define z-r (- 1 x-r y-r)) 155 (define z-g (- 1 x-g y-g)) 156 (define z-b (- 1 x-b y-b)) 157 158 (define x-w 0.313) 159 (define y-w 0.329) 160 (define big-y-w 100.0) 161 162 (define-struct xyz (x y z)) 163 164 (define (xy-big-y->xyz x y big-y) 165 (let ([sigma (/ big-y y)]) 166 (make-xyz 167 (* x sigma) 168 (* y sigma) 169 (* (- 1 x y) sigma)))) 170 171 (define xyz-white (xy-big-y->xyz x-w y-w big-y-w)) 172 173 ;;`((,(xyz-x xyz-white) ,x-r ,x-g ,x-b) 174 ;; (,(xyz-y xyz-white) ,y-r ,y-g ,y-b) 175 ;; (,(xyz-z xyz-white) ,z-r ,z-g ,z-b)) 176 177 ;; sigmas were calculated by soving a set of linear equations based upon ntsc standard phosphors 178 179 (define pre-matrix `((,x-r ,x-g ,x-b) 180 (,y-r ,y-g ,y-b) 181 (,z-r ,z-g ,z-b))) 182 183 (define-values (sigma-r sigma-g sigma-b) 184 (let* ([inversion 185 (matrix-invert pre-matrix)] 186 [sigmas 187 (matrix-multiply inversion `((,(xyz-x xyz-white)) 188 (,(xyz-y xyz-white)) 189 (,(xyz-z xyz-white))))]) 190 (apply values (car (transpose sigmas))))) 191 192 ;; (printf "should be equal to xyz-white: \n~a\n" 193 ;; (matrix-multiply pre-matrix `((,sigma-r) (,sigma-g) (,sigma-b)))) 194 195 (define rgb->xyz-matrix 196 (map (λ (row) 197 (map (λ (row-elt scalar) (* row-elt scalar 1/255)) row `(,sigma-r ,sigma-g ,sigma-b))) 198 pre-matrix)) 199 200 (define xyz->rgb-matrix 201 (matrix-invert rgb->xyz-matrix)) 202 203 ;;(printf "should be identity: \n~a\n" (matrix-multiply rgb->xyz-matrix xyz->rgb-matrix)) 204 205 (define (rgb->xyz r g b) 206 (apply make-xyz (car (transpose (matrix-multiply rgb->xyz-matrix (transpose `((,r ,g ,b)))))))) 207 208 ;;(print-struct #t) 209 ;; (printf "should be xyz-white: \n~a\n" (rgb->xyz 255 255 255)) 210 211 (define (xyz->rgb x y z) 212 (car (transpose (matrix-multiply xyz->rgb-matrix (transpose `((,x ,y ,z))))))) 213 214 ;;l* = 116(y/big-y-n)^1/3 - 16, y/big-y-n > 0.01 215 ;;u* = 13 l*(u-p - u-p-n) 216 ;;v* = 13 l*(v-p - v-p-n) 217 ;; 218 ;;u-p = (4x)/(x+15y+3z) v-p = (9y)/(x+15y+3z) 219 ;;u-p-n = (same but with -n) v-p-n = (same but with -n) 220 221 ;; the following transformation is undefined if the y component 222 ;; is zero. So if it is, we bump it up a little. 223 224 (define (xyz-tweak xyz) 225 (let* ([y (xyz-y xyz)]) 226 (make-xyz (xyz-x xyz) (if (< y 0.01) 0.01 y) (xyz-z xyz)))) 227 228 (define-struct luv (l u v)) 229 230 (define (xyz-denom xyz) 231 (+ (xyz-x xyz) (* 15 (xyz-y xyz)) (* 3 (xyz-z xyz)))) 232 233 (define (xyz-u-p xyz) 234 (/ (* 4 (xyz-x xyz)) (xyz-denom xyz))) 235 236 (define (xyz-v-p xyz) 237 (/ (* 9 (xyz-y xyz)) (xyz-denom xyz))) 238 239 (define (xyz->luv xyz) 240 (let ([xyz (xyz-tweak xyz)]) 241 (let* ([l (- (* 116 (expt (/ (xyz-y xyz) (xyz-y xyz-white)) 242 1/3)) 243 16)] 244 [u-p (xyz-u-p xyz)] 245 [u-p-white (xyz-u-p xyz-white)] 246 [v-p (xyz-v-p xyz)] 247 [v-p-white (xyz-v-p xyz-white)]) 248 (make-luv l (* 13 l (- u-p u-p-white)) (* 13 l (- v-p v-p-white)))))) 249 250 (define (luv-distance a b) 251 (expt (+ (expt (- (luv-l a) (luv-l b)) 2) 252 (expt (- (luv-u a) (luv-u b)) 2) 253 (expt (- (luv-v a) (luv-v b)) 2)) 254 1/2)) 255 256 (define (rgb-color-distance r-a g-a b-a r-b g-b b-b) 257 (let* ([luv-a (xyz->luv (rgb->xyz r-a g-a b-a))] 258 [luv-b (xyz->luv (rgb->xyz r-b g-b b-b))]) 259 (luv-distance luv-a luv-b))) 260 261 ;;(rgb-color-distance 0 0 0 0 0 0) 262 ;; (print-struct #t) 263 ;; (xyz->luv (make-xyz 95.0 100.0 141.0)) 264 ;; (xyz->luv (make-xyz 60.0 80.0 20.0)) 265 266 ;; formulas from https://en.wikipedia.org/wiki/HSL_and_HSV 267 (define (rgb->hsl r255 g255 b255) 268 (define r (/ r255 255)) 269 (define g (/ g255 255)) 270 (define b (/ b255 255)) 271 (define MAX (max r g b)) 272 (define MIN (min r g b)) 273 (define Δ (- MAX MIN)) 274 (define H1 275 (cond 276 [(= Δ 0) 0] 277 [(= MAX r) (* 60 (/ (- g b) Δ))] 278 [(= MAX g) (* 60 (+ 2 (/ (- b r) Δ)))] 279 [(= MAX b) (* 60 (+ 4 (/ (- r g) Δ)))])) 280 (define H (if (H1 . < . 0) (+ H1 360) H1)) 281 (define L (/ (+ MAX MIN) 2)) 282 (define S 283 (cond 284 [(= MAX 0) 0] 285 [(= MIN 1) 0] 286 [else 287 (/ (- MAX L) 288 (min L (- 1 L)))])) 289 (values H S L)) 290 291 (define (hsl->rgb H S L) 292 (define C (* (- 1 (abs (- (* 2 L) 1))) S)) 293 (define H* (/ H 60)) 294 (define X (* C (- 1 (abs (- (MOD H* 2) 1))))) 295 (define-values (R1 G1 B1) 296 (cond 297 [(<= H* 1) (values C X 0)] 298 [(<= H* 2) (values X C 0)] 299 [(<= H* 3) (values 0 C X)] 300 [(<= H* 4) (values 0 X C)] 301 [(<= H* 5) (values X 0 C)] 302 [(<= H* 6) (values C 0 X)])) 303 (define m (- L (/ C 2))) 304 (define-values (R G B) 305 (values (+ R1 m) (+ G1 m) (+ B1 m))) 306 (define (to-255 n) (inexact->exact (round (* n 255)))) 307 (values (to-255 R) (to-255 G) (to-255 B))) 308 309 ;; definition taken from fortran's docs ... 310 (define (MOD A P) 311 (- A (* (floor (/ A P)) P)))