mardi 6 avril 2010

Monde et image libre (didactiel 5/5)

En guise de conclusion, voici quelques remarques sur une des applications spectaculaires de l'effet Droste (qui peut aussi être appliqué à la création de films) : il s'agit, comme dans cet exemple de Julian Turner, du traitement d'une image faisant initialement figurer un porteur de cadre (délimitant une scène quelconque)



et grâce auquel l'image finale est celle du porteur portant un cadre délimitant une scène où le porteur porte un cadre délimitant une scène où ...etc ...
L'effet saisissant réside alors dans le fait que les différentes apparitions du cadre semblent être d'un seul tenant (pas d'inclusion stricte sinon seule une bonne photocopieuse réductrice suffirait).

Dans le principe, l'algorithme reste sensiblement le même : on balaye successivement tous les points de la nouvelle image, pour chacune de leur position on recherche la position correspondante dans l'image d'origine qui va fournir la couleur mais, ici, la recherche est menée dans une image intermédiaire virtuelle correspondant à la mise en abyme que l'on obtiendrait grâce à une photocopieuse.
Quand l'image et le cadre sont centrés sur le même point, ce point est également le point de fuite. Grâce au théorème de Thalès, on peut retrouver le pixel de référence de n'importe quel point de l'image par multiplication ou division du rayon considéré par le coefficient de réduction, autant de fois que nécéssaire, de manière à tomber dans la zone utile (la zone qui n'est pas dans la zone quelconque). En pratique, on sait que le rayon modifié répond au problème quand il est compris entre 2 rayons extrêmes ra et rb déterminés par de simples considérations trigonométriques à partir de l'angle :




Quand l'image et le cadre sont décentrés, il faut d'abord évaluer le point de fuite, qui est l'intersection des deux droites définies par deux paires de sommets correspondants, et la suite du calcul est la même (fourni ici dans sa version vectorisée) :
function [imdroste,mapdroste]= outframe(im,map,cx,cy,dx1,dx2,dy1,dy2)
[dimy, dimx] = size(im);
m=(dx1+dx2)/dimx;
tana=-log(m)/(2*pi);
cosa=sqrt(1/(1+tana^2));
sina=cosa*tana;
imdroste=reshape(1:dimx*dimy,dimx,dimy);

x=(1:dimx)-cx;
y=(1:dimy)-cy;
xx=ones(dimy,1)*x;
yy=y'*ones(1,dimx);
zr=distance2(0,0,xx,yy);
za=atan2(yy,xx);
z1x=log(zr);
z1y=za;

% ici on est passé dans le plan périodisé
z1y/=cosa;
% on y effectue une (anti)rotation pour annuler l'effet spirale
z2x=cosa*z1x+sina*z1y;
z2y=cosa*z1y-sina*z1x;
z3x=dy2*exp(z2x).*cos(z2y);
z3y=dy2*exp(z2x).*sin(z2y);

% ici on est de nouveau dans le plan visible centré sur (0,0) (abyme photocopieuse)
r=distance2(0,0,z3x,z3y);
a=atan2(z3y,z3x);

ra=zeros(dimy,dimx);
dom1=find( -pi<= a & a <= atan2(-dy2,-dx1));
dom2=find(atan2(-dy2,-dx1)<= a & a <= atan2(-dy2,dx2));
dom3=find(atan2(-dy2,dx2) <= a & a <= atan2(dy1,dx2));
dom4=find( atan2(dy1,dx2) <= a & a <= atan2(dy1,-dx1));
dom5=find( atan2(dy1,-dx1)<= a & a <= pi);

ra(dom1) =-dx1./cos(a(dom1));
ra(dom2) =-dy2./sin(a(dom2));
ra(dom3) = dx2./cos(a(dom3));
ra(dom4) = dy1./sin(a(dom4));
ra(dom5) =-dx1./cos(a(dom5));

rb=ra/m;

doma=find(r<ra);
domb=find(r>rb);
r(doma).*=(1/m).^ceil(log(ra(doma)./r(doma))/log(1/m));
r(domb).*=m.^ceil(log(rb(domb)./r(domb))/log(m));

% on a trouvé le rayon du pixel de référence, il ne manque plus que le point :
X=r.*cos(a)+cx;
Y=r.*sin(a)+cy;


% interpolation pour déterminer la couleur la plus vraiscemblable :
i=floor(X);
i(find(X<i))--;
i(find(i+1<X))++;
i(find(i<1))=1;
i(find(i>dimx-2))=dimx-2;
j=floor(Y);
j(find(Y<j))--;
j(find(j+1<Y))++;
j(find(j<1))=1;
j(find(j>dimy-2))=dimy-2;

i0=i;
j0=j;
i1=i+1;
j2=j+1;
idx = dimx*((0:dimy-1)'*ones(1,dimx))+ones(dimy,1)*(1:dimx);
trisup= find(distance2(X,Y,i0,j2) < distance2(X,Y,i1,j0));
j1=j0;
j1(trisup)++;
i2=i1;
i2(trisup)--;
v0=map(dimx*j0+i0+1,:);
v1=map(dimx*j1+i1+1,:);
v2=map(dimx*j2+i2+1,:);
s0=reshape(surface(X,Y,i1,j1,i2,j2),dimx*dimy,1)*ones(1,3);
s1=reshape(surface(i0,j0,X,Y,i2,j2),dimx*dimy,1)*ones(1,3);
s2=reshape(surface(i0,j0,i1,j1,X,Y),dimx*dimy,1)*ones(1,3);
mapdroste(idx,:) = (s0.*v0+s1.*v1+s2.*v2)./(s0+s1+s2);

et voici ce que produit cet appel :

coefreg=2.0;
[imreg,mapreg]=regularize2(im,map,coefreg,coefreg);
[imdroste,mapdroste]=outframe(imreg,mapreg,coefreg*315.0,coefreg*223.4,coefreg*143.0,coefreg*84.0,coefreg*65.6,coefreg*101.4);




où 315.0+223.4i est le point de fuite, 143 et 84 sont les distances de ce point à chacun des 2 segments verticaux du cadre dans l'image initiale, 65.6 et 101.4 celles aux 2 segments horizontaux, l'unité choisie étant le pixel.

Mise en abyme (didactiel 4/5)

La description mathématique de l'effet Droste est plus facile à faire quand l'image d'origine est constituée de zones couronnaires concentriques (d'où le choix de l'image de madnzany ). Dans ce cas, on peut simplement utiliser la classique fonction z ↦ ln(z) et sa réciproque z ↦ exp(z).
Dans le plan complexe en effet, l'image d'un cercle de rayon r centré sur l'origine par z ↦ ln(z/r) est un segment de droite vertical d'extrémités -iπ et +iπ et donc l'image d'une couronne délimitée par deux cercles de rayons r<R centrés sur l'origine est un rectangle de largeur ln(R/r). Cela suppose cependant qu'on a choisi la détermination principale des angles comprise entre -π et +π.
L'effet Droste n'est alors qu'une application possible de la propriété suivante : l'image par la fonction z ↦ exp(z) d'un segment de droite oblique est une portion de spirale logarithmique. L'image par cette fonction d'un rectangle incliné est donc une zone délimitée par 4 portions de spirales logarithmiques et un algorithme de mise en abyme sur une image n'a pour rôle que de gérer les continuités entre ces zones.
En résumé, l'image ln(z/r) d'un pixel z de l'image doit d'abord être recopiée dans tous les rectangles du plan complexe de manière à obtenir une série infinie de rectangles semblables c'est-à-dire de manière à conserver les positions relatives du pixel concerné vis à vis des 4 sommets du rectangle concerné (il suffit de recopier ln(z/r) tous les ln(R/r) sur la composante réelle en conservant sa composante imaginaire grâce à la périodicité des fonctions mises en jeu). Ensuite, le plan complexe est incliné d'un angle à déterminer. Enfin, on obtient l'image (au sens graphique) modifiée en prenant l'image (au sens analytique) par la fonction z ↦ exp(z) du plan complexe incliné (à une nuance près : un coefficient multiplicateur à déterminer doit être préalablement appliqué pour rester cohérent avec le choix de la détermination principale.)
La valeur de l'angle α à déterminer est conditionnée par un critère de continuité : la diagonale d'un rectangle incliné doit devenir verticale après inclinaison pour que les extrémités de cette diagonale ait la meme composante réelle (continuité des rayons). On montre que α vérifie tan α= ln(R/r)/(2π).
La valeur du coefficient k à déterminer est telle que la diagonale d'un rectangle soit égale à 2π (la diagonale devient une hauteur après inclinaison). On montre que k=cos(α).
L'algorithme procède à reculons en cherchant la couleur d'un pixel étant donnée sa position finale : étant donné un pixel z de l'image modifiée, on commence par se ramener au plan complexe incliné par la réciproque de z ↦ exp(z) c'est-à-dire par z ↦ ln(z). Puis on divise par k pour retrouver la taille du rectangle incliné originel. On ramène ce rectangle dans la position où ses cotés sont verticaux et horizontaux par rotation d'un angle -α. Puis on retrouve le rectangle qui lui a servi de modèle et on peut enfin retrouver la position du pixel de l'image d'origine (donc sa couleur) par application de la réciproque de z ↦ ln(z/r) à savoir z ↦ r exp(z) :

function [imdroste,mapdroste]=droste(im,map,r,R,cx,cy)
[dimy, dimx] = size(im);
width=log(R/r);
tana=width/(2*pi);
cosa=sqrt(1/(1+tana^2));
sina=cosa*tana;
imdroste=im;
mapdroste=zeros(dimx*dimy,3);

for x=1:dimx
for y=1:dimy
zr=distance2(0,0,x-cx,y-cy);
za=atan2(y-cy,x-cx);

z1x=log(zr);
z1y=za;

z1y/=cosa;

z2x=cosa*z1x+sina*z1y;
z2y=cosa*z1y-sina*z1x;

n=floor(abs(z2x)/width);
if z2x<0, z2x+=n*width; else z2x-=n*width; end

X=r*exp(z2x)*cos(z2y)+cx;
Y=r*exp(z2x)*sin(z2y)+cy;

if X>0 && X<=dimx && Y>0 && Y<=dimy
idx=dimy*(x-1)+y;
[i,j]=tweak(floor(X),floor(Y),X,Y,1,1);
if i>dimx-2, i=dimx-2;end
if j>dimy-2, j=dimy-2;end
if i<1, i=1;end
if j<1, j=1;end
if distance2(X,Y,i,(j+1)) < distance2(X,Y,(i+1),j)
[i0,i1,i2] = three([i,i+1,i]);
[j0,j1,j2] = three([j,j+1,j+1]);
else
[i0,i1,i2] = three([i,i+1,i+1]);
[j0,j1,j2] = three([j,j,j+1]);
end
v0=map(dimy*i0+j0+1,:);
v1=map(dimy*i1+j1+1,:);
v2=map(dimy*i2+j2+1,:);
s=[surface(X,Y,i1,j1,i2,j2),surface(i0,j0,X,Y,i2,j2),surface(i0,j0,i1,j1,X,Y)];
mapdroste(idx,:) = (s(1)*v0+s(2)*v1+s(3)*v2)/sum(s);
end
end
end

Une fois que la position X+iY du pixel de l'image d'origine à été trouvée, la détermination de la couleur se passe de la meme manière que pour la régularisation de l'image c'est-à-dire par interpolation.

Voici ce que produit cette évaluation :

cx=249;cy=183;
R=65;r=32;
[imcrop,mapcrop]=crop(im,map,cx-R,cx+R,cy-R,cy+R);# !! : le centre n'est plus cx+icy [imdroste,mapdroste]=droste(imcrop,mapcrop,r,R,R,R);



et voici ce qu'il faut calculer pour générer plus de régularité (ici 4x4 soit 16 fois plus long) :

cx=249;cy=183;
R=65;r=32;
[imcrop,mapcrop]=crop(im,map,cx-R,cx+R,cy-R,cy+R);
coefreg=4;
[imreg,mapreg]=regularize2(imcrop,mapcrop,coefreg,coefreg);
[imdroste,mapdroste]=droste(imreg,mapreg,coefreg*r,coefreg*R,coefreg*R,coefreg*R);

Outils de base (didactiel 3/5)

Afin d'illuster mon propos, je vais déjà me munir d'une image dont j'ai bien vérifié les droits (en utilisation et en modification) dans la mesure où je cite le nom de son auteur. Il s'agit donc d'une image de madnzany que l'on peut se procurer ici sous la licence "By" de Creative Commons. Voici de quoi elle a l'air :

Il s'agit de la photo au format 500x375 pixels d'une pièce de mobilier orné d'un disque partagé en 13 secteurs correspondant à 13 des (1600?) Chevaliers de la Table Ronde : le roi Arthur à la base et, dans l'ordre anti-horaire Lancelot,Tristan, Bedwere, Gaheris, Gawain, Kay, Lamorak, Geraint, Gareth, Bors, Galahad et enfin Perceval.
Inétressons-nous justement au secteur attribué à ce dernier et au petit coeur qui orne son blason. Etant données les dimensions de l'image d'origine, on peut faire tenir ce coeur dans une boîte de 11x11 pixels et l'extraction de cette boîte peut se faire gràce au code Octave suivant qui conserve le codage img (liste de couleurs et matrice de positions) en fabriquant ainsi une petite image à part entière :

function [imcrop,mapcrop]=crop (im,map,x1,xn,y1,yn)
dimx = size(im)(2);
dimcropy = yn-y1+1;
xcrop=x1:xn;
ycrop=y1:yn;
imcrop=[];
for x=1:xn-x1+1
row=[];
for y=1:dimcropy,row=[row; dimcropy*(x-1)+y];end
imcrop=[imcrop,row];
end
domcrop=[];
for x = xcrop
for y = ycrop, domcrop=[domcrop,(y-1)*dimx+x]; end
end
mapcrop=[map(:,1)(domcrop),map(:,2)(domcrop),map(:,3)(domcrop)];

On l'utilise par exemple avec
[imcrop,mapcrop]=crop(im,map,225,235,215,225);
ce qui dans notre cas produit une image de 121 pixels que Gnuplot interprète ainsi :




Il semble raisonnable, d'après cette représentation, de concevoir que Gnuplot considère chaque pixel comme le centre d'un carré et, afin d'afficher la petite image, il se charge de remplir chaque carré de la couleur porté par son pixel.
Avant de poursuivre notre chemin vers le 'droste effect', il faut se demander si cette façon de procéder est souhaitable. Comme le principe du 'droste effect' repose sur un repliement de l'image sur elle-même, certaines zones vont se recroqueviller sur elles-mêmes tandis que d'autres vont au contraire s'étaler. Comme la couleur entre deux pixels initialement voisins est inconnue, quand ces deux pixels s'éloignent l'un de l'autre, on augmente la taille d'une zone remplie de pixels de couleurs tout autant inconnues et le problème de coloriage se complique. Remarquons que ce problème n'intervient que quand les couleurs des deux pixels en question sont différentes car, dans le cas contraire, il est raisonnable de considérer que tous les pixels intermédiaires auront la même couleur.
Ainsi, dans le cas général, on ne peut pas se satisfaire d'une représentation ccomme celle de Gnuplot, et on va mettre en oeuvre par interpolation un procédé de régularisation permettant d'approcher une meilleure continuité dans l'image.

Tout d'abord quelques fonctions auxiliaires :

function [a,b] = two (v)
a = v(1);
b = v(2);
function [a,b,c]=three (v)
a=v(1);
b=v(2);
c=v(3);
function d=distance2(x0,y0,x1,y1)
d= sqrt((x1-x0).^2+(y1-y0).^2);
function s=surface(x0,y0,x1,y1,x2,y2)
s=(x0.*(y1-y2)+x1.*(y2-y0)+x2.*(y0-y1))/2;

Une fonction récursive cherchant à partir d'une position approximative le carré qui contient le point considéré :

function [i,j]=tweak(ii,jj,x,y,Dx,Dy)
if ii*Dx<=x
if x<=(ii+1)*Dx
if jj*Dy<=y
if y<=(jj+1)*Dy
i=ii;
j=jj;
else [i,j]=tweak(ii,jj+1,x,y,Dx,Dy);
end
else [i,j]=tweak(ii,jj-1,x,y,Dx,Dy);
end
else [i,j]=tweak(ii+1,jj,x,y,Dx,Dy);
end
else [i,j]=tweak(ii-1,jj,x,y,Dx,Dy);
end


puis une fonction principale pour régulariser une image aux taux donnés par mulx et muly :

function [imint,mapint]=regularize2(im, map,mulx,muly)
[dimy, dimx] = size(im);
[nx,ny]=two(round ([mulx*dimx, muly*dimy]));
[Dx,Dy]=two([1/(dimx-1),1/(dimy-1)]); % dom : [0,1]^2
[dx,dy]=two([1/(nx-1),1/(ny-1)]); % dom : [0,1]^2
mapint=zeros(nx*ny,3);
imint=[];
for xx=1:nx
row=[];
for yy=1:ny
idx = ny*(xx-1)+yy;
row=[row; idx];
x=(xx-1)*dx;
y=(yy-1)*dy;
ii=floor(x/Dx);
jj=floor(y/Dy);
[i,j]=tweak(ii,jj,x,y,Dx,Dy);
if i>dimx-2, i=dimx-2;end
if j>dimy-2, j=dimy-2;end
if distance2(x,y,i*Dx,(j+1)*Dy) < distance2(x,y,(i+1)*Dx,j*Dy)
[i0,i1,i2] = three([i,i+1,i]);
[j0,j1,j2] = three([j,j+1,j+1]);
else
[i0,i1,i2] = three([i,i+1,i+1]);
[j0,j1,j2] = three([j,j,j+1]);
end
v0=map(dimy*i0+j0+1,:);
v1=map(dimy*i1+j1+1,:);
v2=map(dimy*i2+j2+1,:);
x0=i0*Dx;x1=i1*Dx;x2=i2*Dx;
y0=j0*Dy;y1=j1*Dy;y2=j2*Dy;
s=[surface(x,y,x1,y1,x2,y2),surface(x0,y0,x,y,x2,y2),surface(x0,y0,x1,y1,x,y)];
mapint(idx,:) = (s(1)*v0+s(2)*v1+s(3)*v2)/(surface(x0,y0,x1,y1,x2,y2));
end
imint=[imint,row];
end

dont voici l'effet après avoir évalué
[imint,mapint]=regularize2(imcrop,mapcrop,8.0,8.0);

GIMP Immortalizes My Pictures (didactiel 2/5)

Comment parler de libre et d'image sans évoquer GIMP. J'en parlerai d'autant plus dans cet entrée que Scheme, le langage de programmation utilisable dans GIMP, est, comme Haskell, un langage fonctionnel (et pour l'anecdote, c'est en apprenant Scheme que j'ai re-decouvert Python ce qui m'a amené, le moins naturellement du monde, à Haskell!).
Quand on veut manipuler (recadrer, recoloriser, recomposer...) avec GIMP une image comme notre pauvre image créee précédemment avec Octave, le plus simple est de la convertir au format ppm :
saveimage("pourrie.ppm",x,"ppm")
qui a la particularité de générer un fichier binaire (on peut quand même lire le code "P6" en début de fichier confirmant qu'il s'agit bien du format ppm toutefois rien ne garantit de la qualité du fichier produit que, pour ma part, je n'arrive à lire nulle part) .
L'inverse est malheureusement moins direct (ce qui peut se comprendre dans la mesure où Octave n'est pas un outil de traitement d'image a priori) : GIMP n'a pas les moyens de fournir des informations picturales à Octave.
Afin de pouvoir triturer toute sorte d'images (dans la limite de la légalité, bien sûr), le moyen le plus simple auquel j'ai pensé (outre celui de programmer en Scheme mais c'est quand même plutôt pénible) consiste à convertir l'image que l'on a sous la main au format ppm (ppm "P3" cette fois : la version texte de ce format -GIMP nous le propose toujours au moment de l'exportation-) puis, grâce à un petit script, de convertir ce fichier ppm au format img.
Le but de cette entrée est particulièrement de présenter ce script écrit en Haskell et plus généralement de montrer que la syntaxe d'Haskell possède la souplesse et la robustesse permettant d'écrire des codes concis et fiables.
module Main where

import System.IO
import System.Environment (getArgs)

ppm2img inh outh = do
(height, width) <- readHeader inh
hPutStrLn outh "# Created by gimp2octave"
hPutStrLn outh "# name: map"
hPutStrLn outh "# type: matrix"
hPutStrLn outh $ "# rows: " ++ show (height*width)
hPutStrLn outh "# columns: 3"
readWriteLoop inh outh
hPutStrLn outh "# name: img"
hPutStrLn outh "# type: matrix"
hPutStrLn outh $ "# rows: " ++ show height
hPutStrLn outh $ "# columns: " ++ show width
writeLoop outh height width 1 where
readHeader inh = do
inStr <- hGetLine inh -- P3
inStr <- hGetLine inh -- Created by Gimp
inStr <- hGetLine inh -- Height Width
let [hStr,wStr] = words inStr
inStr <- hGetLine inh -- Base
return (read hStr::Int, read wStr::Int)
readWriteLoop inh outh = do
ineof <- hIsEOF inh
if ineof
then return ()
else do
rStr <- hGetLine inh
gStr <- hGetLine inh
bStr <- hGetLine inh
hPutStrLn outh $ " " ++ rStr ++ " " ++ gStr ++ " " ++ bStr
readWriteLoop inh outh
writeLoop outh h w n = do
if n > h
then return ()
else do
let str = foldr (++) "" $ map ((" " ++) . show) [n,n+h..n+(w-1)*h]
hPutStrLn outh str
writeLoop outh h w (n+1)

main = do
args <- getArgs
inh <- openFile (head args) ReadMode
outh <- openFile (head (tail args)) WriteMode
ppm2img inh outh
hClose inh
hClose outh

Image et monde libre (didactiel 1/5)

Après près d'un an de deuil (voir dernières entrées) et à l'approche des vacances, il est temps d'aborder une nouvelle série d'articles.
Cette série portera sur le traitement d'images tout particulièrement du point de vue du monde libre (ou open source comme le disent les moins sérieux.)
La construction ex-nihilo d'un outil générant un 'droste effect' sur une image en sera le fil conducteur.
Les grandes lignes :
* outils de base
* fondements mathématiques et outils sophistiqués
* mise en oeuvre
* exemples
Tout d'abord, je précise que je travaille sur Kubuntu 9.04 et que, sans précisions supplémentaires, les logiciels présentés sont dans les versions standard que me propose aptitude (NB pour les versions antérieures,il peut être également utile d'installer le package xloadimage).
Le langage de choix sera celui d'Octave (et son format d'images RGB pour les fichiers .img) mais je vais, dans la prochaine entrée, surtout parler d'Haskell qui est un langage idéal pour développer de petits utilitaires.
Les images considérées seront des images couleurs et pour créer une image rouge de hauteur 100 et de largeur 50 de type RGB avec Octave, il faut créer un tableau x à 3 dimensions :
base=255;
x(:,:,1)=base*ones(100,50);
x(:,:,2)=zeros(100,50);
x(:,:,3)=zeros(100,50);
et l'affichage avec Gnuplot s'effectue en faisant :
imshow(x)
Quand on demande l'enregistrement de notre image (pourrie) :
saveimage("pourrie.img",x,"img")
Octave crée le fichier au format img (NB les versions antérieures ne semblent pas se tenir à cette présentation et on devra y parenir en externe après avoir sauvegardé en brut avec save "raw" x puis un réencodage maison) qui a l'allure suivante :
# Created by Octave ...
# name: map
# type: matrix
# rows: 5000
# columns : 3
255 0 0
255 0 0
... (5000 lignes identiques)
255 0 0
# name: img
# type: matrix
# rows: 100
# columns: 50
1 101 201 ... 4901
2 102 202 ... 4902
...
100 200 ... 5000
c'est-à-dire une liste de couleurs suivie d'une matrice de positions de pixels.
Pour retrouver notre image le lendemain, cet arrangement nous oblige à effectuer la lecture de la façon suivante :
[im,map]=loadimage("pourrie.img");
et l'affichage :
imshow(ind2rgb(im,map))