% 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 noisefunction freq = images2ts(imgs, do_plot, option, method)close all;index = 1;seq_grup = zeros(1,1);          % return value of empty folder caseseqcoor_grup = zeros(1,1);      % store the value of profile coordinatesfor 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 lengthif 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 imagesC = 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 notChecking_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 pixelnum = 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 neighboursi = 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;