function obj = init_with_ptslist( ptslist, varargin )
% ptslist{i}(:,k) = kth point of ith stroke
% dirlist{i}(:,k) = unit tangent direction of ptslist{i}(:,k)
% curvaturelist{i}(:,k) = curvature at ptslist{i}(:,k)

pp = inputParser;

%switch (x,y) to (y,x) in order to match with the image coordinate
pp.addParamValue('flipxy',false);
pp.addParamValue('dirlist',{});
pp.addParamValue('curvaturelist',{});
pp.parse(varargin{:});
flipxy = pp.Results.flipxy;

dirlist = pp.Results.dirlist;
curvaturelist = pp.Results.curvaturelist;
ttlist = cell(1,length(ptslist));
if isempty(dirlist) || isempty(curvaturelist)
    curvelist = cellfun(@(x)util_curve_from_pts(x),ptslist,'uniformoutput',false);
    for i=1:length(curvelist)
        fc = curvelist{i};
        tmp = sp2pp(fc);
        tt = tmp.breaks;
        
%         tangent directions
        fdir = fnder(fc);
        dirs = fnval(fdir,tt);
        if flipxy
            dirs = flipud(dirs);
        end
        dirs = util_normalize_vectors(dirs);
        if isempty(pp.Results.dirlist)
            dirlist{i} = dirs;
        end
        
%         curvatures
        d1 = fnval(fnder(fc,1),tt);
        d2 = fnval(fnder(fc,2),tt);
        x1 = d1(1,:);
        y1 = d1(2,:);
        x2 = d2(1,:);
        y2 = d2(2,:);
        k = abs(x1.*y2-x2.*y1)./(x1.^2+y1.^2).^(3/2);
        if isempty(pp.Results.curvaturelist)
            curvaturelist{i} = k;
        end
    end
end

BRKTYPE_CORNER = 3;
BRKTYPE_INTERSECT = 4;
BRKTYPE_ARCLENGTH = 5;
BRKTYPE_ENDPOINT = 6;
BRKTYPE_JUNCTION = 7;

NDIM_PTSMARK = 7;

seglist = cell(1,length(ptslist));
ptsmarklist = cell(1,length(ptslist)); %pts(3,i)==1 iff ith point is a keypoint
for i=1:length(ptslist)
    pts = ptslist{i};
    ttlist{i} = linspace(0,1,size(pts,2));
    
    % tangents
    dirs = dirlist{i};
    if flipxy
        dirs = flipud(dirs);
    end
    dirs = util_normalize_vectors(dirs);
    dirlist{i} = dirs;
    
    nseg = size(pts,2)-1;
    
    % segids(:,i) = [idxcurve,idxseg]
    segids = cat(1,repmat(i,1,nseg),1:nseg);
    p1 = pts(:,1:end-1);
    p2 = pts(:,2:end);
    segpts = cat(1,p1,p2);
    seglist{i} = cat(1,segids,segpts);
    
    % marked pts
    ptsmark = pts;
    ptsmark(NDIM_PTSMARK,:) = 0;
    ptsmark(BRKTYPE_ENDPOINT,:) = 0;
    ptsmark(BRKTYPE_ENDPOINT,[1 end]) = true; %end points are marked
    ptsmarklist{i} = ptsmark;
end

%% BREAK POINT CONDITIONS
pts = cat(2,ptslist{:});
boxlens = util_make_aabb_from_pts(pts);

%break polyline per some amount of arclength?
thres_arclength_per_seg= 0.1 * norm(boxlens);

%break polyline at high curvature points?
thres_corner_curvature = 1/(0.05 * norm(boxlens)); %1/R, R is the radius of curvature
thres_corner_angle = 20; %if the accumulated turning has reached this degree, break
linebreak_dist_tol = norm(boxlens) * 0.1;

%% break at intersection
% split each polyline into several subpolylines at intersection points
% find out intersection points
segpts = cat(2,seglist{:});
[adjmat,xmat,ymat] = util_line_intersect(segpts(3:4,:),segpts(5:6,:));

% insert intersection points to polylines
[is,js]=find(adjmat);
for i=1:length(is)
    % segment identification
    idxseg = is(i); %the intersection point occurs at this segment
    curve_id = segpts(1,idxseg);
    seg_id = segpts(2,idxseg);

    % the intersection point
    x = xmat(is(i),js(i));
    y = ymat(is(i),js(i));

    % insert the point to polyline
    pts = ptsmarklist{curve_id};
    inspts = zeros(NDIM_PTSMARK,1);
    inspts(1:2) = [x y];
    inspts(BRKTYPE_INTERSECT) = true;
    pts = cat(2,pts(:,1:seg_id),inspts(:),pts(:,seg_id+1:end));
%     pts(BRKTYPE_INTERSECT,seg_id) = true;
    ptsmarklist{curve_id} = pts;

    % insert tangent
    dirs = dirlist{curve_id};
    dirs = cat(2,dirs(:,1:seg_id),dirs(:,seg_id),dirs(:,seg_id+1:end));
    dirlist{curve_id} = dirs;

    % insert curvature
    ks = curvaturelist{curve_id};
    ks = cat(2,ks(1:seg_id),ks(seg_id),ks(:,seg_id+1:end));
    curvaturelist{curve_id} = ks;
    
    % insert parameter
    tt = ttlist{curve_id};
    tt = cat(2,tt(1:seg_id),tt(seg_id),tt(seg_id+1:end));
    ttlist{curve_id} = tt;
end

%% break at corner
for i=1:length(ptsmarklist)
    pts = ptsmarklist{i};
    ks = curvaturelist{i};
    maskcorner = ks>=thres_corner_curvature;

    % also consider angle between current tangent and starting tangent
    % (accumulated turning)
    dirs = dirlist{i};
    bgdir = dirs(:,1);
    bgpts = pts(1:2,1);
    idxlastbreak = 1;
    for k=1:size(pts,2)
        thisdir = dirs(:,k);
        thispts = pts(1:2,k);
        
        % the first portion of the curve is head
        idxheadpts = idxlastbreak + ceil((k-idxlastbreak)/3);
        dirhead = mean(dirs(:,idxheadpts),2);
        dirhead = dirhead/norm(dirhead);
        
        dir_from_bgpts = thispts - bgpts;
        dir_from_bgpts = dir_from_bgpts/norm(dir_from_bgpts);

%         angle = acosd(dot(bgdir,dir_from_bgpts));
        angle = acosd(dot(dirhead,thisdir));
        if angle > thres_corner_angle
            maskcorner(k) = true;
            bgdir = thisdir;
            bgpts = thispts;
            idxlastbreak = k;
        end
    end
%         curvelen = sum(sqrt(sum(diff(pts,[],2).^2,1)));
%     [~,idxbreak] = dpsimplify(pts(1:2,:)',linebreak_dist_tol);
%     maskcorner(idxbreak) = true;
    if false
        close all
        figure;hold on;
        ptsall = cat(2,ptsmarklist{:});
        ptsall = ptsall(1:2,:);
        util_scatter_pts(ptsall,'bo');
        plot(pts(1,:),pts(2,:),'rx-');
        util_scatter_pts(pts(1:2,idxbreak),'gh');
    end

    % eliminate noise
    maskcorner = imclose(maskcorner,true(1,ceil(length(maskcorner)*0.05)));

    % for consecutive corner points, just take the middle one
    cc = bwconncomp(maskcorner);
    maskcorner(:) = false;
    for k=1:cc.NumObjects
        pixlist = cc.PixelIdxList{k};
        idxcenter = pixlist(ceil(end/2));
        maskcorner(idxcenter) = true;
    end

    pts(BRKTYPE_CORNER,maskcorner) = true;
    ptsmarklist{i} = pts;
end

%% break at arclength
for i=1:length(ptsmarklist)
    pts = ptsmarklist{i};
    p1 = pts(1:2,1:end-1);
    p2 = pts(1:2,2:end);
    dvec = p2-p1;
    seglens = sqrt(dot(dvec,dvec,1));
    cumlens = cumsum(seglens);
    seglbs_by_arclength = ones(1,size(pts,2));
    seglbs_by_arclength(2:end) = ceil(cumlens/thres_arclength_per_seg);
    idxbreak = find(diff(seglbs_by_arclength)~=0)+1;
    pts(BRKTYPE_ARCLENGTH,idxbreak) = true;
    ptsmarklist{i} = pts;
end

%% break at junction
pts_per_stroke = cellfun(@(x)x(1:2,:),ptsmarklist,'uniformoutput',false);
[adjmat,pts_stroke2stroke] = srec_find_stroke2stroke_link(pts_per_stroke);
[is,js] = find(tril(adjmat));
for i=1:length(is)
    ii = is(i);
    jj = js(i);
    
    pi = pts_stroke2stroke(ii,jj,1:2);
    pi = pi(:);
    [~,idxpi] = ismember(pi',pts_per_stroke{ii}','rows');
    assert(idxpi~=0);
    ptsmarklist{ii}(BRKTYPE_JUNCTION,idxpi) = true;
    
    pj = pts_stroke2stroke(ii,jj,3:4);
    pj = pj(:);
    [~,idxpj] = ismember(pj',pts_per_stroke{jj}','rows');
    assert(idxpj~=0);
    ptsmarklist{jj}(BRKTYPE_JUNCTION,idxpj) = true;
end

%% create vecimage object
obj = SBSRVECIMAGE;

ptsmark = cat(2,ptsmarklist{:});
pts = ptsmark(1:2,:);
obj.m_pts = pts;
obj.m_dirs = cat(2,dirlist{:});
obj.m_curvature = cat(2,curvaturelist{:});

segbystroke = zeros(1,size(pts,2));
stklens = cellfun(@(x)size(x,2),ptsmarklist);
cumlens = cumsum(stklens);
cumlens = [0,cumlens(:)'];
for i=1:length(cumlens)-1
    idxbegin = cumlens(i)+1;
    idxend = cumlens(i+1);
    segbystroke(idxbegin:idxend) = i;
end
obj.m_segby_stroke = segbystroke;

obj.m_segby_arclen = break_pts_into_segments(ptsmark,ptsmark(BRKTYPE_ARCLENGTH,:),segbystroke);
obj.m_segby_arclen = cleanup_short_seglbs(obj.m_segby_arclen,segbystroke);

obj.m_segby_intersect = break_pts_into_segments(ptsmark,ptsmark(BRKTYPE_INTERSECT,:),segbystroke);
obj.m_segby_intersect = cleanup_short_seglbs(obj.m_segby_intersect,segbystroke);

obj.m_segby_line = break_pts_into_segments(ptsmark,ptsmark(BRKTYPE_CORNER,:),segbystroke);
obj.m_segby_line = cleanup_short_seglbs(obj.m_segby_line,segbystroke);

obj.m_segby_junction = break_pts_into_segments(ptsmark,ptsmark(BRKTYPE_JUNCTION,:),segbystroke);
obj.m_segby_junction = cleanup_short_seglbs(obj.m_segby_junction,segbystroke);

obj.m_tt = cat(2,ttlist{:});

%% parallel grouping
ptspglist = cellfun(@(x)x([1 2 BRKTYPE_CORNER],:),ptsmarklist,'uniformoutput',false);
pgmarklist = assign_parallel_group(ptspglist);
obj.m_segby_parallel_group = cat(2,pgmarklist{:});
obj.m_segby_parallel_group = cleanup_short_seglbs(obj.m_segby_parallel_group,segbystroke);
end

function segmentlabels = break_pts_into_segments(pts,maskbreak,strokeidx)
% stokeidx(i) = stroke id for point i, points on different stroke cannot
% share the same segment label

pts = pts(1:2,:);
maskstrokebreak = [0,diff(strokeidx(:)')]>0;
maskbreak = maskbreak | maskstrokebreak;

% do not break at the end point, but break at the beginning because a
% segment includes its beginning
maskbreak(1) = true;
maskbreak(end) = false;
maskbreak([1 end]) = false;
idxbreakpts = find(maskbreak);

segmentlabels = zeros(1,size(pts,2));
if isempty(idxbreakpts)
    segmentlabels(:) = 1;
else
    for i=1:length(idxbreakpts)
        idxthis = idxbreakpts(i);
        if i==length(idxbreakpts)
            % this is the last break point
            segmentlabels(idxthis:end) = i;
        else
            %this is the middle break point, count from this break point to
            %all the points ahead
            idxnext = idxbreakpts(i+1);
            segmentlabels(idxthis:idxnext-1) = i;
        end
    end
end

end

function pgmarklist = assign_parallel_group(ptsmarklist)
% ptsmarklist{i}(3,:) = break point markers, 1 if that point is going to
% break the segment
% pgmarklist{i}(k) = parallel group id for ptsmarklist{i}(:,k)

pts = cat(2,ptsmarklist{:});
pts = pts(1:2,:);
boxlens = util_make_aabb_from_pts(pts);
diaglens = norm(boxlens);
% thresholds
thres_max_dist_check = 0.1 * diaglens; %if two subsegments are within this distance, they are checked for parallel grouping
thres_max_angle_parallel = 20; %if two subsegments are within this incident angle, they are considered parallel
thres_max_step_check = 2; %if the temporal distance of two strokes are within this range, they are checked for parallel
thres_max_lens_diff = 0.3; % the lengths of the two segments should not differ too much for them to be considered in a group

% break into segments
seglist = {}; %pts in seglist has the form: pts(3,i)=curveid, pts(4,i)=index in that curve
for i=1:length(ptsmarklist)
    pts = ptsmarklist{i};
    break_markers = pts(3,:)>0;
    pts = pts(1:2,:);
    
    contmask = break_markers == 0;
    cc = bwconncomp(contmask);
    for k=1:cc.NumObjects
        inds = cc.PixelIdxList{k};
        
        % include starting point
        if inds(1)>1
            inds = [inds(1)-1,inds(:)'];
        end
        
        % include ending point
        if inds(end)<size(pts,2)
            inds(end+1) = inds(end)+1;
        end
        
        newpts = pts(:,inds);
        newpts(3,:) = i;
        newpts(4,:) = inds;
        seglist{end+1} = newpts;
    end
end

% %% show segments
% figure;hold on;title('straight segment');
% for i=1:length(seglist)
%     pts = seglist{i}(1:2,:);
%     mk = util_marker(i,GLOBALCONST.PLOT_DIFF_MARKER_LIST);
%     util_scatter_pts(pts,mk);
% end

%%

% compute segment directions
segdirs = zeros(2,length(seglist));
segcenters = zeros(2,length(seglist));
segsteps = zeros(1,length(seglist));
seglens = zeros(1,length(seglist));
for i=1:length(seglist)
    pts = seglist{i};
    segsteps(:,i) = pts(3,1);
    
    pts = pts(1:2,:);
    dirs = pts(:,2:end) - pts(:,1:end-1);
    lens = sqrt(sum(dirs.^2,1));
    seglens(i) = sum(lens);
    
    dirs = util_normalize_vectors(dirs);
    dirs = mean(dirs,2);
    dirs = dirs/norm(dirs);
    segdirs(:,i) = dirs;
    
    center = mean(pts,2);
    segcenters(:,i) = center;
end

% compute parallel grouping
%adjmat(i,j)==1 iff ith and jth segment are in the same paralell group
adjmat_seg2seg = false(length(seglist)); 
for i=1:length(seglist)
    for j=i+1:length(seglist)
        %% show segment pair
%         if i==18
%             tmp_pts = cat(2,seglist{:});
%             tmp_pts = tmp_pts(1:2,:);
%             figure;hold on;title('segment pair');
%             util_scatter_pts(tmp_pts,'.');
%             util_scatter_pts(seglist{i}(1:2,:),'ro');
%             util_scatter_pts(seglist{j}(1:2,:),'gh');
%         end
        
        %%
        c1 = segcenters(:,i);
        c2 = segcenters(:,j);
        dist = norm(c1-c2);
        
        dir1 = segdirs(:,i);
        dir2 = segdirs(:,j);
        angle = acosd(abs(dot(dir1,dir2)));
        
        step1 = segsteps(i);
        step2 = segsteps(j);
        sdist = abs(step1-step2);
        
        len1 = seglens(i);
        len2 = seglens(j);
        lendiff = abs(len1-len2)/mean([len1,len2]);
        
        is_parallel = false;
        if dist<=thres_max_dist_check || sdist<=thres_max_step_check
            is_parallel = angle <= thres_max_angle_parallel && lendiff <= thres_max_lens_diff;
        end
        adjmat_seg2seg(i,j) = is_parallel;
    end
end
adjmat_seg2seg = adjmat_seg2seg | adjmat_seg2seg';

% assign a group number to each segment
[~,pids] = graphconncomp(sparse(adjmat_seg2seg),'weak',true);

% some of the segments are not in any parallel group, their pids are set to
% zero
maskseg_has_group = any(adjmat_seg2seg,2);
pids(~maskseg_has_group) = 0;
[~,~,pids(maskseg_has_group)] = unique(pids(maskseg_has_group)); %re label from 1 to n

pgmarklist = cellfun(@(x)zeros(1,size(x,2)),ptsmarklist,'uniformoutput',false);
for i=1:length(seglist)
    pts = seglist{i};
    idxcurve = pts(3,1);
    inds = pts(4,:);
    thispid = pids(i);
    pgmarklist{idxcurve}(inds) = thispid;
end

end

function seglbs = cleanup_short_seglbs(seglbs,strokeidx)
% clean up short segments
% seglbs(i) = segment label for ith point

seglbs = SBSRVECIMAGE.cleanup_short_seglbs(seglbs,strokeidx);

end

