Outils pour utilisateurs

Outils du site


deflectometry

Différences

Ci-dessous, les différences entre deux révisions de la page.

Lien vers cette vue comparative

Les deux révisions précédentes Révision précédente
Prochaine révision
Révision précédente
deflectometry [2017/10/31 15:56]
bvandepo
deflectometry [2018/01/10 14:46] (Version actuelle)
bvandepo
Ligne 245: Ligne 245:
 {{https://​homepages.laas.fr/​bvandepo/​files/​deflectometry/​simu2.png}} {{https://​homepages.laas.fr/​bvandepo/​files/​deflectometry/​simu2.png}}
  
 +
 +
 +=====2D analysis for uncertainty propagation=====
 +<file matlab simu2D.m>​
 +%B. Vandeportaële & D. Maestro 2017
 +close all
 +clear all
 +%example for the inversion of P1 to P1 scale and offset
 +pixperm=sym('​pixperm','​real'​);​
 +pud=sym('​pud','​real'​);​
 +Kd=[pixperm ​ pud; 0 1]
 +a=sym('​a','​real'​);​
 +b=sym('​b','​real'​);​
 +c=sym('​c','​real'​);​
 +d=sym('​d','​real'​);​
 +Kdinv=[a b;c d ]
 +eqn= reshape(Kdinv*Kd-eye(2),​4,​1)
 +S=solve(eqn,​a,​b,​c,​d)
 +Kdinv=[S.a S.b;S.c S.d]
 +jacobian(reshape(Kdinv,​4,​1),​[pixperm,​pud])
 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +tcx=sym('​tcx','​real'​);​
 +tcy=sym('​tcy','​real'​);​
 +r11=sym('​r11','​real'​);​
 +r12=sym('​r12','​real'​);​
 +r21=sym('​r21','​real'​);​
 +r22=sym('​r22','​real'​);​
 +M=[r11,​r12,​tcx;​r21,​r22,​tcy ;0 0 1]
 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +fuc=sym('​fuc','​real'​);​
 +puc=sym('​puc','​real'​);​
 +K=[fuc ​ puc 0;
 +    0  1 0];
 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +E=[1 0;0 0 ;0 1]
 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +H=K*M*E*Kdinv
 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +tdx=sym('​tdx','​real'​);​
 +tdy=sym('​tdy','​real'​);​
 +thetad=sym('​thetad','​real'​);​
 +Md=[cos(thetad),​-sin(thetad),​ tdx;​sin(thetad),​cos(thetad),​tdy;​ 0, 0, 1]
 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 + 
 +
 +</​file>​
 +
 +
 +
 +=====Fitting of sinus function=====
 +<file matlab phase.m>
 +%@Bertrand Vandeportaele LAAS 2017
 +function phasefunc()
 +
 +close all
 +
 +symbolicComputation=false;​
 +
 +ud=218;
 +vd=153;
 +
 +
 +ud=289;
 +vd=93;
 +
 +of=100;
 +af=73;
 +freq=6/​(1920);​
 +wt=2*pi*freq;​
 +lum=[];
 +lumn=[];
 +n=[];
 +
 +nbPointsPerPeriod=400
 +global gnoiseLevel
 +gnoiseLevel=3
 +
 +%useSyntheticData=true
 +useSyntheticData=false
 +if useSyntheticData
 +    %polynomial for some non linear disturbance
 +    pol=[ 0.003; 0.35; 10]
 +    %polynomial for some no disturbance
 +    %pol=[ 0; 1; 0]
 +    for i=1:​nbPoints*PerPeriod
 +        p=i*2*pi/​nbPointsPerPeriod;​
 +        deltai=ud*wt;​
 +        l=of+af*sin(deltai+p);​
 +        lum=[lum,​l];​
 +        %non linear disturbance
 +        l2=pol(3)+pol(2)*l+pol(1)*l^2;​
 +        %add noise
 +        ln=l2+gnoiseLevel*randn(1);​
 +        lumn=[lumn,​ln];​
 +        n=[n,i-1];
 +    end
 +else
 +    folder_path = '​./​pos_02/​0016/​H/';​
 +    %folder_path = '​./​pos_02/​0032/​H/';​
 +    %folder_path = '​./​pos_02/​0064/​H/';​
 +    %file_extension='​bmp';​
 +    %num_imgs= 11
 +    ​
 +    folder_path = '/​media/​HD500GO/​daniel/​phases/​gray/​m_00/​t4096/​pos_02/​0016/​H/'​
 +    file_extension='​tif';​
 +    num_imgs= 64
 +    ud=128;
 +    vd=1024;
 +    % vd=1
 +    ​
 +    folder_path = '/​media/​HD500GO/​daniel/​phases/​gray/​m_00/​t40/​m_00/​pos_02/​0040/​H/'​
 +    file_extension='​tif';​
 +    num_imgs= 128
 +    ud=1280;
 +    vd=1024;
 +    vd=1000
 +    ​
 +    vd=400
 +
 +    %To avoid reloading the images each time.
 +    % clear I to force it
 +    global I
 +    ​
 +    ​
 +    %if true
 +    if (size(I,​1)==0) || (size(I,​2)==0)
 +        %% Load Images
 +        fprintf('​Start loading data\n'​)
 +        tic;
 +        ​
 +        for indx=0:​num_imgs-1
 +            ima_name ​ = sprintf('​%s%s%03d.%s',​ folder_path,'​img_',​indx,​file_extension);​
 +            if indx==0
 +                I0 = imread(ima_name);​
 +                [ rs, cs] = size(I0);
 +                % Allocate matrix
 +                I = zeros(rs, cs, num_imgs);
 +                I(:,:,1) = I0;
 +            else
 +                %ima_name = fullfile(folder_path,​ image_files(indx).name);​
 +                I(:,:,​indx+1) = imread(ima_name);​
 +            end
 +        end
 +        fprintf([num2str(num_imgs),​ ' Images Loaded. '])
 +        toc
 +    end
 +end
 +n=(0:​num_imgs-1)';​
 +lumn=reshape(I(vd,​ud,:​),​num_imgs,​1)';​
 +lum=lumn;
 +nbPointsPerPeriod=num_imgs
 +%polynomial for no disturbance
 +pol=[ 0; 1; 0]
 +
 +
 +%numberOfSamplesForMeaning=10
 +%mask=(1/​numberOfSamplesForMeaning)*ones(1,​numberOfSamplesForMeaning)
 +%lumn=conv(lumn,​mask,'​same'​)
 +
 +global gpuls
 +gpuls=2*pi/​nbPointsPerPeriod;​
 +
 +
 +%%%%%%%%%%%%%%%%TEST to reduce the number of input data%%%%%%%%%%%%%%%%%%%%
 +% nbPointsPerPeriod=6
 +% n=n(1:​nbPointsPerPeriod)
 +% lum=lum(1:​nbPointsPerPeriod)
 +% lumn=lumn(1:​nbPointsPerPeriod)
 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +%%%%%%%%%%%%%%%TEST pertubate some points to check the robustness
 +%lumn(40)=6000
 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +tic
 +[xLinear covarianceOutLinear]=fitSinusLinear(n,​lumn,​gpuls)
 +toc
 +
 +tic
 +[xLinearOpt covarianceOutLinearOpt,​err]=fitSinusLinearOptimized(n,​lumn,​gpuls)
 +toc
 +
 +[xPerfect covarianceOutPerfect]=fitSinusLinear(n,​lum,​gpuls)
 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +%non linear method
 +if symbolicComputation
 +    ofs = sym('​ofs','​real'​);​
 +    amp = sym('​amp','​real'​);​
 +    gpuls = sym('​puls','​real'​);​
 +    phas = sym('​phas','​real'​);​
 +    t = sym('​t','​real'​);​
 +    y=ofs+amp*sin(gpuls*t+phas)
 +    ymes = sym('​ymes','​real'​);​
 +    res=expand((ymes-y))
 +    %jacobian of res used during optimization
 +    J=jacobian(res,​[ofs,​amp,​phas])'​
 +    %jacobian of res^2, should be 0 at the end of the optimization only
 +    %when there is no noise on the data
 +    %J2=jacobian(res^2,​[ofs,​amp,​phas])'​
 +end
 +%ofs=of+10;
 +%amp=af*0.9;​
 +%phas=0;
 +%x0=[ofs,​amp,​phas];​
 +%x0=x0+[10,​40,​0.3]%add some error
 +[xNonLinear,​resxNonLinear,​JxNonLinear]= sinusNonLinearFit(n,​lumn,​xLinear);​
 +
 +%NON!!!!!!!!!!!!!!!!!
 +%covarianceInNonLinear=diag(resxNonLinear.*resxNonLinear);​
 +%W=inv(covarianceInNonLinear);​
 +
 +%mean of absolute value of errr
 +mean (abs(resxNonLinear),​1)
 +variance=var(resxNonLinear)
 +standardDeviation=variance^0.5
 +W=eye(size(resxNonLinear,​1))*variance^-1;​
 +
 +JxNonLinear
 +
 +InformationOutNonLinear= (JxNonLinear'​*W*JxNonLinear)
 +covarianceOutNonLinear= inv(InformationOutNonLinear);​
 +
 +
 +displayResult(xLinear, ​ covarianceOutLinear,'​xLinear ​     ' )
 +%theese 2 values have comparable variances
 +displayResult(xNonLinear, ​ covarianceOutNonLinear,'​xNonLinear ​  '​ )
 +%return
 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +%non linear method ​ taking into account the non linearity of the
 +%camera/​display
 +if symbolicComputation
 +    ofs = sym('​ofs','​real'​);​
 +    amp = sym('​amp','​real'​);​
 +    gpuls = sym('​puls','​real'​);​
 +    phas = sym('​phas','​real'​);​
 +    t = sym('​t','​real'​);​
 +    yperfect=ofs+amp*sin(gpuls*t+phas);​
 +    ​
 +    nla = sym('​nla','​real'​);​
 +    nlb = sym('​nlb','​real'​);​
 +    nlc = sym('​nlc','​real'​);​
 +    y=nla * yperfect^2 + nlb * yperfect + nlc;
 +    ​
 +    ymes = sym('​ymes','​real'​);​
 +    res=expand((ymes-y))
 +    %jacobian of res used during optimization
 +    J=jacobian(res,​[ofs,​amp,​phas,​nla,​nlb,​nlc])'​
 +    ​
 +end
 +%ofs=of+10;
 +%amp=af*0.9;​
 +%phas=0;
 +%x0=[ofs,​amp,​phas];​
 +%x0=x0+[10,​40,​0.3]%add some error
 +x0_2=[xLinear,​ 0, 1, 0]
 +%x0_2=[1500 1.3 0.5 , 0, 1, 0]
 +[xNonLinear_2,​resxNonLinear_2,​JxNonLinear_2]= sinusNonLinearFit2(n,​lumn,​x0_2)
 +
 +%covarianceOutNonLinear_2=zeros(6,​6);​
 +%covarianceOutNonLinear_2=JxNonLinear_2'​*diag(resxNonLinear_2.*resxNonLinear_2)*JxNonLinear_2
 +%TODO: plot that non linear function and generate the corrected plot
 +%for sinus!
 +
 +%JxNonLinear_2(:,​4)=0
 +
 +pol=xNonLinear_2(4:​6)
 +%!!!!! NON j'ai estimé le polynome inverse!!!!
 +
 +%xNonLinear_2=xNonLinear_2(1:​3);​
 +%covarianceOutNonLinear_2=JxNonLinear_2(:,​1:​3)'​*diag(resxNonLinear_2.*resxNonLinear_2)*JxNonLinear_2(:,​1:​3)
 +%xNonLinear_2(4:​6) =pol
 +%covarianceOutNonLinear_2(6,​6)=0
 +%W=diag(resxNonLinear_2.*resxNonLinear_2);​
 +%covarianceInNonLinear=inv(W)
 +%covarianceInNonLinear=eye(128);​
 +%InformationOutNonLinear_2= (JxNonLinear_2'​*covarianceInNonLinear*JxNonLinear_2)
 +%covarianceOutNonLinear_2= inv(InformationOutNonLinear_2)
 +
 +resxNonLinear_2
 +%mean of absolute value of errr
 +mean (abs(resxNonLinear_2),​1)
 +variance=var(resxNonLinear_2)
 +variancebis=mean(resxNonLinear_2.*resxNonLinear_2) %>give approximately the variance if the mean is close to 0  (no bias)
 +varianceter=(resxNonLinear_2'​*resxNonLinear_2)/​nbPointsPerPeriod
 +standardDeviation=variance^0.5
 + 
 +%
 +
 +
 +W1=eye(size(resxNonLinear_2,​1))*variance^-1;​
 +%W2=diag( ((resxNonLinear_2.*resxNonLinear_2)*nbPointsPerPeriod).^-1 )
 +%W2=diag( ((resxNonLinear_2.^2)*nbPointsPerPeriod).^-1 ) 
 +
 +%the squared error is computed for each data point, with a sample size of
 +%one, so the variance of this 1 set is equal o val^2/1
 +%W2=diag( ((resxNonLinear_2.^2)).^-1 );   ​%<​==>​ W2=diag( (resxNonLinear_2.^-2)) ​  
 +%traceW1=trace(W1)
 +%MtraceW2=trace(W2)
 +%return
 +%[diag(W1),​diag(W2)]'​*ones(128,​2)
 +
 +%return
 +W=W1;
 +
 +
 +    %processing of uncertainty using the whole set of parameters lead to
 +    %non invertible information matrix
 +    InformationOutNonLinear_2= (JxNonLinear_2'​*W*JxNonLinear_2)
 +    covarianceOutNonLinear_2= inv(InformationOutNonLinear_2) %= InformationOutNonLinear_2^-1
 +
 +    ​
 +     ​[u,​s,​v]=svd(covarianceOutNonLinear_2)
 +     
 +covarianceOutNonLinear_2-( ​    ​u*s*v'​)
 +sinv=zeros(6,​6);​
 +for i=1:6
 +    s(i,i) %The sixth term is e-13 and the first is e+16, poor conditionning
 +    if (i~=6)
 +        sinv(i,​i)=1/​s(i,​i);​
 +    end
 +end
 +%sinv=inv(s)
 +sinv
 +inverse=v*sinv*u'​
 +
 +for i=1:6
 +    inverse(i,​i)
 +end
 +    ​
 +rank(covarianceOutNonLinear_2)
 +
 +inverse*covarianceOutNonLinear_2
 +if false
 +end
 +
 +%work on a reduced version focucing on the sinusoidal function parameters,
 +%because the information matrix is of rank 3 only!
 +JxNonLinear_2red=JxNonLinear_2(:,​1:​3);​
 +InformationOutNonLinear_2= (JxNonLinear_2red'​*W*JxNonLinear_2red)
 +
 +rank(InformationOutNonLinear_2)
 +
 +covarianceOutNonLinear_2= inv(InformationOutNonLinear_2)
 +covarianceOutNonLinear_2(1,​1)
 +covarianceOutNonLinear_2(2,​2)
 +covarianceOutNonLinear_2(3,​3)
 +
 +%add 0s to the covariance matrix so it can be displayed later
 +covarianceOutNonLinear_2(6,​6)=0;​
 +
 +
 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +
 +
 +
 +figure
 +subplot(1,​3,​1);​
 +%xx=1:255;
 +xx=1:​(2^12)-1;​
 +yy=polyval(pol,​xx);​
 +plot(xx, yy,'​r-'​)
 +subplot(1,​3,​2);​
 +hold on
 +plot(n, lum,'​r+'​)
 +displaySin(min(n)-1,​max(n)+1,​xPerfect,'​r-'​);​
 +plot(n, lumn,'​b+'​)
 +displaySin(min(n)-1,​max(n)+1,​xLinear,'​k:'​);​
 +displaySin(min(n)-1,​max(n)+1,​xNonLinear,'​g--'​);​
 +
 +displaySinNonLinear(min(n)-1,​max(n)+1,​xNonLinear_2,'​b--'​);​
 +
 +
 +legend('​data without noise','​function without noise','​data with noise','​fitted function using linear method','​fitted function using non linear method','​fitted function using non linear method and non linear correction'​)
 +subplot(1,​3,​3);​
 +
 +%resxNonLinear=randn(size(resxNonLinear,​1),​1);​
 +resToPlot=resxNonLinear_2;​
 +
 +numberOfBinsToPlot=12;​
 +maxErr=round(max(abs(resToPlot)))+1;​
 +binranges = [-maxErr:​maxErr/​(numberOfBinsToPlot-2):​maxErr];​
 +%binranges = [-maxErr:​maxErr/​10:​maxErr];​
 +nbBins=size(binranges,​1);​
 +[nelements,​ind] =histc(resToPlot,​binranges);​
 +plot(binranges,​nelements*nbBins/​nbPointsPerPeriod,'​r+'​);​
 +hold on
 +plot(resToPlot,​zeros(size(resToPlot,​1)),'​g+'​);​
 +meanOfError=mean(resToPlot)
 +varianceOfError=var(resToPlot)
 +standardDeviationOfError=sqrt(varianceOfError)
 +gaussianSamplesx=[-maxErr:​0.1:​maxErr];​
 +gaussianSamplesy=(1/​sqrt(2*pi*varianceOfError))*exp(-( ((gaussianSamplesx-meanOfError).^2)/​(2*varianceOfError)));​
 +plot(gaussianSamplesx,​ gaussianSamplesy ,'​b-'​);​
 +
 +
 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +displayResult(xPerfect,​ covarianceOutPerfect,'​xPerfect ​    '​ )
 +displayResult(xLinear, ​ covarianceOutLinear,'​xLinear ​     ' )
 +displayResult(xNonLinear_2, ​ covarianceOutNonLinear_2,'​xNonLinear_2 ' )
 +
 +%convert the 3sigma uncertainty in phase to pixel position
 +pixelPeriod=16
 +phase3sigma=3*sqrt(covarianceOutLinear(3,​3))
 +pixelLocation3Sigma=phase3sigma*pixelPeriod/​(2*pi)
 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +%compare the distribution of the error with the theroretic one using the
 +%Cumulative distribution function
 +figure
 +resxNonLinear_2sorted=sort(resxNonLinear_2);​
 +plot(resxNonLinear_2sorted,​n/​nbPointsPerPeriod,'​r-'​);​
 +%fonction de répartition:​ https://​en.wikipedia.org/​wiki/​Normal_distribution
 +variance=var(resxNonLinear_2sorted);​
 +x=[-5*sqrt(variance):​5*sqrt(variance)]';​
 +for i=1:​size(x,​1)
 +y(i)=(1/​2)*(1+erf( (x(i)-mean(resxNonLinear_2sorted))/​(sqrt(variance)*sqrt(2))));​
 +end
 +hold on 
 +plot(x,​y,'​b-'​);​
 +legend('​observed','​ideal'​)
 +title('​Cumulative distribution function'​)
 +end
 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +function [x,​resx,​Jx]= sinusNonLinearFit(datax,​datay,​x0)
 +global gdatax
 +global gdatay
 +gdatax=datax;​
 +gdatay=datay;​
 +[res0] = myfun2(x0);
 +options = optimset('​Jacobian','​on','​Display','​iter',​ '​MaxFunEvals',​200,​ '​Algorithm' ​ ,​{'​levenberg-marquardt',​.005});​
 +[x,​resnorm,​residual,​exitflag] = lsqnonlin(@myfun2,​x0,​-Inf,​+Inf,​options); ​   % Invoke optimizer
 +%[resxPerfect] = myfun2(xPerfect);​
 +%norm(resxPerfect,​2)
 +[resx,Jx] = myfun2(x);
 +norm(resx,​2)
 +end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +function [res,J] = myfun2(x)
 +global gpuls
 +global gdatax
 +global gdatay
 +ofs=x(1);
 +amp=x(2);
 +phas=x(3);
 +nbdata=max(size(gdatax));​
 +res=zeros(nbdata,​1);​
 +J=zeros(nbdata,​3);​
 +%J2=zeros(nbdata,​3);​
 +for i=1:nbdata
 +    t=gdatax(i);​
 +    ymes=gdatay(i);​
 +    y=ofs+amp*sin(gpuls*t+phas);​
 +    res(i)=ymes-y;​
 +    %jacobian for residuals
 +    J(i,:​)=[ ​                                        -1,
 +        - cos(gpuls*t)*sin(phas) - sin(gpuls*t)*cos(phas),​
 +        amp*sin(gpuls*t)*sin(phas) - amp*cos(gpuls*t)*cos(phas)];​
 +    %jacobian for residuals^2
 +    %J2(i,:​)=[ ​                                                       2*ofs - 2*ymes + 2*amp*cos(puls*t)*sin(phas) + 2*amp*sin(puls*t)*cos(phas),​
 +    %    2*(cos(puls*t)*sin(phas) + sin(puls*t)*cos(phas))*(ofs - ymes + amp*cos(puls*t)*sin(phas) + amp*sin(puls*t)*cos(phas)),​
 +    %    2*(amp*cos(puls*t)*cos(phas) - amp*sin(puls*t)*sin(phas))*(ofs - ymes + amp*cos(puls*t)*sin(phas) + amp*sin(puls*t)*cos(phas))];​
 +end
 +end
 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +function [res] = displaySin(tmin,​tmax,​x,​style)
 +global gpuls
 +ofs=x(1);
 +amp=x(2);
 +phas=x(3);
 +t=tmin:​0.1:​tmax;​
 +recons=ofs+amp*sin(gpuls.*t+phas);​
 +plot(t,​recons,​style)
 +end
 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +function [x,​resx,​Jx]= sinusNonLinearFit2(datax,​datay,​x0)
 +global gdatax
 +global gdatay
 +gdatax=datax;​
 +gdatay=datay;​
 +%[res0] = myfun3(x0);
 +options = optimset('​Jacobian','​on','​Display','​iter',​ '​MaxFunEvals',​200,​ '​Algorithm' ​ ,​{'​levenberg-marquardt',​.005});​
 +%options = optimset('​Jacobian','​off','​Display','​iter',​ '​MaxFunEvals',​200,​ '​Algorithm' ​ ,​{'​levenberg-marquardt',​.005});​
 +[x,​resnorm,​resx,​exitflag,​output,​lambda,​Jx] = lsqnonlin(@myfun3,​x0,​-Inf,​+Inf,​options); ​   % Invoke optimizer
 +
 +%perturbate to add bias
 +%resx=resx+100*ones(size(resx,​2));​
 +
 +%[resxPerfect] = myfun2(xPerfect);​
 +%norm(resxPerfect,​2)
 +%[resx,Jx] = myfun3(x);
 +norm(resx,​2)
 +end
 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +function [res,J] = myfun3(x)
 +global gpuls
 +global gdatax
 +global gdatay
 +ofs=x(1);
 +amp=x(2);
 +phas=x(3);
 +nla=x(4);
 +nlb=x(5);
 +nlc=x(6);
 +
 +nbdata=max(size(gdatax));​
 +res=zeros(nbdata,​1);​
 +J=zeros(nbdata,​6);​
 +%J2=zeros(nbdata,​3);​
 +for i=1:nbdata
 +    t=gdatax(i);​
 +    ymes=gdatay(i);​
 +    yperfect=ofs+amp*sin(gpuls*t+phas);​
 +    y=nla * yperfect^2 + nlb * yperfect + nlc;
 +    res(i)=ymes-y;​
 +    %jacobian for residuals
 +    J(i,:​)=[ ​                                                                                                                                                                                                                                                         - nlb - 2*nla*ofs - 2*amp*nla*cos(gpuls*t)*sin(phas) - 2*amp*nla*sin(gpuls*t)*cos(phas)
 +        - nlb*cos(gpuls*t)*sin(phas) - nlb*sin(gpuls*t)*cos(phas) - 2*amp*nla*cos(gpuls*t)^2*sin(phas)^2 - 2*amp*nla*sin(gpuls*t)^2*cos(phas)^2 - 2*nla*ofs*cos(gpuls*t)*sin(phas) - 2*nla*ofs*sin(gpuls*t)*cos(phas) - 4*amp*nla*cos(gpuls*t)*sin(gpuls*t)*cos(phas)*sin(phas)
 +        amp*nlb*sin(gpuls*t)*sin(phas) - amp*nlb*cos(gpuls*t)*cos(phas) - 2*amp*nla*ofs*cos(gpuls*t)*cos(phas) - 2*amp^2*nla*cos(gpuls*t)*sin(gpuls*t)*cos(phas)^2 + 2*amp^2*nla*cos(gpuls*t)*sin(gpuls*t)*sin(phas)^2 + 2*amp*nla*ofs*sin(gpuls*t)*sin(phas) - 2*amp^2*nla*cos(gpuls*t)^2*cos(phas)*sin(phas) + 2*amp^2*nla*sin(gpuls*t)^2*cos(phas)*sin(phas)
 +        - ofs^2 - amp^2*cos(gpuls*t)^2*sin(phas)^2 - amp^2*sin(gpuls*t)^2*cos(phas)^2 - 2*amp*ofs*cos(gpuls*t)*sin(phas) - 2*amp*ofs*sin(gpuls*t)*cos(phas) - 2*amp^2*cos(gpuls*t)*sin(gpuls*t)*cos(phas)*sin(phas)
 +        - ofs - amp*cos(gpuls*t)*sin(phas) - amp*sin(gpuls*t)*cos(phas)
 +        -1]';
 +    ​
 +    ​
 +    ​
 +    %    [                                         -1,
 +    %    - cos(gpuls*t)*sin(phas) - sin(gpuls*t)*cos(phas),​
 +    %    amp*sin(gpuls*t)*sin(phas) - amp*cos(gpuls*t)*cos(phas)];​
 +    %jacobian for residuals^2
 +    %J2(i,:​)=[ ​                                                       2*ofs - 2*ymes + 2*amp*cos(puls*t)*sin(phas) + 2*amp*sin(puls*t)*cos(phas),​
 +    %    2*(cos(puls*t)*sin(phas) + sin(puls*t)*cos(phas))*(ofs - ymes + amp*cos(puls*t)*sin(phas) + amp*sin(puls*t)*cos(phas)),​
 +    %    2*(amp*cos(puls*t)*cos(phas) - amp*sin(puls*t)*sin(phas))*(ofs - ymes + amp*cos(puls*t)*sin(phas) + amp*sin(puls*t)*cos(phas))];​
 +end
 +end
 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +function [res] = displaySinNonLinear(tmin,​tmax,​x,​style)
 +global gpuls
 +ofs=x(1);
 +amp=x(2);
 +phas=x(3);
 +nla=x(4);
 +nlb=x(5);
 +nlc=x(6);
 +t=tmin:​0.1:​tmax;​
 +reconsperfect=ofs+amp*sin(gpuls.*t+phas);​
 +recons=nla * reconsperfect.^2 + nlb * reconsperfect + nlc;
 +plot(t,​recons,​style)
 +end
 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +function [X covarianceOut]=fitSinusLinear(datax,​datay,​puls)
 +%linear method for known pulsation
 +%https://​fr.mathworks.com/​matlabcentral/​fileexchange/​41246-sine-function-fit
 +global gnoiseLevel
 +symbolicComputation=false
 +nbdata=max(size(datax));​
 +A=zeros(nbdata,​3);​
 +B=zeros(nbdata,​1);​
 +for i=1:nbdata
 +    t=datax(i);
 +    ymes=datay(i);​
 +    B(i)=ymes;
 +    dx=puls*t;
 +    A(i,:)= [1 cos(dx) sin(dx)] ​   ;
 +end
 +C=(A'​*A)^(-1)*A';​
 +%solve the system by expressing it as a linear product
 +Xa=C*B;
 +%propagate uncertainty from datay to Xa
 +estimateVarianceFromMSE=true
 +if ~estimateVarianceFromMSE
 +    varianceIn=gnoiseLevel^2 * eye(nbdata);​
 +else
 +    %in real case, the input variance should be estimated thanks to the MSE
 +    %varianceIn=zeros(nbdata);​
 +    %for i=1:nbdata
 +    %t=datax(i);​
 +    %ymes=datay(i);​
 +    %dx=puls*t;
 +    %ypred=(Xa(1)*1+ Xa(2)* cos(dx) + Xa(3)* sin(dx));
 +    %err=ymes- ypred;
 +    %varianceIn(i,​i)=err^2;​
 +    %end
 +    err=A*Xa-B;
 +    varianceIn=diag(err.*err);​
 +end
 +%covariance for the internal parameterization
 +covarianceOutIntern=C*varianceIn*C';​
 +%check if there is covariance (diag matrix)
 +covarianceOutIntern-diag([covarianceOutIntern(1,​1),​ covarianceOutIntern(2,​2),​ covarianceOutIntern(3,​3)]);​
 +%this uncertainty has to be transfered to the final parameterization
 +if symbolicComputation
 +    a=sym('​a','​real'​);​
 +    a=a;
 +    b1 = sym('​b1','​real'​);​
 +    b2 = sym('​b2','​real'​);​
 +    b=sqrt(b1*b1+b2*b2);​
 +    %c=wrapTo2Pi(atan2(b1,​b2));​ %bring c between 0 and 2pi
 +    c=atan(b1/​b2);​ %bring c between 0 and 2pi
 +    jacobian(b,​[b1,​b2]) ​ %  [ b1/(b1^2 + b2^2)^(1/​2),​ b2/(b1^2 + b2^2)^(1/​2)]
 +    jacobian(c,​[b1,​b2]) %  [ 1/​(b2*(b*1^2/​b2^2 + 1)), -b1/​(b2^2*(b1^2/​b2^2 + 1))]
 +    jacobian([a,​b,​c],​[a,​b1,​b2])
 +    %[1, 0 , 0;
 +    %0,  b1/(b1^2 + b2^2)^(1/​2),​ b2/(b1^2 + b2^2)^(1/​2);​
 +    %0, 1/​(b2*(b1^2/​b2^2 + 1)), -b1/​(b2^2*(b1^2/​b2^2 + 1))]
 +end
 +a=Xa(1);
 +b1=Xa(2);
 +b2=Xa(3);
 +%compute b & c
 +b=sqrt(b1*b1+b2*b2);​
 +c=wrapTo2Pi(atan2(b1,​b2));​ %bring c between 0 and 2pi
 +%compare c with
 +%wt*ud
 +JJ=[1, 0 , 0;
 +    0,  b1/(b1^2 + b2^2)^(1/​2),​ b2/(b1^2 + b2^2)^(1/​2);​
 +    0, 1/​(b2*(b1^2/​b2^2 + 1)), -b1/​(b2^2*(b1^2/​b2^2 + 1))]
 +covarianceOut=JJ*covarianceOutIntern*JJ'​
 +covarianceOut-diag([covarianceOut(1,​1),​ covarianceOut(2,​2),​ covarianceOut(3,​3)]);​
 +X=[a b c];
 +if symbolicComputation
 +    cdx1 = sym('​cdx1','​real'​);​
 +    sdx1 = sym('​sdx1','​real'​);​
 +    y1=    sym('​y1', ​ '​real'​);​
 +    cdx2 = sym('​cdx2','​real'​);​
 +    sdx2 = sym('​sdx2','​real'​);​
 +    y2=    sym('​y2', ​ '​real'​);​
 +    cdx3 = sym('​cdx3','​real'​);​
 +    sdx3 = sym('​sdx3','​real'​);​
 +    y3=    sym('​y3', ​ '​real'​);​
 +    cdx4 = sym('​cdx4','​real'​);​
 +    sdx4 = sym('​sdx4','​real'​);​
 +    y4=    sym('​y4', ​ '​real'​);​
 +    A=[1 cdx1 sdx1;
 +        1 cdx2 sdx2;
 +        1 cdx3 sdx3]
 +    B=[y1;​y2;​y3]
 +    solution=simplify((A'​*A)^(-1)*A'​*B)
 +    A=[1 cdx1 sdx1;
 +        1 cdx2 sdx2;
 +        1 cdx3 sdx3;
 +        1 cdx4 sdx4]
 +    B=[y1;​y2;​y3;​y4]
 +    C=simplify((A'​*A)^(-1)*A'​)
 +    leastSquaresSolution=simplify(C*B)
 +    %practically,​ an incremental algorithm should keep track of the
 +    %symetric matrix
 +    %(A'​*A)=[ ​                        ​4, ​                    cdx1 + cdx2 + cdx3 + cdx4,                     sdx1 + sdx2 + sdx3 + sdx4]
 +    %       [ cdx1 + cdx2 + cdx3 + cdx4,             ​cdx1^2 + cdx2^2 + cdx3^2 + cdx4^2, cdx1*sdx1 + cdx2*sdx2 + cdx3*sdx3 + cdx4*sdx4]
 +    %       [ sdx1 + sdx2 + sdx3 + sdx4, cdx1*sdx1 + cdx2*sdx2 + cdx3*sdx3 + cdx4*sdx4, ​            ​sdx1^2 + sdx2^2 + sdx3^2 + sdx4^2]
 +    %and the vector
 +    %(A'​*B)= [y1 + y2 + y3 + y4
 +    %         ​cdx1*y1 + cdx2*y2 + cdx3*y3 + cdx4*y4
 +    %         ​sdx1*y1 + sdx2*y2 + sdx3*y3 + sdx4*y4]
 +    var1 = sym('​var1','​real'​);​
 +    var2 = sym('​var2','​real'​);​
 +    var3 = sym('​var3','​real'​);​
 +    var4 = sym('​var4','​real'​);​
 +    %for constant known variance on input data
 +    varIn=diag([var1,​var1,​var1,​var1]);​
 +    simplify(A'​*varIn*A)
 +    %[                           ​4*var1, ​                    ​var1*(cdx1 + cdx2 + cdx3 + cdx4), ​                    ​var1*(sdx1 + sdx2 + sdx3 + sdx4)]
 +    %[ var1*(cdx1 + cdx2 + cdx3 + cdx4), ​            ​var1*(cdx1^2 + cdx2^2 + cdx3^2 + cdx4^2), var1*(cdx1*sdx1 + cdx2*sdx2 + cdx3*sdx3 + cdx4*sdx4)]
 +    %[ var1*(sdx1 + sdx2 + sdx3 + sdx4), var1*(cdx1*sdx1 + cdx2*sdx2 + cdx3*sdx3 + cdx4*sdx4), ​            ​var1*(sdx1^2 + sdx2^2 + sdx3^2 + sdx4^2)]
 +    %for variances obtained from MSE
 +    varIn=diag([var1,​var2,​var3,​var4]);​
 +    simplify(A'​*varIn*A)
 +    %[                     var1 + var2 + var3 + var4,                     ​cdx1*var1 + cdx2*var2 + cdx3*var3 + cdx4*var4, ​                    ​sdx1*var1 + sdx2*var2 + sdx3*var3 + sdx4*var4]
 +    %[ cdx1*var1 + cdx2*var2 + cdx3*var3 + cdx4*var4, ​            ​var1*cdx1^2 + var2*cdx2^2 + var3*cdx3^2 + var4*cdx4^2,​ cdx1*sdx1*var1 + cdx2*sdx2*var2 + cdx3*sdx3*var3 + cdx4*sdx4*var4]
 +    %[ sdx1*var1 + sdx2*var2 + sdx3*var3 + sdx4*var4, cdx1*sdx1*var1 + cdx2*sdx2*var2 + cdx3*sdx3*var3 + cdx4*sdx4*var4, ​            ​var1*sdx1^2 + var2*sdx2^2 + var3*sdx3^2 + var4*sdx4^2]
 +end
 +end
 +% A.X=B
 +% A'​.A.X=A'​.B
 +%p 31 of http://​www.stat.columbia.edu/​~fwood/​Teaching/​w4315/​Fall2009/​lecture_11
 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +function displayResult(res,​covar,​name)
 +%display mean and 3Sigma values
 +nbdata=max(size(res));​
 +s='';​
 +s=sprintf('​%s :',​name);​
 +for n=1:nbdata
 +    s=sprintf('​%s %d +/- %d  ,',​s,​res(n),​3*sqrt(covar(n,​n)));​
 +end
 +disp(s);
 +end
 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 +function test()
 +%computation of the inverse of a symetric 3x3 matrix
 +if false
 +    a1 = sym('​ATAa1','​real'​);​
 +    b1 = sym('​ATAb1','​real'​);​
 +    c1 = sym('​ATAc1','​real'​);​
 +    d1 = sym('​ATAd1','​real'​);​
 +    e1 = sym('​ATAe1','​real'​);​
 +    f1 = sym('​ATAf1','​real'​);​
 +    A1=[a1 b1 c1 ; b1 d1 e1; c1 e1 f1]
 +    ​
 +    %     a2 = sym('​a2','​real'​);​
 +    %     b2 = sym('​b2','​real'​);​
 +    %     c2 = sym('​c2','​real'​);​
 +    %     d2 = sym('​d2','​real'​);​
 +    %     e2 = sym('​e2','​real'​);​
 +    %     f2 = sym('​f2','​real'​);​
 +    %     ​A2=[a2 b2 c2 ; b2 d2 e2; c2 e2 f2]
 +    %     ​prod=A1*A2
 +    %     ​equations=reshape(prod-eye(3),​9,​1)
 +    %     ​solve(equations(1:​6),​[a2 b2 c2  d2 e2  f2])
 +    %     ​solve(equations(1:​6),​a2 )
 +    ​
 +    % http://​www.tutorvista.com/​content/​math/​inverse-of-symmetric-matrix/​
 +    % A-1= adj(A)/​det(A)
 +    A2tmp=adjoint(A1)
 +    determinant=det(A1)
 +    A2=A2tmp/​determinant
 +end
 +end
 +
 +</​file>​
 +
 +<file matlab fitSinusLinearOptimized.m>​
 +function [X covarianceOut,​ err]=fitSinusLinearOptimized(datax,​datay,​puls)
 +nbdata=max(size(datax));​
 +A=zeros(nbdata,​3);​
 +B=zeros(nbdata,​1);​
 +
 +err=zeros(nbdata,​1);​
 +var=zeros(nbdata,​1);​
 +
 +ATA=zeros(3,​3);​
 +ATB=zeros(3,​1);​
 +ATAa=0;
 +ATAb=0;
 +ATAc=0;
 +ATAd=0;
 +ATAe=0;
 +ATAf=0;
 +ATBa=0;
 +ATBb=0;
 +ATBc=0;
 +
 +cdx=zeros(nbdata,​1);​
 +sdx=zeros(nbdata,​1);​
 +cdx2=zeros(nbdata,​1);​
 +sdx2=zeros(nbdata,​1);​
 +cdxsdx=zeros(nbdata,​1);​
 +
 +for i=1:nbdata
 +    x=datax(i);
 +    y=datay(i);
 +    B(i)=y;
 +    dx=puls*x;
 +    cx=cos(dx);
 +    sx=sin(dx);
 +    cx2=cx*cx;
 +    sx2=sx*sx;
 +    cxsx=cx*sx;
 +    cdx(i)=cx;
 +    sdx(i)=sx;
 +    cdx2(i)=cx2;​
 +    sdx2(i)=sx2;​
 +    cdxsdx(i)=cxsx;​
 +    ​
 +    ATAa=ATAa+1;​
 +    ATAb=ATAb+cx;​
 +    ATAc=ATAc+sx;​
 +    ATAd=ATAd+cx2;​
 +    ATAe=ATAe+cxsx;​
 +    ATAf=ATAf+sx2;​
 +    ​
 +    ATBa=ATBa+y;​
 +    ATBb=ATBb+cx*y;​
 +    ATBc=ATBc+sx*y;​
 +end
 +%complete symetric matrices
 +
 +%ATA=[ATAa , ATAb , ATAc ; ATAb , ATAd , ATAe ; ATAc , ATAe , ATAf];
 +%ATB=[ATBa , ATBb , ATBc ]';
 +%leastSquaresSolution=inv(ATA)* ATB
 +%a=leastSquaresSolution(1);​
 +%b1=leastSquaresSolution(2);​
 +%b2=leastSquaresSolution(3);​
 +
 +determinant =- ATAf*ATAb^2 + 2*ATAb*ATAc*ATAe - ATAd*ATAc^2 - ATAa*ATAe^2 + ATAa*ATAd*ATAf;​
 +ATAinva= (ATAd*ATAf - ATAe^2 ​    ​)/​determinant;​
 +ATAinvb= (ATAc*ATAe - ATAb*ATAf ​ )/​determinant;​
 +ATAinvc= (ATAb*ATAe - ATAc*ATAd ​ )/​determinant;​
 +ATAinvd= (ATAa*ATAf - ATAc^2 ​    ​)/​determinant;​
 +ATAinve= (ATAb*ATAc - ATAa*ATAe ​ )/​determinant;​
 +ATAinvf= (ATAa*ATAd - ATAb^2 ​    ​)/​determinant;​
 +
 +ATAinv=[ ATAinva ATAinvb ATAinvc; ATAinvb ATAinvd ATAinve; ATAinvc ATAinve ATAinvf];
 +
 +a =ATAinva * ATBa + ATAinvb * ATBb + ATAinvc * ATBc;
 +b1=ATAinvb * ATBa + ATAinvd * ATBb + ATAinve * ATBc;
 +b2=ATAinvc * ATBa + ATAinve * ATBb + ATAinvf * ATBc;
 +Xa=[a b1 b2]
 +%compute b & c
 +b=sqrt(b1*b1+b2*b2);​
 +c=wrapTo2Pi(atan2(b1,​b2));​ %bring c between 0 and 2pi
 +
 +X=[a b c ];
 +
 +%compute MSE and variance....
 +ATVarAa=0;
 +ATVarAb=0;
 +ATVarAc=0;
 +ATVarAd=0;
 +ATVarAe=0;
 +ATVarAf=0;
 +
 +for i=1:nbdata
 +    er=  1 * a + cdx(i,1) * b1 + sdx(i,1) * b2;
 +    vari= er*er;
 +    err(i)=er;
 +    var(i)=vari;​
 +    ATVarAa= ATVarAa + vari;
 +    ATVarAb= ATVarAb + cdx(i,​1) ​ * vari;
 +    ATVarAc= ATVarAc + sdx(i,​1) ​ * vari;
 +    ATVarAd= ATVarAd + cdx2(i,1) * vari;
 +    ATVarAe= ATVarAe + cdx(i,​1)*sdx(i,​1) ​ * vari;
 +    ATVarAf= ATVarAf + sdx2(i,1) * vari;
 +    ​
 +end
 +ATVarA=[ ATVarAa ATVarAb ATVarAc; ATVarAb ATVarAd ATVarAe; ATVarAc ATVarAe ATVarAf];
 +covarianceOutIntern=ATAinv'​*ATVarA*ATAinv;​
 +JJ=[1, 0 , 0;
 +    0,  b1/(b1^2 + b2^2)^(1/​2),​ b2/(b1^2 + b2^2)^(1/​2);​
 +    0, 1/​(b2*(b1^2/​b2^2 + 1)), -b1/​(b2^2*(b1^2/​b2^2 + 1))];
 +covarianceOut=JJ*covarianceOutIntern*JJ';​
 +end
 +</​file>​
  
deflectometry.1509461793.txt.gz · Dernière modification: 2017/10/31 15:56 par bvandepo