diff --git a/share/grafx2/scripts/samples_2.4/picture/others-8bit/lib/ostro_other.lua b/share/grafx2/scripts/samples_2.4/picture/others-8bit/lib/ostro_other.lua
new file mode 100644
index 00000000..55a459c9
--- /dev/null
+++ b/share/grafx2/scripts/samples_2.4/picture/others-8bit/lib/ostro_other.lua
@@ -0,0 +1,141 @@
+-- ostro_zx.lua : converts a color image into a 
+-- ZX-like image (8+8 fixed colors with color clash) 
+-- using Ostromoukhov's error diffusion algorithm.
+--
+-- Version: 03/21/2017
+--
+-- Copyright 2016-2017 by Samuel Devulder
+--
+-- This program is free software; you can redistribute
+-- it and/or modify it under the terms of the GNU
+-- General Public License as published by the Free
+-- Software Foundation; version 2 of the License.
+-- See 
+
+run('../../thomson/lib/ostromoukhov.lua')
+
+-- get screen size
+local screen_w, screen_h = getpicturesize()
+
+OtherDither = {}
+function OtherDither:new(a)
+	local o = { -- default ZX values
+		-- width of the screen
+		width=a and a.width or 256,
+		-- height of the screen
+		height=a and a.height or 192,
+		-- size of Nx1 clash size
+		clash_size=a and a.clash_size or 8,
+		-- normalize the picture levels (like in imagemagick)
+		normalize=a and a.normalize or 0.005,
+		-- put a pixel
+		pset=a and a.pset or function(self,x,y,c)
+			if c<0 then c=-c-1 end
+			self.screen[x+y*self.width] = c
+		end,
+		-- init gfx data
+		setGfx=a and a.setGfx or function(self)
+			self.screen={}
+		end,
+		-- update gfx to screen
+		updatescreen=a and a.updatescreen or function(self)
+			for i=0,255 do setcolor(i,0,0,0) end
+			for y=0,self.height-1 do
+				for x=0,self.width-1 do
+					putpicturepixel(x,y,self.screen[x+y*self.width] or 0)
+				end
+			end
+			-- refresh palette
+			for i,v in ipairs(self.pal) do
+				local r=v % 16
+				local g=math.floor(v/16)  % 16
+				local b=math.floor(v/256) % 16
+				setcolor(i+thomson._palette.offset-1, 
+						 thomson.levels.pc[r+1], 
+						 thomson.levels.pc[g+1], 
+						 thomson.levels.pc[b+1])
+			end
+			updatescreen()		
+		end,
+		-- palette with thomson ordering (to use thomson's
+		-- lib support)
+		pal= a and a.pal or {
+			0x000,0xF00,0x00F,0xF0F,0x0F0,0xFF0,0x0FF,0xFFF,
+			0x000,0x200,0x002,0x202,0x020,0x220,0x022,0x222
+		}
+	}
+	setmetatable(o, self)
+	self.__index = self
+	return o
+end
+
+-- Converts ZX coordinates (0-255,0-191) into screen coordinates
+function OtherDither:to_screen(x,y)
+	local i,j;
+	if screen_w/screen_h < self.width/self.height then
+		i = x*screen_h/self.height
+		j = y*screen_h/self.height
+	else
+		i = x*screen_w/self.width
+		j = y*screen_w/self.width
+	end
+	return math.floor(i), math.floor(j)
+end
+
+-- return the Color @(x,y) in linear space (0-255)
+-- corresonding to the other platform screen
+OtherDither._getLinearPixel = {} -- cache
+function OtherDither:getLinearPixel(x,y)
+	local k=x+y*self.width
+	local p = self._getLinearPixel and self._getLinearPixel[k]
+	if not p then
+		local x1,y1 = self:to_screen(x,y)
+		local x2,y2 = self:to_screen(x+1,y+1)
+		if x2==x1 then x2=x1+1 end
+		if y2==y1 then y2=y1+1 end
+
+		p = Color:new(0,0,0)
+		for j=y1,y2-1 do
+			for i=x1,x2-1 do
+				p:add(getLinearPictureColor(i,j))
+			end
+		end
+		p:div((y2-y1)*(x2-x1)) --:floor()
+					
+		if self._getLinearPixel then 
+			self._getLinearPixel[k]=p 
+		end
+	end
+	
+	return self._getLinearPixel and p:clone() or p
+end
+
+function OtherDither:ccAcceptCouple(c1,c2)
+	-- bright colors can't mix with dimmed ones
+	return c1~=c2 and ((c1<=8 and c2<=8) or (c1>8 and c2>8))
+end
+
+function OtherDither:dither()
+	local NORMALIZE=Color.NORMALIZE
+	Color.NORMALIZE=self.normalize
+	
+	local dither=OstroDither:new(self.pal)
+	dither.ccAcceptCouple = function(dither,c1,c2) return self:ccAcceptCouple(c1,c2) end
+	dither.clash_size = self.clash_size
+	dither.attenuation = .9
+	
+	self:setGfx()
+	dither:ccDither(self.width,self.height, 
+				  function(x,y) return self:getLinearPixel(x,y) end, 
+				  function(x,y,c) self:pset(x,y,c) end, 
+				  true, 
+				  function(y) 
+					thomson.info("Converting...",
+							math.floor(y*100/self.height),"%") 
+				  end,true)
+	-- refresh screen
+	setpicturesize(self.width,self.height)
+	self:updatescreen()
+	finalizepicture()
+	Color.NORMALIZE=NORMALIZE
+end
diff --git a/share/grafx2/scripts/samples_2.4/picture/others-8bit/ostro_oric.lua b/share/grafx2/scripts/samples_2.4/picture/others-8bit/ostro_oric.lua
new file mode 100644
index 00000000..f5d0ccb8
--- /dev/null
+++ b/share/grafx2/scripts/samples_2.4/picture/others-8bit/ostro_oric.lua
@@ -0,0 +1,27 @@
+-- ostro_zx.lua : converts a color image into a 
+-- Oric image (8+8 fixed colors with color clash) 
+-- using Ostromoukhov's error diffusion algorithm.
+--
+-- Version: 03/21/2017
+--
+-- Copyright 2016-2017 by Samuel Devulder
+--
+-- This program is free software; you can redistribute
+-- it and/or modify it under the terms of the GNU
+-- General Public License as published by the Free
+-- Software Foundation; version 2 of the License.
+-- See 
+
+run('lib/ostro_other.lua')
+
+OtherDither:new{
+	width=240,
+	height=200,
+	clash_size=6,
+	pal={0x000,0x00F,0x0F0,0x0FF,0xF00,0xF0F,0xFF0,0xFFF},
+	pset=function(self,x,y,c)
+		if x<6 then c=0    end
+		if c<0 then c=-c-1 end
+		self.screen[x+y*self.width] = c
+	end
+}:dither()
diff --git a/share/grafx2/scripts/samples_2.4/picture/others-8bit/ostro_zx.lua b/share/grafx2/scripts/samples_2.4/picture/others-8bit/ostro_zx.lua
new file mode 100644
index 00000000..4c0af32a
--- /dev/null
+++ b/share/grafx2/scripts/samples_2.4/picture/others-8bit/ostro_zx.lua
@@ -0,0 +1,17 @@
+-- ostro_zx.lua : converts a color image into a 
+-- ZX image (8+8 fixed colors with color clash) 
+-- using Ostromoukhov's error diffusion algorithm.
+--
+-- Version: 03/21/2017
+--
+-- Copyright 2016-2017 by Samuel Devulder
+--
+-- This program is free software; you can redistribute
+-- it and/or modify it under the terms of the GNU
+-- General Public License as published by the Free
+-- Software Foundation; version 2 of the License.
+-- See 
+
+run('lib/ostro_other.lua')
+
+OtherDither:new{width=256,height=192,clash_size=8}:dither()
diff --git a/share/grafx2/scripts/samples_2.4/picture/thomson/bayer4_mo5.lua b/share/grafx2/scripts/samples_2.4/picture/thomson/bayer4_mo5.lua
new file mode 100644
index 00000000..4128b233
--- /dev/null
+++ b/share/grafx2/scripts/samples_2.4/picture/thomson/bayer4_mo5.lua
@@ -0,0 +1,220 @@
+-- bayer4_mo5.lua : converts an image into TO7/70-MO5
+-- mode for thomson machines (MO6,TO8,TO9,TO9+)
+-- using special bayer matrix that fits well with
+-- color clashes.
+--
+-- Version: 02-jan-2017
+--
+-- Copyright 2016-2017 by Samuel Devulder
+--
+-- This program is free software; you can redistribute
+-- it and/or modify it under the terms of the GNU
+-- General Public License as published by the Free
+-- Software Foundation; version 2 of the License.
+-- See 
+
+-- get screen size
+local screen_w, screen_h = getpicturesize()
+
+run("lib/thomson.lua")
+run("lib/color.lua")
+run("lib/bayer.lua")
+
+-- Converts thomson coordinates (0-319,0-199) into screen coordinates
+local function thom2screen(x,y)
+	local i,j;
+	if screen_w/screen_h < 1.6 then
+		i = x*screen_h/200
+		j = y*screen_h/200
+	else
+		i = x*screen_w/320
+		j = y*screen_w/320
+	end
+	return math.floor(i), math.floor(j)
+end
+
+-- return the pixel @(x,y) in normalized linear space (0-1)
+-- corresonding to the thomson screen (x in 0-319, y in 0-199)
+local function getLinearPixel(x,y)
+	local x1,y1 = thom2screen(x,y)
+	local x2,y2 = thom2screen(x+1,y+1)
+	if x2==x1 then x2=x1+1 end
+	if y2==y1 then y2=y1+1 end
+
+	local p,i,j = Color:new(0,0,0);
+	for i=x1,x2-1 do
+		for j=y1,y2-1 do
+			p:add(getLinearPictureColor(i,j))
+		end
+	end
+	p:div((y2-y1)*(x2-x1)*Color.ONE)
+	
+	return p
+end
+
+local dither = bayer.norm(bayer.double(bayer.double({{1,2},{3,4}})))
+local dx,dy=#dither,#dither[1]
+ 
+-- get thomson palette pixel (linear, 0-1 range)
+local linearPalette = {}
+function linearPalette.get(i)
+	local p = linearPalette[i]
+	if not p then
+		local pal = thomson.palette(i-1)
+		local b=math.floor(pal/256)
+		local g=math.floor(pal/16)%16
+		local r=pal%16
+		p = Color:new(thomson.levels.linear[r+1],
+					  thomson.levels.linear[g+1],
+					  thomson.levels.linear[b+1]):div(Color.ONE)
+		linearPalette[i] = p
+	end
+	return p:clone()
+end
+
+-- distance between two colors
+local distance = {}
+function distance.between(c1,c2)
+	local k = c1..','..c2
+	local d = distance[k]
+	if false and not d then
+		d = linearPalette.get(c1):euclid_dist2(linearPalette.get(c2))
+		distance[k] = d
+	end
+	if not d then
+		local x = linearPalette.get(c1):sub(linearPalette.get(c2))
+		local c,c1,c2,c3=1.8,8,11,8
+		local f = function(c,x) return math.abs(x)*c end
+		d = f(c1,x.r)^c + f(c2,x.g)^c + f(c3,x.b)^c
+		distance[k] = d
+	end
+	return d
+end
+
+-- compute a set of best couples for a given histogram
+local best_couple = {n=0}
+function best_couple.get(h)
+	local k = (((h[1]or 0)*8+(h[2]or 0))*8+(h[3]or 0))*8+(h[4]or 0)
+	.. ',' .. (((h[5]or 0)*8+(h[6]or 0))*8+(h[7]or 0))*8+(h[8]or 0)
+	local best_found = best_couple[k]
+	if not best_found then
+		local dm=1000000
+		for i=1,15 do
+			for j=i+1,16 do
+				local d=0
+				for p,n in pairs(h) do
+					local d1,d2=distance.between(p,i),distance.between(p,j)
+					d = d + n*(d1dm then break; end
+				end
+				if d< dm then dm,best_found=d,{} end
+				if d<=dm then table.insert(best_found, {c1=i,c2=j}) end
+			end
+		end
+	
+		if best_couple.n>10000 then 
+			-- keep memory usage low
+			best_couple = {n=0, get=best_couple.get} 
+		end 
+		best_couple[k] = best_found
+		best_couple.n  = best_couple.n+1
+	end
+	return best_found
+end
+
+-- TO7/70 MO5 mode
+thomson.setMO5()
+
+-- convert picture
+local err1,err2 = {},{}
+local coefs = {0,0.6,0}
+for x=-1,320 do 
+	err1[x] = Color:new(0,0,0) 
+	err2[x] = Color:new(0,0,0) 
+end
+for y = 0,199 do
+	err1,err2 = err2,err1
+	for x=-1,320 do err2[x]:mul(0) end
+	
+	for x = 0,319,8 do
+		local h,q = {},{} -- histo, expected color
+		for z=x,x+7 do
+			local d=dither[1+(y%dx)][1+(z%dx)]
+			local p=getLinearPixel(z,y):add(err1[z])
+			local c=((p.r>d) and 1 or 0) + 
+			        ((p.g>d) and 2 or 0) + 
+					((p.b>d) and 4 or 0) + 1 -- theorical color
+
+			table.insert(q,c)			
+			h[c] = (h[c] or 0)+1
+		end
+		
+		local c1,c2
+		for c,_ in pairs(h) do
+			if c1==nil then c1=c
+			elseif c2==nil then c2=c
+			else c1=nil; break; end
+		end
+		if c1~=nil then
+			c2 = c2 or c1
+		else
+			-- get best possible couples of colors
+			local best_found = best_couple.get(h)
+			if #best_found==1 then
+				c1,c2 = best_found[1].c1,best_found[1].c2
+			else
+				-- keep the best of the best depending on max solvable color clashes
+				function clamp(v) return v<0 and -v or v>1 and v-1 or 0 end
+				local dm=10000000
+				for _,couple in ipairs(best_found) do
+					local d=0
+					for k=1,8 do
+						local q=q[k]
+						local p=distance.between(q,couple.c1)
+
+-- This is my first code in lua, so excuse any bad
+-- coding practice.
+
+-- use a zig zag. If false (recommended value), this gives 
+-- a raster look and feel
+local with_zig_zag = with_zig_zag or false
+
+-- debug: displays histograms
+local debug = false
+
+-- enhance luminosity since our mode divide it by two
+local enhance_lum = enhance_lum or true
+
+-- use fixed levels (default=false, give better result)
+local fixed_levels = fixed_levels or false
+
+-- use void-and-cluster 8x8 matrix (default=false)
+local use_vac = use_vac or false
+
+-- get screen size
+local screen_w, screen_h = getpicturesize()
+
+run("lib/thomson.lua")
+run("lib/color.lua")
+run("lib/bayer.lua")
+
+-- Converts thomson coordinates (0-159,0-99) into screen coordinates
+local function thom2screen(x,y)
+	local i,j;
+	if screen_w/screen_h < 1.6 then
+		i = x*screen_h/100
+		j = y*screen_h/100
+	else
+		i = x*screen_w/160
+		j = y*screen_w/160
+	end
+	return math.floor(i), math.floor(j)
+end
+
+-- return the pixel @(x,y) in linear space corresonding to the thomson screen (x in 0-159, y in 0-99)
+local function getLinearPixel(x,y)
+	local x1,y1 = thom2screen(x,y)
+	local x2,y2 = thom2screen(x+1,y+1)
+	if x2==x1 then x2=x1+1 end
+	if y2==y1 then y2=y1+1 end
+
+	local p,i,j = Color:new(0,0,0);
+	for i=x1,x2-1 do
+		for j=y1,y2-1 do
+			p:add(getLinearPictureColor(i,j))
+		end
+	end
+
+	return p:div((y2-y1)*(x2-x1)):floor()
+end
+
+--[[ make a bayer matrix
+function bayer(matrix)
+	local m,n=#matrix,#matrix[1]
+	local r,i,j = {}
+	for j=1,m*2 do
+		local t = {}
+		for i=1,n*2 do t[i]=0; end
+		r[j] = t;
+	end
+	
+	-- 0 3
+	-- 2 1
+	for j=1,m do
+		for i=1,n do
+			local v = 4*matrix[j][i]
+			r[m*0+j][n*0+i] = v-3
+			r[m*1+j][n*1+i] = v-2
+			r[m*1+j][n*0+i] = v-1
+			r[m*0+j][n*1+i] = v-0
+		end
+	end
+	
+	return r;
+end
+--]]
+
+-- dither matrix
+local dither = bayer.make(4)
+
+if use_vac then
+	-- vac8: looks like FS
+	dither = bayer.norm{
+		{35,57,19,55,7,51,4,21},
+		{29,6,41,27,37,17,59,45},
+		{61,15,53,12,62,25,33,9},
+		{23,39,31,49,2,47,13,43},
+		{3,52,8,22,36,58,20,56},
+		{38,18,60,46,30,5,42,28},
+		{63,26,34,11,64,16,54,10},
+		{14,48,1,44,24,40,32,50}
+	}
+end
+
+-- get color statistics
+local stat = {};
+function stat:clear() 
+	self.r = {}
+	self.g = {}
+	self.b = {}
+	for i=1,16 do self.r[i] = 0; self.g[i] = 0; self.b[i] = 0; end
+end
+function stat:update(px) 
+	local pc2to = thomson.levels.pc2to
+	local r,g,b=pc2to[px.r], pc2to[px.g], pc2to[px.b];
+	self.r[r] = self.r[r] + 1;
+	self.g[g] = self.g[g] + 1;
+	self.b[b] = self.b[b] + 1;
+end
+function stat:coversThr(perc)
+	local function f(stat)
+		local t=-stat[1]
+		for i,n in ipairs(stat) do t=t+n end
+		local thr = t*perc; t=-stat[1]
+		for i,n in ipairs(stat) do 
+			t=t+n 
+			if t>=thr then return i end
+		end
+		return 0
+	end
+	return f(self.r),f(self.g),f(self.b)
+end
+stat:clear();
+for y = 0,99 do
+	for x = 0,159 do
+		stat:update(getLinearPixel(x,y))
+	end
+	thomson.info("Collecting stats...",y,"%")
+end
+
+-- enhance luminosity since our mode divide it by two
+local gain = 1
+if enhance_lum then
+	-- findout level that covers 98% of all non-black pixels
+	local max = math.max(stat:coversThr(.98))
+
+	gain = math.min(2,255/thomson.levels.linear[max])
+
+	if gain>1 then
+		-- redo stat with enhanced levels
+		-- messagebox('gain '..gain..' '..table.concat({stat:coversThr(.98)},','))
+		stat:clear();
+		for y = 0,99 do
+			for x = 0,159 do
+				stat:update(getLinearPixel(x,y):mul(gain):floor())
+			end
+			thomson.info("Enhancing levels..",y,"%")
+		end
+	end
+end
+
+-- find regularly spaced levels in thomson space
+local levels = {}
+function levels.compute(name, stat, num)
+	local tot, max = -stat[1],0;
+	for _,t in ipairs(stat) do
+		max = math.max(t,max)
+		tot = tot + t
+	end	
+	local acc,full=-stat[1],0
+	for i,t in ipairs(stat) do
+		acc = acc + t
+		if acc>tot*.98 then
+			full=thomson.levels.linear[i]
+			break
+		end
+	end	
+	-- sanity
+	if fixed_levels or full==0 then full=255 end
+	local res = {1}; num = num-1
+	for i=1,num do
+		local p = math.floor(full*i/num)
+		local q = thomson.levels.linear2to[p]
+		if q==res[i] and q<16 then q=q+1 end			
+		if not fixed_levels and ires[i]+1 and stat[q-1]>stat[q] then q=q-1 end
+			if q>res[i]+1 and stat[q-1]>stat[q] then q=q-1 end
+			-- 3 corrections? no need...
+			-- if q>res[i]+1 and stat[q-1]>stat[q] then q=q-1 end
+		end
+		res[1+i] = q
+	end
+	
+	-- debug
+	if debug then
+		local txt = ""
+		for _,i in ipairs(res) do
+			txt = txt .. i .. " "
+		end
+		for i,t in ipairs(stat) do
+			txt = txt .. "\n" .. string.format("%s%2d:%3d%% ", name, i, math.floor(100*t/(tot+stat[1]))) .. string.rep('X', math.floor(23*t/max))
+		end
+		messagebox(txt)
+	end
+	
+	return res
+end
+function levels.computeAll(stat)
+	levels.grn = levels.compute("GRN",stat.g,5)
+	levels.red = levels.compute("RED",stat.r,4)
+	levels.blu = levels.compute("BLU",stat.b,3)
+end
+levels.computeAll(stat)
+
+-- put a pixel at (x,y) with dithering
+local function pset(x,y,px)
+	local thr = dither[1+(y % #dither)][1+(x % #dither[1])]
+	local function dither(val,thr,lvls)
+		local i=#lvls
+		local a,b = thomson.levels.linear[lvls[i]],1e30
+		while i>1 and val=thr*(b-a) and 0 or -1)
+	end
+
+	local r = dither(px.r, thr, levels.red);
+	local g = dither(px.g, thr, levels.grn);
+	local b = dither(px.b, thr, levels.blu);
+	
+	local i = r + b*4
+	local j = g==0 and 0 or (11 + g)
+	
+	if with_zig_zag and x%2==1 then
+		thomson.pset(x,y*2+0,j)
+		thomson.pset(x,y*2+1,i)
+	else
+		thomson.pset(x,y*2+0,i)
+		thomson.pset(x,y*2+1,j)
+	end
+end
+
+-- BM16 mode
+thomson.setBM16()
+
+-- define palette
+for i=0,15 do
+	local r,g,b=0,0,0
+	if i<12 then
+		-- r = bit32.band(i,3)
+		-- b = bit32.rshift(i,2)
+		b = math.floor(i/4)
+		r = i-4*b
+	else
+		g = i-11
+	end
+	r,g,b=levels.red[r+1],levels.grn[g+1],levels.blu[b+1]
+	thomson.palette(i,b*256+g*16+r-273)
+end
+
+-- convert picture
+for y = 0,99 do
+	for x = 0,159 do
+		pset(x,y, getLinearPixel(x,y):mul(gain):floor())
+	end
+	thomson.info("Converting...",y,"%")
+end
+
+-- refresh screen
+setpicturesize(320,200)
+thomson.updatescreen()
+finalizepicture()
+
+-- save picture
+do
+	local function exist(file)
+		local f=io.open(file,'rb')
+		if not f then return false else io.close(f); return true; end
+	end
+	local name,path = getfilename()
+	local mapname = string.gsub(name,"%.%w*$","") .. ".map"
+	local fullname = path .. '/' .. mapname
+	-- fullname = 'D:/tmp/toto.map'
+	local ok = not exist(fullname)
+	if not ok then
+		selectbox("Ovr " .. mapname .. "?", "Yes", function() ok = true; end, "No", function() ok = false; end)
+	end
+	if ok then thomson.savep(fullname) end
+end
diff --git a/share/grafx2/scripts/samples_2.4/picture/thomson/lib/bayer.lua b/share/grafx2/scripts/samples_2.4/picture/thomson/lib/bayer.lua
new file mode 100644
index 00000000..55108aab
--- /dev/null
+++ b/share/grafx2/scripts/samples_2.4/picture/thomson/lib/bayer.lua
@@ -0,0 +1,68 @@
+-- bayer.lua : bayer matrix suppport.
+--
+-- Version: 02-jan-2017
+--
+-- Copyright 2016-2017 by Samuel Devulder
+--
+-- This program is free software; you can redistribute
+-- it and/or modify it under the terms of the GNU
+-- General Public License as published by the Free
+-- Software Foundation; version 2 of the License.
+-- See 
+
+if not bayer then
+	bayer = {}
+
+	-- doubles a matrix rows and columns
+	function bayer.double(matrix)
+		local m,n=#matrix,#matrix[1]
+		local r = {}
+		for j=1,m*2 do
+			local t = {}
+			for i=1,n*2 do t[i]=0; end
+			r[j] = t;
+		end
+		
+		-- 0 3
+		-- 2 1
+		for j=1,m do
+			for i=1,n do
+				local v = 4*matrix[j][i]
+				r[m*0+j][n*0+i] = v-3
+				r[m*1+j][n*1+i] = v-2
+				r[m*1+j][n*0+i] = v-1
+				r[m*0+j][n*1+i] = v-0
+			end
+		end
+		
+		return r;
+	end
+
+	-- returns a version of the matrix normalized into  
+	-- the 0-1 range
+	function bayer.norm(matrix)
+		local m,n=#matrix,#matrix[1]
+		local max,ret = 0,{}
+		for j=1,m do
+			for i=1,n do
+				max = math.max(max,matrix[j][i])
+			end
+		end
+		-- max=max+1
+		for j=1,m do
+			ret[j] = {}
+			for i=1,n do
+				ret[j][i]=matrix[j][i]/max
+			end
+		end
+		return ret
+	end
+
+	-- returns a normalized order-n bayer matrix
+	function bayer.make(n)
+		local m = {{1}}
+		while n>1 do n,m = n/2,bayer.double(m) end
+		return bayer.norm(m)
+	end
+
+end -- Bayer
\ No newline at end of file
diff --git a/share/grafx2/scripts/samples_2.4/picture/thomson/lib/color.lua b/share/grafx2/scripts/samples_2.4/picture/thomson/lib/color.lua
new file mode 100644
index 00000000..e5397813
--- /dev/null
+++ b/share/grafx2/scripts/samples_2.4/picture/thomson/lib/color.lua
@@ -0,0 +1,345 @@
+-- color.lua : a color class capable of representing
+-- and manipulating colors in PC-space (gamma=2.2) or
+-- in linear space (gamma=1).
+--
+-- Version: 02-jan-2017
+--
+-- Copyright 2016-2017 by Samuel Devulder
+--
+-- This program is free software; you can redistribute
+-- it and/or modify it under the terms of the GNU
+-- General Public License as published by the Free
+-- Software Foundation; version 2 of the License.
+-- See 
+
+
+if not Color then
+	Color = {ONE=255,NORMALIZE=.005}
+	function Color:new(r,g,b)
+		local o = {};
+		o.r = type(r)=='number' and r or r and r.r or 0;
+		o.g = type(g)=='number' and g or r and r.g or 0;
+		o.b = type(b)=='number' and b or r and r.b or 0;
+		setmetatable(o, self)
+		self.__index = self
+		return o
+	end
+	Color.black = Color:new(0,0,0)
+	
+	function Color.clamp(v,...)
+		if v then 
+			return v<0         and 0 or 
+				   v>Color.ONE and Color.ONE or
+				   v,Color.clamp(...)
+		end
+	end
+	
+	function Color:clone()
+		return Color:new(self.r, self.g, self.b)
+	end
+
+	function Color:tostring()
+		return "(r=" .. self.r .. " g=" .. self.g .. " b=" .. self.b .. ")"
+	end
+	
+	function Color:HSV()
+		local max=math.floor(.5+math.max(self.r,self.g,self.b))
+		local min=math.floor(.5+math.min(self.r,self.g,self.b))
+		
+		local H=(max<=min and 0 or
+			     max<=self.r and (self.g-self.b)/(max-min)+6 or
+			     max<=self.g and (self.b-self.r)/(max-min)+2 or
+			 	 max<=self.b and (self.r-self.g)/(max-min)+4)/6 % 1.0
+		local S=(max==0 or max<=min) and 0 or 1-min/max
+		local V=max/Color.ONE
+		
+		return H,S,V
+	end
+	
+	function Color:intensity()
+		return .3*self.r + .59*self.g + .11*self.b
+	end
+
+	function Color:mul(val)
+		self.r = self.r * val;
+		self.g = self.g * val;
+		self.b = self.b * val;
+		return self;
+	end
+
+	function Color:div(val)
+		return self:mul(1/val);
+	end
+
+	function Color:add(other)
+		self.r = self.r + other.r;
+		self.g = self.g + other.g;
+		self.b = self.b + other.b;
+		return self;
+	end
+	
+	function Color:sub(other)
+		self.r = self.r - other.r;
+		self.g = self.g - other.g;
+		self.b = self.b - other.b;
+		return self;
+	end
+	
+	function Color:dist2(other)
+		return self:euclid_dist2(other)
+		-- return Color.dE2000(self,other)^2
+		-- return Color.dE2fast(self,other)
+	end
+	
+	function Color:euclid_dist2(other)
+		return (self.r - other.r)^2 + 
+		       (self.g - other.g)^2 + 
+			   (self.b - other.b)^2
+	end
+	
+	function Color:floor() 
+		self.r = math.min(math.floor(self.r),Color.ONE);
+		self.g = math.min(math.floor(self.g),Color.ONE);
+		self.b = math.min(math.floor(self.b),Color.ONE);
+		return self;
+	end
+
+	function Color:toPC() 
+		local function f(val) 
+			val = val/Color.ONE
+			-- if val<=0.018 then val = 4.5*val; else val = 1.099*(val ^ (1/2.2))-0.099; end
+			
+			-- works much metter: https://fr.wikipedia.org/wiki/SRGB 
+			if val<=0.0031308 then val=12.92*val else val = 1.055*(val ^ (1/2.4))-0.055 end
+			return val*Color.ONE
+		end;
+		self.r = f(self.r);
+		self.g = f(self.g);
+		self.b = f(self.b);
+		return self;
+	end
+
+	function Color:toLinear() 
+		local function f(val) 
+			val = val/Color.ONE
+			-- if val<=0.081 then val = val/4.5; else val = ((val+0.099)/1.099)^2.2; end
+			
+			-- works much metter: https://fr.wikipedia.org/wiki/SRGB#Transformation_inverse
+			if val<=0.04045 then val = val/12.92 else val = ((val+0.055)/1.055)^2.4 end
+			return val*Color.ONE
+		end;
+		self.r = f(self.r);
+		self.g = f(self.g);
+		self.b = f(self.b);
+		return self;
+	end
+
+	function Color:toRGB()
+		return self.r, self.g, self.b
+	end
+
+	-- return the Color @(x,y) on the original screen in linear space
+	local screen_w, screen_h, _getLinearPictureColor = getpicturesize()
+	function getLinearPictureColor(x,y) 
+		if _getLinearPictureColor==nil then
+			_getLinearPictureColor = {}
+			for i=0,255 do _getLinearPictureColor[i] = Color:new(getbackupcolor(i)):toLinear(); end
+			if Color.NORMALIZE>0 then
+				local histo = {}
+				for i=0,255 do histo[i] = 0 end
+				for y=0,screen_h-1 do
+					for x=0,screen_w-1 do
+						local r,g,b = getbackupcolor(getbackuppixel(x,y))
+						histo[r] = histo[r]+1
+						histo[g] = histo[g]+1
+						histo[b] = histo[b]+1
+					end
+				end
+				local acc,thr=0,Color.NORMALIZE*screen_h*screen_w*3
+				local max
+				for i=255,0,-1 do
+					acc = acc + histo[i]
+					if not max and acc>=thr then
+						max = Color:new(i,i,i):toLinear().r
+					end
+				end
+				for _,c in ipairs(_getLinearPictureColor) do
+					c:mul(Color.ONE/max)
+					c.r,c.g,c.b = Color.clamp(c.r,c.g,c.b)
+				end
+			end
+		end
+		return (x<0 or y<0 or x>=screen_w or y>=screen_h) and Color.black or _getLinearPictureColor[getbackuppixel(x,y)]
+	end
+	
+	-- http://www.brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html
+	function Color.RGBtoXYZ(R,G,B)
+		return 0.4887180*R +0.3106803*G +0.2006017*B,
+			   0.1762044*R +0.8129847*G +0.0108109*B,
+                            0.0102048*G +0.9897952*B
+	end
+	
+	function Color.XYZtoRGB(X,Y,Z)
+		return   2.3706743*X -0.9000405*Y -0.4706338*Z,
+                -0.5138850*X +1.4253036*Y +0.0885814*Z,
+                 0.0052982*X -0.0146949*Y +1.0093968*Z
+	end
+	
+	-- https://fr.wikipedia.org/wiki/CIE_L*a*b*
+	function Color.XYZtoCIELab(X,Y,Z)
+		local function f(t)
+			return t>0.00885645167 and t^(1/3)
+			                        or  7.78703703704*t+0.13793103448
+		end
+		X,Y,Z=X/Color.ONE,Y/Color.ONE,Z/Color.ONE
+		return 116*f(Y)-16,
+			   500*(f(X)-f(Y)),
+			   200*(f(Y)-f(Z))
+	end
+	function Color.CIEALabtoXYZ(L,a,b)
+		local function f(t)
+			return t>0.20689655172 and t^3
+			                        or  0.12841854934*(t-0.13793103448)
+		end
+		local l=(L+16)/116
+		return Color.ONE*f(l),
+		       Color.ONE*f(l+a/500),
+			   Color.ONE*f(l-b/200)
+	end
+	function Color:toLab()
+		return Color.XYZtoCIELab(Color.RGBtoXYZ(self:toRGB()))
+	end
+	
+	-- http://www.brucelindbloom.com/Eqn_DeltaE_CIE2000.html
+	function Color.dE1976(col1,col2)
+		local L1,a1,b1 = col1:toLab()
+		local L2,a2,b2 = col2:toLab()
+		return ((L1-L2)^2+(a1-a2)^2+(b1-b2)^2)^.5
+	end
+	function Color.dE1994(col1,col2)
+		local L1,a1,b1 = col1:toLab()
+		local L2,a2,b2 = col2:toLab()
+		
+		local k1,k2 = 0.045,0.015
+		local kL,kC,kH = 1,1,1
+		
+		local c1 = (a1^2 + b1^2)^.5
+		local c2 = (a2^2 + b2^2)^.5
+		
+		local dA = a1 - a2
+		local dB = b1 - b2
+		local dC = c1 - c2
+		
+		local dH2 = dA^2 + dB^2 - dC^2
+		local dH = dH2>0 and dH2^.5 or 0
+		local dL = L1 - L2
+		
+		local sL = 1
+		local sC = 1 + k1*c1
+		local sH = 1 + k2*c1
+		
+		local vL = dL/(kL*sL)
+		local vC = dC/(kC*sC)
+		local vH = dH/(kH*sH)
+		
+		return (vL^2 + vC^2 + vH^2)^.5
+	end
+	-- http://www.color.org/events/colorimetry/Melgosa_CIEDE2000_Workshop-July4.pdf
+	-- https://en.wikipedia.org/wiki/Color_difference#CIEDE2000
+	function Color.dE2000(col1,col2)
+		local L1,a1,b1 = col1:toLab()
+		local L2,a2,b2 = col2:toLab()
+		
+		local kL,kC,kH = 1,1,1
+
+		local l_p = (L1 + L2)/2
+		
+		function sqrt(x)
+			return x^.5
+		end
+		function norm(x,y)
+			return sqrt(x^2+y^2)
+		end
+		function mean(x,y)
+			return (x+y)/2
+		end
+		local function atan2(a,b)
+			local t=math.atan2(a,b)*180/math.pi
+			return t<0 and t+360 or t
+		end
+		local function sin(x)
+			return math.sin(x*math.pi/180)
+		end
+		local function cos(x)
+			return math.cos(x*math.pi/180)
+		end
+
+		local c1  = norm(a1,b1)
+		local c2  = norm(a2,b2)
+		local c_  = mean(c1,c2)
+		
+		local G   = 0.5*(1-sqrt(c_^7/(c_^7+25^7)))
+		local a1p = a1*(1+G)
+		local a2p = a2*(1+G)
+		
+		local c1p = norm(a1p,b1)
+		local c2p = norm(a2p,b2)
+		local c_p = mean(c1p,c2p)
+		
+		local h1p = atan2(b1,a1p)
+		local h2p = atan2(b2,a2p)
+		
+		local h_p = mean(h1p,h2p) + 
+		            (math.abs(h1p - h2p)<=180 and 0 or
+					  h1p+h2p<360 and 180 or -180)
+		
+		local T   = 1 -
+				    0.17 * cos(    h_p - 30) +
+			        0.24 * cos(2 * h_p     ) +
+			        0.32 * cos(3 * h_p +  6) -
+			        0.20 * cos(4 * h_p - 63)
+
+		local dhp = h2p - h1p + (math.abs(h1p - h2p)<=180 and 0 or
+		                                         h2p<=h1p and 360 or 
+												             -360)
+		local dLp = L2 - L1
+		local dCp = c2p - c1p
+		local dHp = 2*sqrt(c1p*c2p)*sin(dhp/2)
+		
+	
+		local sL = 1 + 0.015*(l_p - 50)^2/sqrt(20+(l_p-50)^2)
+		local sC = 1 + 0.045*c_p
+		local sH = 1 + 0.015*c_p*T
+		
+		local d0 = 30*math.exp(-((h_p-275)/25)^2)
+		
+		local rC = 2*sqrt(c_p^7/(c_p^7+25^7))
+		local rT = -rC * sin(2*d0)
+		
+		return sqrt( (dLp / (kL*sL))^2 + 
+		             (dCp / (kC*sC))^2 +
+			         (dHp / (kH*sH))^2 +
+			         (dCp / (kC*sC))*(dHp / (kH*sH))*rT )
+	end
+	
+	function Color.dE2fast(col1,col2)
+		-- http://www.compuphase.com/cmetric.htm#GAMMA
+		local r1,g1,b1 = Color.clamp(col1:toRGB())
+		local r2,g2,b2 = Color.clamp(col2:toRGB())
+		
+		local rM = (r1+r2)/(Color.ONE*2)
+		
+		return ((r1-r2)^2)*(2+rM) +
+		       ((g1-g2)^2)*(4+1) +
+			   ((b1-b2)^2)*(3-rM)
+	end
+	
+	function Color:hash(M)
+		M=M or 256
+		local m=(M-1)/Color.ONE
+		local function f(x) 
+			return math.floor(.5+(x<0 and 0 or x>Color.ONE and Color.ONE or x)*m) 
+		end
+		return f(self.r)+M*(f(self.g)+M*f(self.b))
+	end
+end -- Color defined
diff --git a/share/grafx2/scripts/samples_2.4/picture/thomson/lib/color_reduction.lua b/share/grafx2/scripts/samples_2.4/picture/thomson/lib/color_reduction.lua
new file mode 100644
index 00000000..1b5604ed
--- /dev/null
+++ b/share/grafx2/scripts/samples_2.4/picture/thomson/lib/color_reduction.lua
@@ -0,0 +1,520 @@
+-- color_reduction.lua : support for reducing the
+-- colors for a thomson image.
+--
+-- Inspire by Xiaolin Wu v2 (Xiaolin Wu 1992).
+-- Greedy orthogonal bipartition of RGB space for
+-- variance minimization aided by inclusion-exclusion
+-- tricks. (Author's description)
+-- http://www.ece.mcmaster.ca/%7Exwu/cq.c
+--
+-- Version: 02-jan-2017
+--
+-- Copyright 2016-2017 by Samuel Devulder
+--
+-- This program is free software; you can redistribute
+-- it and/or modify it under the terms of the GNU
+-- General Public License as published by the Free
+-- Software Foundation; version 2 of the License.
+-- See 
+
+run('color.lua')
+run('bayer.lua')
+run('thomson.lua')
+run('convex_hull.lua')
+
+if not ColorReducer then
+
+-- clamp a value in the 0-255 range
+local function clamp(v) 
+	v=math.floor(v+.5)
+	return v<0 and 0 or v>255 and 255 or v
+end
+
+local Voxel = {}
+
+function Voxel:new() 
+	local o = {m2 = 0, wt=0, mr=0, mg=0, mb=0}
+	setmetatable(o, self)
+	self.__index = self
+	return o
+end
+
+function Voxel:rgb()
+	local n=self.wt; n=n>0 and n or 1
+	return clamp(self.mr/n),
+	       clamp(self.mg/n),
+		   clamp(self.mb/n)
+end
+
+function Voxel:toThomson()
+	local r,g,b=self:rgb()
+	return thomson.levels.linear2to[r]-1,
+		   thomson.levels.linear2to[g]-1,
+		   thomson.levels.linear2to[b]-1
+end
+
+function Voxel:toPal()
+	local r,g,b=self:toThomson()
+	return r+g*16+b*256
+end
+
+function Voxel:tostring()
+	local n=self.wt
+	local r,g,b=self:rgb()
+	return "(n="..math.floor(n*10)/10 .." r=" .. r.. " g="..g .. " b=" .. b.. " rgb=".. table.concat({self:toThomson()},',').. ")"
+end
+
+function Voxel:addColor(color)
+	local r,g,b=color:toRGB()
+	self.wt = self.wt + 1
+	self.mr = self.mr + r
+	self.mg = self.mg + g
+	self.mb = self.mb + b
+	self.m2 = self.m2 + r*r + g*g + b*b
+	return self
+end
+
+function Voxel:add(other,k)
+	k=k or 1
+	self.wt = self.wt + other.wt*k
+	self.mr = self.mr + other.mr*k
+	self.mg = self.mg + other.mg*k
+	self.mb = self.mb + other.mb*k
+	self.m2 = self.m2 + other.m2*k
+	return self
+end
+
+function Voxel:mul(k)
+	return self:add(self,k-1)
+end
+
+function Voxel:module2()
+	return self.mr*self.mr + self.mg*self.mg + self.mb*self.mb
+end
+
+ColorReducer = {}
+
+function ColorReducer:new() 
+	local o = {}
+	setmetatable(o, self)
+	self.__index = self
+	return o
+end
+
+function ColorReducer:v(r,g,b)
+	local i=(r*17+g)*17+b
+	if not self[i] then self[i]=Voxel:new() end
+	return self[i]
+end
+
+function ColorReducer:add(linearColor)
+	local r,g,b=linearColor:toRGB()
+	
+	r,g,b=thomson.levels.linear2to[clamp(r)],
+	      thomson.levels.linear2to[clamp(g)],
+	      thomson.levels.linear2to[clamp(b)]
+	self:v(r,g,b):addColor(linearColor)
+	-- if r==1 and g==1 and b==1 then messagebox(self:v(r,g,b).wt) end
+end
+
+function ColorReducer:M3d()
+	-- convert histogram into moments so that we can
+    -- rapidly calculate the sums of the above quantities
+    -- over any desired box.
+	for r=1,16 do
+		local area={}
+		for i=0,16 do area[i]=Voxel:new() end
+		for g=1,16 do
+			local line=Voxel:new()
+			for b=1,16 do
+				local v = self:v(r,g,b)
+				-- v:mul(0):add(self:v(r-1,g,b)):add(area[b]:add(line:add(v))
+				line:add(v)
+				area[b]:add(line)
+				v:mul(0):add(self:v(r-1,g,b)):add(area[b])
+			end
+		end
+	end
+end
+
+function ColorReducer:Vol(cube)
+	-- Compute sum over a box of all statistics
+	return Voxel:new()
+	       :add(self:v(cube.r1,cube.g1,cube.b1), 1)
+		   :add(self:v(cube.r1,cube.g1,cube.b0),-1)
+		   :add(self:v(cube.r1,cube.g0,cube.b1),-1)
+		   :add(self:v(cube.r1,cube.g0,cube.b0), 1)
+		   :add(self:v(cube.r0,cube.g1,cube.b1),-1)
+		   :add(self:v(cube.r0,cube.g1,cube.b0), 1)
+		   :add(self:v(cube.r0,cube.g0,cube.b1), 1)
+		   :add(self:v(cube.r0,cube.g0,cube.b0),-1)
+end
+
+-- The next two routines allow a slightly more efficient
+-- calculation of Vol() for a proposed subbox of a given
+-- box.  The sum of Top() and Bottom() is the Vol() of a
+-- subbox split in the given direction and with the specified
+-- new upper bound.
+
+function ColorReducer:Bottom(cube,dir)
+	-- Compute part of Vol(cube, mmt) that doesn't
+	-- depend on r1, g1, or b1 (depending on dir)
+	local v=Voxel:new()
+	if dir=="RED" then
+		v:add(self:v(cube.r0,cube.g1,cube.b1),-1)
+		 :add(self:v(cube.r0,cube.g1,cube.b0), 1)
+		 :add(self:v(cube.r0,cube.g0,cube.b1), 1)
+		 :add(self:v(cube.r0,cube.g0,cube.b0),-1)
+	elseif dir=="GREEN" then
+		v:add(self:v(cube.r1,cube.g0,cube.b1),-1)
+		 :add(self:v(cube.r1,cube.g0,cube.b0), 1)
+		 :add(self:v(cube.r0,cube.g0,cube.b1), 1)
+		 :add(self:v(cube.r0,cube.g0,cube.b0),-1)
+	elseif dir=="BLUE" then
+		v:add(self:v(cube.r1,cube.g1,cube.b0),-1)
+		 :add(self:v(cube.r1,cube.g0,cube.b0), 1)
+		 :add(self:v(cube.r0,cube.g1,cube.b0), 1)
+		 :add(self:v(cube.r0,cube.g0,cube.b0),-1)
+	end
+	return v
+end
+
+function ColorReducer:Top(cube,dir,pos)
+	-- Compute remainder of Vol(cube, mmt), substituting
+    -- pos for r1, g1, or b1 (depending on dir)
+	local v=Voxel:new()
+	if dir=="RED" then
+		v:add(self:v(pos,cube.g1,cube.b1), 1)
+		 :add(self:v(pos,cube.g1,cube.b0),-1)
+		 :add(self:v(pos,cube.g0,cube.b1),-1)
+		 :add(self:v(pos,cube.g0,cube.b0), 1)
+	elseif dir=="GREEN" then
+		v:add(self:v(cube.r1,pos,cube.b1), 1)
+		 :add(self:v(cube.r1,pos,cube.b0),-1)
+		 :add(self:v(cube.r0,pos,cube.b1),-1)
+		 :add(self:v(cube.r0,pos,cube.b0), 1)
+	elseif dir=="BLUE" then
+		v:add(self:v(cube.r1,cube.g1,pos), 1)
+		 :add(self:v(cube.r1,cube.g0,pos),-1)
+		 :add(self:v(cube.r0,cube.g1,pos),-1)
+		 :add(self:v(cube.r0,cube.g0,pos), 1)
+	end
+	return v
+end
+
+function ColorReducer:Var(cube)
+	-- Compute the weighted variance of a box
+	-- NB: as with the raw statistics, this is really the variance * size
+	local v = self:Vol(cube)
+	return v.m2 - v:module2()/v.wt
+end
+
+-- We want to minimize the sum of the variances of two subboxes.
+-- The sum(c^2) terms can be ignored since their sum over both subboxes
+-- is the same (the sum for the whole box) no matter where we split.
+-- The remaining terms have a minus sign in the variance formula,
+-- so we drop the minus sign and MAXIMIZE the sum of the two terms.
+
+function ColorReducer:Maximize(cube,dir,first,last,cut,whole)
+	local base = self:Bottom(cube,dir)
+	local max = 0
+	cut[dir] = -1
+	for i=first,last-1 do 
+		local half = Voxel:new():add(base):add(self:Top(cube,dir,i))
+		-- now half is sum over lower half of box, if split at i
+		if half.wt>0 then -- subbox could be empty of pixels!
+			local temp = half:module2()/half.wt
+			half:mul(-1):add(whole)
+			if half.wt>0 then
+				temp = temp + half:module2()/half.wt
+				if temp>max then max=temp; cut[dir] = i end
+			end
+		end
+	end
+	return max
+end
+
+function ColorReducer:Cut(set1,set2)
+	local whole = self:Vol(set1)
+	local cut = {}
+	local maxr = self:Maximize(set1,"RED",  set1.r0+1,set1.r1, cut, whole)
+	local maxg = self:Maximize(set1,"GREEN",set1.g0+1,set1.g1, cut, whole)
+	local maxb = self:Maximize(set1,"BLUE", set1.b0+1,set1.b1, cut, whole)
+	local dir  = "BLUE"
+	if maxr>=maxg and maxr>=maxb then
+		dir = "RED"
+		if cut.RED<0 then return false end -- can't split the box
+	elseif maxg>=maxr and maxg>=maxb then
+		dir = "GREEN"
+	end
+	
+	set2.r1=set1.r1
+	set2.g1=set1.g1
+	set2.b1=set1.b1
+	if dir=="RED" then
+		set1.r1 = cut[dir]
+		set2.r0 = cut[dir]
+		set2.g0 = set1.g0
+		set2.b0 = set1.b0
+	elseif dir=="GREEN" then
+		set1.g1 = cut[dir]
+		set2.g0 = cut[dir]
+		set2.r0 = set1.r0
+		set2.b0 = set1.b0
+	else
+		set1.b1 = cut[dir]
+		set2.b0 = cut[dir]
+		set2.r0 = set1.r0
+		set2.g0 = set1.g0
+	end
+	local function vol(box)
+		local function q(a,b) return (a-b)*(a-b) end
+		return q(box.r1,box.r0) + q(box.g1,box.g0) + q(box.b1,box.b0)
+	end
+	set1.vol = vol(set1)
+	set2.vol = vol(set2)
+	return true
+end
+
+function ColorReducer:boostBorderColors()
+	-- Idea: consider the convex hull of all the colors.
+	-- These color can be mixed to produce any other used
+	-- color, so they are kind of really important. 
+	-- Unfortunately most color-reduction algorithm do not
+	-- retain these color ue to averaging property. The idea
+	-- here is not artifically increase their count so that
+	-- the averaging goes into these colors. 
+
+	-- do return self end -- for testing
+	
+	local hull=ConvexHull:new(function(v)
+		return {v:rgb()}
+	end)
+	
+	-- collect set of points
+	local pts,tot={},0
+	for i=0,17*17*17-1 do
+		local v = self[i]
+		if v then 
+			pts[v] = true
+			tot=tot+v.wt
+		end
+	end
+	
+	-- build convex hull of colors.
+	for v in pairs(pts) do
+		hull:addPoint(v) 
+	end
+	
+	-- collect points near the hull
+	local bdr, hsz, hnb, max = {},0,0,0
+	for v in pairs(pts) do
+		if hull:distToHull(v)>-.1 then
+			bdr[v] = true
+			hnb = hnb+1
+			hsz = hsz+v.wt
+			max = math.max(max,v.wt)
+		end
+	end
+	
+	if tot>hsz then	
+		-- heuristic formula to boost colors of the hull
+		-- not too little, not to much. It might be tuned
+		-- over time, but this version gives satisfying
+		-- result (.51 is important)
+		for v in pairs(bdr) do
+			v:mul(math.min(max,tot-hsz,v.wt*(1+.51*max*hnb/hsz))/v.wt)
+		end
+	end
+	
+	return self
+end
+
+function ColorReducer:buildPalette(max, forceBlack)
+	if self.palette then return self.palette end
+	
+	forceBlack=forceBlack or true
+
+	self:M3d()
+	local function c(r0,g0,b0,r1,g1,b1)
+		return {r0=r0,r1=r1,g0=g0,g1=g1,b0=b0,b1=b1}
+	end
+	local cube = {c(0,0,0,16,16,16)}
+	local n,i = 1,2
+	local vv   = {}
+	while i<=max do
+		cube[i] = c(0,0,0,1,1,1)
+		if forceBlack and i==max then
+			local ko = true;
+			for j=1,max-1 do
+				if self:Vol(cube[j]):toPal()==0 then
+					ko = false
+					break
+				end
+			end
+			if ko then break end -- forcingly add black
+		end
+		if self:Cut(cube[n], cube[i]) then
+			vv[n] = cube[n].vol>1 and self:Var(cube[n]) or 0
+			vv[i] = cube[i].vol>1 and self:Var(cube[i]) or 0
+		else
+			vv[n] = 0
+			cube[i] = nil
+			i=i-1
+		end
+		n = 1; local temp = vv[n]
+		for k=2,i do if vv[k]>temp then temp=vv[k]; n=k; end end
+		if temp<=0 then break end -- not enough color
+		i = i+1
+	end
+
+	-- helper to sort the palette
+	local pal = {}
+	for _,c in ipairs(cube) do
+		local r,g,b=self:Vol(c):toThomson()
+		table.insert(pal, {r=r+1,g=g+1,b=b+1})
+	end
+	-- messagebox(#pal)
+	
+	-- sort the palette in a nice color distribution
+	local function cmp(a,b) 
+		local t=thomson.levels.pc
+		a = Color:new(t[a.r],t[a.g],t[a.b])
+		b = Color:new(t[b.r],t[b.g],t[b.b])
+		local ah,as,av=a:HSV()
+		local bh,bs,bv=b:HSV()
+		as,bs=a:intensity()/255,b:intensity()/255
+		-- function lum(a) return ((.241*a.r + .691*a.g + .068*a.b)/255)^.5 end
+		-- as,bs=lum(a),lum(b)
+		local sat,int=32,256
+		local function quant(ah,as,av)
+			return math.floor(ah*8),
+			       math.floor(as*sat),
+				   math.floor(av*int+.5)
+		end
+		ah,as,av=quant(ah,as,av)
+		bh,bs,bv=quant(bh,bs,bv)
+		-- if true then return ah255 and 255 or x)*m) 
+	end
+	local k=f(linearPixel.r)+M*(f(linearPixel.g)+M*f(linearPixel.b))
+	local c=self[k]
+	if c==nil then
+		local dm=1e30
+		for i,palette in ipairs(self.linear) do
+			local d = palette:dist2(linearColor)
+			if d=b then return v end
+			while v=t and b or a
+		end
+		local t = mat[1+(x%mx)][1+(y%my)]
+		local p=getLinearPixel(x,y)
+		p.r = dith(p.r, t)
+		p.g = dith(p.g, t)
+		p.b = dith(p.b, t)
+		return p
+	end
+	return self:analyze(w,h, function(x,y)
+		return 
+			   dith(x,y)
+		       -- :mul(4):add(getLinearPixel(x,y)):div(5)
+			   -- :add(getLinearPixel(x-1,y))
+			   -- :add(getLinearPixel(x+1,y))
+			   -- :div(3)
+	end, info)
+end
+
+--[[
+function ColorReducer:analyzeBuildWithDither(w,h,max,getLinearPixel,info)
+	if not info then info=function(y) wait(0) end end
+	
+	local dith = ColorReducer:new()
+	dith:analyze(w,h,getLinearPixel,function(y) info(y/2) end)
+	
+	
+	local ostro = OstroDither:new(dith:buildPalette(max),
+	                              thomson.levels.linear, .9)
+	ostro:dither(w,h,
+		function(x,y,xs,err)
+			local p = getLinearPixel(x,y)
+			self:add(err[x]:clone():add(p))
+			self:add(p)
+			return p
+		end, 
+		function(x,y,c)
+			-- self:add(ostro:_linearPalette(c+1))
+		end,true,
+		function(y) info((h+y)/2) end
+	)
+
+	return self:buildPalette(max)
+end
+--]]
+	
+end -- ColorReduction
\ No newline at end of file
diff --git a/share/grafx2/scripts/samples_2.4/picture/thomson/lib/convex_hull.lua b/share/grafx2/scripts/samples_2.4/picture/thomson/lib/convex_hull.lua
new file mode 100644
index 00000000..d4c8e4d8
--- /dev/null
+++ b/share/grafx2/scripts/samples_2.4/picture/thomson/lib/convex_hull.lua
@@ -0,0 +1,219 @@
+-- convxhull.lua : support for computing the convex
+-- hull of a set of points.
+--
+-- inspired from: https://gist.github.com/anonymous/5184ba0bcab21d3dd19781efd3aae543
+--
+-- Version: 02-jan-2017
+--
+-- Copyright 2016-2017 by Samuel Devulder
+--
+-- This program is free software; you can redistribute
+-- it and/or modify it under the terms of the GNU
+-- General Public License as published by the Free
+-- Software Foundation; version 2 of the License.
+-- See 
+
+if not ConvexHull then
+
+local function sub(u,v)
+	return {u[1]-v[1],u[2]-v[2],u[3]-v[3]}
+end
+
+local function mul(k,u)
+	return {k*u[1],k*u[2],k*u[3]}
+end
+
+local function cross(u,v)
+	return {u[2]*v[3] - u[3]*v[2], 
+	        u[3]*v[1] - u[1]*v[3], 
+			u[1]*v[2] - u[2]*v[1]}
+end
+
+local function dot(u,v)
+	return u[1]*v[1] + u[2]*v[2] + u[3]*v[3]
+end
+
+local function unit(u)
+	local d=dot(u,u)
+	return d==0 and u or mul(1/d^.5, u)
+end
+
+ConvexHull = {}
+
+function ConvexHull:new(coordFct) 
+	local o = {
+		points={},
+		coord=coordFct
+	}
+	setmetatable(o, self)
+	self.__index = self
+	return o
+end
+
+function ConvexHull.coord(elt)
+	return {elt[1],elt[2],elt[3]}
+end
+
+function ConvexHull:vect(a,b)
+	return sub(self.coord(b),self.coord(a))
+end
+
+function ConvexHull:normal(face)
+	local u=self:vect(face[1],face[2])
+	local v=self:vect(face[1],face[3])
+	return cross(u,v)
+end
+
+function ConvexHull:printPoint(p)
+	return '('..table.concat(self.coord(p),',')..')'
+end
+
+function ConvexHull:printFace(F)
+	return '['..self:printPoint(F[1])..' '..
+	            self:printPoint(F[2])..' '..
+				self:printPoint(F[3])..']'
+end
+
+function ConvexHull:seen(face,p)
+	local N=self:normal(face)
+	local P=self:vect(face[1],p)
+	return dot(N,P)>=0
+end
+
+function ConvexHull:bdry(faces)
+	local code={n=0}
+	function code.encode(pt,...)
+		if pt then
+			local k = code[pt]
+			if not k then 
+				k = code.n+1
+				code[k]  = pt
+				code[pt] = k
+				code.n   = k
+			end
+			local rest = code.encode(...)
+			return rest and (k..','..rest) or ""..k
+		end
+	end
+	function code.decode(str)
+		local i = str:find(',')
+		if i then
+			local k = str:sub(1,i-1)
+			return code[tonumber(k)],code.decode(str:sub(i+1))
+		else
+			return code[tonumber(str)]
+		end
+	end
+	local set = {}
+	local  function add(...)
+		set[code.encode(...)] = true
+	end
+	local function rem(...)
+		set[code.encode(...)] = nil
+	end
+	local function keys()
+		local r = {}
+		for k in pairs(set) do
+			r[{code.decode(k)}] = true
+		end
+		return r
+	end
+	for F in pairs(faces) do
+		add(F[1],F[2])
+		add(F[2],F[3])
+		add(F[3],F[1])
+	end
+	for F in pairs(faces) do
+		rem(F[1],F[3])
+		rem(F[3],F[2])
+		rem(F[2],F[1])
+	end
+	return keys()
+end
+
+function ConvexHull:addPoint(p)
+	-- first 3 points
+	if self.points then 
+		if p==self.points[1] or p==self.points[2] then return end
+		table.insert(self.points,p)
+		
+		if #self.points==3 then 
+			self.hull={
+				{self.points[1],self.points[2],self.points[3]},
+				{self.points[1],self.points[3],self.points[2]}
+			}
+			self.points=nil
+		end
+	else
+		local seenF,n = {},0
+		for _,F in ipairs(self.hull) do
+			if F[1]==p or F[2]==p or F[3]==p then return end
+			if self:seen(F,p) then seenF[F]=true;n=n+1 end
+		end
+		
+		if n==#self.hull then
+			-- if can see all faces, unsee ones looking "down"
+			local N
+			for F in pairs(seenF) do N=self:normal(F); break; end
+			for F in pairs(seenF) do
+				if dot(self:normal(F),N)<=0 then
+					seenF[F] = nil
+					n=n-1
+				end
+			end
+		end
+		
+		-- remove (old) seen faces
+		local z=#self.hull
+		for i=#self.hull,1,-1 do
+			if seenF[self.hull[i]] then 
+				table.remove(self.hull,i)
+			end
+		end
+
+		-- insert new boundaries with seen faces
+		for E in pairs(self:bdry(seenF)) do
+			table.insert(self.hull,{E[1],E[2],p})
+		end
+	end
+	return self
+end
+
+function ConvexHull:verticesSet()
+	local v = {}
+	if self.hull then
+		for _,F in ipairs(self.hull) do
+			v[F[1]] = true
+			v[F[2]] = true
+			v[F[3]] = true
+		end
+	end
+	return v
+end
+
+function ConvexHull:verticesSize()
+	local n = 0
+	for _ in pairs(self:verticesSet()) do n=n+1 end
+	return n
+end
+
+function ConvexHull:distToFace(F,pt)
+	local N=unit(self:normal(F))
+	local P=self:vect(F[1],pt)
+	return dot(N,P)
+end
+
+function ConvexHull:distToHull(pt)
+	local d
+	for _,F in ipairs(self.hull) do
+		local t = self:distToFace(F,pt)
+		d = d==nil and t or
+			(0<=t and t=t and t>d) and t or
+			d
+		if d==0 then break end
+	end
+	return d
+end
+
+end -- ConvexHull
\ No newline at end of file
diff --git a/share/grafx2/scripts/samples_2.4/picture/thomson/lib/ostromoukhov.lua b/share/grafx2/scripts/samples_2.4/picture/thomson/lib/ostromoukhov.lua
new file mode 100644
index 00000000..bfa6ad4e
--- /dev/null
+++ b/share/grafx2/scripts/samples_2.4/picture/thomson/lib/ostromoukhov.lua
@@ -0,0 +1,629 @@
+-- ostromoukhov.lua : Color dithering using variable
+-- coefficients.
+--
+-- https://liris.cnrs.fr/victor.ostromoukhov/publications/pdf/SIGGRAPH01_varcoeffED.pdf
+--
+-- Version: 02-jan-2017
+--
+-- Copyright 2016-2017 by Samuel Devulder
+--
+-- This program is free software; you can redistribute
+-- it and/or modify it under the terms of the GNU
+-- General Public License as published by the Free
+-- Software Foundation; version 2 of the License.
+-- See 
+
+run('color.lua')
+run('thomson.lua')
+
+if not OstroDither then
+
+OstroDither = {}
+
+local function default_levels() 
+	return {r={0,Color.ONE},g={0,Color.ONE},b={0,Color.ONE}}
+end
+
+function OstroDither:new(palette,attenuation,levels)
+	local o = {
+		attenuation = attenuation or .9, -- works better than 1
+		palette = palette or thomson.default_palette,
+		levels = levels or default_levels(),
+		clash_size = 8 -- for color clash
+	}
+	setmetatable(o, self)
+	self.__index = self
+	return o
+end
+
+function OstroDither:setLevelsFromPalette()
+	local rLevels = {[1]=true,[16]=true}
+	local gLevels = {[1]=true,[16]=true}
+	local bLevels = {[1]=true,[16]=true}
+	local default_palette = true
+	for i,pal in ipairs(self.palette) do
+		local r,g,b=pal%16,math.floor(pal/16)%16,math.floor(pal/256)
+		rLevels[1+r] = true
+		gLevels[1+g] = true
+		bLevels[1+b] = true
+		if pal~=thomson.default_palette[i] then 
+			default_palette = false
+		end
+	end
+	local levels = {r={},g={},b={}}
+	for i,v in ipairs(thomson.levels.linear) do
+		if false then
+			if rLevels[i] and gLevels[i] and bLevels[i] then 
+				table.insert(levels.r, v) 
+				table.insert(levels.g, v) 
+				table.insert(levels.b, v) 
+			end
+		else
+			if rLevels[i] then table.insert(levels.r, v) end
+			if gLevels[i] then table.insert(levels.g, v) end
+			if bLevels[i] then table.insert(levels.b, v) end
+		end
+	end
+	self.levels = levels
+	if default_palette then 
+		self.attenuation = .98
+		self.levels = default_levels()
+	else
+		self.attenuation = .9
+		self.levels = levels
+	end
+end
+
+function OstroDither:_coefs(linearLevel,rgb)
+	if self._ostro==nil then
+		-- original coefs, about to be adapted to the levels
+		local t={ 
+			13,     0,     5,
+			13,     0,     5,
+			21,     0,    10,
+			 7,     0,     4,
+			 8,     0,     5,
+			47,     3,    28,
+			23,     3,    13,
+			15,     3,     8,
+			22,     6,    11,
+			43,    15,    20,
+			 7,     3,     3,
+		   501,   224,   211,
+		   249,   116,   103,
+		   165,    80,    67,
+		   123,    62,    49,
+		   489,   256,   191,
+			81,    44,    31,
+		   483,   272,   181,
+			60,    35,    22,
+			53,    32,    19,
+		   237,   148,    83,
+		   471,   304,   161,
+			 3,     2,     1,
+		   459,   304,   161,
+			38,    25,    14,
+		   453,   296,   175,
+		   225,   146,    91,
+		   149,    96,    63,
+		   111,    71,    49,
+			63,    40,    29,
+			73,    46,    35,
+		   435,   272,   217,
+		   108,    67,    56,
+			13,     8,     7,
+		   213,   130,   119,
+		   423,   256,   245,
+			 5,     3,     3,
+		   281,   173,   162,
+		   141,    89,    78,
+		   283,   183,   150,
+			71,    47,    36,
+		   285,   193,   138,
+			13,     9,     6,
+			41,    29,    18,
+			36,    26,    15,
+		   289,   213,   114,
+		   145,   109,    54,
+		   291,   223,   102,
+			73,    57,    24,
+		   293,   233,    90,
+			21,    17,     6,
+		   295,   243,    78,
+			37,    31,     9,
+			27,    23,     6,
+		   149,   129,    30,
+		   299,   263,    54,
+			75,    67,    12,
+			43,    39,     6,
+		   151,   139,    18,
+		   303,   283,    30,
+			38,    36,     3,
+		   305,   293,    18,
+		   153,   149,     6,
+		   307,   303,     6,
+			 1,     1,     0,
+		   101,   105,     2,
+			49,    53,     2,
+			95,   107,     6,
+			23,    27,     2,
+			89,   109,    10,
+			43,    55,     6,
+			83,   111,    14,
+			 5,     7,     1,
+		   172,   181,    37,
+			97,    76,    22,
+			72,    41,    17,
+		   119,    47,    29,
+			 4,     1,     1,
+			 4,     1,     1,
+			 4,     1,     1,
+			 4,     1,     1,
+			 4,     1,     1,
+			 4,     1,     1,
+			 4,     1,     1,
+			 4,     1,     1,
+			 4,     1,     1,
+			65,    18,    17,
+			95,    29,    26,
+		   185,    62,    53,
+			30,    11,     9,
+			35,    14,    11,
+			85,    37,    28,
+			55,    26,    19,
+			80,    41,    29,
+		   155,    86,    59,
+			 5,     3,     2,
+			 5,     3,     2,
+			 5,     3,     2,
+			 5,     3,     2,
+			 5,     3,     2,
+			 5,     3,     2,
+			 5,     3,     2,
+			 5,     3,     2,
+			 5,     3,     2,
+			 5,     3,     2,
+			 5,     3,     2,
+			 5,     3,     2,
+			 5,     3,     2,
+		   305,   176,   119,
+		   155,    86,    59,
+		   105,    56,    39,
+			80,    41,    29,
+			65,    32,    23,
+			55,    26,    19,
+		   335,   152,   113,
+			85,    37,    28,
+		   115,    48,    37,
+			35,    14,    11,
+		   355,   136,   109,
+			30,    11,     9,
+		   365,   128,   107,
+		   185,    62,    53,
+			25,     8,     7,
+			95,    29,    26,
+		   385,   112,   103,
+			65,    18,    17,
+		   395,   104,   101,
+			 4,     1,     1,
+			 4,     1,     1,
+		   395,   104,   101,
+			65,    18,    17,
+		   385,   112,   103,
+			95,    29,    26,
+			25,     8,     7,
+		   185,    62,    53,
+		   365,   128,   107,
+			30,    11,     9,
+		   355,   136,   109,
+			35,    14,    11,
+		   115,    48,    37,
+			85,    37,    28,
+		   335,   152,   113,
+			55,    26,    19,
+			65,    32,    23,
+			80,    41,    29,
+		   105,    56,    39,
+		   155,    86,    59,
+		   305,   176,   119,
+			 5,     3,     2,
+			 5,     3,     2,
+			 5,     3,     2,
+			 5,     3,     2,
+			 5,     3,     2,
+			 5,     3,     2,
+			 5,     3,     2,
+			 5,     3,     2,
+			 5,     3,     2,
+			 5,     3,     2,
+			 5,     3,     2,
+			 5,     3,     2,
+			 5,     3,     2,
+		   155,    86,    59,
+			80,    41,    29,
+			55,    26,    19,
+			85,    37,    28,
+			35,    14,    11,
+			30,    11,     9,
+		   185,    62,    53,
+			95,    29,    26,
+			65,    18,    17,
+			 4,     1,     1,
+			 4,     1,     1,
+			 4,     1,     1,
+			 4,     1,     1,
+			 4,     1,     1,
+			 4,     1,     1,
+			 4,     1,     1,
+			 4,     1,     1,
+			 4,     1,     1,
+		   119,    47,    29,
+			72,    41,    17,
+			97,    76,    22,
+		   172,   181,    37,
+			 5,     7,     1,
+			83,   111,    14,
+			43,    55,     6,
+			89,   109,    10,
+			23,    27,     2,
+			95,   107,     6,
+			49,    53,     2,
+		   101,   105,     2,
+			 1,     1,     0,
+		   307,   303,     6,
+		   153,   149,     6,
+		   305,   293,    18,
+			38,    36,     3,
+		   303,   283,    30,
+		   151,   139,    18,
+			43,    39,     6,
+			75,    67,    12,
+		   299,   263,    54,
+		   149,   129,    30,
+			27,    23,     6,
+			37,    31,     9,
+		   295,   243,    78,
+			21,    17,     6,
+		   293,   233,    90,
+			73,    57,    24,
+		   291,   223,   102,
+		   145,   109,    54,
+		   289,   213,   114,
+			36,    26,    15,
+			41,    29,    18,
+			13,     9,     6,
+		   285,   193,   138,
+			71,    47,    36,
+		   283,   183,   150,
+		   141,    89,    78,
+		   281,   173,   162,
+			 5,     3,     3,
+		   423,   256,   245,
+		   213,   130,   119,
+			13,     8,     7,
+		   108,    67,    56,
+		   435,   272,   217,
+			73,    46,    35,
+			63,    40,    29,
+		   111,    71,    49,
+		   149,    96,    63,
+		   225,   146,    91,
+		   453,   296,   175,
+			38,    25,    14,
+		   459,   304,   161,
+			 3,     2,     1,
+		   471,   304,   161,
+		   237,   148,    83,
+			53,    32,    19,
+			60,    35,    22,
+		   483,   272,   181,
+			81,    44,    31,
+		   489,   256,   191,
+		   123,    62,    49,
+		   165,    80,    67,
+		   249,   116,   103,
+		   501,   224,   211,
+			 7,     3,     3,
+			43,    15,    20,
+			22,     6,    11,
+			15,     3,     8,
+			23,     3,    13,
+			47,     3,    28,
+			 8,     0,     5,
+			 7,     0,     4,
+			21,     0,    10,
+			13,     0,     5,
+			13,     0,     5}
+		local function process(tab)
+			local tab2={}
+			local function add(i)
+				i=3*math.floor(i+.5)
+				local c0,c1,c2=t[i+1],t[i+2],t[i+3]
+				local norm=self.attenuation/(c0+c1+c2)
+				table.insert(tab2,c0*norm)
+				table.insert(tab2,c1*norm)
+				table.insert(tab2,c2*norm)
+			end
+			local function level(i)
+				return tab[i]*255/Color.ONE
+			end
+			local a,b,j=level(1),level(2),3
+			for i=0,255 do
+				if i>b then a,b,j=b,level(j),j+1; end
+				add(255*(i-a)/(b-a))
+			end
+			return tab2
+		end
+		self._ostro = {r=process(self.levels.r),
+		               g=process(self.levels.g),
+					   b=process(self.levels.b)}
+	end
+	local i = math.floor(linearLevel[rgb]*255/Color.ONE+.5)
+	i = 3*(i<0 and 0 or i>255 and 255 or i)
+	return self._ostro[rgb][i+1],self._ostro[rgb][i+2],self._ostro[rgb][i+3]
+end
+
+function OstroDither:_linearPalette(colorIndex)
+	if self._linear==nil then
+		self._linear = {}
+		local t=thomson.levels.linear
+		for i,pal in ipairs(self.palette) do
+			local r,g,b=pal%16,math.floor(pal/16)%16,math.floor(pal/256)
+			self._linear[i] = Color:new(t[1+r],t[1+g],t[1+b])
+		end
+	end
+	return self._linear[colorIndex]
+end
+
+function OstroDither:getColorIndex(linearPixel)
+	local k=linearPixel:hash(64)
+	local c=self[k]
+	if c==nil then
+		local dm=1e30
+		for i=1,#self.palette do
+			local d = self:_linearPalette(i):dist2(linearPixel)
+			if d M and  M or a
+		end
+		local c0,c1,c2=self:_coefs(linearColor,rgb)
+		if err0 and c0>0 then err0[rgb] = f(err0[rgb],c0) end
+		if err1 and c1>0 then err1[rgb] = f(err1[rgb],c1) end
+		if err2 and c2>0 then err2[rgb] = f(err2[rgb],c2) end
+	end
+	d("r"); d("g"); d("b")
+	
+	return c
+end
+
+function OstroDither:dither(screen_w,screen_h,getLinearPixel,pset,serpentine,info)
+	if not info then info = function(y) thomson.info() end end
+	if not serpentine then serpentine = true end
+	
+	local err1,err2 = {},{}
+	for x=-1,screen_w do 
+		err1[x] = Color:new(0,0,0) 
+		err2[x] = Color:new(0,0,0) 
+	end
+
+	for y=0,screen_h-1 do
+		-- permute error buffers
+		err1,err2 = err2,err1
+		-- clear current-row's buffer
+		for i=-1,screen_w do err2[i]:mul(0) end
+		
+		local x0,x1,xs=0,screen_w-1,1
+		if serpentine and y%2==1 then x0,x1,xs=x1,x0,-xs end
+		
+		for x=x0,x1,xs do
+			local p = getLinearPixel(x,y,xs,err1)
+			local c = self:_diffuse(p,err1[x],err1[x+xs],
+			                          err2[x-xs],err2[x])
+			pset(x,y,c-1)
+		end
+		info(y)
+	end
+end
+
+function OstroDither:ccAcceptCouple(c1,c2)
+	return c1~=c2
+end
+
+function OstroDither:ccDither(screen_w,screen_h,getLinearPixel,pset,serpentine,info) -- dither with color clash
+	local c1,c2	
+	self.getColorIndex = function(self,p)
+		return p:dist2(self:_linearPalette(c1))b.n or a.n==b.n and a.cdm then break end
+					end
+					return d
+				else
+					return dm
+				end
+			end
+			dm=eval()
+		
+			if histo:num(1)>=self.clash_size/2+1 then
+				local z=c2
+				for i=1,#self.palette do c2=i
+					local d=eval()
+					if d0 and 0 or self.clash_size-1) then
+			findC1C2(x,y,xs,err1)
+		end
+		return getLinearPixel(x,y)
+	end
+
+	self:dither(screen_w,screen_h,_getLinearPixel,_pset,serpentine,info)
+end
+
+function OstroDither:dither40cols(getpalette,serpentine) 
+	-- get screen size
+	local screen_w, screen_h = getpicturesize()
+
+	-- Converts thomson coordinates (0-159,0-199) into screen coordinates
+	local function thom2screen(x,y)
+		local i,j;
+		if screen_w/screen_h < 1.6 then
+			i = x*screen_h/200
+			j = y*screen_h/200
+		else
+			i = x*screen_w/320
+			j = y*screen_w/320
+		end
+		return math.floor(i), math.floor(j)
+	end
+
+	-- return the Color @(x,y) in linear space (0-255)
+	-- corresonding to the thomson screen (x in 0-319, 
+	-- y in 0-199)
+	local function getLinearPixel(x,y)
+		local with_cache = true
+		if not self._getLinearPixel then self._getLinearPixel = {} end
+		local k=x+y*thomson.w
+		local p = self._getLinearPixel[k]
+		if not p then
+			local x1,y1 = thom2screen(x,y)
+			local x2,y2 = thom2screen(x+1,y+1)
+			if x2==x1 then x2=x1+1 end
+			if y2==y1 then y2=y1+1 end
+
+			p = Color:new(0,0,0);
+			for j=y1,y2-1 do
+				for i=x1,x2-1 do
+					p:add(getLinearPictureColor(i,j))
+				end
+			end
+			p:div((y2-y1)*(x2-x1)) --:floor()
+						
+			if with_cache then self._getLinearPixel[k]=p end
+		end
+		
+		return with_cache and p:clone() or p
+	end
+		
+	-- MO5 mode
+	thomson.setMO5()
+	self.palette = getpalette(thomson.w,thomson.h,getLinearPixel)
+	
+	-- compute levels from palette
+	self:setLevelsFromPalette()
+	
+	-- convert picture
+	self:ccDither(thomson.w,thomson.h, 
+				  getLinearPixel, thomson.pset, 
+				  serpentine or true, function(y) 
+					thomson.info("Converting...",
+						math.floor(y*100/thomson.h),"%") 
+				  end,true)
+
+	-- refresh screen
+	setpicturesize(thomson.w,thomson.h)
+	thomson.updatescreen()
+	thomson.savep()
+	finalizepicture()
+end
+
+end -- OstroDither
\ No newline at end of file
diff --git a/share/grafx2/scripts/samples_2.4/picture/thomson/lib/thomson.lua b/share/grafx2/scripts/samples_2.4/picture/thomson/lib/thomson.lua
new file mode 100644
index 00000000..a54a4ff4
--- /dev/null
+++ b/share/grafx2/scripts/samples_2.4/picture/thomson/lib/thomson.lua
@@ -0,0 +1,454 @@
+-- thomson.lua : lots of utility for handling
+-- thomson screen.
+--
+-- Version: 02-jan-2017
+--
+-- Copyright 2016-2017 by Samuel Devulder
+--
+-- This program is free software; you can redistribute
+-- it and/or modify it under the terms of the GNU
+-- General Public License as published by the Free
+-- Software Foundation; version 2 of the License.
+-- See 
+
+if not thomson then
+
+run("color.lua") -- optionnal
+
+thomson = {optiMAP=true}
+
+-- RAM banks
+thomson.ramA = {}
+thomson.ramB = {}
+
+function thomson.clear()
+	for i=1,8000 do
+		thomson.ramA[i] = 0
+		thomson.ramB[i] = 0
+	end
+end
+
+-- color levels
+thomson.levels = {
+	-- in pc-space (0-255):
+	pc = {0,100,127,142,163,179,191,203,215,223,231,239,
+		  243,247,251,255}, 
+	-- in linear space (0-255):
+	linear = {},
+	-- maps pc-levels (0-255) to thomson levels (1-16)
+	pc2to={},
+	-- maps linear-levels (0-255) to thomson levels (1-16)
+	linear2to={}
+};
+
+-- pc space to linear space
+local function toLinear(val)
+	-- use the version from Color library
+	if not Color then
+		val = val/255
+		if val<=0.081 then 
+			val = val/4.5; 
+		else 
+			val = ((val+0.099)/1.099)^2.2; 
+		end
+		val = val*255
+		return val;
+	else
+		return Color:new(val,0,0):toLinear().r
+	end
+end
+		
+for i=1,16 do 
+	thomson.levels.linear[i] = toLinear(thomson.levels.pc[i])
+end
+for i=0,255 do 
+	local r,cm,dm;
+	r,cm,dm = toLinear(i),0,1e30
+	for c,v in ipairs(thomson.levels.linear) do
+		local d = math.abs(v-r);
+		if d=0 and i=0 and i 00 02 aa bb
+			if default and num==1 and partial[1]==1 then
+				partial = {0,2,partial[2],car}
+				default = false
+			end
+			-- 00 n xx xx xx 01 bb ==> 00 n+1 xx xx xx bb
+			if default and num==1 and partial[1]==0 and partial[2]<255 then
+				addCarToPartial(car)
+				default = false
+			end
+			-- 00 n xx xx xx 02 bb ==> 00 n+2 xx xx xx bb bb (pas utile mais sert quand combiné ŕ la regle ci-dessus)
+			if default and num==2 and partial[1]==0 and partial[2]<254 then
+				addCarToPartial(car)
+				addCarToPartial(car)
+				default = false
+			end
+		end
+		if default then
+			thomson._append(result, partial)
+			partial = {num,car}
+		end
+		p=p+1
+	end
+	thomson._append(result, partial)
+	return result
+end
+
+-- save a map file corresponging to the current file
+-- if a map file already exist, a confirmation is
+-- prompted to the user
+local function save_current_file()
+	local function exist(file)
+		local f=io.open(file,'rb')
+		if not f then return false else io.close(f); return true; end
+	end
+	local name,path = getfilename()
+	local mapname = string.gsub(name,"%.%w*$","") .. ".map"
+	local fullname = path .. '/' .. mapname
+	local ok = not exist(fullname)
+	if not ok then
+		selectbox("Ovr " .. mapname .. "?", "Yes", function() ok = true; end, "No", function() ok = false; end)
+	end
+	if ok then thomson.savep(fullname) end
+end
+
+-- saves the thomson screen into a MAP file
+function thomson.savep(name)
+	if not name then return save_current_file()	end
+	
+	wait(0) -- allow for key handling
+	local data = thomson._get_map_data()
+	local tmp = {0, math.floor(#data/256), #data%256,0,0}
+	thomson._append(tmp,data,{255,0,0,0,0})
+	local function save(name, buf)
+		local out = io.open(name,"wb")
+		out:write(buf)
+		out:close()
+	end
+	save(name, string.char(unpack(tmp)))
+
+	-- save raw data as well ?
+	local moved, key, mx, my, mb = waitinput(0.01)
+	if key==4123 then -- shift-ESC ==> save raw files as well
+		save(name .. ".rama", string.char(unpack(thomson.ramA)))
+		save(name .. ".ramb", string.char(unpack(thomson.ramB)))
+		local pal = ""
+		for i=0,15 do
+			local val = thomson.palette(i)
+			pal=pal..string.char(math.floor(val/256),val%256)
+		end
+		save(name .. ".pal", pal)
+		messagebox('Saved MAP + RAMA/RAMB/PAL files.')
+	end
+end
+waitbreak(0.01)
+
+function thomson.info(...)
+	local txt = ""
+	for _,t in ipairs({...}) do txt = txt .. t end
+	statusmessage(txt); 
+	if waitbreak(0)==1 then 
+		local ok=false
+		selectbox("Abort ?", "Yes", function() ok = true end, "No", function() ok = false end)
+		if ok then error('Operation aborted') end
+	end
+end
+
+-- copy ramA/B onto GrafX2 screen
+function thomson.updatescreen()
+	-- back out
+	for i=0,255 do
+		setcolor(i,0,0,0)
+	end
+	-- refresh screen content
+	clearpicture(thomson._palette.offset + thomson.border())
+	for y=0,thomson.h-1 do
+		for x=0,thomson.w-1 do
+			local p = thomson.point(x,y)
+			if p<0 then p=-p-1 end
+			thomson._putpixel(x,y,thomson._palette.offset + p)
+		end
+	end
+	-- refresh palette
+	for i=1,thomson._palette.max do
+		local v=thomson._palette[i]
+		local r=v % 16
+		local g=math.floor(v/16)  % 16
+		local b=math.floor(v/256) % 16
+		setcolor(i+thomson._palette.offset-1, 
+				 thomson.levels.pc[r+1], 
+				 thomson.levels.pc[g+1], 
+				 thomson.levels.pc[b+1])
+	end
+	updatescreen()
+end
+
+-- bitmap 16 mode
+function thomson.setBM16() 
+	-- put a pixel onto real screen
+	function thomson._putpixel(x,y,c)
+		putpicturepixel(x*2+0,y,c)
+		putpicturepixel(x*2+1,y,c)
+	end
+	-- put a pixel in thomson screen
+	function thomson.pset(x,y,c)
+		local bank = x%4<2 and thomson.ramA or thomson.ramB
+		local offs = math.floor(x/4)+y*40+1
+		if x%2==0 then
+			bank[offs] = (bank[offs]%16)+c*16
+		else
+			bank[offs] = math.floor(bank[offs]/16)*16+c
+		end
+		-- c=c+thomson._palette.offset
+		-- putpicturepixel(x*2+0,y,c)
+		-- putpicturepixel(x*2+1,y,c)
+	end
+	-- get thomson pixel at (x,y)
+	function thomson.point(x,y)
+		local bank = x%4<2 and thomson.ramA or thomson.ramB
+		local offs = math.floor(x/4)+y*40+1
+		if x%2==0 then
+			return math.floor(bank[offs]/16)
+		else
+			return bank[offs]%16
+		end
+	end
+	-- return internal MAP file
+	function thomson._get_map_data()
+		local tmp = {}
+		for x=1,40 do
+			for y=x,x+7960,40 do
+				table.insert(tmp, thomson.ramA[y])
+			end
+			for y=x,x+7960,40 do
+				table.insert(tmp, thomson.ramB[y])
+			end
+			wait(0) -- allow for key handling
+		end
+		local pal = {}
+		for i=1,16 do 
+			pal[2*i-1] = math.floor(thomson._palette[i]/256)
+			pal[2*i+0] =            thomson._palette[i]%256
+		end
+		-- build data
+		local data={
+			-- BM16
+			0x40,
+			-- ncols-1
+			79,
+			-- nlines-1
+			24
+		};
+		thomson._compress(data, tmp)
+		thomson._append(data,{0,0})
+		-- padd to word
+		if #data%2==1 then table.insert(data,0); end
+		-- tosnap
+		thomson._append(data,{0,128,0,thomson.border(),0,3})
+		thomson._append(data, pal)
+		thomson._append(data,{0xa5,0x5a})
+		return data
+	end
+	
+	thomson.w = 160
+	thomson.h = 200
+	thomson.palette(0,thomson.default_palette)
+	thomson.border(0)
+	thomson.clear()
+end
+
+-- mode MO5 
+function thomson.setMO5() 
+	-- put a pixel onto real screen
+	thomson._putpixel = putpicturepixel
+	-- helpers
+	local function bittst(val,mask) 
+		-- return bit32.btest(val,mask) 
+		return (val % (2*mask))>=mask;
+	end
+	local function bitset(val,mask) 
+		-- return bit32.bor(val, mask)
+		return bittst(val,mask) and val or (val+mask)
+	end
+	local function bitclr(val,mask) 
+		-- return bit32.band(val,255-mask)
+		return bittst(val,mask) and (val-mask) or val
+	end
+	-- put a pixel in thomson screen
+	function thomson.pset(x,y,c)
+		local offs = math.floor(x/8)+y*40+1
+		local mask = 2^(7-(x%8))
+		if c>=0 then
+			thomson.ramB[offs] = (thomson.ramB[offs]%16)+c*16
+			thomson.ramA[offs] = bitset(thomson.ramA[offs],mask)
+		else
+			c=-c-1
+			thomson.ramB[offs] = math.floor(thomson.ramB[offs]/16)*16+c
+			thomson.ramA[offs] = bitclr(thomson.ramA[offs],mask)
+		end
+	end
+	-- get thomson pixel at (x,y)
+	function thomson.point(x,y)
+		local offs = math.floor(x/8)+y*40+1
+		local mask = 2^(7-(x%8))
+		if bittst(thomson.ramA[offs],mask) then
+			return math.floor(thomson.ramB[offs]/16)
+		else
+			return -(thomson.ramB[offs]%16)-1
+		end
+	end
+	-- convert color from MO5 to TO7 (MAP requires TO7 encoding)
+	local function mo5to7(val)
+		-- MO5: DCBA 4321
+		--      __    
+		-- TO7: 4DCB A321
+		local t=((val%16)>=8) and 0 or 128
+		val = math.floor(val/16)*8 + (val%8)
+		val = (val>=64 and val-64 or val+64) + t
+		return val
+	end
+	-- return internal MAP file
+	function thomson._get_map_data()
+		-- create columnwise data
+		local tmpA,tmpB={},{};
+		for x=1,40 do
+			for y=x,x+7960,40 do
+				table.insert(tmpA, thomson.ramA[y])
+				table.insert(tmpB, thomson.ramB[y])
+			end
+			wait(0) -- allow for key handling
+		end
+		if thomson.optiMAP then
+			-- optimize
+			for i=2,8000 do
+				local c1,c2 = math.floor(tmpB[i-0]/16),tmpB[i-0]%16
+				local d1,d2 = math.floor(tmpB[i-1]/16),tmpB[i-1]%16
+				
+				if tmpA[i-1]==255-tmpA[i] or c1==d2 and c2==c1 then
+					tmpA[i] = 255-tmpA[i]
+					tmpB[i] = c2*16+c1
+				elseif tmpA[i]==255 and c1==d1 or tmpA[i]==0 and c2==d2 then
+					tmpB[i] = tmpB[i-1]
+				end
+			end
+		else
+			for i=1,8000 do
+				local c1,c2 = math.floor(tmpB[i]/16),tmpB[i]%16
+				
+				if tmpA[i]==255 or c1
+
+run('lib/ostromoukhov.lua')
+
+local dith=OstroDither:new()
+local tmp=dith.setLevelsFromPalette
+dith.setLevelsFromPalette = function(self)
+	tmp(self)
+	self.attenuation=0
+end
+dith:dither40cols(function(w,h,getLinearPixel)
+	local pal={}
+	for i=0,15 do pal[i+1] = thomson.palette(i) end
+	return pal
+end)
\ No newline at end of file
diff --git a/share/grafx2/scripts/samples_2.4/picture/thomson/none_to8.lua b/share/grafx2/scripts/samples_2.4/picture/thomson/none_to8.lua
new file mode 100644
index 00000000..a02f8346
--- /dev/null
+++ b/share/grafx2/scripts/samples_2.4/picture/thomson/none_to8.lua
@@ -0,0 +1,82 @@
+-- ostro_to8.lua : convert a color image to a BM16
+-- (160x200x16) thomson image using the Ostromoukhov's
+-- error diffusion algorithm.
+--
+-- Version: 02-jan-2017
+--
+-- Copyright 2016-2017 by Samuel Devulder
+--
+-- This program is free software; you can redistribute
+-- it and/or modify it under the terms of the GNU
+-- General Public License as published by the Free
+-- Software Foundation; version 2 of the License.
+-- See 
+
+run('lib/thomson.lua')
+run('lib/ostromoukhov.lua')
+run('lib/color_reduction.lua')
+-- run('lib/zzz.lua')
+	
+-- get screen size
+local screen_w, screen_h = getpicturesize()
+
+-- Converts thomson coordinates (0-159,0-199) into screen coordinates
+local function thom2screen(x,y)
+	local i,j;
+	if screen_w/screen_h < 1.6 then
+		i = x*screen_h/200
+		j = y*screen_h/200
+	else
+		i = x*screen_w/320
+		j = y*screen_w/320
+	end
+	return math.floor(i*2), math.floor(j)
+end
+
+-- return the Color @(x,y) in normalized linear space (0-1)
+-- corresonding to the thomson screen (x in 0-319, y in 0-199)
+local function getLinearPixel(x,y)
+	local x1,y1 = thom2screen(x,y)
+	local x2,y2 = thom2screen(x+1,y+1)
+	if x2==x1 then x2=x1+1 end
+	if y2==y1 then y2=y1+1 end
+
+	local p = Color:new(0,0,0);
+	for j=y1,y2-1 do
+		for i=x1,x2-1 do
+			p:add(getLinearPictureColor(i,j))
+		end
+	end
+	p:div((y2-y1)*(x2-x1)) --:floor()
+	
+	return p
+end
+
+local red = ColorReducer:new():analyzeWithDither(160,200,
+	getLinearPixel,
+    function(y)
+		thomson.info("Collecting stats...",math.floor(y/2),"%")
+	end)
+	
+-- BM16 mode
+thomson.setBM16()
+
+-- define palette
+local palette = red:boostBorderColors():buildPalette(16)
+thomson.palette(0, palette)
+
+-- convert picture
+OstroDither:new(palette, 0)
+           :dither(thomson.h,thomson.w,
+             function(y,x) return getLinearPixel(x,y) end,
+			 function(y,x,c) thomson.pset(x,y,c) end, 
+			 true,
+			 function(x) thomson.info("Converting...",math.floor(x*100/160),"%") end)
+
+-- refresh screen
+setpicturesize(320,200)
+thomson.updatescreen()
+finalizepicture()
+
+-- save picture
+thomson.savep()
diff --git a/share/grafx2/scripts/samples_2.4/picture/thomson/none_to9.lua b/share/grafx2/scripts/samples_2.4/picture/thomson/none_to9.lua
new file mode 100644
index 00000000..9d7f40f6
--- /dev/null
+++ b/share/grafx2/scripts/samples_2.4/picture/thomson/none_to9.lua
@@ -0,0 +1,55 @@
+-- ostro_mo5.lua : converts a color image into a 
+-- TO9 image (320x200x16 with color clashes) 
+-- using Ostromoukhov's error diffusion algorithm.
+--
+-- Version: 02-jan-2017
+--
+-- Copyright 2016-2017 by Samuel Devulder
+--
+-- This program is free software; you can redistribute
+-- it and/or modify it under the terms of the GNU
+-- General Public License as published by the Free
+-- Software Foundation; version 2 of the License.
+-- See 
+
+run('lib/thomson.lua')
+thomson.optiMAP = false
+
+run('lib/ostromoukhov.lua')
+run('lib/color_reduction.lua')
+
+local dith=OstroDither:new()
+local tmp=dith.setLevelsFromPalette
+dith.setLevelsFromPalette = function(self)
+	tmp(self)
+	self.attenuation=0
+end
+dith:dither40cols(function(w,h,getLinearPixel)
+	local c16 = true
+	for y=0,h-1 do
+		for x=0,w-1 do
+			if getbackuppixel(x,y)>15 then c16 = false end
+		end
+	end
+	
+	local pal
+	if c16 then
+		pal = {}
+		for i=0,15 do
+			local r,g,b=getbackupcolor(i)
+			r = thomson.levels.pc2to[r]
+			g = thomson.levels.pc2to[g]
+			b = thomson.levels.pc2to[b]
+			pal[i+1] = r+g*16+b*256-273
+		end
+	else
+		pal=ColorReducer:new():analyzeWithDither(w,h,
+			getLinearPixel,
+			function(y)
+				thomson.info("Building palette...",math.floor(y*100/h),"%")
+			end):buildPalette(16)
+	end
+	thomson.palette(0, pal)
+
+	return pal
+end)
\ No newline at end of file
diff --git a/share/grafx2/scripts/samples_2.4/picture/thomson/ostro_mo5.lua b/share/grafx2/scripts/samples_2.4/picture/thomson/ostro_mo5.lua
new file mode 100644
index 00000000..e3160ba2
--- /dev/null
+++ b/share/grafx2/scripts/samples_2.4/picture/thomson/ostro_mo5.lua
@@ -0,0 +1,21 @@
+-- ostro_mo5.lua : converts a color image into a 
+-- MO5 image (16 fixed colors with color clash) 
+-- using Ostromoukhov's error diffusion algorithm.
+--
+-- Version: 02-jan-2017
+--
+-- Copyright 2016-2017 by Samuel Devulder
+--
+-- This program is free software; you can redistribute
+-- it and/or modify it under the terms of the GNU
+-- General Public License as published by the Free
+-- Software Foundation; version 2 of the License.
+-- See 
+
+run('lib/ostromoukhov.lua')
+
+OstroDither:new():dither40cols(function(w,h,getLinearPixel)
+	local pal={}
+	for i=0,15 do pal[i+1] = thomson.palette(i) end
+	return pal
+end)
\ No newline at end of file
diff --git a/share/grafx2/scripts/samples_2.4/picture/thomson/ostro_to7.lua b/share/grafx2/scripts/samples_2.4/picture/thomson/ostro_to7.lua
new file mode 100644
index 00000000..6e35fa51
--- /dev/null
+++ b/share/grafx2/scripts/samples_2.4/picture/thomson/ostro_to7.lua
@@ -0,0 +1,21 @@
+-- ostro_mo5.lua : converts a color image into a 
+-- TO7 image (8 fixed colors with color clash) 
+-- using Ostromoukhov's error diffusion algorithm.
+--
+-- Version: 02-jan-2017
+--
+-- Copyright 2016-2017 by Samuel Devulder
+--
+-- This program is free software; you can redistribute
+-- it and/or modify it under the terms of the GNU
+-- General Public License as published by the Free
+-- Software Foundation; version 2 of the License.
+-- See 
+
+run('lib/ostromoukhov.lua')
+
+OstroDither:new():dither40cols(function(w,h,getLinearPixel)
+	local pal={}
+	for i=0,7 do pal[i+1] = thomson.palette(i) end
+	return pal
+end)
\ No newline at end of file
diff --git a/share/grafx2/scripts/samples_2.4/picture/thomson/ostro_to8.lua b/share/grafx2/scripts/samples_2.4/picture/thomson/ostro_to8.lua
new file mode 100644
index 00000000..d2212961
--- /dev/null
+++ b/share/grafx2/scripts/samples_2.4/picture/thomson/ostro_to8.lua
@@ -0,0 +1,82 @@
+-- ostro_to8.lua : convert a color image to a BM16
+-- (160x200x16) thomson image using the Ostromoukhov's
+-- error diffusion algorithm.
+--
+-- Version: 02-jan-2017
+--
+-- Copyright 2016-2017 by Samuel Devulder
+--
+-- This program is free software; you can redistribute
+-- it and/or modify it under the terms of the GNU
+-- General Public License as published by the Free
+-- Software Foundation; version 2 of the License.
+-- See 
+
+run('lib/thomson.lua')
+run('lib/ostromoukhov.lua')
+run('lib/color_reduction.lua')
+-- run('lib/zzz.lua')
+	
+-- get screen size
+local screen_w, screen_h = getpicturesize()
+
+-- Converts thomson coordinates (0-159,0-199) into screen coordinates
+local function thom2screen(x,y)
+	local i,j;
+	if screen_w/screen_h < 1.6 then
+		i = x*screen_h/200
+		j = y*screen_h/200
+	else
+		i = x*screen_w/320
+		j = y*screen_w/320
+	end
+	return math.floor(i*2), math.floor(j)
+end
+
+-- return the Color @(x,y) in normalized linear space (0-1)
+-- corresonding to the thomson screen (x in 0-319, y in 0-199)
+local function getLinearPixel(x,y)
+	local x1,y1 = thom2screen(x,y)
+	local x2,y2 = thom2screen(x+1,y+1)
+	if x2==x1 then x2=x1+1 end
+	if y2==y1 then y2=y1+1 end
+
+	local p = Color:new(0,0,0);
+	for j=y1,y2-1 do
+		for i=x1,x2-1 do
+			p:add(getLinearPictureColor(i,j))
+		end
+	end
+	p:div((y2-y1)*(x2-x1)) --:floor()
+	
+	return p
+end
+
+local red = ColorReducer:new():analyzeWithDither(160,200,
+	getLinearPixel,
+    function(y)
+		thomson.info("Collecting stats...",math.floor(y/2),"%")
+	end)
+	
+-- BM16 mode
+thomson.setBM16()
+
+-- define palette
+local palette = red:boostBorderColors():buildPalette(16)
+thomson.palette(0, palette)
+
+-- convert picture
+OstroDither:new(palette)
+           :dither(thomson.h,thomson.w,
+             function(y,x) return getLinearPixel(x,y) end,
+			 function(y,x,c) thomson.pset(x,y,c) end, 
+			 true,
+			 function(x) thomson.info("Converting...",math.floor(x*100/160),"%") end)
+
+-- refresh screen
+setpicturesize(320,200)
+thomson.updatescreen()
+finalizepicture()
+
+-- save picture
+thomson.savep()
diff --git a/share/grafx2/scripts/samples_2.4/picture/thomson/ostro_to9.lua b/share/grafx2/scripts/samples_2.4/picture/thomson/ostro_to9.lua
new file mode 100644
index 00000000..01d1f0a8
--- /dev/null
+++ b/share/grafx2/scripts/samples_2.4/picture/thomson/ostro_to9.lua
@@ -0,0 +1,47 @@
+-- ostro_mo5.lua : converts a color image into a 
+-- TO9 image (320x200x16 with color clashes) 
+-- using Ostromoukhov's error diffusion algorithm.
+--
+-- Version: 02-jan-2017
+--
+-- Copyright 2016-2017 by Samuel Devulder
+--
+-- This program is free software; you can redistribute
+-- it and/or modify it under the terms of the GNU
+-- General Public License as published by the Free
+-- Software Foundation; version 2 of the License.
+-- See 
+
+run('lib/ostromoukhov.lua')
+run('lib/color_reduction.lua')
+
+OstroDither:new():dither40cols(function(w,h,getLinearPixel)
+	local c16 = h==200 and w==320
+	for y=0,h-1 do
+		for x=0,w-1 do
+			if getbackuppixel(x,y)>15 then c16 = false end
+		end
+	end
+	
+	local pal
+	if c16 then
+		pal = {}
+		for i=0,15 do
+			local r,g,b=getbackupcolor(i)
+			r = thomson.levels.pc2to[r]
+			g = thomson.levels.pc2to[g]
+			b = thomson.levels.pc2to[b]
+			pal[i+1] = r+g*16+b*256-273
+		end
+	else
+		pal=ColorReducer:new():analyzeWithDither(w,h,
+			getLinearPixel,
+			function(y)
+				thomson.info("Building palette...",math.floor(y*100/h),"%")
+			end):boostBorderColors():buildPalette(16)
+	end
+	
+	thomson.palette(0, pal)
+
+	return pal
+end)
\ No newline at end of file