function obj = init_with_svg( svgfile,varargin )

pp = inputParser;

% resize according to database
pp.addParamValue('resizebydb',[],@(x)isempty(x) || isa(x,class(SBSRDATABASE)));
% pp.addParamValue('imgsize',[]); %the size of final image
% pp.addParamValue('aabb',[]); %[x,y,width,height] that specifies where the svg image will be in the output image
pp.addParamValue('opt',SBSRSKETCHOPT.default_for_vector);
pp.parse(varargin{:});

skopt = pp.Results.opt;
[~,curvelist] = util_read_svg(svgfile);

% take points
npts_per_curve = 100;
ptslist = cell(1,length(curvelist));
dirlist = cell(1,length(curvelist));
curvaturelist = cell(1,length(curvelist)); %curvatures

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

seglist = cell(1,length(curvelist));
ptsmarklist = cell(1,length(curvelist)); %pts(3,i)==1 iff ith point is a keypoint
for i=1:length(curvelist)
    fc = curvelist{i};
    tt = linspace(0,1,npts_per_curve);
    pts = fnval(fc,tt);
    pts = flipud(pts); %sync with image coordinate
    ptslist{i} = pts;
    
    % tangents
    fdir = fnder(fc);
    dirs = fnval(fdir,tt);
    dirs = flipud(dirs);
    dirs = util_normalize_vectors(dirs);
    dirlist{i} = dirs;
    
    % 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);
    curvaturelist{i} = k;
    
    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(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 at intersection points?
break_at_intersection = skopt.m_breakseg_by_intersection; 

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

%break polyline at high curvature points?
break_at_corner = skopt.m_breakseg_by_corner; 
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;

if 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(6,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;
    end
end

pargroupidslist = cell(1,length(ptsmarklist));
if break_at_corner || skopt.m_use_parallel_group
    ptsmarklist_for_parallel = cell(1,length(ptsmarklist));
    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);
%         for k=1:size(pts,2)
%             thisdir = dirs(:,k);
%             thispts = pts(1:2,k);
%             dir_from_bgpts = thispts - bgpts;
%             
%             angle = acosd(dot(bgdir,dir_from_bgpts));
%             if angle > thres_corner_angle
%                 maskcorner(k) = true;
%                 bgdir = thisdir;
%                 bgpts = thispts;
%             end
%         end
%         curvelen = sum(sqrt(sum(diff(pts,[],2).^2,1)));
        [~,idxbreak] = dpsimplify(pts(1:2,:)',linebreak_dist_tol);
%         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');
        maskcorner(idxbreak) = true;
        
        % 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
        
        if break_at_corner
            pts(BRKTYPE_CORNER,maskcorner) = true;
            ptsmarklist{i} = pts;
        end
        
        if skopt.m_use_parallel_group
            newpts = pts;
            newpts(3,:) = 0;
            newpts(3,maskcorner) = true;
            ptsmarklist_for_parallel{i} = newpts;
        end
    end
    
    if skopt.m_use_parallel_group
        pargroupidslist = assign_parallel_group(ptsmarklist_for_parallel);

%         %% show parallel grouping
%         tmp_pgids = cat(2,pargroupidslist{:});
%         tmp_npargroup = max(tmp_pgids);
%         tmp_pts = cat(2,ptsmarklist{:});
%         tmp_pts = tmp_pts(1:2,:);
%         figure;hold on;title('pargroup');
%         util_scatter_pts(tmp_pts,'.');
%         for i=1:tmp_npargroup
%             mk = util_marker(i,GLOBALCONST.PLOT_COLOR_MARKER_LIST);
%             mask = tmp_pgids == i;
%             util_scatter_pts(tmp_pts(:,mask),mk);
%         end
    end
end
show_ptsmarklist(ptsmarklist);

if 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;
%         idxbreak = find(cumlens>=thres_arclength_per_seg)+1;
        pts(BRKTYPE_ARCLENGTH,idxbreak) = true;
        ptsmarklist{i} = pts;
    end
end
show_ptsmarklist(ptsmarklist);

% scale to proper size
pts = cat(2,ptsmarklist{:});
pts = pts(1:2,:);
if isempty(pp.Results.resizebydb)
    [boxlens,mincorner] = util_make_aabb_from_pts(pts);
    imgsize = ceil(boxlens)+1;
    imgsize = imgsize(:)';
    for i=1:length(ptsmarklist)
        ptsmarklist{i}(1:2) = bsxfun(@plus,ptsmarklist{i}(1:2),-mincorner(:)+1);
    end
else
    boxlens_svg = util_make_aabb_from_pts(pts);
    ar = boxlens_svg(1)/boxlens_svg(2);
    
    % find example in database with closest h/w ratio
    K = 5;
    database = pp.Results.resizebydb;
    dbobjlist = database.m_dbobjlist;
    ars_db = zeros(1,length(dbobjlist));
    boxlens_db = zeros(2,length(dbobjlist));
    for i=1:length(dbobjlist)
        dbobj = dbobjlist{i};
        dbsketch = dbobj.get_sketch;
        [is,js] = find(dbsketch);
        dbaabb = util_make_aabb_from_pts([is(:) js(:)]');
        ars_db(i) = dbaabb(1)/dbaabb(2);
        boxlens_db(:,i) = dbaabb;
    end
    idxhit = knnsearch(ars_db(:),ar,'K',K);
    boxlens_expected = mean(boxlens_db(:,idxhit),2);
    
    scalef = min(boxlens_expected./boxlens_svg(:));
    pts = pts * scalef;
    [boxlens_svg_new,corner_svg_new] = util_make_aabb_from_pts(pts);
    boxcenter_svg_new = corner_svg_new(:) + boxlens_svg_new(:)/2;
    
    % translate the pts to the center of image
    imgsize = database.m_dbobjlist{1}.imgsize;
    imgcenter = floor((imgsize+1)/2);
    transvec = imgcenter(:) - boxcenter_svg_new(:);
    
    % although useless, add here for debug purpose
    pts = bsxfun(@plus,pts,transvec(:));
    
    % transform the pts
    for i=1:length(ptsmarklist)
        ptsmarklist{i}(1:2,:) = bsxfun(@plus,ptsmarklist{i}(1:2,:)*scalef,transvec(:));
    end
end
show_ptsmarklist(ptsmarklist);

% break up polylines
newptslist = {};
newdirlist = {};

% pargroupid(i)=k iff ith segment is in the jth parallel group
pargroupid_per_subseg = [];

% segcontgroup(i)=j iff ith subsegment belongs to jth curve
segcontgroup_per_subseg = [];

% segtiming(i)=t iff ith subsegment is drawn at time t
segtiming_per_subseg = [];

for i=1:length(ptsmarklist)
    pts = ptsmarklist{i};
    dirs = dirlist{i};
    
    markers = any(pts(3:end,:),1);
    idxkeypts = find(markers);
    nsubseg = length(idxkeypts)-1;
    for k=1:nsubseg
        idx = idxkeypts(k):idxkeypts(k+1);
        newptslist{end+1} = pts(1:2,idx);
        newdirlist{end+1} = dirs(1:2,idx);
        if ~isempty(pargroupidslist{i})
            pgids = pargroupidslist{i}(idx);
            pargroupid_per_subseg(end+1) = mode(pgids);
        end
    end
    segcontgroup_per_subseg(end+1:end+nsubseg) = i;
    
    dt = 0.01;
    segtiming_per_subseg(end+1:end+nsubseg) = i+(dt:dt:dt*nsubseg);
end
seglbs_per_subseg = 1:length(newptslist);

% rasterize curves in image
sketchimg = false(imgsize);
seglbsmap = zeros(imgsize);
gdx = zeros(imgsize);
gdy = zeros(imgsize);
imgpts_per_subseg = cell(1,length(newptslist));
for i=1:length(newptslist)
    ptsample = newptslist{i};
    dirsample = newdirlist{i};
    mask = util_rasterize_polyline(ptsample,imgsize);
    mask = srec_skeletonize_sketch(mask);
    
    % for each rasterized point, find the closest sample and use its
    % tangent as the sample tangent
    [is,js] = find(mask);
    ptsuse = [is(:),js(:)]';
    idxsample = knnsearch(ptsample',ptsuse');
    diruse = dirsample(:,idxsample);
    
    [~,idxsort] = sort(idxsample);
    imgpts_per_subseg{i} = ptsuse(:,idxsort);
    
    idx = find(mask);
    sketchimg(idx) = true;
    seglbsmap(idx) = seglbs_per_subseg(i);
    gdx(idx) = diruse(1,:);
    gdy(idx) = diruse(2,:);
end

% filter imgpts_per_subseg because some points may have been covered by
% other seglbs
for i=1:length(imgpts_per_subseg)
    pts = imgpts_per_subseg{i};
    inds = sub2ind(size(seglbsmap),pts(1,:),pts(2,:));
    seglbs = seglbsmap(inds);
    pts = pts(:,seglbs == i);
    imgpts_per_subseg{i} = pts;
end

obj = SBSRSKETCH;
obj.m_vecpts_per_seglb = newptslist;
obj.m_vecdirs_per_seglb = newdirlist;
obj.m_imgpts_per_seglb = imgpts_per_subseg;

[is,js] = find(sketchimg);
ptsketch = [is(:) js(:)]';
obj.m_cxpts = ptsketch;
obj.m_datafile = svgfile;
obj.m_gdx = gdx;
obj.m_gdy = gdy;
obj.m_original_sketchimg = sketchimg;

% continuity grouping
nseglb = max(seglbs_per_subseg);
segcontmat = zeros(nseglb,nseglb);
ncontgroup = max(segcontgroup_per_subseg);
for i=1:ncontgroup
    maskingroup = segcontgroup_per_subseg == i;
    seglbs = seglbs_per_subseg(maskingroup);
    segcontmat(seglbs,seglbs) = i;
end
obj.m_segcontgroup_seglb2seglb = segcontmat;

% parallel grouping
if skopt.m_use_parallel_group
    nseglb = max(seglbs_per_subseg);
    pargroupadjmat = zeros(nseglb,nseglb);
    npargroup = max(pargroupid_per_subseg);
    for i=1:npargroup
        maskingroup = pargroupid_per_subseg == i;
        seglbs = seglbs_per_subseg(maskingroup);
        pargroupadjmat(seglbs,seglbs) = i;
    end
    obj.m_parallelgroup_seg2seg = pargroupadjmat;
end

% straight line decomposition
show_ptsmarklist(ptsmarklist);
pts_per_line = {};
for i=1:length(ptsmarklist)
    pts = ptsmarklist{i};
    breakmask = pts(BRKTYPE_CORNER,:)>0;
    pts = pts(1:2,:);
    linelist = break_pts_into_segments(pts,breakmask);
    pts_per_line = cat(2,pts_per_line,linelist);
end
obj.m_vecpts_per_line = pts_per_line;

lineidx = cell(1,length(pts_per_line));
for i=1:length(pts_per_line)
    lineidx{i} = repmat(i,1,size(pts_per_line{i},2));
end
lineidx = cat(2,lineidx{:});
pts = cat(2,pts_per_line{:});
imgpts = util_pts_from_bwimg(sketchimg);
idx = knnsearch(pts',imgpts');
lineidx_per_imgpts = lineidx(idx);
imgpts_per_line = cell(1,length(pts_per_line));
for i=1:length(pts_per_line)
    thisidx = lineidx_per_imgpts == i;
    thisimgpts = imgpts(:,thisidx);
    
    % sort the points
    thislinepts = pts_per_line{i};
    idx = knnsearch(thislinepts',thisimgpts');
    [~,idxsort] = sort(idx);
    thisimgpts = thisimgpts(:,idxsort);
    
    imgpts_per_line{i} = thisimgpts;
end
obj.m_imgpts_per_line = imgpts_per_line;

seglbs_per_sketchpts = seglbsmap(sketchimg);
obj.m_seglbs_per_sketchpts = seglbs_per_sketchpts(:)';
obj.m_segtimes = segtiming_per_subseg;
obj.m_sketchimg = sketchimg;

% % find parallel group
% if skopt.m_use_parallel_group
%     fprintf(1,'SKETCH: grouping by parallel lines...\n');
%     skimg = sketchimg;
%     gdxy = cat(3,obj.m_gdx,obj.m_gdy);
%     idxpix = find(skimg);
%     pixpairs = srec_pixgroup_parallel_from_sketch(skimg,'gdxy',gdxy);
%     [~,pixpairs(1,:)] = ismember(pixpairs(1,:),idxpix);
%     [~,pixpairs(2,:)] = ismember(pixpairs(2,:),idxpix);
%     adjmat = sparse(pixpairs(1,:),pixpairs(2,:),true,nnz(skimg),nnz(skimg))>0;
%     adjmat = adjmat | adjmat';
%     obj.m_parallelgroup_sketchpts2sketchpts = adjmat;
% end

obj.m_sketchopt = skopt;

% sequence stroke points
obj = srec_sequence_stroke_pts(obj);

end

function [seglist,idxlist] = break_pts_into_segments(pts,maskbreak)

seglist = {}; %seglist{i} = pts in ith segment
idxlist = {}; %idxlist{i} = index of pts in ith segment
break_markers = maskbreak;
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);
    seglist{end+1} = newpts;
    idxlist{end+1} = inds;
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 show_ptsmarklist(ptsmarklist,breakid)
return;
if nargin == 1
    breakid = 5;
end
seglist = {};
for i=1:length(ptsmarklist)
    pts = ptsmarklist{i};
    tt = break_pts_into_segments(pts(1:2,:),pts(breakid,:)>0);
    seglist = cat(2,seglist,tt);
end

figure;hold on;
for i=1:length(seglist)
    hold on;
    pts = seglist{i};
    util_scatter_pts(pts,util_marker(i,GLOBALCONST.PLOT_DIFF_MARKER_LIST));
end
end





