clear all; 
l = imread('D:\\ear\\earpic\\8-4.bmp');
m = 1;
n = 1;
sizel = size(l);
%l_temp = zeros(sizel(1)/4, sizel(2)/4);
for i =1 :4:sizel(1)    
    for j = 1:4:sizel(2)
        l_temp(m, n) = l(i,j);
        n = n+1;
    end
    m = m+1;
n =1;
end
compressear = uint8(l_temp);
[count,x] = find(compressear<100);
countsize = size(count);
compressear_temp = compressear;
for i = 1:countsize(1)
pixx = x(i);
pixy = count(i);
compressear_temp(pixy,pixx) = 0;
end

compressear_edge = edge(compressear_temp,'sobel');
figure, imshow(compressear_edge);
% BW = imread('blobs.png');
       [B,L,N,A] = bwboundaries(compressear_edge);
       imshow(compressear_edge); hold on;
       for k=1:length(B),
         if(~sum(A(k,:)))
           boundary = B{k};
           plot(boundary(:,2), boundary(:,1), 'r','LineWidth',2);
           for l=find(A(:,k))'
             boundary = B{l};
             plot(boundary(:,2), boundary(:,1), 'g','LineWidth',2);
           end
         end
 
附錄二  特征提取PCA程序

close all;
clear all;
class_num=10;
sample_num = 6;
row=163;
line=123;
rowline=row*line;
load earwave2;
 
disp('ear eigenear now begin');
right_num = [];
rr_0 = [];
for count =2         %count=1,2,3
    tic;
    count
    train_num = count;
    train_total = train_num * class_num;
    test_num = sample_num - train_num;
    test_total = test_num * class_num;
    
    %% -------- 1.calculate Sb Sw -------- %%
    ear_train = zeros(rowline, train_total);
    ear_m_all = zeros(rowline,1);
    ear_m = zeros(rowline,class_num);
    for i = 1:class_num
        for j = 1:train_num
            kk = (i-1)*train_num + j;
           %ear_temp = double(ear(:,:,i,j));
            ear_temp = double(compressear(:,:,i,j));
            for p = 1:row
                for q = 1:line
                      site = (p-1)*line + q;  
                      ear_train(site, kk) = ear_temp(p,q);
                end
            end
            ear_m_all = ear_m_all + ear_train(:,kk);
            ear_m(:,i) = ear_m(:,i) + ear_train(:,kk);
        end
        ear_m(:,i) = ear_m(:,i)/train_num;
    end
    ear_m_all = ear_m_all / train_total;
    
   ear_test = zeros(rowline, test_total);
    for i = 1:class_num
        for j = 1:test_num
            kk = (i-1)*test_num + j;
            %ear_temp = double(ear(:,:,i,j+train_num));
            ear_temp = double(compressear(:,:,i,j+train_num));
            for p = 1:row
                for q = 1:line
                    site = (p-1)*line + q;  
                    ear_test(site, kk) =ear_temp(p,q);
                end
            end
        end
    end        
    
    Sb = zeros(rowline,class_num);
    for i = 1:class_num
        Sb(:,i) = ear_m(:,i) - ear_m_all;
    end
%    Sb0 = (Sb * Sb') / class_num;  %% [rowline, rowline]      
 
    St = zeros(rowline,train_total);
    for i = 1:class_num
        for j = 1:train_num
            kk = (i-1)*train_num + j;
            St(:,kk) =ear_train(:,kk) -ear_m_all;
        end    
    end
%    St0 = (St * St') / train_total;  %% [rowline, rowline]
    
    clear ear_m_all ear_m;
    
    %% -------- 2.calculate transform of St -------- %%
    ss = (St' * St) / train_total;   %% [train_total, train_total]
    rr = rank(ss);
    rr_0 = [rr_0; rr];
    
    [V_St, D_St] = eig(ss);
    dd = abs(eig(ss));
    [dd_value, dd_site] = sort(dd);
    temp2 = rank(ss);
    temp1 = train_total - temp2 + 1;
    yy = [];
    Dt = zeros(rr);
    jj = 0;
 
    dd_demo = zeros(temp2,1);
    dd_count = 1;
    
    for i = train_total:-1:temp1
        yy = [yy; V_St(:,dd_site(i))'];
        jj = jj + 1;
        Dt(jj,jj) = dd_value(i)^(-0.5);   
        
        dd_demo(dd_count) = dd_value(i);
        dd_count = dd_count + 1; 
    end
    yy = yy';
    zz = St * yy * Dt;
    ww = zz';
    
    figure, plot(dd_demo);
    for ii = 1:6
        ww_temp = ww(ii,:);
        ww_temp = ww_temp - min(ww_temp);
        ww_temp = ww_temp / max(ww_temp);
        ww_temp = floor(ww_temp * 255);
        ear_temp = zeros(row, line);
        for p = 1:row
           for q = 1:line
                 site = (p-1)*line + q;  
                 ear_temp(p,q) = ww_temp(site);
           end
           end
        ear_temp = uint8(ear_temp);
        ss = sprintf('eigenear_%d.bmp',ii);
        imwrite(ear_temp, ss,'bmp');
        end    
    disp('ear_eigenear is over');
   
    %% -------- 4.calculate transform of test and train -------- %%
    ear_train_new = ww * ear_train;
    ear_test_new = ww *ear_test;
 
    error = [];
    for i = 1:class_num
        for j = 1:test_num
            %i,j
            kk = (i-1)*test_num + j;
            res_y =ear_test_new(:,kk);
        
            dis = [];
            for ii = 1:class_num
                for jj = 1:train_num
                    res_x =ear_train_new(:,(ii-1)*train_num+jj);
                    res = (res_y - res_x).^2;%最小距離
                    temp = sqrt(sum(res));
                    %res = abs(res_y - res_x);%最近鄰
                    %temp = sum(res);
                    
                    dis = [dis; temp];
                end
            end
            [value, site] = min(dis);
            site_new = floor((site-1)/train_num) + 1;
            
            if site_new ~= i
                error = [error; i*test_num+j];
            end
        end
    end
    right_num = [right_num; (1.0-length(error)/(class_num*test_num))*100 ];
    right_num'
    toc;
end