%%
%マスク付き2次元フィッティングの基本コード
%
%製作者：堀越宗一
%最終更新日：2019/08/10

%%
%m_dataという２次元テストデータを準備

x=1:300;
y=1:250;
[X,Y]=meshgrid(x,y);

A=1+0.2*(rand-0.5);
cx=150+20*(rand-0.5);
cy=125+20*(rand-0.5);
wx=20+5*(rand-0.5);
wy=20+5*(rand-0.5);
y0=0.1*(rand-0.5)+0.5*(rand(length(y),length(x))-0.5);

m_data=A*exp(-0.5*(X-cx).^2/wx^2).*exp(-0.5*(Y-cy).^2/wy^2)+y0;
m_noise=2.5*exp(-0.5*(X-70).^2/15^2).*exp(-0.5*(Y-40).^2/5^2);

m_data=m_data+m_noise;

%%
%ここからマスク付き2次元フィッティング
%スポットノイズを持つm_dataという２次元データを解析したいとする

[ysize,xsize]=size(m_data);%２次元データのサイズ取得
x=1:xsize;
y=1:ysize;
[X,Y]=meshgrid(x,y);
XY=[X(:),Y(:)];

%マスクの設定
Mask_cx=70;%マスク中心x座標
Mask_cy=40;%マスク中心y座標
Mask_w=120;%マスク幅
Mask_h=40;%マスク高さ
Mask_xmin=Mask_cx-Mask_w/2;
Mask_xmax=Mask_cx+Mask_w/2;
Mask_ymin=Mask_cy-Mask_h/2;
Mask_ymax=Mask_cy+Mask_h/2;

%マスク外側のデータ座標
k=find(not(XY(:,1)>Mask_xmin & XY(:,1)<Mask_xmax & XY(:,2)>Mask_ymin & XY(:,2)<Mask_ymax));

%フィッティング関数（ここでは２次元ガウス分布）
%p(1):振幅
%p(2):x方向中心座標
%p(3):x方向幅
%p(4):y方向中心座標
%p(5):y方向幅
%p(6):オフセット
f_fit = @(p,XY) p(1)*exp(-0.5*(XY(:,1)-p(2)).^2/p(3)^2).*exp(-0.5*(XY(:,2)-p(4)).^2/p(5)^2)+p(6);

p_ini=[1,150,20,125,20,0];%フィッティング初期値（データに合わせて入力）
p_lb=[0,0,1,0,1,-1];%フィッティング下限値（データに合わせて入力）
p_ub=[2,300,100,250,100,1];%%フィッティング上限値（データに合わせて入力）

[p_fit,resnorm,resid,exitflag,output,lambda,J] = lsqcurvefit(f_fit,p_ini,XY(k,:),m_data(k),p_lb,p_ub);%マスク外側のデータでフィッティング
p_err = nlparci(p_fit,resid,'jacobian',J);%95% の信頼区間を計算

p_out=[p_fit',p_err];%フィッティング結果
m_fit=reshape(f_fit(p_fit,XY),[ysize,xsize]);%２次元フィッティングデータ

%%
%プロット

subplot(2,2,1)
imagesc(m_data)
hold on;
rectangle('Position',[Mask_cx-Mask_w/2 Mask_cy-Mask_h/2 Mask_w Mask_h],'LineStyle','--','EdgeColor','g')
hold off;
axis xy;
axis equal;
title 'Data';

subplot(2,2,2)
imagesc(m_fit)
axis xy;
axis equal;
title 'Fitting';

subplot(2,2,3)
plot(x,sum(m_data,1),'k-')
hold on;
plot(x,sum(m_fit,1),'r-')
plot([Mask_xmin Mask_xmin],[-20 100],'g--')
plot([Mask_xmax Mask_xmax],[-20 100],'g--')
hold off;
xlabel('x')
legend ('Data','Fitthing','Mask')
xlim([0,300])
ylim([-20,100])

subplot(2,2,4)
plot(y,sum(m_data,2),'k-')
hold on;
plot(y,sum(m_fit,2),'r-')
plot([Mask_ymin Mask_ymin],[-20 100],'g--')
plot([Mask_ymax Mask_ymax],[-20 100],'g--')
hold off;
xlabel('y')
legend ('Data','Fitthing','Mask')
xlim([0,250])
ylim([-20,100])