Home » Featured

# Category Archives: Featured

## Distance Transform Path Planning Algorithm

clc;close all;clear all;

Input Map

```RGB = imread('Mapn.png');
figure('Position',[600 0 600 1000],'color','k');
image(RGB);
hold on;
axis off
axis image
hold off;
% Convert to binary image matrix and inverse matrix values
I = rgb2gray(RGB);
binaryImage = im2bw(RGB, 0.3);
binaryImage = 1-binaryImage;
% Inital pose
xs = 540; ys = 750;
% Goal Pose
xg = 150; yg = 150;```

Distance to Goal

```sz  = size(binaryImage);
for i = 1:sz(1)
for j = 1:sz(2)
rg(i,j) = sqrt((i-yg)^2 +(j - xg)^2);
end
end
% Pose Matrix
pose = zeros(sz(1),sz(2));
pose(ys,xs) = 1;
pose(yg,xg) = 1;```
```Di = 15; % Distance of influence of the obstacles
% Compute for every pixel the distance to the nearest obstacle(non zero
% elements which after inversion corespondent to the obstacles)
[D,IDX] = bwdist(binaryImage);
for i = 1:sz(1)
for j = 1:sz(2)
rd = D(i,j);
if (rd <= Di)
rg(i,j)  = 1e6;
end
end
end```

Distance Transform Path Planning Algorithm

```figure('Position',[600 0 600 1000],'color','k');
hold on;
image(RGB);
colormap gray
contour(rg,15);
spy(pose,'*r');
axis off
axis image
title('Distance Transform Path Planning Algorithm','color','w');
x = xs; y = ys;
last = rg(y-1,x-1);
while (x ~= xg)||(y ~= yg)
dis =[ rg(y-1,x-1), rg(y-1,x),rg(y-1,x+1);
rg(y,x-1), rg(y,x) ,rg(y,x+1);
rg(y+1,x-1), rg(y+1,x),rg(y+1,x+1)];
m = min(dis(:));
[r,c] = find(dis == m);
rg(y,x) = 1E6;
y = y-2+r;
x = x-2+c;
pose(y,x) = 1;
end
spy(pose,'.r');
set(gcf, 'InvertHardCopy', 'off');
hold off;```

## Voronoi Road Map, Path Planing

clc;close all;clear all;

Input Map

```RGB = imread('Map.png');
figure('Position',[100 0 600 700],'color','k');
imagesc(RGB);
axis image
title('Input Map','color','w');
axis off
set(gcf, 'InvertHardCopy', 'off');
hold off;```

Convert to binary image matrix and inverse matrix values

```I = rgb2gray(RGB);
binaryImage = im2bw(RGB, 0.3);
sz  = size(binaryImage);```

Inital pose

```xs = 225; ys = 355;
% Goal Pose
xg = 50; yg = 50;
% Pose Matrix
pose = zeros(sz(1),sz(2));
pose(ys,xs) = 1;
pose(yg,xg) = 1;```

Voronoi Road Map

```figure('Position',[100 0 600 700],'color','k');
hold on;
set(gcf, 'InvertHardCopy', 'off');
sk = 1;
for i = 1:sz(1)
for j = 1:sz(2)
if(binaryImage(i,j) == 1 )
xv(sk) = j;
yv(sk) = i;
sk = sk+1;
end
end
end
[vrx,vry]  = voronoi(xv,yv);
set(gca,'YDir','Reverse');
plot(xv,yv,'y.',vrx,vry,'b-');
axis tight
%axis([1 sz(2) 1 sz(1) ]) ;
axis image
title({'Voronoi Road Map';
'* - Represents Starting and Goal Position '},'color','w');
axis off
spy(pose,'*r');```

Robot Path

```x = xs; y = ys;
indk = 1;
while (abs(x - xg) > 5)||(abs(y - yg)> 5)
path(:,indk) = [x y];
dist = sqrt((vrx - x).^2+(vry - y).^2);
goal = sqrt((xg - x).^2+(xg - y).^2);
% Find closest vertices of the Voronoi edges
[mn ind] = min(dist(:));
xt = vrx(ind);
yt = vry(ind);
goalj = sqrt((xg - xt).^2+(xg - yt).^2);
if (goalj < goal )
x = xt;
y = yt;
vrx(ind) = 1E6;
vry(ind) = 1E6;
indk = indk +1;
else
vrx(ind) = 1E6;
vry(ind) = 1E6;
end
end
path(:,indk) = [xg yg];
plot(path(1,:),path(2,:),'-r');
title({'Voronoi Road Map and Robot Path';
'* - Represents Starting and Goal Position '},'color','w');```

## Potential Field Path Planning

clc;close all;clear all;

Input Map

```RGB = imread('Mapn.png');
figure('Position',[600 0 600 1000],'color','k');
image(RGB);
hold on;
axis off
axis image```

Convert to binary image matrix and inverse matrix values

```I = rgb2gray(RGB);
binaryImage = im2bw(RGB, 0.3);
binaryImage = 1-binaryImage;
% Inital pose
xs = 540; ys = 750;
% Goal Pose
xg = 150; yg = 150;```

Attractive potential

```sz  = size(binaryImage);
for i = 1:sz(1)
for j = 1:sz(2)
rg = sqrt((i-yg)^2 +(j - xg)^2);
U_att(i,j) = 0.5*rg^2;
end
end
contour(U_att,40);

set(gcf, 'InvertHardCopy', 'off');
% Pose Matrix
pose = zeros(sz(1),sz(2));
pose(ys,xs) = 1;
pose(yg,xg) = 1;
spy(pose,'*r');
title({'Attractive Potential Field';...
'Circles Center Corespondents to the Goal Position '},'color','w');
hold off;```

Repulsive potential

```Di = 35; % Distance of influence of the object
% Compute for every pixel the distance to the nearest obstacle(non zero
% elements which after inversion corespondent to the obstacles)
[D,IDX] = bwdist(binaryImage);
for i = 1:sz(1)
for j = 1:sz(2)
rd = D(i,j);
if (rd <= Di)
U_rep(i,j)  = 0.5*(rd - Di)^2;
else
U_rep(i,j)  = 0;
end
end
end
figure('Position',[600 0 600 1000],'color','k');
hold on;
colormap jet;
contourf(U_rep,5);
spy(pose,'*r');
axis off
axis image
set(gcf, 'InvertHardCopy', 'off');
title('Repulsive Potential Field','color','w');
hold off;```

Summary Potential Field

```k_s  = 2; % Scaling factor
U_sum = U_rep/max(U_rep(:))+k_s*U_att/max(U_att(:));
figure('Position',[600 0 600 1000],'color','k');
hold on;
contourf(U_sum,15);
spy(pose,'*r');
axis off
axis image
title('Summary Potential Field','color','w');
set(gcf, 'InvertHardCopy', 'off');
hold off;```

Potential Field Path Planning

```figure('Position',[600 0 600 1000],'color','k');
hold on;
contourf(U_sum,15);
spy(pose,'*r');
axis off
axis image
title('Potential Field Path Planning,*Note the Potential Trap','color','w');
x = xs; y = ys;
last = U_sum(y-1,x-1);
while (x ~= xg)||(y ~= yg)
dis =[ U_sum(y-1,x-1), U_sum(y-1,x),U_sum(y-1,x+1);
U_sum(y,x-1), U_sum(y,x) ,U_sum(y,x+1);
U_sum(y+1,x-1), U_sum(y+1,x),U_sum(y+1,x+1)];
m = min(dis(:));
[r,c] = find(dis == m);
U_sum(y,x) = inf;
y = y-2+r;
x = x-2+c;
pose(y,x) = 1;
end
spy(pose,'.r');
set(gcf, 'InvertHardCopy', 'off');
hold off;```

## Computing a Skeleton of a Map for Mobile Robots,Path Planning

clc;close all;clear all;

Input Map

```RGB = imread('Map1.png');
figure('Position',[0 0 800 800]);
imshow(RGB);```

Convert to binary image matrix and inverse matrix values

```I = rgb2gray(RGB);
binaryImage = im2bw(RGB, 0.3);
binaryImage = 1-binaryImage;
figure('Position',[0 0 800 800]);
imshow(binaryImage);```
```Warning: Image is too big to fit on screen; displaying at
50%```

Compute for every pixel the distance to the nearest obstacle(non zero elements which after inversion corespondent to the obstacles)

```[D,IDX] = bwdist(binaryImage);
figure('Position',[0 0 800 800],'color','k')
colormap bone
image(D);
axis off          % Remove axis ticks and numbers
axis image
set(gcf, 'InvertHardCopy', 'off');```

Taking the Laplacian of the distance transform. Play with the treshold value 0.15 to remove or add more futures in skeleton

```lap = del2(D);
sz  = size(lap);
for i = 1:sz(1)
for j = 1:sz(2)
if( abs(lap(i,j)) < 0.15)
lap(i,j) = 0;
end
end
end
lap = flipdim(lap ,1);           %# vertical flip
figure('Position',[0 0 800 800],'color','k');
colormap winter;
contour(lap,'DisplayName','lap');
axis off
axis image
set(gcf, 'InvertHardCopy', 'off');```

## Mesosphere-Stratosphere-Troposphere(MST) Radar Data Analysis

In this example we  analyse data  from MST (mesosphere-stratosphere-troposphere) radar observations. MST radars used to do observation of the dynamics of the lower and middle atmosphere to study winds, waves, turbulence and instabilities generate irregularities in the atmosphere. The reflected radar signals from the random irregularities are collected by the receiver antenna. The return signal strength is highly depending on the refractive index which is a function of atmospheric parameters such as humidity, temperature, and pressure and electron density. Hence those parameters will highly affect the signal to noise ratio (SNR). Input data file structure {UT, Altitude,Signal amplitude (linear),Signal-to-noise ratio (SNR), dB, Zonal wind, m/s, Meridional wind, m/s}

```clear all; clc; close all;
% Data files name ID
fl = {'150 m height resolutions';
'1200 m height resolutions';
'Barker Coding';
'Complementary Coding';
'Uncoded Data';
};
% Files for three different observation days
fid  = { 'TXT_20061211_test1.fca','TXT_20071010_test1.fca',...
'TXT_20081014_test1.fca';   % 150 m height resolutions
'TXT_20061211_test2.fca','TXT_20071010_test2.fca',...
'TXT_20081014_test2.fca';   % 1200 m height resolutions
'TXT_20061211_test3.fca','TXT_20071010_test3.fca',...
'TXT_20081014_test3.fca';   % Barker coding
'TXT_20061211_test4.fca','TXT_20071010_test4.fca',...
'TXT_20081014_test4.fca';   % Complementary coding
'TXT_20061211_test5.fca','TXT_20071010_test5.fca',...
'TXT_20081014_test5.fca';   % Uncoded data
};
dates = ['2006 Dec 11'; '2007 Oct 10'; '2008 Dec 14' ];```

SNR value as a function of universal time and altitude for three different observation days

```fs = size(fid,1);
fd = size(fid,2);
for jd = 1:fd
for id = 1:fs
clear SNR;
fname = fid{id,jd};
data    = load(fname);
time    = unique(data(:,1));
alt     = unique(data(:,2));
size_t  = size(time,1);
size_a  = size(alt,1);
% Signal-to-noise ratio (SNR), dB
for i = 1:size_a
for j = 1:size_t
SNR(i,j)= data(i+((j-1)*size_a),4);
end
end
% Plot
if(id > 2 )
FigHandle = figure(2);
hold on;
ii = id - 2;
subplot(3,3,ii+(jd-1)*3);
colormap(jet);
pcolor(time,alt,SNR);
shading flat;
caxis([-25,30]);
if(id == 5)
colorbar;
end
if(id == 3)
ylabel('Altitude [km]');
end
if((id == 4))
xlabel(['UT/Date: ',dates(jd,:) ]);
end
if((jd == 1)&&(id == 4) )
title(['Signal-to-Noise Ratio(SNR) [dB]', fl(id)]);
else
if(jd == 1)
title(fl(id));
end
end
set(FigHandle, 'Position', [100, 0, 800, 800]);
else
FigHandle = figure(1);
hold on;
subplot(3,2,id+(jd-1)*2);
colormap(jet);
pcolor(time,alt,SNR);
shading flat;
caxis([-25,30]);
if(id == 2)
colorbar;
end
if(id == 1)
ylabel('Altitude [km]');
end
xlabel(['UT/Date: ',dates(jd,:) ]);
if((jd == 1))
title(['Signal-to-Noise Ratio(SNR) [dB]', fl(id)]);
end
set(FigHandle, 'Position', [100, 0, 600, 800]);
end
end
end ```

## Magnitude and direction of horizontal winds

The figures below provide the horizontal wind variation over the altitude and universal time for three different days. The direction of the wind is mostly from west to east. The results are expected as the wind in this layer of atmosphere moves west to east because of the Coriolis acceleration due to force caused by the rotation of the earth.

```clear all; clc; close all;
dates = ['2006 Dec 11'; '2007 Oct 10'; '2008 Dec 14' ];
% Data files name ID
fid = {'TXT_20061211_test4.fca','TXT_20071010_test4.fca',...
'TXT_20081014_test4.fca' };     % Complementary coding
% Radar with Complementary coding
for id = 1:3
data   = load(fid{id});
time   = unique(data(:,1));
alt    = unique(data(:,2));
size_t = size(time,1);
size_a = size(alt,1);
% Signal-to-noise ratio (SNR), dB
for i = 1:size_a
for j = 1:size_t
zWind(i,j)= data(i+((j-1)*size_a),5); % Zonal Wind, East
mWind(i,j)= data(i+((j-1)*size_a),6); % Meridional Wind, North
end
end
hWind = sqrt(zWind.^2 + zWind.^2);       % Horizontal Wind [m/s]
%dWind=  atan2(zWind,mWind)*180/pi + 180; % Wind Direction angle, Zero in North direction
a_max = 60;
a_min = 30;
FigHandle = figure(2 + id);
set(FigHandle, 'Position', [100, 0, 800, 400]);
subplot(2,2,[1 3]);
colormap(cool);
pcolor(time,alt(1:a_max),hWind(1:a_max,:));
shading flat;
colorbar;
xlabel(['UT/Date:',dates(id,:)]);
ylabel('Altitude [km]');
title('Horizontal wind speed [m/s], Complementary coded signal');
subplot(2,2,[2 4]);
hold on;
whitebg([0.0 .0 .2]);
quiver(time,alt(1:a_max),zWind(1:a_max,:),mWind(1:a_max,:),'r');
%contour(time,alt(1:a_max),hWind(1:a_max,:));
xlabel(['UT/Date:',dates(id,:)]);
title('Horizontal wind direction(top - North,right - East)');
axis([min(time),max(time),min(alt(1:a_max)),max(alt(1:a_max))]);
hold off;
end```