External Function and Program Repository (F-J)

---

| HOME |   BACK   |

fitline.fun

/ 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

---

| HOME |   BACK   |

fitparabola.fun

/ 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

---

| HOME |   BACK   |

fitplane.fun

/ 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

---

| HOME |   BACK   |

gcp2enh.fun

/ 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

---

| HOME |   BACK   |

gen_test_image.prg

// This Program is to generate an image consisting of a regular grid of cross symbols

// set path "d:\doc_ku_class\photo_203312\lab_material\sony\"

set path "where your data are"

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

// Please change the following variable values

// BEFORE running the program

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

// n_row = no of rows of the resulting image

// n_col = no of columns of the resulting image

// spacing = spacing distance between each cross in horizontal and vertical

// tick_length = length of the cross symbol

// test_image = TIF file name of the output image

n_row = 900

n_col = 1600

spacing = 150

tick_length = 81

border = 50

a = Matrix_uch(n_row, n_col,255)

b = Matrix_uch(tick_length, tick_length, 255)

for ii = 0, tick_length-1

b(int(tick_length/2) , ii) = 0

end

for ii = 0, tick_length-1

b(ii, int(tick_length/2)) = 0

end

start_row = border-1

start_col = border-1

no_tick_hor = int((n_col-2*border)/spacing)

no_tick_ver = int((n_row-2*border)/spacing)

for i=start_row, start_row+no_tick_ver*spacing, spacing

for j=start_col, start_col+no_tick_hor*spacing, spacing

a.setsymb(i, j, b)

end

end

a.savetif("test_image")

---

| HOME |   BACK   |

Image_rgb2uch.fun

/ This function is to generate an Image_uch object from an Image_RGB object

/ The resulting image is a grey scale image whose brightness value is the average

/ of those of color image.

 / The class function "RGB" of Image_rgb almost does the same as this function does,
 / except, instead of generating an Image_uch object, it generates an Matrix_uch object.

/ Author : anonymous

Function Image_uch img_out = Image_rgb2uch(Image_rgb img)

no_row = img.nrow()

no_col = img.ncol()

R = img.r().float()

G = img.g().float()

B = img.b().float()

M = (R + G + B)/ 3

M = M.uchar()

pt1 = img.lower_left()

pt2 = img.upper_right()

img_out = Image_uch(no_row, no_col)

img_out.matrix() = M

img_out.lower_left() = pt1

img_out.upper_right() = pt2

return

---

| HOME |   BACK   |