st4orient = zeros(length(st4idx_X50), 1);

for i=1:length(st4idx_X50), 
    
    fprintf('%d (%d): ', i, st4idx_X50(i));
    
    % exclude ubiquitous patterns
    j = length(find(st4dmesh_X50(:, i) ~= 0));
    if j > 160,
        fprintf('ubiquitous (%d triangles stained)\n', j);
        continue;
    end
    
    % find the major domain
% - first sort the vectorized array (this way you don't matter about
%    its size:
% s = sort(array(:))
% - compute running differences and prepend 1 for the first value :
% d = [ 1 diff(s) ]
% - find where d is greater than 0 :
% f = find(d>0)
% - the values are s(f)
% - the counts are diff([0 f])
    
    
    if max(st4label_X50_2(:, i)) > 1,
        [ct_freq, ct_label] = hist(st4label_X50_2(st4label_X50_2(:,i) > 0, i), unique(st4label_X50_2(st4label_X50_2(:,i) > 0, i)));
        [ct_max, ct_idx] = max(ct_freq(ct_label));
        j = ct_label(ct_idx);
    else
        j = 1;
    end
    fprintf('label=%d, ', j);
    
    [VV, DD] = eig(cov([t_xy(1,st4label_X50_2(:,i) == j)', t_xy(2,st4label_X50_2(:,i) == j)' ]));
    [DDmval, DDmidx] = max(diag(DD));
    
    j = abs(atan(VV(1,DDmidx)/VV(2,DDmidx))*180/pi);
    fprintf('angle=%f ', j);
    
    
    if j < 20,
        st4orient(i) = 1;
    elseif j > 75
        st4orient(i) = 2;
    end
    
    fprintf('-> orientation=%d\n', st4orient(i));
end

clear ct_freq ct_label ct_max ct_idx;
clear VV DD DDmval DDmidx;