images2ts.m

彩色图像纹理提取的一种算法

                findpixel = true;
                break;
            end;
        end;    
    end;
    if(findpixel == false)
        disp('could not find a valid pixel in the near area, please check it!');
        return;
    end;
end;
if do_plot == 1
    subplot(2,1,1);
    plot(j,i, 'b+');
end;

prior_i = -1;
prior_j = -2;
start_i = i;
start_j = j;
a_order = 1;
back_trace = 1;

while a_order <= n(1)
  
    %%% look careful for this backtrace, need to modify wisely
    if (prior_i == i) & (prior_j == j) % need to expand the neighboring scanning.
        step_len = step_len + 1; %  neighbour distance
        if step_len > 3
            if (a_order - back_trace) < 1
                break;
            end;
            i = imgmatrix_x + 1 - seq_coordinate(2,a_order-back_trace);
            j = seq_coordinate(1, a_order-back_trace);
            back_trace = back_trace + 1;
            prior_i = i;
            prior_j = j;
            step_len =  1;
            %break;
        end;
    else 
        step_len = 1;
        back_trace = 1;
        if a_order > n(1)/2 && abs(start_i - i) < 3 && abs(j - start_j) < 3
            break;% get the round trip!
        end;
        seq_coordinate(1,a_order) = j;   % keep such info into matrix seq_coordinate
        seq_coordinate(2,a_order) = imgmatrix_x + 1 - i;
        if prior_i ~= -1 & prior_j ~= -2 % not the first time to enter this loop
            if j > prior_j | (i < prior_i & j == prior_j) % go right direction, need to change the direction of search
                N1 = N;
            elseif j < prior_j | (i > prior_i & j == prior_j)
                N1 = N11;
            end;
        end;
        prior_i = i;
        prior_j = j;
        if a_order~=1
            Checking_matrix(i,j) = 1;% specify this pixel has been visited
        end;
        a_order = a_order + 1;
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Draw the sketch based on the scanning
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% process, attention : matrix x y & img
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% coordinate x y
        if do_plot == 1      
            subplot(2,1,1);
             plot(j,i, 'r');
       end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
    end 
% scan clockwise to find the next leaf pixel
    %disp(step_len);
    
    switch (step_len)
    case 1, total_nb_cand = size(N1,2);
            for nb_cand = 1 : total_nb_cand
                nb_x = min(i + N1(1, nb_cand),imgmatrix_x);
                nb_y = min(j + N1(2, nb_cand),imgmatrix_y);
                if(nb_x <= 0)
                    nb_x = 1;
                end
                if(nb_y <= 0)
                    nb_y = 1;
                end
                if (Leaf_edge(nb_x, nb_y)~=0) & (Checking_matrix(nb_x, nb_y) == 0) % leaf and not been selected out before
                    i = nb_x;
                    j = nb_y;
                    break;% jump out the neighbour checking if this is the neighbouring leaf pixel
                end;    
            end; 
    case 2, total_nb_cand = size(N2,2);
            for nb_cand = 1 : total_nb_cand
                nb_x = min(i + N2(1, nb_cand),imgmatrix_x);
                nb_y = min(j + N2(2, nb_cand),imgmatrix_y);
                if(nb_x <= 0)
                    nb_x = 1;
                end
                if(nb_y <= 0)
                    nb_y = 1;
                end
                if (Leaf_edge(nb_x, nb_y)~=0) & (Checking_matrix(nb_x, nb_y) == 0)% leaf and not been selected out before
                    i = nb_x;
                    j = nb_y;
                    break;% jump out the neighbour checking if this is the neighbouring leaf pixel
                end;    
            end; 
    case 3, total_nb_cand = size(N3,2);
            for nb_cand = 1 : total_nb_cand
                nb_x = min(i + N3(1, nb_cand),imgmatrix_x);
                nb_y = min(j + N3(2, nb_cand),imgmatrix_y);
                if(nb_x <= 0)
                    nb_x = 1;
                end
                if(nb_y <= 0)
                    nb_y = 1;
                end
                if (Leaf_edge(nb_x, nb_y)~=0) & (Checking_matrix(nb_x, nb_y) == 0)% leaf and not been selected out before
                    i = nb_x;
                    j = nb_y;
                    break;% jump out the neighbour checking if this is the neighbouring leaf pixel
                end;    
            end; 
    end; % end switch
end% while

%output the difference between the original edge image with the one we scaned.

S = nonzeros(seq_coordinate);
seqcoor = S';
seq_num = size(seqcoor,2)/2;
TimeSeriesSeq = zeros(1,seq_num);

x_coor = seqcoor(1:2:end-1);
y_coor = seqcoor(2:2:end);
center_j = round(mean(x_coor));
center_i = imgmatrix_x + 1 - round(mean(y_coor));

if do_plot == 1
    subplot(2,1,1);
    plot(center_j, center_i, 'w+');
end;


if method == 1
  shift = 60;
  for Angle_index = 1 : seq_num
     LeftIndex = Angle_index - shift;
     if (LeftIndex < 1)
         LeftIndex = LeftIndex + seq_num;
     end;
    RightIndex = Angle_index + shift;
     if (RightIndex > seq_num)
         RightIndex = RightIndex - seq_num;
     end;
     Left = [x_coor(1,LeftIndex); y_coor(1,LeftIndex)];
     Middle = [x_coor(1,Angle_index); y_coor(1,Angle_index)];
     Right = [x_coor(1,RightIndex); y_coor(1,RightIndex)];
     TimeSeriesSeq(Angle_index) = angle_cal(Left, Middle, Right) - pi;% calculate the angle given 3 points
 end;
elseif method == 0
  for num_pt = 1 : seq_num
    temp_j = x_coor(1,num_pt);
    temp_i = imgmatrix_x + 1 - y_coor(1,num_pt);
    TimeSeriesSeq(num_pt) = sqrt((center_j - temp_j)*(center_j - temp_j) + (center_i - temp_i) * (center_i - temp_i));
  end;
end;

%%%%% plot the outline and its timeseries into one axes and adding lines to connect the angles up.
if do_plot == 1
    
    [x,y] = size(TimeSeriesSeq);
    TSProcess = TimeSeriesSeq;
    TSProcess(1,:) = (TSProcess(1,:) - mean(TSProcess(1,:))) / std(TSProcess(1,:));
    
    for xindex = 1:y
         xaxes1(xindex) = xindex;
     end;
    xaxes = (xaxes1 - mean(xaxes1))/std(xaxes1);

    
    %figure('Name','Object & TimeSeires'),
    subplot(2,1,2);
    
    
    plot(xaxes, TSProcess);
    hold on,
    imgx_shrink = 10;
    imgy_shrink = 2;
     
     x_coordinate = (seq_coordinate(1,1:seq_num) - mean(seq_coordinate(1,1:seq_num)))...
         /std(seq_coordinate(1,1:seq_num))/imgx_shrink ;
     y_coordinate = (seq_coordinate(2,1:seq_num) - mean(seq_coordinate(2,1:seq_num)))...
         /std(seq_coordinate(2,1:seq_num))/imgy_shrink + 5;
     for draw=1:seq_num  
         plot(x_coordinate(draw), y_coordinate(draw));
     end;
     h = gca;
     set(h, 'ylim', [-4, 8], 'xlim', [-2, 2]);
     
     %%% find out the corresponding points and draw lines between them in one axes
    
     width = 50; %window size
    % 
     pos_threshold = mean(TSProcess) + (max(TSProcess) - mean(TSProcess))/3;
     neg_threshold = mean(TSProcess) - (mean(TSProcess) - min(TSProcess))/3;
    % 
     TS = [0;0] ; % format the TimeSeries coordinate array
     LF = [0;0] ; % Leaf coordinate array
     count = 0;
     index = 2;
     while index < seq_num - width -1
         for run_index = 1:2
             if run_index == 1
                 [c, i] = max(TSProcess(index:index+width));
             else 
                 [c, i] = min(TSProcess(index:index+width));
             end;
             i = i+index-1;
             if (c > pos_threshold && TSProcess(i) > TSProcess(i-1)...
                 && TSProcess(i) > TSProcess(i+1))|| (c < neg_threshold &&...
                 TSProcess(i) < TSProcess(i-1)...
                 && TSProcess(i) < TSProcess(i+1))
                 count = count+1;
                 lf_x = seq_coordinate(1,i);
                 lf_y = seq_coordinate(2,i);
                 lf_x1 = (lf_x - mean(seq_coordinate(1,1:seq_num)))...
                         /std(seq_coordinate(1,1:seq_num))/imgx_shrink ;
                 lf_y1 = (lf_y - mean(seq_coordinate(2,1:seq_num)))...
                         /std(seq_coordinate(2,1:seq_num))/imgy_shrink + 5;
                 LF(:,count) = [lf_x1;lf_y1];
            
                 i1 = (i - mean(xaxes1))/std(xaxes1);
                 c1 = (c - mean(TSProcess))/std(TSProcess);
                 TS(:,count) = [i1;c1];
            
                 hold on,
                 line('XData', [lf_x1 i1], 'YData', [lf_y1 c]);
             end;
         end;
         index = index + width;
     end;
end;


% =========================================
function N = resampleAndNorm(T)

[row col] = size(T);
shortest = inf;
longest = 0;
for i = 1 : row
    if T(i,col) < shortest
        shortest = T(i,col);
    end
    if T(i,col) > longest
        longest = T(i,col);
    end
end
    
% resample to average length
length = round((shortest + longest) / 2);
N = zeros(row,length);
for i = 1 : row
    oldlen = T(i, end);
    seq = T(i,1:oldlen);
    N(i,:) = resample(seq, length, oldlen);
    N(i,:) = (N(i,:) - mean(N(i,:)))/std(N(i,:));   % normalize the data
end