############################################################################
#
#	File:     colrmodl.icn
#
#	Subject:  Procedures to convert between color models
#
#	Author:   Ralph E. Griswold
#
#	Date:     December 5, 1995
#
############################################################################
#
#   This file is in the public domain.
#
############################################################################
#
#     These procedures convert between various color models. A color
#  value is represented by a record (see the declarations below).
#
#  Color values are normalized to a maximum of 1.0.
#
############################################################################
#
#  Acknowledgement:  Some of the procedures here are based on information
#  given in Computer Graphics; Principles and Practice, second edition;
#  James D. Foley, Andries van Dam, Steven K. Feiner, and John F. Hughes;
#  Addison-Wesley Publishing Company; 1990.
#
############################################################################
#
#  Note:  These procedures have not been extensively tested. Those related
#  to the YIQ model are particularly in question.
#
############################################################################
#
#  Links:  matrix, numbers
#
############################################################################

link matrix
link numbers

record rgb(r, g, b)
record cmy(c, m, y)
record cmyk(c, m, y, k)
record yiq(y, i, q)
record hsv(h, s, v)
record hls(h, l, s)

procedure rgb2cmy(color)

   return cmy(1.0 - color.r, 1.0 - color.g, 1.0 - color.b)

end

procedure cmy2rgb(color)

   return rgb(1.0 - color.c, 1.0 - color.m, 1.0 - color.y)

end

#  Note: The following procedure illustrates the principle of
#  undercolor removal, but for pragmatic reasons, it does not
#  produce acceptable results in process printing.

procedure cmy2cmyk(color)
   local k

   k := min(color.c, color.m, color.y)

   return cmyk(color.c - k, color.m - k, color.y - k, k)

end

procedure cmyk2cmy(color)
   local kdelta

   kdelta := color.k / 3

   return cmy(color.c + kdelta, color.m + kdelta, color.y + kdelta)

end

#
#  Note: The RGB specification is assumed to be based on the standard
#  NTSC phosphors.  See the reference cited above.

procedure rgb2yiq(color)
   static M, R, Y

   initial {
      M := create_matrix(3, 3)
      M[1, 1] := 0.299
      M[1, 2] := 0.587
      M[1, 3] := 0.114
      M[2, 1] := 0.596
      M[2, 2] := -0.275
      M[2, 3] := -0.321
      M[3, 1] := 0.212
      M[3, 2] := -0.528
      M[3, 3] := 0.311
      }

   R := create_matrix(3, 1)
   R[1][1] := color.r
   R[2][1] := color.g
   R[3][1] := color.b

   Y := mult_matrix(M, R)

   return yiq(Y[1][1], Y[2][1], Y[3][1])

end

procedure yiq2rgb(color)
   static M, R, Y

   initial {
      M := create_matrix(3, 3)
      M[1, 1] := 1.0031
      M[1, 2] := 0.9548
      M[1, 3] := 0.6179
      M[2, 1] := 0.9968
      M[2, 2] := -0.2707
      M[2, 3] := -0.6448
      M[3, 1] := 1.0084
      M[3, 2] := -1.1005
      M[3, 3] := 1.6996
      }

   Y := create_matrix(3, 1)
   Y[1][1] := color.y
   Y[2][1] := color.i
   Y[3][1] := color.q

   R := mult_matrix(M, Y)

   return rgb(R[1][1], R[2][1], R[3][1])

end

procedure rgb2hsv(color)
   local maximum, minimum, delta, h, s, v

   maximum := max(color.r, color.g, color.b)
   minimum := min(color.r, color.g, color.b)
   delta := maximum - minimum

   v := maximum

   if maximum ~= 0 then s := delta / maximum
   else s := 0

   if s = 0 then h := -1.0		# undefined
   else {
      if color.r = maximum then {
         h := (color.g - color.b) / delta
         }
      else if color.g = maximum then {
         h := 2 + (color.b - color.r) / delta
         }
      else if color.b = maximum then {
         h := 4 + (color.r - color.g) / delta
         }
      h := h * 60
      if h < 0 then h +:= 360.0		# make sure hue is nonnegative
      }

   return hsv(h, s, v)

end

procedure hsv2rgb(color)

   local h, i, f, p, q, t, s, v

   if color.s = 0 then {
      if color.h = -1 then {
         return rgb(color.v, color.v, color.v)
         }
      else stop("*** error in HSV to RGB conversion")
      }
   else {
      h := color.h
      v := color.v
      s := color.s
      if h = 360.0 then h := 0.0
      h /:= 60
      i := floor(h)
      f := h - i
      p := v * (1.0 - s)
      q := v * (1.0 - s * f)
      t := v * (1.0 - (s * (1.0 - f)))
      return case i of {
         0:  rgb(v, t, p)
         1:  rgb(q, v, p)
         2:  rgb(p, v, t)
         3:  rgb(p, q, v)
         4:  rgb(t, p, v)
         5:  rgb(v, p, q)
         default:  stop("*** error in HSV to RGB conversion")
         }
      }

end

procedure rgb2hls(color)
   local maximum, minimum, delta, sum, h, s, l

   maximum := max(color.r, color.b, color.g)
   minimum := min(color.r, color.b, color.g)

   delta := maximum - minimum
   sum := maximum + minimum
   l := sum / 2				# lightness

   if maximum = minimum then {		# achromatic case
      s := 0.0
      h := -1.0
      }
   else {
      if l <= 0.5 then
         s := delta / sum
      else s := delta / (2 - sum)

      if color.r = maximum then
         h := (color.g - color.r) / delta
      else if color.g = maximum then
         h := 2 + (color.b - color.r) / delta
      else if color.b = maximum then
         h := 4 + (color.r - color.g) / delta
      h *:= 60				# convert to degrees
      if h < 0.0 then h +:= 360.0	# make positive

      return hls(h, l, s)
      }

end

procedure hls2rgb(color)
   local h, l, s, m1, m2

   h := color.h
   l := color.l
   s := color.s

   if l <= 0.5 then m2 := l * (1 + s)
   else m2 := l + s - l * s
   m1 := 2 * l - m2
   if s = 0 then {			# achromatic case
      if h = -1.0 then return rgb(l, l, l)
      else stop("*** error in HLS specification")
      }
   else {
      return rgb(
         color_value(m1, m2, h + 120.0),
         color_value(m1, m2, h),
         color_value(m1, m2, h - 120.0)
         )
      }

end

procedure color_value(m1, m2, h)

   if h > 360.0 then h -:= 360.0
   else if h < 0.0 then h +:= 360.0
   if h < 60.0 then return m1 + (m2 - m1) * h / 60.0
   else if h < 180.0 then return m2
   else if h < 240.0 then return m1 + (m2 - m1) * (240.0 - h) / 60.0
   else return m1

end