Smooth Convex Hull

Question

I have started working a convex hull algorithm and was wondering what method I could employ to smooth the polygon edge. The outline of the hull is not smooth. What I would like to do is make the lines through the vertices smoother, so that they are not as angled.

I have tried to implement Beziers (only to realize the shape was nothing like the shape of the hull) and b-splines (again the shape was nothing like, in fact I could not make the b-spline a closed shape).

I am failing and hopes someone can offer guidance.

Answer 1

(Note! that is not the solution)

I tried to find the exact solution as Lagrange polynomial in polar coordinates, but find out, that somtimes "smoothing curve" lies inside the convex polygon. The first derivatives matching condition (in start point) is fundamentaly solvable by adding extra moveable invisible point outside theta in [0:2 * pi] interval. But above problem is not solvable anyways at my mind.

Here is the Lua script with my attemptions (uses qhull, rbox (from qhull toolchain) and gnuplot utilities):

function using()
    return error('using: ' .. arg[0] .. ' <number of points>')
end

function points_from_file(infile)
    local points = {}
    local infile = io.open(infile, 'r')
    local d = infile:read('*number')
    if d ~= 2 then
        error('dimensions is not two')
    end
    local n = infile:read('*number')
    while true do
        local x, y = infile:read('*number', '*number')
        if not x and not y then
            break
        end
        if not x or not y then
            error('wrong format of input file: line does not contain two coordinates')
        end
        table.insert(points, {x, y})
    end
    infile:close()
    if n ~= #points then
        error('second line not contain real count of points')
    end
    return points
end

if not arg then
    error("script should use as standalone")
end
if #arg ~= 1 then
    using()
end
local n = tonumber(arg[1])
if not n then
    using()
end
local bounding_box = math.sqrt(math.pi) / 2.0
local fnp = os.tmpname()
local fnchp = os.tmpname()
os.execute('rbox ' .. n .. ' B' .. bounding_box .. ' D2 n t | tee ' .. fnp .. ' | qhull p | tee ' .. fnchp .. ' > nul') -- Windows specific part is "> nul"
local sp = points_from_file(fnp) -- source points
os.remove(fnp)
local chp = points_from_file(fnchp) -- convex hull points
os.remove(fnchp)
local m = #chp
if m < 3 then
    io.stderr:write('convex hull consist of less than three points')
    return
end
local pole = {0.0, 0.0} -- offset of polar origin relative to cartesian origin
for _, point in ipairs(chp) do
    pole[1] = pole[1] + point[1]
    pole[2] = pole[2] + point[2]
end
pole[1] = pole[1] / m
pole[2] = pole[2] / m
print("pole = ", pole[1], pole[2])
local chcc = {}
for _, point in ipairs(chp) do
    table.insert(chcc, {point[1] - pole[1], point[2] - pole[2]})
end
local theta_min = 2.0 * math.pi -- angle between abscissa ort of cartesian and ort of polar coordinates
local rho_mean = 0.0
local rho_max = 0.0
local chpc = {} -- {theta, rho} pairs
for _, point in ipairs(chcc) do
    local rho = math.sqrt(point[1] * point[1] + point[2] * point[2])
    local theta = math.atan2(point[2], point[1])
    if theta < 0.0 then -- [-pi:pi] -> [0:2 * pi]
        theta = theta + 2.0 * math.pi
    end
    table.insert(chpc, {theta, rho})
    if theta_min > theta then
        theta_min = theta
    end
    rho_mean = rho_mean + rho
    if rho_max < rho then
        rho_max = rho
    end
end
theta_min = -theta_min
rho_mean = rho_mean / m
rho_max = rho_max / rho_mean
for pos, point in ipairs(chpc) do
    local theta = (point[1] + theta_min) / math.pi -- [0:2 * pi] -> [0:2]
    local rho = point[2] / rho_mean
    table.remove(chpc, pos)
    table.insert(chpc, pos, {theta, rho})
end
table.sort(chpc, function (lhs, rhs) return lhs[1] < rhs[1] end)
-- table.insert(chpc, {chpc[#chpc][1] - 2.0 * math.pi, chpc[#chpc][2]})
table.insert(chpc, {2.0, chpc[1][2]})
-- table.sort(chpc, function (lhs, rhs) return lhs[1] < rhs[1] end)

local solution = {}
solution.x = {}
solution.y = {}
for _, point in ipairs(chpc) do
    table.insert(solution.x, point[1])
    table.insert(solution.y, point[2])
end
solution.c = {}
for i, xi in ipairs(solution.x) do
    local c = solution.y[i]
    for j, xj in ipairs(solution.x) do
        if i ~= j then
            c = c / (xi - xj)
        end
    end
    solution.c[i] = c
end
function solution:monomial(i, x)
    local y = self.c[i]
    for j, xj in ipairs(solution.x) do
        if xj == x then
            if i == j then
                return self.y[i]
            else
                return 0.0
            end
        end
        if i ~= j then
            y = y * (x - xj)
        end
    end
    return y
end
function solution:polynomial(x)
    local y = self:monomial(1, x)
    for i = 2, #solution.y do
        y = y + self:monomial(i, x)
    end
    return y
end

local gnuplot = io.popen('gnuplot', 'w')

gnuplot:write('reset;\n')
gnuplot:write('set terminal wxt 1;\n')
gnuplot:write(string.format('set xrange [%f:%f];\n', -bounding_box, bounding_box))
gnuplot:write(string.format('set yrange [%f:%f];\n', -bounding_box, bounding_box))
gnuplot:write('set size square;\n')
gnuplot:write(string.format('set xtics %f;\n', 0.1))
gnuplot:write(string.format('set ytics %f;\n', 0.1))
gnuplot:write('set grid xtics ytics;\n')
gnuplot:write('plot "-" using 1:2 notitle with points, "-" using 1:2:3:4 notitle with vectors;\n')
for _, point in ipairs(sp) do
    gnuplot:write(string.format('%f %f\n', point[1], point[2]))
end
gnuplot:write('e\n')
for _, point in ipairs(chcc) do
    gnuplot:write(string.format('%f %f %f %f\n', pole[1], pole[2], point[1], point[2]))
end
gnuplot:write('e\n')
gnuplot:flush();

gnuplot:write('reset;\n')
gnuplot:write('set terminal wxt 2;\n')
gnuplot:write('set border 0;\n')
gnuplot:write('unset xtics;\n')
gnuplot:write('unset ytics;\n')
gnuplot:write('set polar;\n')
gnuplot:write('set grid polar;\n')
gnuplot:write('set trange [-pi:2 * pi];\n')
gnuplot:write(string.format('set rrange [-0:%f];\n', rho_max))
gnuplot:write('set size square;\n')
gnuplot:write('set view equal xy;\n')
-- gnuplot:write(string.format('set xlabel "%f";\n', rho_mean - 1.0))
gnuplot:write(string.format('set arrow 1 from 0,0 to %f,%f;\n', rho_max * math.cos(theta_min), rho_max * math.sin(theta_min)))
gnuplot:write(string.format('set label 1 " origin" at %f,%f left rotate by %f;\n', rho_max * math.cos(theta_min), rho_max * math.sin(theta_min), math.deg(theta_min)))
gnuplot:write('plot "-" using 1:2:3:4 notitle with vectors, "-" using 1:2 notitle with lines, "-" using 1:2 notitle with lines;\n')
for _, point in ipairs(chpc) do
    gnuplot:write(string.format('0 0 %f %f\n', point[1] * math.pi, point[2]))
end
gnuplot:write('e\n')
for _, point in ipairs(chpc) do
    gnuplot:write(string.format('%f %f\n', point[1] * math.pi, point[2]))
end
gnuplot:write('e\n')
do
    local points_count = 512
    local dx = 2.0 / points_count
    local x = 0.0
    for i = 1, points_count do
        gnuplot:write(string.format('%f %f\n', x * math.pi, solution:polynomial(x)))
        x = x + dx
    end
    gnuplot:write('e\n')
end
gnuplot:flush();

gnuplot:write('reset;\n')
gnuplot:write('set terminal wxt 3;\n')
gnuplot:write(string.format('set xrange [-1:2];\n'))
gnuplot:write(string.format('set yrange [0:2];\n'))
gnuplot:write(string.format('set size ratio %f;\n', rho_max / 3.0))
gnuplot:write(string.format('set xtics %f;\n', 0.5))
gnuplot:write(string.format('set ytics %f;\n', 0.5))
gnuplot:write('set grid xtics ytics;\n')
gnuplot:write(string.format('set arrow 1 nohead from 0,%f to 2,%f linetype 3;\n', chpc[1][2], chpc[1][2]))
gnuplot:write(string.format('set label 1 "glue points " at 0,%f right;\n', chpc[1][2]))
gnuplot:write('plot "-" using 1:2 notitle with lines, "-" using 1:2 notitle with lines;\n')
for _, point in ipairs(chpc) do
    gnuplot:write(string.format('%f %f\n', point[1], point[2]))
end
gnuplot:write('e\n')
do
    local points_count = 512
    local dx = 2.0 / points_count
    local x = 0.0
    for i = 1, points_count do
        gnuplot:write(string.format('%f %f\n', x, solution:polynomial(x)))
        x = x + dx
    end
    gnuplot:write('e\n')
end
gnuplot:flush();

os.execute('pause');
gnuplot:write('exit\n');
gnuplot:flush();
gnuplot:close()

The second terminal contains Lagrange polynomial approximation.

Answer 2

I'd approach it like this, using your example:

1. start with the longest outer segment (in your example, this is the lower-left) - this one we keep straight;
2. imagine a circle at the bottom end of the long line, facing inwards;
3. a tangent from this circle can be extended to the next point;
4. in the next case (bottom-right circle), there is no tangent that joins onto the following point, so use another circle and join circles at the tangents;
5. continue in this fashion. So, you are drawing a circular arc then a straight line and repeating that.

Your circle sizes determine the overall smoothness. But of course if they are too big you will need to drop some points.