let oil_painting matrix (width, height) =
  let radius = 2 and intensityLevels = 20 and
      new_matrix = Array.make_matrix width height (0, 0, 0) in
  let averageR = Array.make intensityLevels 0 and
      averageG = Array.make intensityLevels 0 and
      averageB = Array.make intensityLevels 0 and
      intensityCount = Array.make intensityLevels 0 in
  for i = 0 to width - 1 do
    let left = max 0 (i - radius) and right = min (i + radius) (width - 1) in
    for j = 0 to height - 1 do
      let top = max 0 (j - radius) and bottom = min (j + radius) (height - 1) and maxIndex = ref (-1) in
      for l = top to bottom do
        for k = left to right do
          if l >= 0 && l <= height - 1 && k >= 0 && k <= width - 1 then
            let (r, g, b) = matrix.(k).(l) in
            let r_ = foi r and g_ = foi g and b_ = foi b and
                intensityLevels_ = foi intensityLevels in
            let intensityIndex = iof ((((r_+.g_+.b_)/.3.)*.intensityLevels_)/.256.) in
            intensityCount.(intensityIndex) <- intensityCount.(intensityIndex) + 1;
            averageR.(intensityIndex) <- averageR.(intensityIndex) + r;
	    averageG.(intensityIndex) <- averageG.(intensityIndex) + g;
	    averageB.(intensityIndex) <- averageB.(intensityIndex) + b;
            if !maxIndex= -1 || intensityCount.(!maxIndex) < intensityCount.(intensityIndex) then
              maxIndex := intensityIndex;
        done;
       done;
       let curMax = intensityCount.(!maxIndex) in
       let new_r = averageR.(!maxIndex) / curMax and
           new_g = averageG.(!maxIndex) / curMax and
           new_b = averageB.(!maxIndex) / curMax in
       new_matrix.(i).(j) <- (fix_rgb new_r, fix_rgb new_g, fix_rgb new_b);
     done;
  done;
  new_matrix