%{
This file demonstrates how to create a component graph in order to formulate the global interpretation MIP problem.

The graph is about a set of candidate components, distributed at different locations on the sketch image, 
where each component correspondes to some portion of the sketch.

The essential information in the graph includes:
    * The sketch image, and its low-level segment structure. 
    Each sketch segment is a near-straight curve, which can be obtained
    by breaking the strokes at high-curvature points, X junctions and T
    junctions.

    * Locations and labels of the candidate components

    * Possible connections of the candidate components. The "quality" of
    each connection is also included, which represents the likelihood of
    the connection.

    * Fitness score of each candidate component, which quantifies how well
    it matches to the sketch.

    * Each candidate component covers which set of sketch segments

    * The conflict relationship among the set of the sketch segments 
    covered by a single candidate component. Two sketch segments conflict
    if we are sure that the candidate component will cover only one of them, 
    but currently we do not know which one, and the resolution of the
    coverage is left for the MIP to solve.

After the graph is created, it can be used to formulate a MIP problem, and
solve for the component selection. For detail, please read the paper.
%}

%% load example data
x=load('sample.mat');
fulldata = x.fulldata;
old_data = fulldata.m_samplegraph;

%% create a new component graph step-by-step
samplegraph = SBSRSAMPLEGRAPH;

%% set the sketch image
samplegraph.m_sketchimg = old_data.m_sketchimg;
imshow(samplegraph.m_sketchimg);

%% sequence the sketch points
[is,js] = find(samplegraph.m_sketchimg);
pts = [is,js]';
samplegraph.m_samplepts = pts;

%% break the strokes into sketch segments
% you can break it in any way you like, the output is a vector seglbs,
% where seglbs(i) is the segment id of samplegraph.m_samplepts(:,i)
% In the paper, the sketch strokes are broken at high curvature point, T
% junction and X junctions
seglbs = old_data.m_seglbs_per_sample;
samplegraph.m_seglbs_per_sample = seglbs;

% show the sketch segments in different colors
seglbmap = zeros(size(samplegraph.m_sketchimg));
seglbmap(samplegraph.m_sketchimg) = seglbs;
cm = distinguishable_colors(max(seglbs),'w');
cm = cat(1,[1 1 1],cm);
imshow(seglbmap+1,cm);

%% position of each candidate component
% samplegraph.m_compts(:,i) = the position (in image coordinate) of the ith candidate component 
samplegraph.m_compts = old_data.m_compts;

% show it
figure;
pts = samplegraph.m_compts;
imshow(~samplegraph.m_sketchimg);
hold on;
scatter(pts(2,:),pts(1,:));

%% label of the each component
% samplegraph.m_complabels(i) = the label of ith candidate component
samplegraph.m_complabels = old_data.m_complabels;

% show all candidate components that belong to the 'chair leg' class
mask = samplegraph.m_complabels == GLOBALCONST.CHAIR_LEG;
pts = samplegraph.m_compts(:,mask);
figure;
imshow(~samplegraph.m_sketchimg);
hold on;
scatter(pts(2,:),pts(1,:));

%% fitness of each component
% the higher the score, the better the component matches the sketch locally
% By default, the scores will be normalized to [0,1] when formulating 
% the MIP problem, though you can turn off the normalization if desired.
% see SBSRMIPOPT "cost quantization".
samplegraph.m_compscores = old_data.m_compscores;

%% coverage of each component
% A point of the sketch is said to be covered by a component, if it
% matches to some point of the component.
% samplegraph.m_adjmat_comp2sample(i,j)=1 iff sketch point j 
% (i.e.,samplegraph.m_samplepts(:,j)) is covered by component i
samplegraph.m_adjmat_comp2sample = old_data.m_adjmat_comp2sample;

% compute overlapping of the components
% samplegraph.m_overlapcost_comp2comp(i,j)=c if the overlap cost of
% component i and component j is c
samplegraph = samplegraph.update_overlapcost_comp2comp();

%% coverage per sketch segment
% samplegraph.m_coverlen_seg2comp(i,j)=k iff sketch segment i covers k
% percentage of component j.
% suppose complen(i)=the stroke length of component i, you can compute the
% samplegraph.m_coverlen_seg2comp as follows
%{
    nseg = max(samplegraph.m_seglbs_per_sample);
    ncomp =size(samplegraph.m_compts,2);
    samplegraph.m_coverlen_seg2comp = zeros(nseg,ncomp);
    for i=1:nseg
        mask = samplegraph.m_seglbs_per_sample == i;
        comp2sample = samplegraph.m_adjmat_comp2sample(:,mask);
        samplegraph.m_coverlen_seg2comp(i,:) = sum(comp2sample,2)./complen(:);
    end
%}
samplegraph.m_coverlen_seg2comp = old_data.m_coverlen_seg2comp;

%% connection quality
% samplegraph.m_structscore_comp2comp_raw(i,j) is the connection quality of
% connecting component i and component j. If it is 0, then component i and
% component j cannot connect
samplegraph.m_structscore_comp2comp_raw = old_data.m_structscore_comp2comp_raw;

%% impossible connection mask
% samplegraph.m_adjmat_invalid_comp2comp(i,j)=true iff component i and
% component j CANNOT connect. This may happen if these components cover
% nearly the same portion of the sketch, or their connection quality is too
% low, or their labels are incompatible, for example, a human head cannot
% connect to a human leg, etc. This condition will be formulated into the
% MIP problem.
samplegraph.m_adjmat_invalid_comp2comp = old_data.m_adjmat_invalid_comp2comp;

%% conflicting sketch segments
% samplegraph.m_conflict_seg2seg2comp(i,j,k)=true if sketch segment i and j
% CANNOT be simultaneously assigned to component k. The detail is explained
% in the "sketch segment assignment" constraint in the paper.
samplegraph.m_conflict_seg2seg2comp = old_data.m_conflict_seg2seg2comp;

%% symmetry enforcement
% If some candidate components are symmetric to each other (e.g., the left
% and right arms of a chair), you can force them to be selected together
% when solving the MIP. In order to use symmetry enforcement, you need to
% detect possibly symmetric components in some way, and record it in
% samplegraph.
% samplegraph.m_symask_comp2comp(i,j)=true if candidate components i and j
% are symmetric for sure.
% samplegraph.m_symscore_comp2comp(i,j)=c if the possibility of 
% candidate components i and j being symmetric is c. You can just copy
% m_symask_comp2comp here.
samplegraph.m_symask_comp2comp = old_data.m_symask_comp2comp;
samplegraph.m_symscore_comp2comp = double(samplegraph.m_symask_comp2comp);

%% all possible labels and their connectivity degree
% samplegraph.m_dblbs = all possible labels
% samplegraph.m_maxdgmat(i,j) = connection degree of label i and label j
%{  
    for example

    if:
    m_dblbs(3) = GLOBALCONST.CHAIR_ARM
    m_dblbs(6) = GLOBALCONST.CHAIR_SEAT
    m_maxdgmat(6,3) = 2
    m_maxdgmat(3,6) = 1

    then:
    1 CHAIR_SEAT can connect to at most 2 CHAIR_ARM
    1 CHAIR_ARM can connect to at most 1 CHAIR_SEAT
%}
samplegraph.m_dblbs = old_data.m_dblbs;
samplegraph.m_maxdgmat = old_data.m_maxdgmat;

%% formulate MIP
mipopt = SBSRMIPOPT; %there are a lot of options inside

%{
% sketch segment grouping, see the paper for detail.
% using this may lead to better or worse result.

sketchdata = fulldata.m_sketchdata;

%continuity grouping, segcontgroup{i} = sketch segments that belong to ith
%continuity group, because they belong to the same stroke
segcontgroup = sketchdata.segcontgroups;

% parallel grouping, segpargroups{i} = sketch segments that are grouped
% because they are nearly parallel
segpargroup = sketchdata.segpargroups;

% temporal grouping, segtimegroup{i} = sketch segments that are grouped
% because they are near in temporal space (that is, they belong to
% consecutive strokes)
segtimegroup = sketchdata.segtimegroups;

% put the grouping information into mipopt
seglbsgroup = cat(2,segcontgroup,segtimegroup,segpargroup);
mipopt.m_segcontgroup = segcontgroup;
mipopt.m_seglbsgroup = seglbsgroup;
%}

[gmodel,linmodel] = srec_formulate_mip_new(samplegraph,mipopt);

%% solve MIP
params.TimeLimit = 600; %time limit for solving the MIP problem, or you can use ctrl-c to stop it while solving and still get a suboptimal result
mip_solution = gurobi(gmodel,params);

%% explain the solution
y = linmodel.explain_solution(mip_solution.x, linmodel);

%% the selected components and their connections
% the selected components
idxcompts = y.idxcompts;

% the selected connections
% adjmat(i,j)=true means component idxcompts(i) and idxcompts(j) are
% determined to be connected
adjmat = y.adjmat(idxcompts,idxcompts);

% plot it
pts = samplegraph.m_compts(:,idxcompts);
lbs = samplegraph.m_complabels(idxcompts);
ulbs = unique(lbs);
figure;
imshow(~samplegraph.m_sketchimg);
hold on;
for i=1:length(ulbs)
	thislabel = ulbs(i);
	thispts = pts(:,lbs==thislabel);
	scatter(thispts(2,:),thispts(1,:));
end
namelist = util_get_name_by_label(ulbs);
legend(namelist);

% show selected connections
hold on;
[is,js] = find(tril(adjmat));
for i=1:length(is)
	p0 = pts(:,is(i));
	p1 = pts(:,js(i));
	plot([p0(2),p1(2)],[p0(1),p1(1)],'-.');
end

%% the sketch segment assignment
% seglbs2comp(i,j)=true iff sketch segment i is assigned to component j
seglbs2comp = y.maskassign_seglbs2comp;
