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
Small Satellites 