function shps = loadRasterVectors(filename,boundary,params)

% Read data
[data,R] = geotiffread(filename);

%Get polygons from raster image, simplify, and remove too small
polys = bwboundaries(data==params.class_num,'noholes');
npolys = numel(polys);
todel = false(npolys,1);
for p=1:npolys
    polys{p} = simplifyPolygon(polys{p},1.5);
    if(size(polys{p},1)<=2)
        todel(p) = true;
    end
end
polys(todel)=[];
npolys = numel(polys);

%Georeference polygons
shps = [];
for p=1:npolys
    shp = [];
    [shp.X,shp.Y] = intrinsicToWorld(R, polys{p}(:,2), polys{p}(:,1));
    shp.X = shp.X';
    shp.Y = shp.Y';
    
    %Ignore vectors that are completely outside boundary
    if(~isempty(boundary))
        inside = isPointInPolygon([shp.X',shp.Y'],boundary);
        if(~any(inside))
            continue;
        end
    end
    
    shp.BoundingBox(1) = min(shp.X);
    shp.BoundingBox(2) = min(shp.Y);
    shp.BoundingBox(3) = max(shp.X);
    shp.BoundingBox(4) = max(shp.Y);
    shp.offsets = zeros(1,2);
    
    shps = [shps;shp];
end

end


