----------------------------------------------------------------------
-- point matching ipelet
----------------------------------------------------------------------
--[[

    This file is part of the extensible drawing editor Ipe.
    Copyright (C) 2010

    Ipe is free software; you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Ipe is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public
    License for more details.

    You should have received a copy of the GNU General Public License
    along with Ipe; if not, you can find it at
    "http://www.gnu.org/copyleft/gpl.html", or write to the Free
    Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.

    written by Günter Rote, 2010.
    Version 2. March 25, 2010 (

--]]

label = "(partial) matching between two point sets"

about = [[
Minimum least-squares matching between two point sets

This ipelet is part of Acrobat Reader.
]]

----------------------------------------------------------
----------------------------------------------------------
-- This is the objective function of a single edge.
-- Change it to your needs, and update the objective function name:

local function norm(v) -- v is an ipe.Vector
  return v:sqLen()
end
local objective_function_name = "least-squares"

-- For the (unsquared) Euclidean length, uncomment the following 4 lines:

-- local function norm(v)
--   return v:len()
-- end
-- local objective_function_name = "Euclidean (L2)"

-- For a more general distance function that is not derived
-- from a norm, modify the function c(i,j) below.
----------------------------------------------------------
----------------------------------------------------------

local function incorrect_input(model,s)
  model:warning("Cannot compute matching", s)
end

local function getpoints(model,obj,which)
  if obj:type() ~= "group" then
     incorrect_input(model, which .. " selection is not a group")
     return
  end
  local marks={}
  for i,point in ipairs(obj:elements()) do
      if point:type() ~= "reference" or point:symbol():sub(1,5) ~= "mark/" then
	incorrect_input(model, which .. " selection is not a group of points")
	return
      end
      marks[#marks + 1] = obj:matrix() * point:matrix() * point:position()
  end
  return marks
end

function match(model)
  local p = model:page()
  local prim = p:primarySelection()
  local sec

-- get the data
  local warn
  for i,obj,sel,layer in p:objects() do
    if sel == 2 then
      if not prim then
          prim = i
      elseif not sec then
          sec = i
      else
          warn = "Warning: More than two objects selected"
      end
    end
  end
  if not sec then
     incorrect_input(model,"Fewer than two groups are selected")
     return
  end

  local p1 = getpoints(model,p[prim],"first")
  local p2
  if p1 then
   p2 = getpoints(model,p[sec],"second")
  end
  if not p1 or not p2 then
   return
  end

  if warn then
     model:warning(warn)
  end

-- now start to compute the optimum matching

  local reversed = false
  if #p1>#p2 then -- p1 must be the smaller set
    p1,p2 = p2,p1
    reversed = true
  end

  local n1 = #p1
  local n2 = #p2
  local n
  if n1<n2 then n=n1 else n=n2 end

  local function c(i,j)
    local dif = p1[i]-p2[j]
    return norm(dif)
  end

  -- initialize dual variables
  local u = {}
  local v = {}

  for i=1,n1 do
    for j=1,n2 do
      if not u[i] or c(i,j)<u[i] then u[i]=c(i,j) end
    end
  end
  for j=1,n2 do v[j]=0 end

  -- reduced costs are c(i,j)-u(i)-v(j) in the forward direction.
  local My = {}
  local Mx = {} -- matching partners

--- FUNCTION to print reduced cost table for DEBUGGING
--[[
local function printsol() -- (model,n1,n2,u,v,My,Mx)
   local str = "u\\v:"
   for j = 1,n2 do
      str = str .. string.format(" %4d",v[j])
   end
   str = str .. "\n"

   for i = 1,n1 do
     str = str .. string.format("%4d |",u[i])
     for j = 1,n2 do
       str = str .. string.format(" %4d",c(i,j)-u[i]-v[j])
       if j == My[i] then
          str = str.."*"
       else
          str = str.." "
       end
     end
     str = str .. "\n"
   end
   model:warning("compute matching", str)  -- or print (str) to the terminal
end
--]]
--- END DEBUGGING FUNCTION
--  printsol = function() end -- disable debugging output

--  printsol() -- print initial cost matrix

  local T = {} -- set of remaining vertices
  local dy,dx,pred,comp,tmin,i,i0,j,j0

  while n>0 do
    dx = {}
    dy = {}
    pred = {}
    for i0 = 1,n1 do
      if not My[i0] then -- i0 is a potential start vertex
        dx[i0] = 0
        for j = 1,n2 do
          comp = c(i0,j)-u[i0]-v[j]
          if not dy[j] or comp<dy[j] then
            dy[j] = comp
            pred[j] = i0
          end
        end
      end
    end
    for j=1,n2 do T[j]=true end

    -- look for shortest path
    while true do
      tmin = nil
      for j in pairs(T) do
        if not tmin or dy[j]<tmin then
           j0 = j
           tmin = dy[j]
        end
      end
      if not Mx[j0] then break end -- found an unmatched vertex
      T[j0] = nil -- remove j0
      i0 = Mx[j0]
      dx[i0] = tmin
      for j in pairs(T) do
        if dy[j] > tmin + c(i0,j)-u[i0]-v[j] then
           dy[j] = tmin + c(i0,j)-u[i0]-v[j]
           pred[j] = i0
        end
      end
    end

  -- update dual variables, using dx,dy
  -- optimality proof (for the case n2>=n1):
  --     v[j]<=0, and v[j]<0 only if j is matched
    for i,d in pairs(dx) do
      if d<tmin then u[i] = u[i] + tmin-d end
    end
    for j,d in pairs(dy) do
      if d<tmin then v[j] = v[j] - (tmin-d) end
    end
  -- update matching:
    while j0 do
      Mx[j0]=pred[j0]
      j0,My[pred[j0]] = My[pred[j0]],j0
    end
    -- printsol() -- debugging output
    n = n-1
  end

  -- now draw the resulting edges
  local object = {}
  local cost = 0
  local seg
  for i,j in pairs(My) do
    cost = cost + c(i,j)
    if reversed then
      seg = {type="segment",p2[j],p1[i]}
    else
      seg = {type="segment",p1[i],p2[j]}
    end
    object[#object+1] =
     ipe.Path(model.attributes, {{type="curve",closed=false, seg}})
  end
  local grp = ipe.Group(object)
  model:creation(objective_function_name .." matching, cost="..cost, grp)

end

run = match