% 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 <sblusseau@gmail.com>, 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
% <http://www.gnu.org/licenses/>.
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;
