images2ts.m

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

% Taken from Dr. Keogh's image2bmp.m, modified by Adam Meadows

% function: images2ts
% description: convert a cell array of images into a collection of time
% series
% author: Li Wei & Xiaopeng Xi (modified by Adam Meadows)
% last updated: 06/05/2006
%
% - add an input parameter 'do_plot'
%                           do_plot = 1, draw the extraction process
%                           otherwise,   don't draw
% - add two input parameters 'method'
%                           method =0, centroid extraction
%                           method =1, angle extraction
%                           'option'
%                           option = 1, nomalize and interpolate to the same length
%                           otherwise, don't do it
% - update the profile_extraction function (using the distance to centroid
% now) and fix a backtrace bug caused by noise

function freq = images2ts(imgs, do_plot, option, method)

close all;

index = 1;
seq_grup = zeros(1,1);          % return value of empty folder case
seqcoor_grup = zeros(1,1);      % store the value of profile coordinates

for i = 1 : length(imgs)
    
    old = imgs{i};

    % ##### rescale the images ####
    sizevector = size(old);
    num_of_dim = ndims(old);
    newim = old;
    if (num_of_dim == 3)
        if (sizevector(1) ~= sizevector(2))
            if (sizevector(1) > sizevector(2))                        
                newim(:, sizevector(2)+1:sizevector(1), :) = 0;                        
            else                        
                newim(sizevector(1)+1:sizevector(2), :, :) = 0;
            end
        end
    end
    if (num_of_dim == 2)
        if (sizevector(1) ~= sizevector(2))
            if (sizevector(1) > sizevector(2))                        
                new(:, sizevector(2)+1:sizevector(1)) = 0;                        
            else                        
                newim(sizevector(1)+1:sizevector(2), :) = 0;
            end
        end
    end   
    index = index + 1;

    [seq, seqcoor] = profile_extraction(old, do_plot, method);

    %keep track of the size of Timeseries
    sizes = size(seq);
    seq(sizes(2)+1) = sizes(2); 

    %keep track of the size of coordinate
    sizes2 = size(seqcoor);
    seqcoor(sizes2(2)+1) = sizes2(2); 


    %====== Adds the newly generated time series -- seq into the 
    %existing time series -- seq_grup ===================
    [newseq_num, max_size] = size(seq);
    if newseq_num == 0
        continue;
    end
    oldsize = size(seq_grup);
    if max_size > oldsize(2)    % all the previous stored seq need to be restored
        if oldsize(2) == 1      % seq_grup contains no useful data yet
            seq_grup = seq;     % reallocate the seq_grup
        else
            tmp = zeros(newseq_num+oldsize(1), max_size);   % already has values inside seq_grup
            padnum = max_size - oldsize(2);                 % get the pad array for the new storation
            pad = zeros(1,padnum);
            for i = 1:oldsize(1)
                tmp(i,:) = [seq_grup(i,:),pad];
                tmp(i,oldsize(2)) = [0];                    % clear old size info
                tmp(i,max_size) = seq_grup(i,oldsize(2));   % move the size info backwards
            end
            tmp(oldsize(1)+1:oldsize(1)+newseq_num,:) = seq;
            seq_grup = tmp;
        end                 
    else                    %no need to reallocate the original seq_grup, just copy the new seq into the new rows of it.
        padnum = oldsize(2) - max_size; % the new seq need be expanded by pad
        pad = zeros(1,padnum);
        for padind = 1:padnum
            pad(padind) = 0;
        end;
        tmp = zeros(newseq_num, oldsize(2));
        for i = 1:newseq_num
            tmp(i,:) = [seq(i,:),pad];
            tmp(i,max_size) = [0];%clear old size info
            tmp(i,oldsize(2)) = seq(i,max_size);%move the size info backwards
        end;
        seq_grup(oldsize(1)+1:oldsize(1) + newseq_num,:) = tmp;
    end;

    % ====== Adds the newly generated coordinates -- seqcoor into the 
    % existing time series -- seqcoor_grup ===================
    [newseq_num, max_size] = size(seqcoor);
    if newseq_num == 0
        continue;
    end
    oldsize = size(seqcoor_grup);
    if max_size > oldsize(2)  % all the previous stored seq need to be restored
        if oldsize(2) == 1% seq_grup contains no useful data yet
            seqcoor_grup = seqcoor;%reallocate the seq_grup
        else
            tmp = zeros(newseq_num+oldsize(1), max_size);% already has values inside seq_grup
            padnum = max_size - oldsize(2);    % get the pad array for the new storation
            pad = zeros(1,padnum);
            for i = 1:oldsize(1)
                tmp(i,:) = [seqcoor_grup(i,:),pad];
                tmp(i,oldsize(2)) = [0];%clear old size info
                tmp(i,max_size) = seqcoor_grup(i,oldsize(2));%move the size info backwards
            end
            tmp(oldsize(1)+1:oldsize(1)+newseq_num,:) = seqcoor;
            seqcoor_grup = tmp;
        end;%ending up the case which need to expand the original seq_grup
    else %no need to reallocate the original seq_grup, just copy the new seq into the new rows of it.
        padnum = oldsize(2) - max_size; % the new seq need be expanded by pad
        pad = zeros(1,padnum);
        for padind = 1:padnum
            pad(padind) = 0;
        end
        tmp = zeros(newseq_num, oldsize(2));
        for i = 1:newseq_num
            tmp(i,:) = [seqcoor(i,:),pad];
            tmp(i,max_size) = [0];%clear old size info
            tmp(i,oldsize(2)) = seqcoor(i,max_size);%move the size info backwards
        end
        seqcoor_grup(oldsize(1)+1:oldsize(1) + newseq_num,:) = tmp;
    end              
end

% normalize seq_grup and resample them to be same length
if option == 1
    newseq_grup = resampleAndNorm(seq_grup);
    %save(freqfile, 'newseq_grup', '-ascii');
    freq = newseq_grup;
else
    %save(freqfile, 'seq_grup', '-ascii');    
    freq = seq_grup;
end


% =========== local function ===============
function [TimeSeriesSeq, seqcoor] = profile_extraction(image, do_plot, method)

I = image;

row = size(I,1);
col = size(I,2);

%TBD: for now assume color images
C = 1;
if C == 1
     J = rgb2gray(I);
     se0 = strel('line', 2, 90);
     se1 = strel('line', 2, 0);
     J = imdilate(J, [se0,se1]);
     J = imfill(J,'hole');
     seD = strel('diamond',1);
     J = imerode(J,seD);
     J = im2bw(J,0.2);
     J = edge(J, 'sobel',(graythresh(I) * .1));
    
     k1 = find(J==0);
     k2 = find(J==1);
     if size(k1,1) == 0 | size(k2,1) == 0
        J = rgb2gray(I);
        J = im2bw(J,0.6);
        se0 = strel('line', 1, 90);
        se1 = strel('line', 1, 0);
        J = imdilate(J, [se0,se1]);
        seD = strel('diamond',1);
        J = imerode(J,seD);
        J = medfilt2(J,[6 6]);
        J = edge(J, 'sobel',(graythresh(I) * .1));
    end;
else
    temp = double(I);
    J = edge(temp, 'sobel',(graythresh(temp) * .1));
end;

if do_plot == 1
    figure, subplot(2,1,1);
    imshow(J);
    hold on;
end;

%    Get the coordinate of edge:			
%   Scan according to the order up-down, left-right;
%   Treat the matrix as a sparse matrix; 'find' returns the x/y of non-0 pixels]
%   c = [i, j, v] to display i, j, v in three columns.

Leaf_edge = double(J);
[imgmatrix_x,imgmatrix_y] = size(Leaf_edge);

%allocate the all 0 matrix as same size as the 0/1 image, to later on check if the neighbour are selected out or not
Checking_matrix = zeros(imgmatrix_x,imgmatrix_y);
[i_coordinate,j_coordinate,v] = find(Leaf_edge); 
%plot(j_coordinate(20), i_coordinate(20), 'g+');

%    Start from the first non-0 value clock-wise to find the sequence of angles
%   (adjacent two pixels).eg : Check (i, j)'s neighbour. Order should  be: 
%   (i-1,j), (i-1,j+1), (i,j+1), (i+1,j+1), (i+1,j), (i+1, j-1), (i,j-1), (i-1,j-1). 
%   The neighbor that has checked before should be omitted from the 8 possible neighbors checking. 
%   Only 1/7 is selected out.
c = [i_coordinate, j_coordinate].'; %two rows, each column are the coordinate of the leaf pixel
num = nnz(Leaf_edge);
N =   [1 1 0 -1 -1 -1 0 1; 0 1 1 1 0 -1 -1 -1];
N1 =  [1 1 0 -1 -1 -1 0 1; 0 1 1 1 0 -1 -1 -1];
N11 = [-1 -1 0 1 1 1 0 -1; 0 -1 -1 -1 0 1 1 1];
N2 = [-2 -2 -2 -1  0  1  2  2  2  2  2  1  0  -1  -2 -2; 
       0 1  2  2  2  2  2  1  0 -1 -2 -2 -2  -2  -2 -1];
N3 = [-3 -3 -3 -3 -2 -1  0  1  2  3  3  3  3  3  3  3  2  1  0 -1 -2 -3 -3 -3;  
       0  1  2  3  3  3  3  3  3  3  2  1  0 -1 -2 -3 -3 -3 -3 -3 -3 -3 -2 -1];

%%%%%%   seq = timeseries(i_coordinate, j_coordinate);
% convert the coordinate of the edge to the sequence of edge shape.
%%%%%%   function seq = timeseries(x, y)
n = size(i_coordinate);
seq_coordinate = zeros(2,n(1));  % Preallocate matrix for clockwise ordered leaf pixels


%find non-zero pixel in the 8 neighbours
i = i_coordinate(1);
j = j_coordinate(1);
%plot(j,i,'b+');

findpixel = false;
while Leaf_edge(i,j) == 0
    total_nb_cand = size(N1,2);
    for nb_cand = 1 : total_nb_cand
        i1 = i + N1(1, nb_cand);
        j1 = j + N1(2, nb_cand);
        if (Leaf_edge(i1,j1)~=0)
            i = i1;
            j = j1;
            findpixel = true;
            break;
        end;
    end;
    if(findpixel == false)   % find non-zero pixels in the farther neighbours
        total_nb_cand = size(N2,2);
        for nb_cand = 1 : total_nb_cand
            i1 = i + N2(1,nb_cand);
            j1 = j + N2(2, nb_cand);
           if (Leaf_edge(i1,j1)~=0)
                i = i1;
                j = j1;
                findpixel = true;
                break;
            end;
        end;
    end;
    if(findpixel == false)
        total_nb_cand = size(N3,2);
        for nb_cand = 1 : total_nb_cand
            i1 = i + N3(1,nb_cand);
            j1 = j + N3(2, nb_cand);
           if (Leaf_edge(i1,j1)~=0)
                i = i1;
                j = j1;