MATLAB处理自旋极化电子显微镜数据自动化（2）- compound color image

单峰及双峰高斯拟合后得到Asymmetry图像—SPLEEM001(ip1,ip2,pp)

function [SPLEEM] = SPLEEM001(ip1,ip2,pp)
% Get SPLEEM asymmetry fitting results

% fitting
ip1 = double(ip1);
ip2 = double(ip2);
pp  = double(pp);
% binwidth; m; sigma
ip1_fit = Gfit_singlepeak(ip1,4,32890,0.4);
ip2_fit = Gfit_singlepeak(ip2,6,32850,0.4);
% binwidth; m; sigma1; sigma2, peak1; peak2
pp_fit  = Gfit_doublepeak(pp,10,32880,0.2,0.2,0.2,0.2);

% Calculate asymmetry results
ip1 = (ip1 - ip1_fit.c)/(2^15-1);
ip2 = (ip2 - ip2_fit.c)/(2^15-1);
pp = (pp - pp_fit.c)/(2^15-1);
% ip1 = (ip1 - 30231)/(2^15-1);
% ip2 = (ip2 - 30188)/(2^15-1);
% pp = (pp - 31850)/(2^15-1);

SPLEEM.ip1 = ip1;
SPLEEM.ip2 = ip2;
SPLEEM.pp  = pp;
% % denoise
% f = fspecial('average',[2,2]);
% SPLEEM.ip1 = imfilter(ip1,f);
% SPLEEM.ip2 = imfilter(ip2,f);
% SPLEEM.pp  = imfilter(pp,f);

% ip1 fitting result plot
cl = [1 1 1]*0.7;
fsize = 14;
figure(1)
clf
hold on
scatter(ip1_fit.xBins,ip1_fit.yHist,'linewidth',.75)
% plot(xBins1,yHist1,'linewidth',2,'color','k')
plot(ip1_fit.xBins,ip1_fit.yFit,'linewidth',2,'color','r')
txt = ['\leftarrow',num2str(round(ip1_fit.c))];
text(ip1_fit.c,ip1_fit.intMax/5,txt)
line([0 0]+ 2^15-1,[0 1]*ip1_fit.intMax,'linewidth',1,'color',cl,'linestyle','--')
line([0 0]+ ip1_fit.c,[0 1]*ip1_fit.intMax,'linewidth',1,'color'...
    ,'r','linestyle','--')
box on
% set(gca,'XTick',-90:90:270);
xlim([ip1_fit.xBins(1) ip1_fit.xBins(end)]);
ylim([0 ip1_fit.intMax]);
xlabel('Asymmetry','interpreter','latex','fontsize',fsize)
ylabel('Counts','interpreter','latex','fontsize',fsize)
% set(gca,'DataAspectRatio',[1 1e3 1])
hold off

% ip2 fitting result plot
figure(2)
clf
hold on
scatter(ip2_fit.xBins,ip2_fit.yHist,'linewidth',.75)
% plot(xBins1,yHist1,'linewidth',2,'color','k')
plot(ip2_fit.xBins,ip2_fit.yFit,'linewidth',2,'color','r')
txt = ['\leftarrow',num2str(round(ip2_fit.c))];
text(ip2_fit.c,ip2_fit.intMax/5,txt)
line([0 0]+ 2^15-1,[0 1]*ip2_fit.intMax,'linewidth',1,'color',cl,'linestyle','--')
line([0 0]+ ip2_fit.c,[0 1]*ip2_fit.intMax,'linewidth',1,'color'...
    ,'r','linestyle','--')
box on
% set(gca,'XTick',-90:90:270);
xlim([ip2_fit.xBins(1) ip2_fit.xBins(end)]);
ylim([0 ip2_fit.intMax]);
xlabel('Asymmetry','interpreter','latex','fontsize',fsize)
ylabel('Counts','interpreter','latex','fontsize',fsize)
% set(gca,'DataAspectRatio',[1 1e3 1])
hold off

% pp fitting result plot
figure(3)
clf
hold on
scatter(pp_fit.xBins,pp_fit.yHist,'linewidth',.75)
% plot(xBins1,yHist1,'linewidth',2,'color','k')
plot(pp_fit.xBins,pp_fit.yFit,'linewidth',2,'color','r')
txt = ['\leftarrow',num2str(round(pp_fit.c))];
text(pp_fit.c,pp_fit.intMax/5,txt)
line([0 0]+ 2^15-1,[0 1]*pp_fit.intMax,'linewidth',1,'color',cl,'linestyle','--')
line([0 0]+ pp_fit.c,[0 1]*pp_fit.intMax,'linewidth',1,'color'...
    ,'r','linestyle','--')
box on
% set(gca,'XTick',-90:90:270);
xlim([pp_fit.xBins(1) pp_fit.xBins(end)]);
ylim([0 pp_fit.intMax]);
xlabel('Asymmetry','interpreter','latex','fontsize',fsize)
ylabel('Counts','interpreter','latex','fontsize',fsize)
% set(gca,'DataAspectRatio',[1 1e3 1])
hold off

end

function [fit] = Gfit_singlepeak(im,binwidth,m,sigma)
r = sort(im(:));
d = max(max(r))- min(min(r));
numBins = round(d/binwidth);
% single_peak Gaussian function
fun = @(c,x) ...
    c(1) * exp((x(:,1) - c(2)).^2/(-2*c(3)^2));
options = optimset('TolX',1e-6,'TolFun',1e-6,...
    'MaxFunEvals',10^4,'display','off');
[yHist,xHist] = histcounts(im,numBins);
xBins = (xHist(1:end-1) + xHist(2:end))/2;
intMax = max(yHist(:));
width = xBins(end) - xBins(1);
coefs = [intMax m sigma*width];
coefs = lsqcurvefit(fun,coefs,xBins(:),yHist(:),[],[],options);
yFit = fun(coefs,xBins(:));
c = coefs(1,2);
% package fitting results
fit.intMax = intMax;
fit.xBins = xBins;
fit.yHist = yHist;
fit.yFit = yFit;
fit.c = c;
fit.A = 0.00;
end

function [fit] = Gfit_doublepeak(im,binwidth,m,sigma1,sigma2,peak1,peak2)
r = sort(im);
d = max(max(r))- min(min(r));
numBins = round(d/binwidth);
% double_peak Gaussian function
fun = @(c,x) ...
    c(1) * exp((x(:,1) - c(2)).^2/(-2*c(3)^2)) + ...
    c(4) * exp((x(:,1) - c(5)).^2/(-2*c(6)^2));
options = optimset('TolX',1e-6,'TolFun',1e-6,...
    'MaxFunEvals',10^4,'display','off');
[yHist,xHist] = histcounts(im,numBins);
xBins = (xHist(1:end-1) + xHist(2:end))/2;
intMax = max(yHist(:));
width = xBins(end) - xBins(1);
coefs = [ ...
    intMax m+peak2*width sigma2*width...
    intMax m-peak1*width sigma1*width];
coefs = lsqcurvefit(fun,coefs,xBins(:),yHist(:),[],[],options);
yFit = fun(coefs,xBins(:));
c = (coefs(1,2)+coefs(1,5))/2;
% package fitting results
fit.intMax = intMax;
fit.xBins = xBins;
fit.yHist = yHist;
fit.yFit = yFit;
fit.c = c;
fit.A = (coefs(1,5)-coefs(1,2))/c/2;
end

不同自旋极化条件下图像合成彩图—SPLEEM002(SPLEEM)

function [m] = SPLEEM002(SPLEEM)
% Convert SPLEEM data to m
% m(:,:,1) point right
% m(:,:,2) point down
% m(:,:,3) point up perpedicular to screen

% Coordinate in spin-space (FOV = 25um)
% |------- y (ip158)
% |o z
% |
% |
% x (ip68)

% Increase the DW signal
seSize = 2; % DW width multiplication
A = 2; % IP-DW signal enhance
B = 0.9; % OOP-Domain signal decrease
f = fspecial('average',[2,2]); % mean filter

% I1 is the smaller ip direction
I1 = 90;
I2 = 180;
t1 = (I1 - 157)/180*pi;
t2 = (I2 - 247)/180*pi;
proj1 = [cos(t1) -sin(t2) 0];
proj2 = [sin(t1)  cos(t2) 0];
proj3 = [   0       0    -1];
T = [proj1;proj2;proj3];
N = size(SPLEEM.ip1);
xyz = zeros(N(1),N(2),3);

for i = 1:3
    xyz(:,:,i) = T(i,1)*SPLEEM.ip1 + T(i,2)*SPLEEM.ip2 ...
        + T(i,3)*SPLEEM.pp;
end

% Normalized vector direction
% |------- alpha
% |o gamma
% |
% |
% beta
alpha = xyz(:,:,2)./sqrt(xyz(:,:,1).^2 + xyz(:,:,2).^2 + xyz(:,:,3).^2);
beta  = xyz(:,:,1)./sqrt(xyz(:,:,1).^2 + xyz(:,:,2).^2 + xyz(:,:,3).^2);
gamma = xyz(:,:,3)./sqrt(xyz(:,:,1).^2 + xyz(:,:,2).^2 + xyz(:,:,3).^2);
alpha(isnan(alpha)) = 0;
beta(isnan(beta)) = 0;
gamma(isnan(gamma)) = 0;

% % increase the in-plane signal through increasing DW width in pp
% threshold = 0.95;
% gamma(abs(gamma)<=threshold) = gamma(abs(gamma)<=threshold)*0;

% Detect edges
thresh = [0.2 0.45];
sigmaEdge = 1.0;
Iedge = edge(gamma,'canny',thresh,sigmaEdge);
% figure(1)
% imagesc(Iedge)
% colormap(gray(256))
% axis equal off

% Dilate edges
se = strel('disk', seSize, 0);
Iedge_dilated = imdilate(Iedge, se);
% figure(2)
% imagesc(Iedge_dilated)
% colormap(gray(256))
% axis equal off

% % Dilate domain wall
% se = strel('disk',1,0);
% alpha = imdilate(abs(alpha),se).*sign(alpha);
% beta  = imdilate(abs(beta),se).*sign(alpha);

% adjust 3DSPLEEM contrast
alpha(Iedge_dilated) = alpha(Iedge_dilated)*A;
beta(Iedge_dilated) = beta(Iedge_dilated)*A;
alpha(~Iedge_dilated) = alpha(~Iedge_dilated)*B;
beta(~Iedge_dilated) = beta(~Iedge_dilated)*B;
gamma(~Iedge_dilated) = gamma(~Iedge_dilated)*A;
alpha = imfilter(alpha,f);
beta  = imfilter(beta,f);
gamma = imfilter(gamma,f);

% SPLEEMgrayplot
figure(1);
set(gca,'position',[0 0 1 1])
imagesc([SPLEEM.ip1 SPLEEM.ip2 SPLEEM.pp])
colormap(gray)
axis equal off

figure(2);
set(gca,'position',[0 0 1 1])
imagesc([alpha beta gamma])
colormap(gray)
axis equal off

%package
m(:,:,1) = alpha;
m(:,:,2) = beta;
m(:,:,3) = gamma;

% SPLEEMcolorplot
% contrast = 0.65;
theta = 0;
HSI = spin2hsi2(m,theta);
imcolorplot = hsi2rgb(HSI);
h4 = figure(4);
hold on
set(h4,'outerposition',[300 300 500 500])
imshow(imcolorplot);
imwrite(imcolorplot,'SPLEEM_3D.tiff');

% SPLEEMcolorplot-in-plane
theta = 0;
HSI = spin2hsi3(m,theta);
imcolorplot = hsi2rgb(HSI);
h5 = figure(5);
hold on
set(h5,'outerposition',[300 300 500 500])
imshow(imcolorplot);
imwrite(imcolorplot,'SPLEEM_3D_IP.tiff');

% plot colorwheel
h3 = figure(3);
set(h3,'outerposition',[300 300 300 300])
Colorwheel2(theta);

end

function [hsi] = spin2hsi(spin,theta,contrast)
% spin space to hsi space
alpha = spin(:,:,1);
beta  = spin(:,:,2);
gamma = spin(:,:,3);
H = atan2(-alpha,-beta) + pi/2;
H(H<0) = H(H<0) + 2*pi;
H = H + theta/180*pi;
H(H<0) = H(H<0) + 2*pi;
H(H>2*pi) = H(H>2*pi) - 2*pi;
S = min(1,contrast*(1-gamma));
I = min(1,contrast*(1+gamma));
m = alpha.^2 + beta.^2 + gamma.^2;
H(m==0) = 0;
S(m==0) = 0;
I(m==0) = 1;
hsi = cat(3,H,S,I);
end

function [hsi] = spin2hsi2(spin,theta)
% spin space to hsi space
alpha = spin(:,:,1);
beta  = spin(:,:,2);
gamma = spin(:,:,3);
H = atan2(-alpha,-beta) + pi/2;
H(H<0) = H(H<0) + 2*pi;
H = H + theta/180*pi;
H(H<0) = H(H<0) + 2*pi;
H(H>2*pi) = H(H>2*pi) - 2*pi;
% S = ones(size(H));
% I = ones(size(H))/2;
S = sqrt(alpha.^2 + beta.^2);
I = (gamma+1)/2;
m = alpha.^2 + beta.^2 + gamma.^2;
H(m==0) = 0;
S(m==0) = 0;
I(m==0) = 1;
hsi = cat(3,H,S,I);
end

function [hsi] = spin2hsi3(spin,theta)
% spin space to hsi space
alpha = spin(:,:,1);
beta  = spin(:,:,2);
gamma = spin(:,:,3);
H = atan2(-alpha,-beta) + pi/2;
H(H<0) = H(H<0) + 2*pi;
H = H + theta/180*pi;
H(H<0) = H(H<0) + 2*pi;
H(H>2*pi) = H(H>2*pi) - 2*pi;
S = ones(size(H));
I = ones(size(H))/2;
% S = sqrt(alpha.^2 + beta.^2);
% I = (gamma+1)/2;
m = alpha.^2 + beta.^2 + gamma.^2;
H(m==0) = 0;
S(m==0) = 0;
I(m==0) = 1;
hsi = cat(3,H,S,I);
end

function [rgb] = hsi2rgb(hsi)
% hsi space to rgb space
H = hsi(:,:,1);
S = hsi(:,:,2);
I = hsi(:,:,3);
R = zeros(size(hsi,1),size(hsi,2));
G = zeros(size(hsi,1),size(hsi,2));
B = zeros(size(hsi,1),size(hsi,2));
idx = find((0 <= H)&(H < 2*pi/3));
B(idx) = I(idx).*(1-S(idx));
R(idx) = I(idx).*(1+S(idx).*cos(H(idx))./cos(pi/3-H(idx)));
G(idx) = 3*I(idx)-(B(idx)+R(idx));
idx = find((2*pi/3 <= H)&(H < 4*pi/3));
R(idx) = I(idx).*(1-S(idx));
G(idx) = I(idx).*(1+S(idx).*cos(H(idx)-2*pi/3)...
    ./cos(pi-H(idx)));
B(idx) = 3*I(idx)-(G(idx)+R(idx));
idx = find((4*pi/3 <= H)&(H <= 2*pi));
G(idx) = I(idx).*(1-S(idx));
B(idx) = I(idx).*(1+S(idx).*cos(H(idx)-4*pi/3)...
    ./cos(5*pi/3-H(idx)));
R(idx) = 3*I(idx)-(G(idx)+B(idx));
rgb = cat(3,R,G,B);
rgb = max(min(rgb,1),0);
end

function [colorwheelRGBImage] = Colorwheel(r,azimuth)
x = -r:r;
x = x/r;
y = x;
[X,Y] = meshgrid(x,y);
phi = atan2(Y,X);
phi(phi<0) = phi(phi<0)+2*pi;
phi = phi + azimuth/180*pi;
phi(phi<0) = phi(phi<0) + 2*pi;
phi(phi>2*pi) = phi(phi>2*pi) - 2*pi;
hue = phi/2/pi;
saturation = sqrt(X.^2+Y.^2);
saturation(saturation>1) = 0;
value = ones(size(X));
imHSV = cat(3,hue,saturation,value);
colorwheelRGBImage = hsv2rgb(imHSV);
colorwheelRGBImage = flipud(colorwheelRGBImage);
end

function [] = Colorwheel2(theta)
x0 = 0;
y0 = 0;
t = (0:.02:2)*pi;
r1 = 0:.02:.3;   % interior white
r2 = .3:.02:.8;  % colorwheel
r3 = .8:.02:1;   % exterior black
xp1 = x0+cos(t)'*r1;
yp1 = y0+sin(t)'*r1;
xp2 = x0+cos(t)'*r2;
yp2 = y0+sin(t)'*r2;
xp3 = x0+cos(t)'*r3;
yp3 = y0+sin(t)'*r3;
surface(xp2,yp2,t'*(r2 == r2),'EdgeColor','none');
surface(xp3,yp3,t'*(r3 == r3),'FaceColor',[0 0 0]);
surface(xp1,yp1,t'*(r1 == r1),'FaceColor',[1 1 1],'EdgeColor','none');
colormap(hsv(256))
axis equal off
view(theta,90)
end

function [] = Colorwheel3(azimuth)
r = linspace(0,1,100);
theta = linspace(0,2*pi,1000);
[rg,thg] = meshgrid(r,theta);
[x,y] = pol2cart(thg,rg);
pcolor(x,y,thg);
colormap(hsv)
view(azimuth,90)
shading flat
axis equal off
end