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.