%B. Vandeportaele October 2017 %Cubic spline interpolation in Homography space for deflectometry close all clear all %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %set the size of the spline control point grid sx=15; sy=10; %how many sampled pixels between 2 control points in each direction for %the generation of the image samplingfactor=100; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %choose the type of interpolation for homography coefficients through %function pointer InterpolationFunction=@bicubicInterpolationMulti; %InterpolationFunction=@bilinearInterpolationMulti; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% HomographyGrid=zeros(sy,sx,9); %identity matrix by default HomographyDefault=reshape(eye(3),9,1); %if the default value has to be set differently %HomographyDefault=reshape([1,0,65; 0,1,0; 0,0,1],9,1); %fill the grid with the default value for y=1:sy for x=1:sx HomographyGrid(y,x,:)=HomographyDefault; end end %perturbate some control points of the spline HomographyGrid(5,4,:)=reshape([1,0,0; 0,1,0; 0,0,0.9],9,1); HomographyGrid(4,4,:)=reshape([1.3,0,0; 0,1.3,0; 0,0,1],9,1); % HomographyGrid(5,4,:)=reshape([1,0,0; 0,1,0; 0,0,0.22],9,1); % HomographyGrid(5,5,:)=reshape([1,0,0; 0,1,0; 0,0,1.3],9,1); % HomographyGrid(6,5,:)=reshape([1,0,0; 0,1,0; 0,0,1.5],9,1); % HomographyGrid(6,4,:)=reshape([1,0,0; 0,1,0; 0,0,0.3],9,1); % HomographyGrid(6,4,:)=reshape([1.2,0,0; 0,0.9,0; 0,0,1],9,1); HomographyGrid(6,5,:)=reshape([1,0,0; 0,1,0; 0,0,1],9,1); HomographyGrid(6,8,:)=reshape([1,0,60; 0,1,-30; 0,0,1],9,1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %interpolate the homography in a small region at low resolution for %visualization purpose %very unefficient code because the size of the destination matrix are not %known at the begining=> reallocation indy=1; for yy=3:0.1:sy-2 indx=1; for xx=3:0.1:sx-2 %val = InterpolationFunction( xx,yy,HomographyGrid); val = feval(InterpolationFunction, xx,yy,HomographyGrid); indx=indx+1; res(indy,indx,1)=xx; res(indy,indx,2)=yy; res(indy,indx,3)=val(9); InterpolatedFunction(indy,indx,:)=val; end indy=indy+1; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure axis equal plot3(res(:,:,1),res(:,:,2),res(:,:,3),'r+') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure for channel=1:9 subplot(3,3,channel) imagesc(InterpolatedFunction(:,:,channel)) end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %compute correspondance map correspondancemapu=zeros(sy*samplingfactor,sx*samplingfactor); correspondancemapv=zeros(sy*samplingfactor,sx*samplingfactor); for yy=2*samplingfactor:sy*samplingfactor-3*samplingfactor for xx=2*samplingfactor:sx*samplingfactor-3*samplingfactor %if xx==1099 % display 'ici' %end H =reshape( InterpolationFunction( xx/samplingfactor,yy/samplingfactor,HomographyGrid),3,3); p2=H*[xx,yy,1]'; %transformation homographique inverse %ud=round(p2(1)/p2(3)); %vd=round(p2(2)/p2(3)); ud=p2(1)/p2(3); vd=p2(2)/p2(3); correspondancemapu(yy,xx)=ud; correspondancemapv(yy,xx)=vd; end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %save to a correspondance map using an openCV compatible format %file saved with .map extension map=[correspondancemapv;correspondancemapu]; imagedirectory='./'; mapname=sprintf('%s/correspondance.map',imagedirectory); saveRawFloat(map, mapname); display('correspondance.map file saved in imagedirectory'); %future use: process the images using python bindings of opencv remap %function: see /media/HD500GO/saveHDDgarossos/menage_ocamcalib/remap %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Apply the interpolated homography to a synthetic image. %%% it visualizes the equivalent of the control grid but in the source %%% image plane imoutsinu=zeros(sy*samplingfactor,sx*samplingfactor); imoutsinv=zeros(sy*samplingfactor,sx*samplingfactor); imoutsinuv=zeros(sy*samplingfactor,sx*samplingfactor); imoutsinvu=zeros(sy*samplingfactor,sx*samplingfactor); imoutcircle=zeros(sy*samplingfactor,sx*samplingfactor); imoutradial=zeros(sy*samplingfactor,sx*samplingfactor); imoutgrid=zeros(sy*samplingfactor,sx*samplingfactor); for yy=2*samplingfactor:size(imoutsinu,1)-3*samplingfactor for xx=2*samplingfactor:size(imoutsinu,2)-3*samplingfactor ud=correspondancemapu(yy,xx); vd=correspondancemapv(yy,xx); %for a square grid valpix=127*mod(floor(ud/samplingfactor),2) + 127*mod(floor(vd/samplingfactor),2) ; imoutgrid(yy,xx)=valpix; %on attribue la couleur du pixel concerné %for a sinus wave in u direction valpix=uint8(127+127*sin(ud/2.5)); imoutsinu(yy,xx)=valpix; %on attribue la couleur du pixel concerné %for a sinus wave in v direction valpix=uint8(127+127*sin(vd/2.5)); imoutsinv(yy,xx)=valpix; %on attribue la couleur du pixel concerné %for a sinus wave in uv diagonal direction valpix=uint8(127+127*sin((vd+ud)/(2*2.5))); imoutsinuv(yy,xx)=valpix; %on attribue la couleur du pixel concerné %for a sinus wave in vu diagonal direction valpix=uint8(127+127*sin((vd-ud)/(2*2.5))); imoutsinvu(yy,xx)=valpix; %on attribue la couleur du pixel concerné %for a sinus wave in circle direction valpix=uint8(127+127*sin((norm([(ud-750);(vd-500)],2) /2.5))); imoutcircle(yy,xx)=valpix; %on attribue la couleur du pixel concerné %for a sinus wave in radial direction angle=atan2((vd-500),(ud-750))*20; valpix=uint8(127+127*sin(angle)); imoutradial(yy,xx)=valpix; %on attribue la couleur du pixel concerné end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %display images figure subplot(3,1,1); image(imoutgrid); axis equal;colormap(gray(256)); subplot(3,1,2); image(imoutsinu); axis equal;colormap(gray(256)); subplot(3,1,3); image(imoutsinv); axis equal;colormap(gray(256)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %save images imwrite(uint8(imoutgrid),'imoutgrid.png'); imwrite(uint8(imoutsinu),'imoutu.png'); imwrite(uint8(imoutsinv),'imoutv.png'); imwrite(uint8(imoutsinuv),'imoutuv.png'); imwrite(uint8(imoutsinvu),'imoutvu.png'); imwrite(uint8(imoutcircle),'imoutcircle.png'); imwrite(uint8(imoutradial),'imoutradial.png'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %generate sequences of images nn=0; for angle=0:2*pi/60:2*pi if angle~=0 nn=nn+1 end for yy=2*samplingfactor:size(imoutsinu,1)-3*samplingfactor for xx=2*samplingfactor:size(imoutsinu,2)-3*samplingfactor ud=correspondancemapu(yy,xx); vd=correspondancemapv(yy,xx); %for a sinus wave in u direction valpix=uint8(127+127*sin(angle+(ud/2.5))); imoutsinu(yy,xx)=valpix; %on attribue la couleur du pixel concerné %for a sinus wave in v direction valpix=uint8(127+127*sin(angle+(vd/2.5))); imoutsinv(yy,xx)=valpix; %on attribue la couleur du pixel concerné end end imgname=sprintf('%s/sinu%02i.png',imagedirectory,nn); imwrite(uint8(imoutsinu),imgname); imgname=sprintf('%s/sinv%02i.png',imagedirectory,nn); imwrite(uint8(imoutsinv),imgname); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %compute the mpg videos... does not work from within matlab, invoke it %through a linux shell %system('avconv -f image2 -i sinu%02d.png -f avi -vcodec mpeg4 -b 4000k -g 300 -bf 2 -y ./videosinu.mpg') %system('avconv -f image2 -i sinv%02d.png -f avi -vcodec mpeg4 -b 4000k -g 300 -bf 2 -y ./videosinv.mpg');