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

/ PROGRAM AT_rc2xpyp

/ written by : Parallax
/ date : 12/29/2001

/ This program prepare data for Aerial Triangulation Bundle Adjustment
/ It reads row and column measurement of Each photo from an input file
/ converts them to photo coordinate and write to another file

/ Photo DOC file(s) must exist, otherwise an error is reported.
/ Output file must not exist, if exist new data will be appended to the file.

/ Photo ID must be an integer number
/ Photo, including IO data, must exist in the hard disk with the same name
/ as the photo ID (e.g. if ID = 12 then photo files are 12.raw, and 12.txt)

/ INPUT file has the following format
/ Photo ID
/ Pt row col
/ Pt row col
/ Pt row col
/ Pt row col
/ ....

/ Photo ID
/ Pt row col
/ Pt row col
/ Pt row col
/ Pt row col
/ ....


/ OUTPUT file has the following format
/ Photo ID
/ Pt x y
/ Pt x y
/ Pt x y
/ ....

/ Photo ID
/ Pt x y
/ Pt x y
/ Pt x y
/ ....

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

fname_inp = String()
fname = String()

input "input file name : " fname_inp
input "output file name : " fname

MAX_NO_PHOTO = 100

VecID = VecInt(MAX_NO_PHOTO)

VecData = VecStr()
VecData.load(fname_inp)

/ Search for Photo ID and store in VecID

n = VecData.size()

for i=0, n-1
line = VecData(i).ltrim().rtrim()

/ Neglect comment lines and empty lines
if line.len()==0 | line.left(1) == "/"
continue
end

no = line.no_word()
if no==1
ID = line.extract(1)(0).val()
VecID.pushback(ID)
end
end

/ remember the old fout name
old_fout = getfout()

set fout fname

n = VecID.size()

a = Photo()
v_rc = VecIdPt2d()

/ Load measuring data of each Photo into v_rc, using function "loadsch"
/ Then convert to photo coordinates and store in v_xpyp

for i=0, n-1


/ Load each Photo (ID is file name)
ID = VecID(i)
a.vload(int2str(ID))

v_rc.loadsch(fname_inp, ID)
v_xpyp = a.vidrc2xpyp(v_rc)

/ write v_xpyp to the output file
fprint ID
n1 = v_xpyp.size()
set width 15
set precision 3
for j=0, n1-1
fprint v_xpyp(j).ID();
fprint v_xpyp(j).x();
fprint v_xpyp(j).y()
end
fprint
fprint
print "read" ID

end

/ set fout back to the old name
set fout old_fout
 

/ This is to create a buffer image
/ result is an image of 0 and 1.
/ Area with 1 is those who have distance <= the distance "buf" from the object
/ the object is defined by pixel with value > 0, in the original image A

/ Author : anonymous

/ *****************************************************************************
  Function Image_uch result = buffer(Image_uch A, double buf)

/ A      = original image, containing object (pixel value > 0)
/ buf    = buffer distance
/ result = buffer zone image (contain value of 0 or 1)
/          Area with 1 is those who have distance <= "buf" from the object
/          the object is defined by pixel with value > 0, in the image A
/ *****************************************************************************

  result = A.blank()

  B = A.dist()

  vecxy = findimg(B<=buf)

  result.set(vecxy,1)

  return
 

/ This function creates a colorbar image, using the current color map

/ The return result is an Image_rgb object

/ Author : anonymous

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

  Function Image_rgb Img_out = colorbar(int m, int n)

/
/ m = no of row of each color
/ n = no of col of each color

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

  colormap = getcolormap()

  no_color = colormap.nrow()

  tmp      = matrix_flt()

  A        = matrix_flt()

  for i=0,no_color-1
    tmp.init(m,n,i)
    A = A.concatdown(tmp)
  end

  Img_out = A.rgb()

  return
 

/ This function is calculate grid convergence or grid declination of a   point (given by latitude and longitude).

/ Usage : call a = convergence(lat, lon, proj) ; where proj is Map Projection object

/ Author : chainman


  Function double conv = convergence(double lat, double lon, Projection proj)

  lat2 = lat + 0.001/3600

/ 0.001/3600 degree = 0.03 meters on ground

  xy1 = proj.geo2xy(lat, lon)
  xy2 = proj.geo2xy(lat2, lon)

  dx = xy2.x() - xy1.x()
  dy = xy2.y() - xy1.y()

  conv = atan2(-dx, dy)


  return
 

/ This is to create an AUX file for read into GeoMatic FreeView as a
/ Generic Binary (*.RAW)
/
/ On return, a file name with an extension .AUX will be created
/ Do make sure that the file fname.aux DOES NOT exist before calling this function
/ or else this function will append the content of AUX file to the existing file

/ Author : anonymous

/ *****************************************************************************
  Function doc2aux(String fname)

/ fname = file name of the Matrix object without an extension

/ data type can be uchar, short or float ONLY

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


  V_doc = VecStr()
  V_doc.load(fname)

  n = V_doc.size()

  V_doc_new = VecStr(n)

/ Take away comment lines and blank lines
/ and store in another VecStr

  for i = 0, n-1
    text = V_doc(i).ltrim().rtrim()
    if(text.left(1)=="/" | text.len()==0)
      continue
    else
      V_doc_new.pushback(text)
    end
  end

  word1 = V_doc_new(0).upper()

/ Check if this is the Document file of a Matrix object
  if(word1<>"MATRIX")
     print "error in function doc2aux : not a Matrix object data file"
     return
  end

  V_aux = VecStr(10)

/ Raw file in Geomatix FreeView only handle 8U,16S, 16U and 32R

  no_byte = round(V_doc_new(1).val())
  if(no_byte<>1 & no_byte<>2 & no_byte<>4)
    print "error in function doc2aux : no of byte of data must be 1 2 or 4 "
    return
  end

  if(no_byte==1)
    data_str = "8U"
  else
    if(no_byte==2)
      data_str = "16S"
    else
      data_str = "32R"
    end
  end

  no_row = round(V_doc_new(2).val())
  no_col = round(V_doc_new(3).val())

  V_aux.pushback("AuxilaryTarget: " + fname + ".raw")
  V_aux.pushback("RawDefinition: " + int2str(no_col) + " " + int2str(no_row) + " 1")
  V_aux.pushback("ChanDefinition-1: " + data_str + " 0 " + int2str(no_byte) + " " +   int2str(no_byte*no_col) + " Swapped")

/ memorize the previous FOUT file name
  fout_old = getfout()

  fout_new = fname + ".aux"
  set fout fout_new

  fprint V_aux(0)
  fprint V_aux(1)
  fprint V_aux(2)

/ recover the previous FOUT file name
  set fout fout_old

  return