Patch title: Release 94 bulk changes
Abstract:
File: /pliant/graphic/sample/color.pli
Key:
    Removed line
    Added line
   
module "/pliant/math/functions.pli"
module "/pliant/language/compiler.pli"
module "/pliant/language/stream.pli"
module "/pliant/graphic/color/spectrum.pli"
module "/pliant/graphic/color/color.pli"

if true
  constant grid_file "file:/pliant_data/pliant/graphic/color/epson4000/cmyk"
  # constant grid_file "file:/tmp/cmyk"
  constant deaden 0.04
  constant dot -0.04
  constant speed 2
  constant limit 0.01
eif false
  constant grid_file "file:/pliant_data/pliant/graphic/color/cromalin/cmyk"
  constant deaden 0.12
  constant dot -0.30
else
  constant grid_file "file:/pliant_data/pliant/graphic/color/helio/cmyk"
  constant deaden 0.04
  constant dot -0.10


constant steps 5
gvar Array:ColorSpectrum32 grid

constant trace true
constant spectral false
constant km true


#---------------------------------------------------------------
#   sample math function


function sinh x -> y
  arg Float x y
  y := 0.5*((exp x)-(exp -x))

function cosh x -> y
  arg Float x y
  y := 0.5*((exp x)+(exp -x))

function tanh x -> y
  arg Float x y
  y := sinh:x/cosh:x

function coth x -> y
  arg Float x y
  y := cosh:x/sinh:x

function tanh s1 -> s
  arg ColorSpectrum32 s1 s
  for (var Int i) 0 ColorSpectrum32:size\Float32:size-1
    addressof:s map Float32 i := tanh (addressof:s1 map Float32 i)


function kubelka_munk K S Rg X -> R
  arg Float K S Rg X R
  var Float a := (S+K)/S
  var Float b := (a^2-1)^0.5
  var Float cot := coth b*S*X
  R := ( 1-Rg*(a-b*cot) ) / (a-Rg+b*cot)

type KubelkaMunkSample
  field Float Rg X R

constant K_divided_by_S_maxi 500 # Rinfini cannot be less than 0.001

function kubelka_munk_setup samples Kmini Kmaxi Kstep Smini Smaxi Sstep K S
  arg Array:KubelkaMunkSample samples ; arg Float Kmini Kmaxi Kstep Smini Smaxi Sstep ; arg_w Float K S
  var Float best := 1e9
  var Float Ktest := Kmini
  while Ktest<Kmaxi+Kstep/2
    var Float Stest := Smini
    while Stest<Smaxi+Sstep/2
      if Ktest/Stest<=K_divided_by_S_maxi
        var Float err := 0
        for (var Int i) 0 samples:size-1
          var Float delta := (kubelka_munk (max Ktest 1e-5) (max Stest 1e-6) samples:i:Rg samples:i:X) - samples:i:R
          err += delta*delta
        if err<best
          K := Ktest ; S := Stest
          best := err
      Stest += Sstep
    Ktest += Kstep

function kubelka_munk_setup samples K S
  arg Array:KubelkaMunkSample samples ; arg_w Float K S
  kubelka_munk_setup samples 0 10 0.1 0 10 0.1 K S
  kubelka_munk_setup samples (max K-1 0) K+1 0.01 (max S-1 0) S+1 0.01 K S
  kubelka_munk_setup samples (max K-0.1 0) K+0.1 0.001 (max S-0.1 0) S+0.1 0.001 K S
  if trace
    console "K = " K " , S = " S
    var Float x := K/S
    var Float Rinf := 1+x-(x^2+2*x)^0.5
    console " , K/S = " (string K/S "fixed 3") " , quotient = " (string samples:0:R/Rinf "fixed 3") eol
  if K>=9.9
    console "K overflow !" eol
  if S>=9.9
    console "S overflow !" eol


function spectral_distance s1 s2 -> d
  arg ColorSpectrum32 s1 s2 ; arg Float d
  d := 0
  for (var Int i) 0 ColorSpectrum32:size\Float32:size-1
    # d := max d (abs (addressof:s1 map Float32 i)-(addressof:s2 map Float32 i))
    d := max d (abs (addressof:Y_spectrum map Float32 i)*((addressof:s1 map Float32 i)-(addressof:s2 map Float32 i)))
  
method s is_ok -> c
  arg ColorSpectrum32 s ; arg CBool c
  for (var Int i) 0 ColorSpectrum32:size\Float32:size-1
    if (addressof:s map Float32 i)=undefined
      return false  
  c := true


#---------------------------------------------------------------
#   test simulation on measured grid


function load
  grid size := steps^4  
  (var Stream s) open grid_file in
  while not s:atend
    if (s:readline parse (var Int cyan) (var Int magenta) (var Int yellow) (var Int black) ":" (var ColorSpectrum spectrum))
      grid (cyan+1)\(256\(steps-1))+(magenta+1)\(256\(steps-1))*steps+(yellow+1)\(256\(steps-1))*steps^2+(black+1)\(256\(steps-1))*steps^3 := spectrum

gvar Array:ColorSpectrum32 Kink
gvar Array:ColorSpectrum32 Sink

function setup
  Kink size := 4 ; Sink size := 4
  var Array:KubelkaMunkSample samples
  for (var Int i) 0 3
    for (var Int j) 0 ColorSpectrum32:size\Float32:size-1
      if not trace
        console "  " i*ColorSpectrum32:size\Float32:size+j "/" 4*ColorSpectrum32:size\Float32:size "   [cr]"
      if false
        samples size := steps-1
        for (var Int k) 0 steps-2
          samples:k Rg := (addressof grid:0) map Float32 j
          samples:k X := exposure (k+1)/(steps-1) dot
          samples:k R := addressof:(grid (k+1)*steps^i) map Float32 j
      else
        samples size := 2
        samples:0 Rg := (addressof grid:0) map Float32 j
        samples:0 X := 1
        samples:0 R := addressof:(grid (steps-1)*steps^i) map Float32 j
        var Int black := shunt i<>3 (steps-1)*steps^3 (steps-1)*(1+steps+steps^2)
        samples:1 Rg := (addressof grid:black) map Float32 j
        samples:1 X := 1
        samples:1 R := addressof:(grid black+(steps-1)*steps^i) map Float32 j
        # var Int black := steps^4-1
        # samples:1 Rg := addressof:(grid black-(steps-1)*steps^i) map Float32 j
        # samples:1 X := 1
        # samples:1 R := addressof:(grid black) map Float32 j
      if trace
        console "ink " i " wavelength " 400+10*j " is " (addressof:(grid (steps-1)*steps^i) map Float32 j) " : "
      kubelka_munk_setup samples (var Float K) (var Float S)
      (addressof Kink:i) map Float32 j := K
      (addressof Sink:i) map Float32 j := S
  if trace
    for (var Int i) 0 3
      console "K" i " = " (cast Kink:i ColorSpectrum) eol
      console "S" i " = " (cast Sink:i ColorSpectrum) eol
  else
    console "[cr]            [cr]"

# lap 0 is naive (exp -density)
# lap 1 is same with deaden (adding logs on a non linear scale)
# lap 2 is test
# lap 3 is Kubelka Munk

function test
  for (var Int lap) 0 (shunt km 3 2)
    var Float davr := 0 ; var Float dmax := 0 ; var Int worse := 0
    var ColorLCh measured computed
    var Float s_davr := 0 ; var Float s_dmax := 0 ; var Int s_worse := 0
    var ColorLCh s_measured s_computed
    for (var Int i) 0 grid:size-1
      var ColorSpectrum32 f Kfinal Sfinal
      var Float total_density := 0
      var ColorSpectrum32 total_diffusion := cast 0 ColorSpectrum32
      for (var Int dim) 0 3
        var Float density := (i\steps^dim)%steps/(steps-1)
        var ColorSpectrum32 ink := (grid (steps-1)*steps^dim)/grid:0
        total_density += density
        total_diffusion += density*ink
      var Float deaden2 := 0.025*(tanh total_density)
      if lap=0
        f := cast 1 ColorSpectrum32
      eif lap=1
        f := cast 0 ColorSpectrum32
      eif lap=2
        if true
          f := cast 1 ColorSpectrum32
        else
          f := cast 0 ColorSpectrum32
      eif lap=3
        Kfinal := cast 0 ColorSpectrum32
        Sfinal := cast 0 ColorSpectrum32
      for (var Int dim) 0 3
        var Float density := (i\steps^dim)%steps/(steps-1)
        var ColorSpectrum32 ink := (grid (steps-1)*steps^dim)/grid:0
        if lap=0
          f *= ink^density
        eif lap=1
          f += exposure (-1)*(log ink^density) -deaden
        eif lap=2
          if true
            f *= exp density*log:ink
          else
            f += exposure (-1)*(log ink^density) -deaden2
        eif lap=3
          density := exposure density dot
          Kfinal += density*Kink:dim
          Sfinal += density*Sink:dim
      if lap=0
        void
      eif lap=1
        f := exp (-1)*(unexposure f -deaden)
      eif lap=2
        if true
          var Float z := speed*total_density 
          z := tanh z # 1-(exp -z)
          f += limit*z/total_density*total_diffusion
          # f += limit*z*(cast 1 ColorSpectrum32)
        else
          f := exp (-1)*(unexposure f -deaden2)
      eif lap=3
        for (var Int j) 0 ColorSpectrum32:size\Float32:size-1
          var Float K := addressof:Kfinal map Float32 j
          var Float S := addressof:Sfinal map Float32 j
          var Float Rg := (addressof grid:0) map Float32 j
          var Float R := kubelka_munk K S Rg 1
          addressof:f map Float32 j := R/Rg
      if f:is_ok
        var Float dist := cmc_distance (filter_XYZ grid:i/grid:0) filter_XYZ:f
        davr += dist/grid:size
        if dist>dmax
          worse := i
          dmax := dist
          measured := cast (filter_XYZ grid:i/grid:0) ColorLCh
          computed := cast filter_XYZ:f ColorLCh
        var Float s_dist := spectral_distance grid:i/grid:0 f
        s_davr += s_dist/grid:size
        if s_dist>s_dmax
          s_worse := i
          s_dmax := s_dist
          s_measured := cast (filter_XYZ grid:i/grid:0) ColorLCh
          s_computed := cast filter_XYZ:f ColorLCh
      else
        console "spectral formulation oops for sample " i eol
    if lap=0
      console eol
      console "MULTIPLY" eol
    eif lap=1
      console "DEADEN" eol
    eif lap=2
      console "ADVANCED DEADEN" eol
    eif lap=3
      console "KUBELKA MUNK" eol
    console "average distance is " (string davr "fixed 1") eol
    console "maximum distance is " (string dmax "fixed 1") " for sample"
    for (var Int dim) 0 3
      console " " (cast worse\steps^dim%steps/(steps-1)*100 Int) "%"
    console eol
    console "computed LCh " (string computed:L "fixed 1") " " (string computed:C "fixed 1") " " (string computed:h "fixed 0") eol
    console "measured LCh " (string measured:L "fixed 1") " " (string measured:C "fixed 1") " " (string measured:h "fixed 0") eol
    if spectral
      console "average spectral distance is " (string s_davr "fixed 3") eol
      console "maximum spectral distance is " (string s_dmax "fixed 3") " for sample "
      for (var Int dim) 0 3
        console " " (cast s_worse\steps^dim%steps/(steps-1)*100 Int) "%"
      console eol
      console "computed LCh " (string s_computed:L "fixed 1") " " (string s_computed:C "fixed 1") " " (string s_computed:h "fixed 0") eol
      console "measured LCh " (string s_measured:L "fixed 1") " " (string s_measured:C "fixed 1") " " (string s_measured:h "fixed 0") eol
    console eol  

      
load
if km
  setup
test