/ This function fits a line to a set of x y coordinate
/ using Least Squares technique

/ Given a list of x y coordinates in a VecPt2D object
/ The line function is
/                y = ax + b
/ Return a matrix containing values of a, b

/ Author : anonymous

/ ************************************************

  Function Matrix Coef = fitline(VecPt2D Vecxy)

/ ************************************************

  Coef = Matrix()

  tol = 0.0001

  x_min = Vecxy.min().x()
  x_max = Vecxy.max().x()

  if (x_max - x_min < tol)
    print
    print "Error in function (fitline) : vertical line found..."
    print
    return
  end

  n = Vecxy.size()

  Mrc = Vecxy.matrix()
  Lb = Mrc.getcol(1)
  A = Mrc.getcol(0).concat(Matrix(n,1,1))

  N = A.tsp()*A
  Y = A.tsp()*Lb
  X = N.inv()*Y

  Coef = X

  return

/ This function fits a parabola function to a set of x y coordinate
/ using Least Squares technique

/ Given a list of x y coordinates in a VecPt2D object
/ The parabola function is
/           y = ax^2 + bx + c
/ Return a matrix containing values of a, b, c

/ Author : anonymous

/ ************************************************

  Function Matrix Coef = fitparabola(VecPt2D Vecxy)

/ ************************************************

  n = Vecxy.size()

  mx = Vecxy.mean().x()

  A  = Matrix(n,3)
  Lb = Matrix(n,1)

  for i=0, n-1
    x = Vecxy(i).x()
    y = Vecxy(i).y()

/ below technique help increase accuracy when x is large

    A(i,0) = x*x/(mx*mx)
    A(i,1) = x/mx
    A(i,2) = 1

    Lb(i,0) = y
  end

  N = A.tsp()*A
  Y = A.tsp()*Lb
  X = N.inv()*Y

  X(0,0) = X(0,0)/(mx*mx)
  X(1,0) = X(1,0)/mx

  Coef = X

  return

 

/ This function is to determine plane parameter
/ using Least Squares technique

/ Given a list of x y z coordinates in a VecPt3D object
/ The plane equation is
/       z = Ax + By + C
/ Return a matrix containing values of A, B, and C

/ Author : anonymous


  Function Matrix pABC = fitplane(VecPt3D vecpt)

  Mxyz = vecpt.matrix()

  xyz_mean = vecpt.mean()
  mx = xyz_mean.x()
  my = xyz_mean.y()

  n = Mxyz.nrow()

/ mx and my help increase accuracy when x y are large

  y = Mxyz.getcol(2)

  A_x = Mxyz.getcol(0)/mx
  A_y = Mxyz.getcol(1)/my
  A_1 = Matrix(n,1,1)

  A = A_x.concat(A_y).concat(A_1)

  N = A.tsp()*A
  c = A.tsp()*y
  X = N.inv()*c

/ need to recover the original parameters A and B
  X(0,0) = X(0,0)/mx
  X(1,0) = X(1,0)/my

  pABC = X

  return

 

/ This function convert a set of GCP points (XYZ)
/     where XY is in a certain map projection system and Z is elevation
/ To a local 3D cartesian system (ENH), origin at lat0, lon0 and h0

/ input elevation are assumed to be ellipsoid height
/ P is the input map projection

/ Author : anonymous

/**********************************************************************************************************
/ function prototype:

  Function VecIdPt3D gcp_xyz = gcp2enh(VecIdPt3D gcp_utm, Projection P, double lat0, double lon0, double h0)

/ gcp_xyz = an output vector of IdPt3D in ENH system
/ gcp_utm = an input vector of IdPt3D in x y z (x y = map proj coord and z = elevation)
/ P       = input map projection
/ lat0    = latitude origin of the local ENH system
/ lon0    = longitude origin of the local ENH system
/ h0      = elevation origin of the local ENH system

/**********************************************************************************************************

  E = P.ellipsoid()

  n = gcp_utm.size()

/ convert input coordinate to a matrix
  mat = gcp_utm.matrix()

/ extract x y coordinate from matrix
  xy = mat.window(0,n-1,1,2)

  gcp_xy = xy.vecpt2d()

/ convert xy to lat lon coordinate
  gcp_geo = P.vecxy2geo(gcp_xy)

/ convert lat lon coordinate to a matrix
  mat_geo = gcp_geo.matrix()

/ put lat lon coord back to the original matrix
/ Now we have ID lat lon height
  mat.setwindow(0,1, mat_geo)

/ convert the matrix to a VecPt3D object
  gcp_geoh = mat.window(0,n-1,1,3).vecpt3d()

  gcp_XYZ = E.vecgeo2xyz(gcp_geoh)

  gcp_xyz = E.vecxyz2enh(gcp_XYZ,lat0, lon0, h0)

  mat_xyz = gcp_xyz.matrix()

  mat.setwindow(0,1,mat_xyz)

/ convert the matrix back to an VecIdPt3D obje3ct
  gcp_xyz = mat.vecidpt3d()

  return