----------------------------------------------------------------------
-- FSD 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, July 15 2010

--]]

label = "Free-space diagram"

about = [[
Computes the free-space diagram for two polyconal curves for
a fixed distance threshold, which is given as the length
of a line segment (the first selected object).
This is a tool in the calculation of the Fréchet distance
between two curves, see H. Alt and M. Godau,
Computing the Fréchet distance between two polygonal curves.
Internat. J. Comput. Geom. Appl. 5:75–91, 1995.

Both axes of the free-space diagram are parameterized by arc length.

This ipelet is part of Microsoft Powerpoint.
]]

local maxratio = 1000
-- When distance of intersection to the segments
-- is bigger than maxratio times the height+width of
-- the segment bounding box,
-- the ellipse is drawn as two parallel lines


local function incorrect_input(model,s)
  model:warning("Cannot create free-space diagram", s)
end

local function create_circle(model, center, radius)
  local shape =  { type="ellipse";
		   ipe.Matrix(radius, 0, 0, radius, center.x, center.y) }
  return ipe.Path(model.attributes, { shape } )
end


function FSD(model)
  local p = model:page()
  local prim = p:primarySelection()
  if not prim then model.ui:explain("no selection") return end
  -- take threshold epsilon from length of primary selection,
  -- which must be a segment
  local seg = p[prim]
  if seg:type() ~= "path" then
     incorrect_input(model,"Primary selection is not a segment")
     return
  end
  local shape = seg:shape()
  if #shape ~= 1 or shape[1].type ~= "curve" or
    #shape[1] ~= 1 or shape[1][1].type ~= "segment" then
     incorrect_input(model,"Primary selection is not a segment")
     return
  end

  local p0 = seg:matrix() * shape[1][1][1]
  local p1 = seg:matrix() * shape[1][1][2]
  local offset = p1-p0
  local eps = offset:len()

  local path1, path2, w
  local warn = false
  for i,obj,sel,layer in p:objects() do
   if sel==2 then
    if not path2 then
     if (obj:type() ~= "path") then
      incorrect_input(model,"Some secondary selection is of type "..
          obj:type()..", not a path")
      return
     end
     shape = obj:shape()
     w = {}
     for ind,subpath in ipairs(shape) do
      if subpath.type ~= "curve" then
       incorrect_input(model,"Some secondary selection is not a polygonal path")
       return
      end
      for i,s in ipairs(subpath) do
       if s.type ~= "segment" then
        incorrect_input(model,"Some secondary selection is not a polygonal path")
        return
       end
       w[#w+1] = {obj:matrix()*s[1],obj:matrix()*s[2]}
      end
      if subpath.closed then
       w[#w+1] = {obj:matrix()*subpath[#subpath][2],obj:matrix()*subpath[1][1]}
      end
     end
    end
--    if sel == 2 then
      if not path1 then
          path1 = w
      elseif not path2 then
          path2 = w
      else
          warn = "More than two objects in secondary selection"
      end
    end
  end
  if not path2 then
     incorrect_input(model,"Less than two objects in secondary selection")
     return
  end
  if warn then
     model:warning(warn)
  end

-------- Now we have gathered the input data

--- Prepare for further processing:

  local off1 = {0}
  local len1 = {}
  for i,s in ipairs(path1) do
    offset = s[2]-s[1]
    len1[i] = offset:len()
    off1[i+1] = off1[i]+len1[i]
  end
  local extent1 = off1[#off1]

  local off2 = {0}
  local len2 = {}
  for i,s in ipairs(path2) do
    offset = s[2]-s[1]
    len2[i] = offset:len()
    off2[i+1] = off2[i]+len2[i]
  end
  local extent2 = off2[#off2]

--  draw elipses
  local seg1, seg2, line1, line2, dir1, dir2, p0, rectbbox
  local origx,origy,ang,mu,mv,mat,elli, ll,lr,ul,ur
  local freeshape, freeobject, free, clipshape
  local xanchor, offsetvector
  local shift, distance, halfwidth
  local freespace = {}
          
  for i,s in ipairs(path1) do
    seg1  = ipe.Segment(s[1],s[2])
    line1 = ipe.LineThrough(s[1],s[2])
    dir1 = (s[2]-s[1]):normalized()

    for j,t in ipairs(path2) do
      seg2 = ipe.Segment(t[1],t[2])
      if seg1:intersects(seg2) or
         seg1:distance(t[1])<=eps or
         seg1:distance(t[2])<=eps or
         seg2:distance(s[1])<=eps or
         seg2:distance(s[2])<=eps then -- otherwise: distance >eps
        line2 = ipe.LineThrough(t[1],t[2])
        dir2 = (t[2]-t[1]):normalized()
        p0 = line1:intersects(line2)
        rectbbox = ipe.Rect()
        for i,p in ipairs{s[1],s[2],t[1],t[2]} do rectbbox:add(p) end
        if not p0 or seg1:distance(p0) >
                  maxratio*(rectbbox:height()+rectbbox:width()) then
          -- draw two parallel lines, -> i.e., a rectangle
          if dir1..dir2>0 then
            mat = ipe.Matrix(0.5,-0.5,0.5,0.5,  off1[i],off2[j])
            xanchor = off1[i]
            offsetvector = s[1]-t[1]
            dir = (dir1+dir2)*0.5 
          else
            mat = ipe.Matrix(-0.5,-0.5,-0.5,0.5,  off1[i+1],off2[j])
            xanchor = off1[i+1]
            offsetvector = s[2]-t[1]
            dir = (dir2-dir1)*0.5 
          end
          shift = -offsetvector .. dir
          distance = offsetvector .. dir:orthogonal()
          halfwidth = math.sqrt(eps^2-distance^2)
	  ll = ipe.Vector(shift-halfwidth,0)
	  ul = ipe.Vector(shift+halfwidth,0)
	  lr = ipe.Vector(shift-halfwidth,len1[i]+len2[j])
	  ur = ipe.Vector(shift+halfwidth,len1[i]+len2[j])
          freeshape = {{type="curve",closed=true,
                       {type="segment",ll,lr},
                       {type="segment",lr,ur},
                       {type="segment",ur,ul}}}
          freeobject = ipe.Path(model.attributes,freeshape)
          freeobject:setMatrix(mat)
        else
          -- draw ellipse
          origx = off1[i] - ((s[1]-p0)..dir1)
          origy = off2[j] - ((t[1]-p0)..dir2) -- center of ellipse
          ang = dir1 .. dir2
          mu = eps/math.sqrt(2*(1-ang))
          mv = eps/math.sqrt(2*(1+ang))
          mat = ipe.Matrix( mu, mu, -mv, mv, origx, origy)
          freeshape = {{type="ellipse",mat}}
          freeobject = ipe.Path(model.attributes,freeshape)
        end 
        freeobject:set("pathmode","strokedfilled","black","white")

        -- now clip to rectangle
        free = ipe.Group{freeobject}
        ll = ipe.Vector(off1[i],off2[j])
        lr = ipe.Vector(off1[i+1],off2[j])
        ul = ipe.Vector(off1[i],off2[j+1])
        ur = ipe.Vector(off1[i+1],off2[j+1])
        clipshape = {{type="curve",closed=true,
                      {type="segment",ll,lr},
                      {type="segment",lr,ur},
                      {type="segment",ur,ul}}}
        free:setClip(clipshape)
        freespace[#freespace+1] = free

      end
    end
  end

--  draw grid

  local gridshape = {}
  for i,x in ipairs(off1) do
    local gridline = {type="curve",closed=false,
       {type="segment",ipe.Vector(x,0),ipe.Vector(x,extent2)}}
    gridshape[#gridshape+1] = gridline
  end
  for j,y in ipairs(off2) do
    local gridline = {type="curve",closed=false,
       {type="segment",ipe.Vector(0,y),ipe.Vector(extent1,y)}}
    gridshape[#gridshape+1] = gridline
  end
  local grid = ipe.Path(model.attributes, gridshape )
  grid:set("pathmode","stroked")

  local ll = ipe.Vector(0,0)
  local lr = ipe.Vector(extent1,0)
  local ul = ipe.Vector(0,extent2)
  local ur = ipe.Vector(extent1,extent2)
  local rectshape = {{type="curve",closed=true,
       {type="segment",ll,lr},
       {type="segment",lr,ur},
       {type="segment",ur,ul}}}
  local boundingrectancle = ipe.Path(model.attributes,rectshape)
  boundingrectancle:set("pathmode","filled",nil,"gray")

  local fsgroup = ipe.Group(freespace)
  local obj = ipe.Group{boundingrectancle,fsgroup,grid}
  model:creation("create free-space diagram", obj)

end

run = FSD