Code covered by the BSD License  

Highlights from
immiscible LB

immiscible LB

by

 

23 Jul 2009 (Updated )

Implements Immiscible Lattice Boltzmann (ILB, D2Q9) method for two phase flows

image_curvature.m
close all

% I= single(imread('eight.tif')); 
 
 
% I= zeros(30,30); I(10:20,15:25)=1; I(15:18,5:15)=1; I(20:25,25:28)=1;

 if(1)
mask=zeros(20,20); N=100; radius=40;
xax = (1:N)-ceil(N/2); 
x = repmat(xax, N, 1);  y= repmat(xax', 1, N);  
mask= ((x-0).^2+(y-0).^2) <= radius.^2 ; % mask the inner circle of diam N
imshow (mask,[]);
I=double(mask);

 end
 
 
 figure, imshow(I,[]);
 [n,m]=size(I);
 
 X=repmat((1:m),n,1);
 Y=repmat((1:n)',1,m);
  
[K,H,Pmax,Pmin] = surfature(X,Y,I);
% K AND H ARE THE GAUSSIAN AND MEAN CURVATURES, RESPECTIVELY.
figure, imshow(K,[]); title('gaussian')

figure, imshow(H,[]); title(' mean')



 [gradPhiX gradPhiY]=gradient(I);
 % magnitude of gradient of I 
 absGradPhi=sqrt(gradPhiX.^2+gradPhiY.^2);
 % normalized gradient of I - eliminating singularities
 normGradPhiX=gradPhiX./(absGradPhi+(absGradPhi==0));
 normGradPhiY=gradPhiY./(absGradPhi+(absGradPhi==0));

 [divXnormGradPhiX divYnormGradPhiX]=gradient(normGradPhiX);
 [divXnormGradPhiY divYnormGradPhiY]=gradient(normGradPhiY);
 % curvature is the divergence of normalized gradient of I
 
K = -(divXnormGradPhiX + divYnormGradPhiY);
 
figure, imshow(K,[]); title(' div norm grad');


K_Lish = gradPhiX.*gradPhiY.*(divXnormGradPhiY + divYnormGradPhiX) ...
    -(gradPhiY.^2).*divXnormGradPhiX - (gradPhiX.^2).*divYnormGradPhiY ;

figure, imshow(K_Lish,[]); title('lishchuk');


if(0)
Cx_2d=zeros(3,3); Cy_2d=zeros(3,3); 
% STENCILS for Gradient calculatiob
Cx_2d=[[-ones(3,1)],[zeros(3,1)],[ones(3,1)]]; % for Colour Gradient x calculation
Cy_2d=[[ones(1,3)];[zeros(1,3)];[-ones(1,3)]]; % For CGy calculation

else
Cy_2d=zeros(3,3); Cy_2d(1,2)=1; Cy_2d(3,2)=-1;
Cx_2d=zeros(3,3); Cx_2d(2,1)=-1; Cx_2d(2,3)=1;
end

 CGx=imfilter(I, Cx_2d,'circular','same'); % Grad x
 CGy=imfilter(I,-Cy_2d,'circular','same'); % Grad y 
 
 M_Grad=sqrt(CGx.^2+CGy.^2); % Mod CG (Magnitude )
 %M_Grad=abs(CGx) + abs(CGy); % Mod CG (fast Magnitude as abs value)
 CGy=CGy./(M_Grad+(M_Grad==0));   CGx=CGx./(M_Grad+(M_Grad==0)); 
% CGy=CGy./(M_Grad);   CGx=CGx./(M_Grad); 
  
 CGxx=imfilter(CGx, Cx_2d,'circular','same'); % Grad x^2
 CGyy=imfilter(CGy,-Cy_2d,'circular','same'); % Grad y^2
 
 CGxy=imfilter(CGx, -Cy_2d,'circular','same'); % Grad y Grad x
 CGyx=imfilter(CGy, Cx_2d,'circular','same');  % Grad x Grad y 

KK_Lish = CGx.*CGy.*(CGyx + CGxy) ...
    -(CGy.^2).*CGxx - (CGx.^2).*CGyy ;


figure, imshow(KK_Lish,[]); title('imfilter lishchuk');

KL_simple= -(CGxx + CGyy) ;

figure, imshow(KL_simple,[]); title('K smple');


[HH KK]=HK(I);

figure
subplot(2,2,1)
imshow(KK,[]); xlabel(' K Mean Curvature of surface') 
subplot(2,2,2);
imshow(HH,[]); xlabel(' H Gaussian Curvature of surface') 


Contact us