Matlab code to a heart and display happy valentine day.

close all; % Clear all windows
clear all; % Clear all variables
clc; % Clear console

h=figure('color','k');axis off;
set(h,'menubar','none','toolbar','none');
hold on;
syms x y
s=uicontrol('style','pushbutton','units','normal','position',[.05 .05 .12 .05],'string','Close','callback','close all');
i=0.25;x=0;
a=[-1.5-(x*i) 1.5+(x*i) -1.5-(x*i) 1.5+(x*i)];
text(a(1)+0.01,a(3)+0.05,'Happy Valentine''s Day','color','b','fontsize',2*get(0,'defaulttextfontsize'));
a=randn(1,75);c=randn(1,75);
h1=plot(a,c);
set(h1,'linestyle','none','marker','*','markersize',10,'color','r');axis tight;axis off;axis equal;
iter=0;
while(iter<10)
x=3;i=0.25;
while((x*i)<=2.00 && (x*i)>0)
f1=inline('x.^6+y.^6+3*(x.^4)*(y.^2)+3*(x.^2)*(y.^4)-(x.^2)*(y.^3)-3*(x.^4)-6*(x.^2)*(y.^2)-3*(y.^4)+3*(x.^2)+3*(y.^2)-1');
f2=vectorize(f1);
h=ezplot(f2);axis([-1.5-(x*i) 1.5+(x*i) -1.5-(x*i) 1.5+(x*i)]);
set(h,'LineWidth',6);
set(h,'LineStyle','-.');
x=x-1;
drawnow;
end
pause(0.3);
x=3;i=0.25;
while((x*i)<=0.90)
f1=inline('x.^6+y.^6+3*(x.^4)*(y.^2)+3*(x.^2)*(y.^4)-(x.^2)*(y.^3)-3*(x.^4)-6*(x.^2)*(y.^2)-3*(y.^4)+3*(x.^2)+3*(y.^2)-1');
f2=vectorize(f1);
h=ezplot(f2);axis([-1.5-(x*i) 1.5+(x*i) -1.5-(x*i) 1.5+(x*i)]);
set(h,'LineWidth',6);
set(h,'LineStyle','-.');
x=x+1;
drawnow;
end
iter=iter+1;
end



x=-2:0.001:2;
f=sqrt(1-((abs(x)-1).^2));
g=acos(1-abs(x))-pi;
figure
plot(x,f,x,g)
title ('drawing a heart ..have fun');

Published with MATLAB® R2015a

sensitivity and specificity test using matlab

Matlab code example to generate some data a hypothetical scenario

Applications: sensitivity (true positive rate) and specificity (true negative rate) analysis In medical diagnosis and equipment test, sensitivity is a test of correctly identify those with the disease. specificity is a test of correctly identify those without the disease.

Contents

• close
• perform linear classification
• Apply classifier to the data

close

close all; % Clear all windows
clear; % Clear all variables
clc; % Clear console
clf;

subplot(211);
% this is the training data set
% group 1 (centered at 0.2, 0.2)
M1 = 20;
c1 = [0.2,0.2];
x1 = randn(M1,1)+c1(1);
y1 = randn(M1,1)+c1(2);

h1 = line(x1,y1,'linestyle','none','marker','o','DisplayName','normal');
line(c1(1),c1(2),'marker','x','color','b','markersize',12);


% group 2 (centered at 0.4,0.3)
M2 = 20;
c2 = [2, 3];
x2 = randn(M2,1)+c2(1);
y2 = randn(M2,1)+c2(2);

h2 = line(x2,y2,'linestyle','none','marker','^','color','r','DisplayName','atrial tach');
line(c2(1),c2(2),'marker','x','color','r','markersize',12);

xlim([-3,5]);
ylim([-3,7]);

legend([h1, h2]);
xlabel('\propto Heart rate');
ylabel('\propto P-R interval');
title('Hypothetical classification data');

perform linear classification

% concatenate the data
X = [x1, y1, ones(M1,1); x2, y2, ones(M2,1)];

% ground truth (desired classifier output)
d = [repmat(0,[M1,1]); repmat(1,[M2,1])];

% scatter plot
subplot(212);
plot(X(:,1),X(:,2),'o');
xlim([-3,5]);
ylim([-3,7]);

% solve for classifier weights
w = inv(X.' * X) * X.' * d;

% plot decision line
xx = linspace(min(xlim()), max(xlim()), 64);
yy = -w(1)/w(2)*xx - (w(3)-0.5)/w(2);
line(xx,yy,'color','k');

% note: decision boundary depends on the realization of the training data!

Apply classifier to the data

d_out = X*w>0.5;

figure;
cla
line(X((d_out==0),1), X((d_out==0),2),'marker','o','color','b','linestyle','none');
line(X((d_out==1),1), X((d_out==1),2),'marker','^','color','r','linestyle','none');

line(xx,yy,'color','k');

xlim([-3,5]);
ylim([-3,7]);
title('classifier output');

Published with MATLAB® R2015a

Calculate magnetic field of solenoid

Contents

• Biot-Savart integration on a generic curve
• Domain discretization
• Induction curve
• Declaration of variables
• Numerical integration of Biot-Savart law
• Graphic

Biot-Savart integration on a generic curve

close all windows:

clc;
clear all;
close all;

Domain discretization

ND = 7;
Dom = [-1.1 1;
       -1.1 1;
        0.1 6];
% Induction constant
gamma = 1;

% Integration step size
ds = 0.1;

Induction curve

        theta = linspace(0, 15*pi, 70);
        L = [cos(theta') sin(theta') theta'/10];

Nl = numel(L)/3; % Number of points of the curve

Declaration of variables

Induction vector components B = (U, V, W);

U = zeros(ND, ND, ND);
V = zeros(ND, ND, ND);
W = zeros(ND, ND, ND);

% Volume Mesh
[X, Y, Z] = meshgrid(linspace(Dom(1,1), Dom(1,2), ND), ...
                     linspace(Dom(2,1), Dom(2,2), ND), ...
                     linspace(Dom(3,1), Dom(3,2), ND));

Numerical integration of Biot-Savart law

Wait = waitbar(0, 'Integrating, please wait...');
for i = 1:ND
    for j = 1:ND
        for k = 1:ND
            waitbar(sub2ind([ND ND ND],k,j,i)/ND/ND/ND, Wait)
            % Ptest is the point of the field where we calculate induction
            pTest = [X(i,j,k) Y(i,j,k) Z(i,j,k)];
            % The curve is discretized in Nl points, we iterate on the Nl-1
            % segments. Each segment is discretized with a "ds" length step
            % to evaluate a "dB" increment of the induction "B".
            for pCurv = 1:Nl-1
                % Length of the curve element
                len = norm(L(pCurv,:) - L(pCurv+1,:));
                % Number of points for the curve-element discretization
                Npi = ceil(len/ds);
                if Npi < 3
                    close(Wait);
                    error('Integration step is too big!!')
                end
                % Curve-element discretization
                Lx = linspace(L(pCurv,1), L(pCurv+1,1), Npi);
                Ly = linspace(L(pCurv,2), L(pCurv+1,2), Npi);
                Lz = linspace(L(pCurv,3), L(pCurv+1,3), Npi);
                % Integration
                for s = 1:Npi-1
                    % Vector connecting the infinitesimal curve-element
                    % point and field point "pTest"
                    Rx = Lx(s) - pTest(1);
                    Ry = Ly(s) - pTest(2);
                    Rz = Lz(s) - pTest(3);
                    % Infinitesimal curve-element components
                    dLx = Lx(s+1) - Lx(s);
                    dLy = Ly(s+1) - Ly(s);
                    dLz = Lz(s+1) - Lz(s);
                    % Modules
                    dL = sqrt(dLx^2 + dLy^2 + dLz^2);
                    R = sqrt(Rx^2 + Ry^2 + Rz^2);
                    % Biot-Savart
                    dU = gamma/4/pi*(dLy*Rz - dLz*Ry)/R/R/R;
                    dV = gamma/4/pi*(dLz*Rx - dLx*Rz)/R/R/R;
                    dW = gamma/4/pi*(dLx*Ry - dLy*Rx)/R/R/R;
                    % Add increment to the main field
                    U(i,j,k) = U(i,j,k) + dU;
                    V(i,j,k) = V(i,j,k) + dV;
                    W(i,j,k) = W(i,j,k) + dW;
                end
            end
        end
    end
end
close(Wait);

Graphic

figure(1)
M=sqrt(U.^2+V.^2+W.^2);
subplot 121

quiver3(X,Y,Z,U./M,V./M,W./M), hold on, axis equal, grid off
plot3(L(:,1),L(:,2),L(:,3),'r-o','linewidth',3)
title('Normalized field')
view(3)

subplot 122

quiver3(X,Y,Z,U,V,W), hold on, axis equal, grid off
plot3(L(:,1),L(:,2),L(:,3),'r-o','linewidth',3)
title('Magnitude field')
view(3)

Published with MATLAB® R2015a

Calculate the magnetic field of square and circle loop

Contents

• calculate R-vector from the loop(X-Y plane)to Y-Z plane.

%Calculation of Square loop Magnetic Field ..

% close all windows:
clc;
clear all;
close all;

%Square loop is in the X-Y plane
% every point in the Y-Z plane(X=0)
% Magnetic Field is Evaluated

Nz=51;  % No. of grids in Z-axis
Ny=51;  % No. of grids in Y-axis
S=12;   % Length of the Square Loop
N=(S+1)*4;   % No of grids in the loop ( X-Y plane)
N4=N/4;
a=S/2;
I=3;    % current in the loop
u0=1;   % for simplicity, u0 is taken as 1 (permitivity)


% input for x plane
Xc(1:N4)=-a:1:a;
Xc(N4+1:2*N4)=a;
Xc(2*N4+1:3*N4)=a:-1:-a;
Xc(3*N4+1:N)=-a;

% input for y plane
Yc(1:13)=-a;
Yc(14:26)=-a:1:a;
Yc(27:39)=a;
Yc(40:52)=a:-1:-a;

% input for z plane

yp(1:51)=-25:1:25; % Y-coordinates of the plane where we are interested
zp(1:51)=0:1:50;% Z-coordinates of the plane where we are interested

% input for y-z coordinates
Y(1:Ny,1:Nz)=0; % This array is for 1-d to 2-d conversion of coordinates
Z(1:Ny,1:Nz)=0;

for i=1:Ny
    Y(i,:)=yp(i); % all y-coordinates value in 2-d form
end
for i=1:Nz
    Z(:,i)=zp(i);% all z-coordinates value in 2-d form
end

calculate R-vector from the loop(X-Y plane)to Y-Z plane.

%1. find the magnetic field and also the dl-vector along the current loop
%  R is the position vector pointing from loop (X-Y plane) to the
% magnetic field in Y-Z plane
% dl is the current element vector which will make up the square loop

for A=1:Ny  % for all points in the Y-Z plane, find the magnetic field
for B=1:Nz  % from the loop situated in the X-Y plane

for i=1:N-1
Rx(i)=-0.5*(Xc(i)+Xc(i+1));
Ry(i)=(Y(A,B)-(0.5*(Yc(i)+Yc(i+1))));
Rz(i)=Z(A,B);
dlx(i)=Xc(i+1)-Xc(i);
dly(i)=Yc(i+1)-Yc(i);
end
Rx(N)=-0.5*(Xc(N)+Xc(1));
Ry(N)=(Y(A,B)-(0.5*(Yc(N)+Yc(1))));
Rz(N)=Z(A,B);
dlx(N)=-Xc(N)+Xc(1);
dly(N)=-Yc(N)+Yc(1);


% calculate dl cross R is the curl of vector dl and R
for i=1:N
Xcross(i)=dly(i).*Rz(i);
Ycross(i)=-dlx(i).*Rz(i);
Zcross(i)=(dlx(i).*Ry(i))-(dly(i).*Rx(i));
R(i)=sqrt(Rx(i).^2+Ry(i).^2+Rz(i).^2);
end

%The biot savarts law equation.
Bx1=(I*u0./(4*pi*(R.^3))).*Xcross;
By1=(I*u0./(4*pi*(R.^3))).*Ycross;
Bz1=(I*u0./(4*pi*(R.^3))).*Zcross;

% Initialize sum magnetic field to be zero first
BX(A,B)=0;
BY(A,B)=0;
BZ(A,B)=0;

% create a loop by adding all magnetic field from different current elements

for i=1:N   % loop over all current elements along loop
    BX(A,B)=BX(A,B)+Bx1(i); %BX is a null field
    BY(A,B)=BY(A,B)+By1(i);
    BZ(A,B)=BZ(A,B)+Bz1(i);
end


end
end

figure(1);
plot(Xc,Yc,'linewidth',3);
axis([-20 20 -20 20]);
xlabel('X-axis','fontsize',14);
ylabel('Y-axis','fontsize',14);
title('square loop co-ordinates','fontsize',14);
h=gca;
get(h,'FontSize');
set(h,'FontSize',14);
h = get(gca, 'ylabel');
fh = figure(1);
set(fh, 'color', 'white');
grid on;

figure(2);
lim1=min(min(BZ));
lim2=max(max(BZ));
steps=(lim2-lim1)/100;
contour(zp,yp,BZ,lim1:steps:lim2);
axis([1 50 -25 25]);
xlabel('Z-axis','fontsize',14);
ylabel('Y-axis','fontsize',14);
title('BZ component','fontsize',14);
colorbar('location','eastoutside','fontsize',14);
h=gca;
get(h,'FontSize');
set(h,'FontSize',14);
h = get(gca, 'ylabel');
fh = figure(2);
set(fh, 'color', 'white');


figure(3);
lim1=min(min(BY));
lim2=max(max(BY));
steps=(lim2-lim1)/100;
contour(zp,yp,BY,lim1:steps:lim2)
axis([1 50 -25 25])
xlabel('Z-axis','fontsize',14)
ylabel('Y-axis','fontsize',14)
title('BY component','fontsize',14)
colorbar('location','eastoutside','fontsize',14);
h=gca;
get(h,'FontSize')
set(h,'FontSize',14)
h = get(gca, 'ylabel');
fh = figure(3);
set(fh, 'color', 'white');


figure(4);
quiver(zp,yp,BZ,BY,2);
axis([1 50 -25 25]);
xlabel('Z-axis','fontsize',14)
ylabel('Y-axis','fontsize',14)
title('B-field Vector flow','fontsize',14)
h=gca;
get(h,'FontSize')
set(h,'FontSize',14)
h = get(gca, 'ylabel');
fh = figure(4);
set(fh, 'color', 'white');


%%Simulation of Cicular Coil Magnetic Field

%Coil is in the X-Y plane and every point in the Y-Z plane(X=0)
%calculate Magnetic Field

Nz=51;  % No. of grids in Z-axis
Ny=51;  % No. of grids in Y-axis
N=25;   % No of grids in the coil ( X-Y plane)
Ra=6;    % Radius of the coil in the X-Y plane
I=3;    % current in the coil
u0=1;   % for simplicity, u0 is taken as 1 (permitivity)
phi=-pi/2:2*pi/(N-1):3*pi/2; % For describing a circle (coil)

Xc=Ra*cos(phi); % X-coordinates of the coil
Yc=Ra*sin(phi); % Y-coordinates of the coil

yp(1:51)=-25:1:25; % Y-coordinates of the plane
zp(1:51)=0:1:50;% Z-coordinates of the plane

Y(1:Ny,1:Nz)=0; % This array is for 1-d to 2-d conversion of coordinates
Z(1:Ny,1:Nz)=0;

for i=1:Ny
    Y(i,:)=yp(i); % all y-coordinates value in 2-d form
end
for i=1:Nz
    Z(:,i)=zp(i);% all z-coordinates value in 2-d form
end

% calculate R-vector from the coil(X-Y plane)to Y-Z plane
%find the magnetic field and also the dl-vector along the coil current
% R is the position vector pointing from  (X-Y plane)coil to
%  the magnetic field in Y-Z plane


for a=1:Ny  %in for loop is along Y
for b=1:Nz  %is along Z-axis-

for i=1:N-1
Rx(i)=-0.5*(Xc(i)+Xc(i+1));
Ry(i)=(Y(a,b)-(0.5*(Yc(i)+Yc(i+1))));
Rz(i)=Z(a,b);
dlx(i)=Xc(i+1)-Xc(i);
dly(i)=Yc(i+1)-Yc(i);
end
Rx(N)=-0.5*(Xc(N)+Xc(1));
Ry(N)=(Y(a,b)-(0.5*(Yc(N)+Yc(1))));
Rz(N)=Z(a,b);
dlx(N)=-Xc(N)+Xc(1);
dly(N)=-Yc(N)+Yc(1);

%dl is the current element vector in the coil
%dl cross R is the curl of vector dl and R
%calculate dl cross R

for i=1:N
Xcross(i)=dly(i).*Rz(i);
Ycross(i)=-dlx(i).*Rz(i);
Zcross(i)=(dlx(i).*Ry(i))-(dly(i).*Rx(i));
R(i)=sqrt(Rx(i).^2+Ry(i).^2+Rz(i).^2);
end


% the biot savarts law equation

Bx1=(I*u0./(4*pi*(R.^3))).*Xcross;
By1=(I*u0./(4*pi*(R.^3))).*Ycross;
Bz1=(I*u0./(4*pi*(R.^3))).*Zcross;

% magnetic field from all current
BX(a,b)=0;       % Initialize sum magnetic field to be zero first
BY(a,b)=0;
BZ(a,b)=0;

for i=1:N   % loop over all current elements along coil
    BX(a,b)=BX(a,b)+Bx1(i);
    BY(a,b)=BY(a,b)+By1(i);
    BZ(a,b)=BZ(a,b)+Bz1(i);
end

end
end

figure(5)
plot(Xc,Yc,'linewidth',3)
axis([-20 20 -20 20])
xlabel('X-axis','fontsize',14)
ylabel('Y-axis','fontsize',14)
title('square loop co-ordinates','fontsize',14)
h=gca;
get(h,'FontSize')
set(h,'FontSize',14)
h = get(gca, 'ylabel');
fh = figure(5);
set(fh, 'color', 'white');
grid on

figure(6)
lim1=min(min(BZ));
lim2=max(max(BZ));
steps=(lim2-lim1)/100;
contour(zp,yp,BZ,lim1:steps:lim2)
axis([1 50 -25 25])
xlabel('Z-axis','fontsize',14)
ylabel('Y-axis','fontsize',14)
title('BZ component','fontsize',14)
colorbar('location','eastoutside','fontsize',14);
h=gca;
get(h,'FontSize')
set(h,'FontSize',14)
h = get(gca, 'ylabel');
fh = figure(6);
set(fh, 'color', 'white');


figure(7)
lim1=min(min(BY));
lim2=max(max(BY));
steps=(lim2-lim1)/100;
contour(zp,yp,BY,lim1:steps:lim2)
axis([1 50 -25 25])
xlabel('Z-axis','fontsize',14)
ylabel('Y-axis','fontsize',14)
title('BY component','fontsize',14)
colorbar('location','eastoutside','fontsize',14);
h=gca;
get(h,'FontSize')
set(h,'FontSize',14)
h = get(gca, 'ylabel');
fh = figure(7);
set(fh, 'color', 'white');


figure(8)
quiver(zp,yp,BZ,BY,2)
axis([1 50 -25 25])
xlabel('Z-axis','fontsize',14)
ylabel('Y-axis','fontsize',14)
title('B-field Vector flow','fontsize',14)
h=gca;
get(h,'FontSize')
set(h,'FontSize',14)
h = get(gca, 'ylabel');
fh = figure(8);
set(fh, 'color', 'white');

ans =

    10


ans =

    10


ans =

    10


ans =

    10


ans =

    10


ans =

    10

Published with MATLAB® R2015a

calculate the elegtromagnetic of cylindrical wires

Contents

• Calculate additonial parameter
• Plot the Current density

% This code is used to calculate the elegtromagnetic of cylindrical wires.

% close all windows:
clc;
clear all;
close all;

%%PART I - Calculate integral parameter of wire for different frequencies
%In PART I the specific resistance and inductance parameters are calculated
% for a couple of frequencies for a copper wire. As well the skindepth and
% the ratios of different resistancemodels are also calculated. The values
% are written in a cellarray, which is the result [tab] of this function.
% Input data for a copper wire of diameter 3.3 mm
% Gauge: 8 AWG. Average Wire Diameter: 0.1285 in (3.264) mm
%copper wire info: https://en.wikipedia.org/wiki/American_wire_gauge
%Define global variables
global u0
u0=4*pi*1e-7; %permeability of free space constant in [Vs/Am]
ur          = 1;           % relative Permeability of Material   e.g. Copper
%Define coil parameters
rho         = 1.72e-8;     % Resistivity in [Ohm*m],             e.g. Copper
D_wire_m    = 3.3e-3;        % Wire-Diameter 0.002677 m
f           = [1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 13.5e6, 1e9, 2e9, 3e9, 4e9, 10e9, 20e9]';       % frequency vector


Calculate additonial parameter
R_m         = D_wire_m/2;  % Wire-Radius in [m²]
area        = R_m^2*pi;    % Circle_area in [m²]
sigma       = 1/rho;       % Conductivity, [S/m]

R_DC = 1/(sigma*area);     % DC-Resistance load per unit length in [Ohm/m]
omega  =  2*pi*f;                                  % angular frequency in [Hz]
delta  = 1./sqrt(omega.*sigma.*u0.*ur./2);         % skindepth in [m]


% calculate impedance of this wire;
z_math = Z_wire (omega, R_m,  sigma, u0, ur);

% R_DC_Ratio by an tubemodel-calculation (DC-Resistance in depence of skindepth and wire-radius)
R_dc_tube_ratio = tube_model_func(delta ,R_m ,sigma);


% *Collect the data for the table*
tab_header = {'frequency', 'skindepth', 'L lpul', 'R_DC lpul', 'R_DC-ratio ', 'R_DC-tubemodel_ratio'; ...
              'Hz',        'm',          'H/m',   'Ohm/m' ,    '[1]',          '[1]' };
tab2       = num2cell([ f, ...                                              % 1. column frequency
                        delta, ...                                          % 2. column skindepth
                        imag(z_math)./omega,  ...                           % 3. column impedance  (H/m)
                        real(z_math), ...                                   % 4. column resistance (Ohm/m)
                        real(z_math)/R_DC,...                               % 5. column resistance-Ratio [1]
                        R_dc_tube_ratio   ]);                               % 6. column resistance tube_model ratio [1]
tab       = [tab_header; cellfun(@num2str,(tab2), 'UniformOutput',0)];
                                                                             % lpul: load per unit length

%%PART II - Plot currentdensity on wire for 13.5e6 Hz
% In PART II the currentdensity is calculated as a function of radius of the
% wire for a frequency of 13.56 MHz and a current of 1 Amps. The current
% distribution is plotted as 3D-surface plot as well as over the radius.
%
f_circle_Hz = 13.5e6;     % frequency in Hz
current_A   = 1;      % current in wire in Amps
dr          = 200;    % number of steps (dicretization

% create a grid for that circle
[X,Y] = meshgrid(-R_m: R_m/(dr-1) :R_m , -R_m: R_m/(dr-1) :R_m  );

% Calculate radius for each Point:
r = sqrt(X.^2+Y.^2);

% calculation of currentdensity:
A      =  sqrt(-1j*2*pi*f_circle_Hz*sigma*u0*ur);       % Nomalized argument for Bessel function
[J0]   =  besselj ( 0 , A.* r );                        % first bessel-function of zeroth order
[J1]   =  besselj ( 1 , A.*R_m );                       % first bessel-function of first order
J_vec  =  A.*current_A.*1./(2.*pi.*R_m).*J0./J1;

% Recalc only the wire with radius
J_vec((X.^2+Y.^2)>=R_m^2)   = NaN;
X(isnan(J_vec))             = NaN;
Y(isnan(J_vec))             = NaN;

% Plot results


Plot the Current density
figure (1);
surf(X,Y,real(J_vec),'LineStyle','none');
colorbar;
axis square;
view([0 0 90]);
title (['Current density in wire (R=', num2str(R_m),' m, I=',num2str(current_A),' A, f=', num2str(f_circle_Hz), ' Hz)'  ]);

%%Plot the Current density over radius
figure (2);
plot(Y(X==0), J_vec(X==0));
grid on;
xlabel('radius in m');
ylabel('J in A/m²');
title (['Current density distribution over radius in wire (R=', num2str(R_m)])

%%PART III - Calculate the electrical field strength and the voltage drop for the parameters of PART II
% In PART III the electric field strength for the wire of a length of 1000m
% and the voltage drop are calculated. The parameter for resistance and
% inductance are also computed.


length_m    = 1000;  % lenghth of wire [m];

E_Vpm       = current_A .* A./(2 .* pi.* R_m.* sigma)  .* besselj ( 0 , A.* R_m )./ besselj ( 1 , A.*R_m )   % electrical field_strength in Volt per m
deltaU_V    = E_Vpm.*length_m                              % complex voltage drop in Volt
absU_V      = abs(deltaU_V)                                % absolute voltage drop in Volt
R_Ohm       = real(deltaU_V/current_A)                     % Resistance in Ohm
Li_H        = imag(deltaU_V/(2*pi*f_circle_Hz*current_A))  % inductance in Henr



I0=current_A; %Coil current in Amps
a=.15; %Coil radius in m

%Define coordinates of coil center point
x_p=0; y_p=0; z_p=0;

%%calculate the magnetic field at a single point in space
%Input test point
x=0; y=0; z=.1;

[Bx,By,Bz] = magnetic_field_current_loop(x,y,z,x_p,y_p,z_p,a,I0)

%Input vector of points

x=0; y=0; z=linspace(0,.25,100); %These default coordinates calculates the magnetic field along the z axis through the center of the coil

[Bx,By,Bz] = magnetic_field_current_loop(x,y,z,x_p,y_p,z_p,a,I0);

figure(3);
plot(z,Bz);
xlabel('z [m]');
ylabel('Bz [T]');
title('1D magnetic field tests');


%%input mesh of points in 2D plane

x=0; [y,z]=meshgrid(linspace(-.05,.05,25),linspace(0,.1,25)); %this is a 2d plane over the x=0 plane that extends away from the coils in the yz plane.

[Bx,By,Bz] = magnetic_field_current_loop(x,y,z,x_p,y_p,z_p,a,I0);

figure (5);
surf(y,z,Bz);
xlabel('y [m]');
ylabel('z [m]');
zlabel('Bz [T]');
title('2D magnetic field tests');
colorbar; %add colorbar
shading flat; %Removes black lines from the mesh

Warning: Imaginary parts of complex X and/or Y arguments ignored 

E_Vpm =

   0.0929 + 0.0924i


deltaU_V =

  92.8570 +92.3501i


absU_V =

  130.9617


R_Ohm =

   92.8570


Li_H =

   1.0887e-06


Bx =

     0


By =

     0


Bz =

   2.4129e-06


    





































Published with MATLAB® R2015a

matlab code to identify organs of a CT image

Contents

• this matlab code provides some techniques and methods to view and identify organs of a CT image
• plot histogram of the image
• perform the color map of a CT image
• HSV space
• filtering demo on the CT data
• average filter by hand
• high-pass by hand
• average filter
• gaussian filter
• unsharp mask
• laplacian
• laplacian of gaussian
• threshold
• prewitt edge filter
• sobel edge filter
• standard deviation edge filter
• skew filter
• kurtosis filter

this matlab code provides some techniques and methods to view and identify organs of a CT image

The first method is to use histogram to segment the second method is to use color a CT image third method is filter.

% read a CT image
I=imread('./Abdominal wall normal anat (8).png');

plot histogram of the image

subplot(311); % plot to ensure a correct image is loaded.
imhist(I)

subplot(3,1,2:3);
imshow(I);
colorbar

perform the color map of a CT image

colormap(jet);
caxis
caxis([0,255]);

ans =

     0   255

colormap( [repmat([0,0,0],[50,1]); ... % 0-50
          repmat([0,1,0],[50,1]);...  % 51-100
          repmat([1,0,0],[75,1]);...  % 100-175
          repmat([0,0,1],[50,1]);...  % 176-225
          repmat([1,1,1],[56-25,1])]);
colorbar

HSV space

convert to single-precision floating point

I = single(I);
[ny,nx] = size(I);
Ihsv = zeros(ny,nx,3);
Ihsv(:,:,2) = 0.5;
Ihsv(:,:,3) = I/max(I(:));

% add color
th1 = 50;
th2 = 100;
th3 = 170;
Ihsv(:,:,1) = Ihsv(:,:,1) + (I<th1) * 0;
Ihsv(:,:,1) = Ihsv(:,:,1) + (I>th1 & I<th2) * 0.1;
Ihsv(:,:,1) = Ihsv(:,:,1) + (I>th2 & I<th3) * 0.2;
Ihsv(:,:,1) = Ihsv(:,:,1) + (I>th3 & I<225) * 0.5;
Ihsv(:,:,1) = Ihsv(:,:,1) + (I>225) * 0.8;

imshow(hsv2rgb(Ihsv));

filtering demo on the CT data

clf;
imagesc(I);colormap(gray);
title('unfiltered');

average filter by hand

n = 5;
b = ones(n,n)/ (n*n);
I2 = conv2(I,b,'same');
imagesc(I2);

high-pass by hand

b = [-1, -1, -1; -1, 7.9866, -1; -1, -1, -1];
I2 = conv2(I,b,'same');
imagesc(I2);
colorbar
mean(I2(:))
caxis([-200,300])

ans =

   -0.0026

average filter

I=single(I);
b = fspecial('average',[9,9]);
I2 = conv2(I,b,'same');
imagesc(I2);

gaussian filter

note: pick filter size appropriately imagausfilt avoids this problem

I=single(I);
b = fspecial('gaussian',32,2);
I2 = conv2(I,b,'same');
imagesc(I2);

unsharp mask

b = fspecial('unsharp');
subplot(221);
imagesc(b);
colorbar;
subplot(222);
imagesc(I);
subplot(224);
I2 = conv2(I,b,'same');
imagesc(I2);

laplacian

clf
b = fspecial('laplacian');
I2 = conv2(I,b,'same');
imagesc(I2);

laplacian of gaussian

b = fspecial('log');
I2 = conv2(I,b,'same');
imagesc(I2);
colorbar

threshold

imagesc(abs(I2)>60);

prewitt edge filter

b = fspecial('prewitt');
I2 = conv2(I,b.','same');
imagesc(abs(I2));
colorbar
caxis([0,300])

imagesc(abs(I2)>150)

sobel edge filter

b = fspecial('sobel');
I2 = conv2(I,b,'same');
imagesc(I2);
colorbar

imagesc(abs(I2)>200)

standard deviation edge filter

I2 = nlfilter(I, [3,3], 'std2');
imagesc(I2);
colorbar;

imagesc(I2>10);

skew filter

f = @(X) skewness(X(:));
I2 = nlfilter(I, [3,3], f);
imagesc(I2);
colorbar;

kurtosis filter

f = @(X) kurtosis(X(:));
I2 = nlfilter(I, [16,16], f);
imagesc(I2);
colorbar;

Published with MATLAB® R2016a