#################################
### COPYRIGHT: Anders Lyhagen ###
#################################


module CAUL_FAQS
module QuadGrid
  
  #########################
  ########## TRIS #########
  #########################

  class Tris
    attr_reader :kind

    def initialize(vs, fs, es, enforce_kind)
      
      #ENFORCE = 0, 1, 2, nil <--default
      is_flat = flat?(vs)
      
      if(is_flat && (enforce_kind == 2 || enforce_kind == nil))   
        @kind = 2
        @epsilon = 0.00000001
        @ps = vs.map { |v| v.position.clone }
        @tr_to = getTransformation(fs[0].normal)
        @tr_from = @tr_to.inverse
        @v0 = (@ps[1] - @ps[0]) 
        @v1 = (@ps[2] - @ps[3])
      else
        @ts = []
        e = fs.length > 1 ? fs[0].edges & fs[1].edges : []
        @kind = e.length > 0 ? (e[0].used_by?(vs[0]) ? 1 : 0) : 0
        
        if is_flat && enforce_kind != nil && enforce_kind != @kind
          @kind = enforce_kind
        end
        
        f0 = fs[0]
        f1 = fs.length > 1 ? fs[1] : fs[0]
    
        if @kind == 0
          @ts << Tri.new(vs[0], es[0], 1, es[3], 1, f0)
          @ts << Tri.new(vs[2], es[2], -1, es[1], -1, f1)
        end
        if @kind == 1
          @ts << Tri.new(vs[1], es[0], -1, es[1], 1, f0)
          @ts << Tri.new(vs[3], es[2], 1, es[3], -1, f1)
        end
      end  
    end

    def flat?(vs)
      ps = vs.map { |v| v.position }
      pl = Geom::fit_plane_to_points ps
      return ps.all? { |p| p.on_plane? pl }
    end
  
    #Compute transformation for flat quad
    def getTransformation(n)
      #1) align with z axis
      tr1 = Matrix.getAlignmentMatrix(n, Geom::Vector3d.new(0, 0, 1))
      @ps.each { |p| p.transform!(tr1) }
    
      #2) translate to origo
      tr2 = Geom::Transformation.translation(Geom::Vector3d.new(-@ps[0].x, -@ps[0].y, -@ps[0].z))
      @ps.each { |p| p.transform!(tr2) }
    
      #3) rotate so that u coincides with x
      v0 = @ps[1] - @ps[0]
      v1 = Geom::Vector3d.new(1, 0, 0)
      v2 = Geom::Vector3d.new(0, 1, 0)
      a = v0.angle_between(v1)
      a = v2.dot(v0) < 0 ? a : - a
      tr3 = Geom::Transformation.rotation(@ps[0], [0, 0, 1], a)
      @ps.each { |p| p.transform!(tr3) }
      return tr3 * tr2 * tr1
    end 
    
    ###### FLAT COMPUTATION ####
    def getPointFlat(u, v)
      x = @v0.x * u + (@ps[3].x + u * (@v1.x - @v0.x))  * v
      y = (@ps[3].y + @v1.y * u) * v
      pp = Geom::Point3d.new(x, y, 0)
      return pp.transform!(@tr_from)
    end
    
    def getUVFlat(p)
      pp = p.transform(@tr_to)
      a = @v1.y * @v0.x
      b = ((pp.y * (@v1.x - @v0.x)) + @ps[3].y * @v0.x - @v1.y * pp.x)
      c = @ps[3].x * pp.y - @ps[3].y * pp.x 
       
      d1 = -b
      d2 = b * b - 4 * a * c
      d3 = 2 * a
      
      u = 0
      if(d2 >= 0 && d3.abs > 0.0000001)
        d2 = Math::sqrt(d2)
        u = (d1 - d2) / d3
      elsif d3.abs <= 0.0000001 #linear system
        u = pp.x / @v0.length
      end
      
      #compute the v coordinate
      v2x = @ps[3].x + @v1.x * u - @v0.x * u; v2y = @ps[3].y + @v1.y * u
      v3x = pp.x - @v0.x * u; v3y = pp.y
      l0 = Math::sqrt(v2x * v2x + v2y * v2y)
      l1 = Math::sqrt(v3x * v3x + v3y * v3y)
      v = l0 <= 0 ? 0 : l1 / l0
      return truncate(u), truncate(v)
    end
    
    def truncate(x)
      x = 1.0 if (x - 1).abs < @epsilon
      x = 0.0 if x.abs < @epsilon
      return x
    end
  
  ###########################
  #### PUBLIC INTERFACE #####
  ###########################
  
    def getUV(p)
      if @kind == 2
        getUVFlat(p)
      else
        u0, v0 = @ts[0].pointOnTri?(p) ? @ts[0].getUV(p) : @ts[1].getUV(p)
        return u0, v0
      end
    end
    
    def getPoint(u, v)
      if @kind == 2
        getPointFlat(u, v)
      else
        p0 = @ts[0].onTri?(u, v) ? @ts[0].getPoint(u, v) : @ts[1].getPoint(u, v)
        return p0
      end
    end
  end

  class Tri
  
    def initialize(v, eu, u_dir, ev, v_dir, f)
      @v = v
      @u_dir = u_dir
      @v_dir = v_dir
      @u_ve = eu.other_vertex(@v).position - @v.position
      @v_ve = ev.other_vertex(@v).position - @v.position
      @u_len_inv = 1.0 / eu.other_vertex(@v).position.distance(@v.position)
      @v_len_inv = 1.0 / ev.other_vertex(@v).position.distance(@v.position) 
      @f = f
      @fold_ve = ev.other_vertex(@v).position - eu.other_vertex(@v).position 
      @up = eu.other_vertex(@v).position
    end
    
    def onTri?(u, v)
      return (@u_dir > 0 ? u : 1 - u) + (@v_dir > 0 ? v : 1 - v) <= 1
    end
    
    #determine whether a given point is on the tri (provided it's on the face)
    def pointOnTri?(p)
      return @fold_ve.cross(p - @up).dot(@f.normal) * @u_dir <= 0
    end
    
    #the point is assumed to be on the tri
    def getPoint(u, v)
      t0 = Geom::Point3d.new(@v.position.x, @v.position.y, @v.position.z)
      a = @u_dir > 0 ? u : 1 - u
      b = @v_dir > 0 ? v : 1 - v
      t0.offset!([@u_ve.x * a + @v_ve.x * b, @u_ve.y * a + @v_ve.y * b, @u_ve.z * a + @v_ve.z * b])
      return t0
    end
    
    #the [x, y, z] point is assumed to be on the tri
    def getUV(p)
      pu = Geom::intersect_line_line([p, @v_ve], [@v.position, @u_ve])
      pv = Geom::intersect_line_line([p, @u_ve], [@v.position, @v_ve])
      return nil,nil if pu == nil || pv == nil
      u = pu.distance(@v.position) * @u_len_inv
      v = pv.distance(@v.position) * @v_len_inv
      #correct directions if necesarry
      u = @u_dir > 0 ? u : 1 - u
      v = @v_dir > 0 ? v : 1 - v
      return u, v
    end
  end

  
  #########################
  ######### QUAD ##########
  #########################
  
  class Quad
    
    attr_reader :fs, :es, :vs, :ns, :kind
    
    def initialize(fs, loop, enforce_kind)
      #settings
      @normal_mode = 1 #0 = MWA, 1 = MWE
      @adjust_mode = 0 #0 = off, 1 = on
      #variables
      @fs = fs
      @es = []
      @vs = []
      @ns = [] #vertex normals...

      loop.each { |a| @vs << a[0]; @es << a[1] }
      @tris = Tris.new(@vs, @fs, @es, enforce_kind)
      @kind = @tris.kind
    end
    
  
    def computeNormals
      @vs.each { |v| @ns << getVertexNormal(v) }   
    end
    
    def getVertexNormal(v)
      n = Geom::Vector3d.new(0, 0, 0)
      v.faces.each { |f| 
        a =  @normal_mode == 0 ? getVertexAngle(v, f) : 1.0
        n.x= n.x + a * f.normal.x; n.y= n.y + a * f.normal.y; n.z= n.z + a * f.normal.z 
      }
      n.normalize!
      if @adjust_mode == 1
        div = 0
        v.faces.each { |f| div += n.dot(f.normal) }
        div = div / v.faces.length
        n.length = n.length / div
      end
      return n
    end
    
    def getVertexAngle(v, f)
      es = []
      f.edges.each { |e| es << e if e.used_by?(v) }
      v1 = es[0].vertices[0] == v ? es[0].vertices[1] : es[0].vertices[0]
      v2 = es[1].vertices[0] == v ? es[1].vertices[1] : es[1].vertices[0]
      return (v1.position - v.position).angle_between(v2.position - v.position)
    end
   
    def setGlobalUV(min_u, max_u, min_v, max_v)
      @min_u = min_u
      @max_u = max_u
      @min_v = min_v
      @max_v = max_v
      @u_span_inv = 1.0 / (max_u - min_u)
      @v_span_inv = 1.0 / (max_v - min_v)
    end    
    
    def getOffset(p, n, d)
      p0 = p.clone
      p0.x = p0.x + n.x * d
      p0.y = p0.y + n.y * d
      p0.z = p0.z + n.z * d
      return p0
    end  

    #### INTERFACE ####
    
    def getPoint(u_global, v_global)
      ul = (u_global - @min_u) * @u_span_inv
      vl = (v_global - @min_v) * @v_span_inv   
      return @tris.getPoint(ul, vl)
    end
    
    def getLocalUV(u_global, v_global)
      ul = (u_global - @min_u) * @u_span_inv
      vl = (v_global - @min_v) * @v_span_inv
      return ul, vl
    end
    
    def getGlobalUV2(u_local, v_local)
      ug = (u_local / @u_span_inv) + @min_u
      vg = (v_local / @v_span_inv) + @min_v
      return ug, vg
    end
    
    def getGlobalUV(p, f)
      ul, vl = @tris.getUV(p)
      return nil, nil if ul == nil
      u_global = (ul / @u_span_inv) + @min_u 
      v_global = (vl / @v_span_inv) + @min_v 
      return u_global, v_global
    end
    
    def getNormal(u_global, v_global)
      ul = (u_global - @min_u) * @u_span_inv
      vl = (v_global - @min_v) * @v_span_inv 
      return getLocalNormal(ul, vl)
    end
    
    def getLocalNormal(u, v)
      return getLocalNormalInter(u, v)
    end
    
    #observe, the returned normal should not be expected to be normalized..
    def getLocalNormalInter(u, v)
      ui = 1 - u;
      vi = 1 - v;
      n0 = [@ns[0].x * ui + @ns[1].x * u, @ns[0].y * ui + @ns[1].y * u, @ns[0].z * ui + @ns[1].z * u]  
      n1 = [@ns[1].x * vi + @ns[2].x * v, @ns[1].y * vi + @ns[2].y * v, @ns[1].z * vi + @ns[2].z * v]
      n2 = [@ns[3].x * ui + @ns[2].x * u, @ns[3].y * ui + @ns[2].y * u, @ns[3].z * ui + @ns[2].z * u]
      n3 = [@ns[0].x * vi + @ns[3].x * v, @ns[0].y * vi + @ns[3].y * v, @ns[0].z * vi + @ns[3].z * v]
      nx = (n0[0] * vi  + n1[0] * u + n2[0] * v + n3[0] * ui) * 0.5
      ny = (n0[1] * vi  + n1[1] * u + n2[1] * v + n3[1] * ui) * 0.5
      nz = (n0[2] * vi  + n1[2] * u + n2[2] * v + n3[2] * ui) * 0.5
      return Geom::Vector3d.new(nx, ny, nz)
    end
    
    def getOffsetCell(d, tr)
      ps = []
      (0..3).each { |i| ps << getOffset(@vs[i].position, @ns[i], d).transform!(tr) } 
      
      fs = []
      if @fs.length == 1 
        fs << [ps[0], ps[3], ps[2], ps[1]]
      elsif @kind == 0 # / -case
        fs << [ps[0], ps[3], ps[1]]
        fs << [ps[1], ps[3], ps[2]]
      else
        fs << [ps[0], ps[3], ps[2]]
        fs << [ps[0], ps[2], ps[1]]
      end
      fold = @fs.length == 2 ? (@kind == 0 ? [ps[1], ps[3]] : [ps[0], ps[2]]) : nil
      return ps, fs, fold
    end
    
    def getFoldPoints
      return nil if @fs.length == 1
      return @kind == 0 ? [@vs[1].position.clone, @vs[3].position.clone] : 
                          [@vs[0].position.clone, @vs[2].position.clone]
    end
    
    def getFoldPointsUV
      return nil if @fs.length == 1
      return @kind == 0 ? [@min_u, @max_v, @max_u, @min_v] : 
                          [@min_u, @min_v, @max_u, @max_v]    
    end
           
  end

  #########################
  ######### GRID ##########
  #########################

  class Grid
    
    attr_reader :u_map, :v_map, :grid
  
    #grid is a 2D array of cells
    def initialize(grid)
      @grid = grid
      u_path = getUPath()
      v_path = getVPath()
      @du = 1.0 / u_path.length
      @dv = 1.0 / v_path.length
      @u_map = getEquiUVMap(u_path)
      @v_map = getEquiUVMap(v_path)
      
      (0..@grid.length-1).each { |u|
        (0..@grid[0].length-1).each { |v|
          @grid[u][v].setGlobalUV(@u_map[u], @u_map[u + 1], @v_map[v], @v_map[v + 1])
      }}
    
      #compute the grid normals
      (0..@grid.length-1).each { |u| 
          (0..@grid[0].length-1).each { |v| 
            @grid[u][v].computeNormals
      }}
      
      #Create a face hash for ray -> xyz -> uv computations
      @face_hash = FaceHash.new(getFaceHash())
      @d_factor = 1
    end
    
    def getUPath()
      return (0..@grid.length - 1).inject([]) { |p, i| p << [@grid[i][0].vs[0], @grid[i][0].es[3]]}
    end
    
    def getVPath()
      return @grid[0].map { |quad| [quad.vs[0], quad.es[0]] }
    end
       
    def getEquiUVMap(p)
      dx = 1.0 / p.length
      m = (0..p.length - 1).inject([]) { |a, i| a << dx * i }
      m << 1.0
      return m
    end
     
    def getQuad(a, mode)
      map = mode == 0 ? @u_map : @v_map
      return (map.length-2) if a == 1
      ind = (a * (map.length - 1)).floor
      return -1 if ind < 0 || ind > map.length - 2
      return ind
    end

    ####################
    ##### INTERFACE ####
    ####################
    ## UV -> XYZ
    def getPointOnGrid(u, v)
      return @grid[getQuad(u, 0)][getQuad(v, 1)].getPoint(u, v)
    end
    
    def getPoint(u, v, d)
      nd = @d_factor * d
      q = @grid[getQuad(u, 0)][getQuad(v, 1)]
      p = q.getPoint(u, v)
      n = q.getNormal(u, v)
      n.set!(n.x * nd, n.y * nd, n.z * nd)
      p.offset!(n)
      return p
    end
    
    def getNormal(u, v)
      q = @grid[getQuad(u, 0)][getQuad(v, 1)]
      n = q.getNormal(u, v)
      return n
    end
    
    #Ray -> XYZ -> UV-d
    def getUVD(ray)
      #get point and face from faceHash
      p, f, u_ind, v_ind = @face_hash.getPoint(ray)
      return nil if p == nil
      u, v = @grid[u_ind][v_ind].getGlobalUV(p, f)
      d = (p - ray[0]).dot(f.normal) < 0 ? ray[0].distance(p) : -ray[0].distance(p)
      return [u, v, d]
    end
        
    def getXYZ(ray)
      p, f, u_ind, v_ind = @face_hash.getPoint(ray)
      return p
    end
    
    def getFolds
      folds = []
      (0..@grid.length-1).each { |u| 
          (0..@grid[0].length-1).each { |v| 
            folds << @grid[u][v].getFoldPoints unless @grid[u][v].fs.length == 1 
      }}
      return folds
    end
    
    def getFoldsUV
      folds = []
      (0..@grid.length-1).each { |u| 
          (0..@grid[0].length-1).each { |v| 
            folds << @grid[u][v].getFoldPointsUV unless @grid[u][v].fs.length == 1 
      }}
      return folds
    end
    
    def getFaceHash()
      fh = {}
      (0..@grid.length-1).each { |u|
        (0..@grid[0].length-1).each { |v|
          fs = @grid[u][v].fs
          fs.each { |f| fh[f] = [u, v] }
      }}
      return fh
    end
   
    #create an offset mesh at distance d
    def getOffsetMesh(d, tr)
      mesh = Geom::PolygonMesh.new
      folds = []
      (0..@grid.length-1).each { |u|
        (0..@grid[0].length-1).each { |v|
          ps, fs, fold = @grid[u][v].getOffsetCell(d, tr)
          folds << fold if fold != nil
          ps.each { |p| mesh.add_point p }
          fs.each { |f| mesh.add_polygon f }
      }}
     
      p0 = getPoint(0, 0, d).transform(tr)
      p1 = getPoint(1, 0, d).transform(tr)
      return mesh, p0, p1, folds
    end
    
    def to_s
      return 'u: '+grid.length.to_s+'  '+'v: '+grid[0].length.to_s
    end
    
  end
  
end
end
