% Generation and detection of alignments in Gabor patterns % % This code is part of the following publication and was subject % to peer review: % % "Generation and detection of alignments in Gabor patterns" by % Samy Blusseau and Rafael Grompone von Gioi. % % Copyright (c) 2016 Samy Blusseau , Rafael Grompone % von Gioi % % This program is free software: you can redistribute it and/or modify % it under the terms of the GNU Affero General Public License as % published by the Free Software Foundation, either version 3 of the % License, or (at your option) any later version. % % This program is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the % GNU Affero General Public License for more details. % % You should have received a copy of the GNU Affero General Public % License along with this program. If not, see % . function aligned_gabors(fn_out_im, fn_out_im_al, fn_out_im_algo, ... fn_out_pts, fn_out_ors, fn_out_al, fn_out_nfa_detec, ... fn_out_nfa_target, fn_out_report, N_in, n_in, jitter_level_in) rng('shuffle'); N = str2double(N_in); n = str2double(n_in); jitter_level = str2double(jitter_level_in); orientation_jitter = jitter_level*pi/2;% Jitter level must belong to (0, 1] %% GENERATION ALGORITHM (Algorithm 1 in the manuscript) % Correspondance of variable names (when non identical): % ------------------------------------- % Manuscript <--> Source code % ------------------------------------- % \alpha <--> orientation_jitter % P <--> all_points % \Theta <--> all_orientations % idx_a <--> idx_alignment % \mathcal{I} <--> IMG_stimulus [all_points, all_orientations, p, idx_alignment, IMG_stimulus] ... = generation(N, n, orientation_jitter); % This is not part of the algorithm but useful for the demo % Same parameters for the stimulus image I = 500; % Size of the image freq = 0.12; % spatial frequence for the gabor patches, in cycles per pixel patch_side = 15; % size of each patch's side, in pixels sigma = 1/(4*freq);% space constant, in pixels IMG_alignment = build_stimulus(I, patch_side, ... all_points(idx_alignment,:),... all_orientations(idx_alignment), sigma, freq); %% Save stimulus files % write final image imwrite(IMG_stimulus, fn_out_im); % write image of alignment imwrite(IMG_alignment, fn_out_im_al); % save all points' cordinates and orientations, and aligned points' indexes dlmwrite(fn_out_pts, all_points, '\t'); dlmwrite(fn_out_ors, all_orientations, '\t'); dlmwrite(fn_out_al, idx_alignment, '\t'); %% DETECTION ALGORITHM (Algorithm 2 in the manuscript) % % Correspondance of variable names (when non identical): % ------------------------------------- % Manuscript <--> Source code % ------------------------------------- % x <--> points % \Theta <--> all_orientations % NFA^* <--> best_nfa % x_i, i\in C^* <--> best_line_pos % theta_i, i\in C^* <--> best_line_ors % C^* <--> best_chain % [best_nfa, best_line_pos, best_line_ors] = ... detection(all_points, all_orientations); dlmwrite(fn_out_nfa_detec, best_nfa); IMG_out = build_stimulus(I, patch_side, best_line_pos, best_line_ors, ... sigma, freq); imwrite(IMG_out, fn_out_im_algo); %% NFA of the target alignment % This is not part of the algorithm but useful for the demo nfa_target = nfa(idx_alignment, all_points, all_orientations); dlmwrite(fn_out_nfa_target, nfa_target); %% Report writing on generation, density check and detection % This is not part of the algorithm but useful for the demo report_stimulus = sprintf(['The resulting stimulus contains %d elements'... ', the target alignment %d elements'], ... length(all_points), length(idx_alignment)); if length(idx_alignment) < n, report_stimulus = sprintf(['%s (instead of %d:\nindeed, for a'... ' stimulus containing %d elements, we'... ' set the maximum number of aligned'... ' elements to %d.)'], report_stimulus,... n, N, length(idx_alignment)); else report_stimulus = sprintf('%s.', report_stimulus); end if p == 0, report_density_cue = sprintf(['The density cue could not be checked'... ' because not enough Voronoi cells'... ' are reliable in the alignment.\n']); else report_density_cue = sprintf(['The hypothesis of equal average' ... ' Voronoi yields a p-value p = '... '%0.2g.\n(A small p-value, typically'... ' p < 0.05, indicates a likely'... ' density cue.)\n'], p); end if nfa_target == best_nfa, report_detection = sprintf(['The detected structure is the target'... ' alignment and its NFA is %0.2g.'],... best_nfa); else report_detection = sprintf(['The NFA of the target alignment is'... ' %0.2g, and the NFA of the detected'... ' structure is %0.2g'], nfa_target, ... best_nfa); end report = sprintf('%s \n\n%s\n%s', report_stimulus, report_density_cue, ... report_detection); fid = fopen(fn_out_report, 'w'); fwrite(fid, report); fclose(fid); function [all_points, all_orientations, p, idx_alignment, IMG_stimulus] ... = generation(N, n, orientation_jitter) % [all_points, all_orientations, p, idx_alignment, IMG_stimulus] ... % = generation(N, n, orientation_jitter) % N : integer, total number of elements in the simulus. % n : integer, number of aligned (target) elements % orientation_jitter : float in [0, 1], level of angular jitter % all_points : N x 2 matrix containing the coordinates of the N elements % all_orientations : N x 1 vector containing the N elements' orientations % p : p-value of the density check statistical test % idx_alignment : 1 x n vector of integers; the indexes of the target % elements % IMG_stimulus : 500x500 pixels image; the actual stimulus % % % Correspondance of variable names (when non identical): % ------------------------------------- % Manuscript <--> Source code % ------------------------------------- % \alpha <--> orientation_jitter % h <--> patch_side % f <--> freq % P <--> all_points % \Theta <--> all_orientations % idx_a <--> idx_alignment % \mathcal{I} <--> IMG_stimulus % Parameters for the image I = 500; % Size of the image freq = 0.12; % spatial frequence for the gabor patches, in cycles per pixel patch_side = 15; % size of each patch's side, in pixels sigma = 1/(4*freq);% space constant, in pixels %% Set the coordinates of the patches [all_points, idx_alignment, theta_a] = set_coordinates(I, N, n); %% Check density cue p = check_density_cue(all_points, I, idx_alignment); %% Set orientations actual_N = length(all_points); all_orientations = set_orientations(theta_a, orientation_jitter, ... idx_alignment, actual_N); %% Stimulus rendering % Build images IMG_stimulus = build_stimulus(I, patch_side, ... all_points, all_orientations, sigma, freq); function [P, idx_alignment, theta_a] = set_coordinates(I, N, n_input) % % Correspondance of variable names (when non identical): % ------------------------------------- % Manuscript <--> Source code % ------------------------------------- % n_{in} <--> n_input % \rho <--> resolution % idx_a <--> idx_alignment % S <--> aligned_pts %% Set minimal distance between points lambda = 1.5; d_a = 2*I/(sqrt(2*sqrt(3)*N)+4);% between aligned points d_b = lambda*I/(sqrt(2*sqrt(3)*N)+2*lambda);% between background points %% Place the aligned points [aligned_pts, theta_a] = place_alignment(I, n_input, d_a, d_b); n = length(aligned_pts); idx_alignment = 1:n; %% Place background points resolution = 0.5; P = placeBackground(aligned_pts, N, resolution, I, d_b, d_b); function [aligned_pts, theta_a] = place_alignment(I, n_input, d_a, d_marg) % % Correspondance of variable names (when non identical): % ------------------------------------- % Manuscript <--> Source code % ------------------------------------- % n_{in} <--> n_input % S <--> aligned_pts % mu_x and mu_y define the position of the alignment's center % in the authorized rectangle [x_c_min x_c_max] x [y_c_min y_c_max] mu_x = rand; % in [0, 1]; the closer to 0 (resp. 1), the closer x_c %to x_c_min (resp. x_c_max) mu_y = rand; % in [0, 1]; the closer to 0 (resp. 1), the closer y_c %to y_c_min (resp. y_c_max) theta_a = 2*pi*rand; % orientation of the alignment, in [0, 2*pi] n_max = floor((I-2*d_marg)/d_a + 1);%Max. authorized num. of aligned points n = min(n_input, n_max);% Decreasing num. of aligned points if necessary L = (n - 1)*d_a; % length of the alignment % Set the authorized rectangle [x_c_min x_c_max] x [y_c_min y_c_max] % from which the alignment center is picked x_c_min = d_marg + 0.5*L*abs(cos(theta_a)); x_c_max = I - d_marg - 0.5*L*abs(cos(theta_a)); y_c_min = d_marg + 0.5*L*abs(sin(theta_a)); y_c_max = I - d_marg - 0.5*L*abs(sin(theta_a)); %% Place alignment % set the alignment center x_c = x_c_min + mu_x*( x_c_max - x_c_min ); y_c = y_c_min + mu_y*( y_c_max - y_c_min ); % set the aligned points' coordinates v_a = [cos(theta_a) sin(theta_a)];% alignment's directing vector pt1 = [x_c y_c] - 0.5*L*v_a;% first aligned point k_s = (0 : n-1)'; aligned_pts = ones(n, 1)*pt1 + d_a*k_s*v_a;% all aligned points function all_points = placeBackground(fixed_points, N, resolution, I,... d_marg, d_min) % % Correspondance of variable names (when non identical): % ------------------------------------- % Manuscript <--> Source code % ------------------------------------- % S <--> fixed_pts % r <--> resolution (replaced by r in the first line) % P <--> all_points % I' <--> I2 % d_min' <--> d_min2 % d_marg' <--> d_marg2 % S' <--> fp_2 % \mathcal{I}' <--> image2 % placed <--> n_placed_pts % B' and B <--> bg_points % \mathcal{A} <--> remaining_idx r = resolution; % pixelian precision % Side, in pixels, of the transformed image "image2" I2 = ceil(I/r); % Convert the minimum distance between elements d_min2 = ceil(d_min/r); % Convert the minimum distance from the image's edges d_marg2 = ceil(d_marg/r); % Convert the fixed element coordinates fp_2 = ceil(fixed_points/r); fp_x = fp_2(:, 1); fp_y = fp_2(:, 2); % Starting pool of possible positions: only margins are excluded % image2 is the binary image noted \mathcal{I}^{prime} in the article image2 = false(I2, I2); % Central square is turned to true (loop for, lines 6-10 in manuscript) image2(d_marg2 + 1 : I2 - d_marg2, d_marg2 + 1 : I2 - d_marg2) = true; % Create a disk with a radius equal to d_min % to dilate white pixels in a binary image with dilate_disk = strel('disk',d_min2,0); pool = false(I2, I2); %% Place existing elements % Mark the fixed elements with single white pixels in a sea of % black pool(sub2ind([I2 I2],fp_x,fp_y)) = true; % Dilate them using the disk with the d_min radius pool = imdilate(pool, dilate_disk); % Subtract them from the overall pool of possible positions image2 = image2 & ~pool; %% Start placing background points n_placed_pts = length(fixed_points); bg_points = []; % Indexes of the points in the set noted \mathcal{A} in the article % Determine which positions are still possible remaining_idx = find(image2); while n_placed_pts < N && ~isempty(remaining_idx), % Draw a batch of candidate points, keep only the unique ones new_point = remaining_idx(randi(length(remaining_idx))); [new_x,new_y] = ind2sub([I2 I2],new_point); bg_points = [bg_points; new_x, new_y]; pool(:) = false; pool(new_point) = true; pool = imdilate(pool,dilate_disk); image2 = image2 & ~pool; remaining_idx = find(image2);% Positions that are still possible n_placed_pts = n_placed_pts+1; end % Convert points coordinates back to fit in an image with dimensions dims bg_points = (bg_points - 1)*r + rand(size(bg_points))*r; all_points = [fixed_points; bg_points]; function orientations = set_orientations(theta_a, orientation_jitter, ... idx_alignment, N) % % Correspondance of variable names (when non identical): % ------------------------------------- % Manuscript <--> Source code % ------------------------------------- % \alpha <--> orientation_jitter % idx_a <--> idx_alignment % \Theta <--> orientations n = length(idx_alignment);% Number of aligned points orientations = 2*pi*rand(N, 1);% Initialize all orientations randomly jitter = orientation_jitter*(2*rand(n,1) - 1);% orientation noise to add orientations_alignment = theta_a*ones(n, 1) + jitter; orientations(idx_alignment) = orientations_alignment; function im = build_stimulus(img_side, patch_side, ... points, orientations, sigma, f) % % Correspondance of variable names (when non identical): % ------------------------------------- % Manuscript <--> Source code % ------------------------------------- % I <--> img_side % h <--> patch_side % P <--> points % \Theta <--> orientations %\mathcal{I} <--> im % Here we work with square images and square patches img_height = img_side; img_width = img_side; patch_h = patch_side; patch_w = patch_side; N = size(points, 1); im = 0.5*ones(img_height, img_width); sc_w= 1 : patch_w ; sc_h= patch_h : -1 : 1 ; [x y] = meshgrid(sc_w,sc_h); cw = (1 + patch_w)/2; ch = (1 + patch_h)/2; for i = 1:N, gab_patch = G((x-ch), (y-cw), sigma, f, orientations(i)); xc = points(i,1); yc = img_width - points(i,2); x_start = round(xc - patch_w/2); y_start = round(yc - patch_h/2); im(y_start : y_start + patch_h -1, x_start : x_start + patch_w -1)... = gab_patch; end function val = G(x, y, sigma, f, theta) sin_fun = cos( 2*pi*f*(-sin(theta)*x + cos(theta)*y)); exp_fun = exp(-0.5*(x.^2 + y.^2)/sigma^2); val = 0.5*( 1 + sin_fun.*exp_fun ); function p = check_density_cue(all_points, I, idx_alignment) % % Correspondance of variable names (when non identical): % ------------------------------------- % Manuscript <--> Source code % ------------------------------------- % P <--> all_points % idx_a <--> idx_alignment % \mathcal{S}_a <--> locsurf_al % \mathcal{S}_b <--> locsurf_bg % Construct the Voronoi diagram [V,C] = voronoin(all_points); % Retrieve the vertices of each cell, remove those cells that have % vertices outside the display dimensions cells_inside = false(1,size(C,1)); for i = 1:size(C,1) Vx = V(C{i},1); Vy = V(C{i},2); if all(inpolygon(Vx,Vy,[0 0 I I],[0 I I 0]) ) cells_inside(i) = true; end end idx_in = find(cells_inside); % indexes of the Voronoi cells included in the image idx_in_align = idx_alignment(ismember(idx_alignment, idx_in)); % indexes of the the aligned points' Voronoi cells included in the image idx_in_bg = idx_in(~ismember(idx_in, idx_alignment)); % indexes of the the background points' Voronoi cells included in the image if length(idx_in_align) < 0.8*length(idx_alignment), p = 0; return end % Compute the surface of each cell locsurf = zeros(1,size(C,1)); for i = 1:size(C,1) Vx = V(C{i},1); Vy = V(C{i},2); locsurf(i) = polyarea(Vx,Vy); end % Surface of the Voronoi cells locsurf_al = locsurf(idx_in_align);% Cells containing aligned elements locsurf_bg = locsurf(idx_in_bg);% Cells containing background elements % The Matlab function ttest2 requires a rejection threshold alhpa, but we % prefer to provide only the p-value and let the user choose the threshold. p_min = 0.05; % Rejection threshold, not used here % but is necessary to use ttest2 [~, p] = ttest2(locsurf_al, locsurf_bg, p_min,'both','equal'); function [best_nfa, best_line_pos, best_line_ors] = detection(points,... orientations) % [best_nfa best_line_pos best_line_ors] = detection(points, orientations) % % - points : a N x 2 matrix; each line contains the coordinates of a point % - orientations : a vector of length N, containing the orientations of the % patches % % - best_nfa : the smallest found NFA % - best_line_pos : the coordinates of the elements of the detected % alignment % - best_line_ors : the orientations of the elements of the detected % alignment % % Correspondance of variable names (when non identical): % ------------------------------------- % Manuscript <--> Source code % ------------------------------------- % x <--> points % \Theta <--> orientations % NFA^* <--> best_nfa % x_i, i\in C^* <--> best_line_pos % theta_i, i\in C^* <--> best_line_ors % C^* <--> best_chain % lastElt <--> last_point N = length(points);% total number of elements gamma = compute_gamma(points); %% Explore the family of tests and compute the NFA for each tested chain best_nfa = inf; best_chain = []; for i = 1:N, % Testing chains starting at point i neighbours_of_i = find(gamma(i,:)); % indexes of the neighbours of % point i for v = 1: length(neighbours_of_i), j = neighbours_of_i(v); chain = [i, j]; % initialize a chain last_point = chain(2); while length(chain) < floor(sqrt(N)) && ~isempty(last_point), last_point = find_next_point(points, gamma, chain); chain = [chain, last_point]; nfa_chain = nfa(chain, points, orientations); if nfa_chain < best_nfa, best_nfa = nfa_chain; best_chain = chain; end end end end best_line_pos = points(best_chain, :); % coordinates of the elements in the detected alignment best_line_ors = orientations(best_chain, :); % orientations of the elements in the detected alignment function gamma = compute_gamma(points) % gamma = compute_gamma(points) % % Compute the 6 nearest nearest neighbors oriented graph and removes % arrows longer than 2avg_d % % - points : a N x 2 matrix; each line contains the coordinates of a point % - gamma : a N x N adjacence matrix representing an oriented graph. % gamma(i, j) = 1 if there is an arrow from point i to point j. % gamma(i, j) = 0 otherwise. dist_matrix = dist2(points, points); % matrix of square distances between points nn = 6; % number of nearest neighbours to consider % sort distances in ascending order, row by row; % variable idx contains the new order of the original distance values [D_sorted, idx] = sort(dist_matrix,2); % Collect all minimum distances to nearest neighbor, D_sorted(:, 2); % D_sorted(:, 1) contains only zeros: distances of a point to itself min_distances = D_sorted(:, 2); avg_d = mean(sqrt(min_distances)); gamma = zeros(size(dist_matrix)); for i = 1: size(dist_matrix, 1), d_i = D_sorted(i, :);% all distances to point i, in ascending order idx_i = idx(i,:);% the new order of the original distance values for j = 2:nn+1, % check the nn nearest neighbors of point i if d_i(j) <= 4*avg_d^2, gamma(i, idx_i(j)) = 1; end end end function n2 = dist2(x, c) %DIST2 Calculates squared distance between two sets of points. % % Description % D = DIST2(X, C) takes two matrices of vectors and calculates the % squared Euclidean distance between them. Both matrices must be of % the same column dimension. If X has M rows and N columns, and C has % L rows and N columns, then the result has M rows and L columns. The % I, Jth entry is the squared distance from the Ith row of X to the % Jth row of C. % % See also % GMMACTIV, KMEANS, RBFFWD % % Copyright (c) Christopher M Bishop, Ian T Nabney (1996, 1997) [ndata, dimx] = size(x); [ncentres, dimc] = size(c); if dimx ~= dimc error('Data dimension does not match dimension of centres') end n2 = (ones(ncentres, 1) * sum((x.^2)', 1))' + ... ones(ndata, 1) * sum((c.^2)',1) - ... 2.*(x*(c')); function index_next_point = find_next_point(points, adj_matrix, chain) % index_next_point = find_next_point(points, adj_matrix, chain) % % Finds, if any, the point that expands the already existing chain with % minimum change in direction. This new point must be one of the neighbors % of the chain's last point, in the graph defined by the adjacency matrix. % % - points : a N x 2 matrix; each line contains the coordinates of a point % - adj_matrix : a N x N adjacence matrix representing an oriented graph. % adj_matrix(i, j) = 1 if there is an arrow from point i to point j. % adj_matrix(i, j) = 0 otherwise. % - chain : a set of indexes defining the current chain as points(chain,:). % % - index_next_point : the index of the point supposed to expand the % current chain. i1 = chain(end-1); p1 = points(i1, :); i2 = chain(end); p2 = points(i2, :); max_cos_theta = -inf; neighbors_i2 = find(adj_matrix(i2, :)); i_star = [];% will receive the index of the point achieving the smallest % direction change for i3 = neighbors_i2, % neighbors of the last point in the current chain p3 = points(i3, :); % candidate to be added to the chain % compute cos_theta: cosine of the angle bewteen vectors p1p2 and p2p3 scalar_prod = (p2-p1)*(p3-p2)'; cos_theta = scalar_prod/sqrt(((p2-p1)*(p2-p1)')*((p3-p2)*(p3-p2)')); if cos_theta > max_cos_theta, max_cos_theta = cos_theta; i_star = i3; end end if ~ismember(i_star, chain), index_next_point = i_star; else index_next_point = []; end function [nfa] = nfa(chain, points, orientations) n = length(chain); N = length(points); %% compute angles omega omegas = zeros(size(chain)); omegas(1) = abs( angular_diff(points(chain(1), :),... points(chain(2), :),... orientations(chain(1))) ); omegas(end) = abs( angular_diff(points(chain(end-1), :),... points(chain(end), :),... orientations(chain(end))) ); for i = 2: n-1, alpha_iL = abs( angular_diff(points(chain(i-1), :),... points(chain(i), :),... orientations(chain(i))) ); alpha_iR = abs( angular_diff(points(chain(i), :),... points(chain(i+1), :),... orientations(chain(i))) ); omegas(i) = max(alpha_iL, alpha_iR); end %% omega_star = max(omegas); p = (2*omega_star/pi)^2*F(omega_star)^(n-2); Ntests = 6*N*sqrt(N); nfa = Ntests*p; function [delta_ang] = angular_diff(pt1, pt2, theta) % [angular_errors] = angular_diff(pt1, pt2, ors) % % - pt1 and pt2 are two 1 x 2 vectors, each representing a point's % coordinates. % - theta is an angle in [0, 2*pi) % - delta_ang is the angular difference, in (-pi/2, pi/2), between theta % and theta_12 % % Returns an angular difference in [-pi/2, pi/2], % between the orientation theta and the direction defined % by the line (pt1, pt2). theta_12 = arg(pt1, pt2); % the direction defined by the line (pt1, pt2) delta_ang = atan( tan( theta - theta_12 ) ); function [ang] = arg(pt1, pt2) % Returns the angle, in ]-pi/2, pi/2] btw line (pt1, pt2) and the x-axis a = pt1; b = pt2; if abs(b(1)-a(1))> 0; ang = atan( (b(2)-a(2)) / (b(1)-a(1))); else ang = pi/2 ; end function u = F(t) u1 = (12/pi^2)*(t >= 0).*(t < pi/12).*t.^2; u2 = (12/pi^2)*(t >= pi/12).*(t < 5*pi/12).*(pi*t/6 - (pi/12)^2); u3 = (12/pi^2)*(t >= 5*pi/12).*(t < pi/2).*((pi-2*t).*(3*t/2-pi/4)... + (2*t-pi/2).*(2*t-5*pi/6)); u4 = (t >= pi/2); u = u1 + u2 + u3 + u4;