Developed Algorithms

Gaussian Smoothing

Often when images are recorded by a camera, high frequency pixels are seen in the image caused by stuck pixels on the CCD. These pixels can cause problems for feature detectors because of their high intensity value.  In order to filter this noise a 2D Gaussian kernel was convolved over the images generated in this project. Figure 1 displays a noisy image, Figure 2 shows Figure 1 after smoothing.

Corner Detection

To detect features the images a variety of image processing techniques were applied to analyze the color channels of the images[3][4][8]. Corner detection is one type of feature detection which uses the slopes of the image pixel intensity’s in X and Y to find the pixels that define a good corner. Figure 4 shows the corners detected in Figure 3.

 

                 Figure 3. Original Image[11].                                                         Figure 4. Corners Detected[11].            

Thresholding

Thresholding[2][7] involves examining each pixel’s color intensity in an image. If the color intensity is above a certain threshold then that pixel’s color intensity is given a value of 1; likewise if it is below, it is given a value of 0. This technique was used to remove background noise and distortion in the images. Figure 5 shows Figure 4 after thresholding.

        
                    Figure 5. Thresholded Image[11].        

Image Dilation

Dilation is an operation that moves a pixel window over an image calculating the highest pixel intensity in that window and then assigning it to the pixel in the center of that region. The outcome is an image making the pixels with the highest intensities very easy to extract. Figure 6 shows Figure 3 after the dilation operation has been used.

Figure 6. Dilated Image[11].

Feature Extraction

Corner detection outputs an image where the highest intensity points represent the best corners. To extract the coordinates of these corners a technique using thresholding and dilation was applied to the images. Figure 7 shows the corners extracted from Figure 4 represented as green dots.

Figure 7. Extracted Features[11].

Angle Extraction

Figure 8. CCD Design[11]

However because the lens system on digital cameras are not as straight forward as the pinhole design, some distortion occurs.  An experiment was conducted to find the range of the error. The camera which recorded the images was placed on the turn table which recorded a video file as the motor rotated. The rate at which the turn table was rotated was designed so each column in the video at the resolution of 480 x 640 would represent a degree of rotation. These conversion values were computed by doing a calibrated test where the camera was rotated ten degrees and the correlated points could be used to find the per column rotation constant. For the case of this camera the constant was 0.0636 degrees per column. An LED was mounted on the far wall and  used to allow the tracking to be accomplished easily. By analyzing the location of the LED in the frames the lens distortion could be investigated- see figures 9a and 9b. Through looking at the position of the LED it was seen that there was minimal noise caused by lens distortion. This was proven by taking the curve that represented the center of the LED at each frame and fitting a second degree polynomial to the curve. The motor was calibrated so that the frames that the camera recorded were displaced by a column.  Therefore if there was no distortion a line could be used to represent this data and the slope of the line would be 1. The actual curve that was fit was a parabola which could be modeled with f(x) = 0.0000x^2 + 0.9674x^2 + 3.9578 where f(x) is the position of the LED and x is the frame. The coefficient of x^2 which is 0.0000 was so small that it was rounded to zero. The error was then computed from this polynomial which resulted in an range between -3 and 1.5 columns which in degrees means between -0.1908 and 0.0954 degrees. Therefore this shows that it is possible to extract the angles directly from the camera.

Figure 9a. LED Trajectory[11]                                                               Figure 9b. Lens Distortion[11]

Code

 MatLab code routines written specifically to perform image processing operations for this project

function g =  Guass_Pulse(size,width)

e = 2.71828183;

x = zeros(size);

y = zeros(size);

[x,y] = meshgrid(-size/2:size/2);

 g = (1/(2*pi*width))*e.^(-(x.^2+y.^2)/(2*width^2));

function [cim] = Corner(im,g,k)

if nargin == 2

    k = 0.04;

end

% Derivative masks

dx = [-1 0 1; -1 0 1; -1 0 1];  

dy = dx';

% Image derivatives 

Ix = conv2(im, dx, 'same');    

Iy = conv2(im, dy, 'same');

% Smoothed squared image derivatives

Ix2 = conv2(Ix.^2, g, 'same');

Iy2 = conv2(Iy.^2, g, 'same');

Ixy = conv2(Ix.*Iy, g, 'same');

% Compute Corners

cim = (Ix2.*Iy2 - Ixy.^2) - k*(Ix2 + Iy2).^2;

function [r,c, mx] = FeatureExtraction(cim,im,radius,thresh)

if nargin < 3

    radius = 3;

    thresh = 100;

end

% Size of dilation mask.

sze = 2*radius+1;          

% Grey-scale dilate.

mx = ordfilt2(cim,sze^2,ones(sze));

% Find Border of image   

bordermask = zeros(length(im(:,1)),length(im(1,:)));

bordermask(radius+1:end-radius, radius+1:end-radius) = 1;

% Find maxima, threshold, and apply bordermask

cimmx = (cim==mx) & (cim>thresh*10^6) & bordermask;

% Find row,col coords.

[r,c] = find(cimmx);

function FData = Track(R,cor)

D_Thresh = 30;

Sw = 4;

C_Thresh = 1500;

for j = 1:length(cor(1).r)

     feat(j).r(1) = cor(1).r(j);

     feat(j).c(1) = cor(1).c(j);

    for i = 1:length(R)-1

        w = double(R(i).cdata(feat(j).r(i)-Sw/2:feat(j).r(i)+Sw/2,feat(j).c(i)-Sw/2:feat(j).c(i)+Sw/2));

        for h = 1:length(cor(i+1).r)

     S(h) = sum(sum(abs(w-R(i+1).cdata(cor(i+1).r(h)-Sw/2:cor(i+1).r(h)+Sw/2,cor(i+1).c(h)-Sw/2:cor(i+1).c(h)+Sw/2))));

        end

%         %calculate closest feature

         p = find(sqrt((cor(i+1).c-feat(j).c(i)).^2+  (cor(i+1).r-feat(j).r(i)).^2)<D_Thresh)';

         Cv = find(S<C_Thresh);

         c = 0;

         for h = 1:length(p)

             if ~isempty(find(p(h) == Cv))

                 c = c+1;

        Fp(c) = find(p(h) == Cv);

             end

         end

         if ~isempty(Fp)

             if length(feat(j).r) ~= 1

              Sri = polyfit(1:length(feat(j).r),feat(j).r,1);

             Sci = polyfit(1:length(feat(j).c),feat(j).c,1);

             for q = 1:length(Fp)

                  Srf = polyfit(1:length(feat(j).r)+1,[feat(j).r,cor(i+1).r(Cv(Fp(q)))],1);

                  Scf = polyfit(1:length(feat(j).c)+1,[feat(j).c,cor(i+1).c(Cv(Fp(q)))],1);

                  Sd(q) = abs(Srf(1)-Sri(1))+abs(Scf(1)-Sci(1));

             end

              feat(j).r(i+1) = cor(i+1).r(Cv(Fp( find(min(Sd)))));

         feat(j).c(i+1) = cor(i+1).c(Cv(Fp( find(min(Sd)))));

             else

              feat(j).r(i+1) = cor(i+1).r(Cv(Fp(1)));

         feat(j).c(i+1) = cor(i+1).c(Cv(Fp( 1)));

             end

         else

             break

         end

        Fp = [];

    end

end

dat_thresh = 5;

f = 0;

for i = 1:length(feat)

    if length(feat(i).c) == length(R)

        f = f+1;

        Data(f).c = feat(i).c;

        Data(f).r = feat(i).r;

    end

end

h = 0;

for i = 1:length(Data)-1

     for j = 1:length(Data(i).c)-1

        ar(j) = abs(Data(i).r(j)-Data(i).r(j+1));

        ac(j) = abs(Data(i).c(j)-Data(i).c(j+1));

     end

     if isempty(find(ar>dat_thresh+sum(ar)/length(ar))) && isempty(find(ac>dat_thresh+sum(ac)/length(ac)))

    h = h+1

    FData(h).r = Data(i).r;

    FData(h).c = Data(i).c;

    end  

end

for j = 1:length(FData)

    cH(j) = FData(j).c(length(R));

    rH(j) = FData(j).r(length(R));

end

h = 0;

for i = 1:length(FData)

   g = find(cH(i) == cH & rH(i) == rH)

   if length(g) > 1

       h = h+1;

       Ef(h) = g(1);

       cH(g) = [-2*i,-i*100];

   end

end

    FData(Ef) = [];

Crista Device

The Crista device is an inertial measurement unit(IMU) which has 3 accelerometers and 3 gyros internally integrated see figure 10. The data is sampled and read into the computer through the use of the serial port. The sampling process is controlled through C++ drivers. Once the data has been sampled it is stored in text files where it can easily be read into MatLab.

Figure 10. Crista Device[11]

MatLab

MatLab[1] is a technical computing software package that has a wide range of features including image processing and algorithm development. Data in MatLab is stored in n-dimensional matrices or arrays, this allows for images to be handled easily. Also there are many functions designed for image processing which come with MatLab. The majority of the programs designed for this project were written in MatLab.

Maya

Maya[9] is a commercial 3D application that can be used for all purposes of 3D modeling and most applications in 3D production. Maya was used to create and render the environments that were used for the simulation testing.