include <tbset.h>

# t_vecfea-- compute the necessary parameters to register a 1-D  image  by
#            using common features of image data
#
# R.L. Williamson II, 29-Apr-92 Task derived from REGIST_SX in old SDAS
#
#   This (quite general) routine enables two vectors to be registered.
#   It needs two features-found-tables for each of the input science files.
#   As long as there is no reflection or distortion, the registration will
#   return the offset and magnification necessary to register the files.
#   A Linear Leastsquares algorithm is used.
#
#   This routine is based on the following:  if we have a list of
#   coordinates in the original vector in x and the positions of the
#   corresponding points in the secondary vector in x', then the two are
#   related by:
#
#                   x' = x * z + dx,
#
#   where dx is the offset and z the magnification factor between the two
#   vectors, which is simply a linear equation with z as the slope and dx as
#   the intercept.  Thus, these can be got by a straightforward linear
#   least squares fit to the x and x' positions which are input via a
#   Located Features File.

procedure t_vecfea()

char reftable[SZ_FNAME]                # reference input table name
char sectable[SZ_FNAME]                # second input table name
char outtable[SZ_FNAME]                # output table name
char colname[SZ_COLNAME]               # input column name
double parmval[4]                      # output parameters of registration
int mode, nterms                       # mode, number of terms for polyq
int rnfeat, snfeat,mxfea               # number of features in each table, max
int nfeat                              # number of features comparing
int i                                  # Counter
int cnt                                # Counter
double chisq                           # chi squared from result of fit
pointer sp
pointer rp                             # reference table descriptor
pointer secp                           # secondary table descriptor
pointer otp                            # output table pointer
pointer rcolptr, scolptr               # pointers to columns of input tables
pointer nullflag1                      # null pointer
pointer nullflag2                      # null pointer
pointer tempx1,tempx2                  # pointers to temporary arrays
pointer refx,secx                      # pointers to arrays of input coordinates
pointer coeff                          # scratch for fit coefficients
pointer sigmay                         # scratch array for sigmas
pointer out[4]                         # pointer for columns of output table
int tbpsta()                           #Functions
pointer tbtopn()

begin
        call smark (sp)

        #Get size of the reference group

        call clgstr ("reftable", reftable, SZ_FNAME)
        rp = tbtopn (reftable, READ_ONLY, 0)
        rnfeat = tbpsta(rp, TBL_NROWS)

        #Now get the secondary

        call clgstr ("sectable", sectable, SZ_FNAME)
        secp = tbtopn (sectable, READ_ONLY, 0)
        snfeat = tbpsta(secp, TBL_NROWS)

        nfeat = min (rnfeat, snfeat)
        mxfea = max (rnfeat,snfeat)

        # Read the name of the table column from the par file

        call clgstr ("column", colname, SZ_COLNAME) 

        # Read the data from the tables

        call salloc (nullflag1, mxfea, TY_BOOL)
        call salloc (nullflag2, mxfea, TY_BOOL)
        call salloc (refx, rnfeat, TY_DOUBLE)
        call salloc (secx, snfeat, TY_DOUBLE)
        call salloc (tempx1, rnfeat, TY_DOUBLE)
        call salloc (tempx2, snfeat, TY_DOUBLE)
        
        call tbcfnd (rp, colname, rcolptr, 1)
        call tbcgtd (rp, rcolptr, Memd[tempx1], Memb[nullflag1], 1, rnfeat)

        # Weed out nulls
        cnt=0
        do i=0,rnfeat-1 {
            if(!Memb[nullflag1+i]) {
            Memd[refx+cnt]=Memd[tempx1+i]
            cnt=cnt+1
            }
        }
        rnfeat=cnt

        call tbcfnd (secp, colname, scolptr, 1)
        call tbcgtd (secp, scolptr, Memd[tempx2], Memb[nullflag2], 1, snfeat)

        # Weed out nulls
        cnt=0
        do i=0,rnfeat-1 {
            if(!Memb[nullflag2+i]) {
            Memd[secx+cnt]=Memd[tempx2+i]
            cnt=cnt+1
            }
        }
        snfeat=cnt

        # Close input tables
        call tbtclo (rp)
        call tbtclo (secp)

        # Allocate some memory for fit coefficients, sigma array

        call salloc (coeff, 2, TY_DOUBLE)
        call salloc (sigmay,mxfea,TY_REAL)

        # Set mode and number of terms for polylq
        mode = 0
        nterms = 2

        # Call the POLYLQ to do the calculations

        call polylq( Memd[secx], Memd[refx],
                                Memr[sigmay], nfeat, nterms, mode,
                                Memd[coeff], chisq)

           
        # Place Values into array
        parmval[1] = Memd[coeff]
        parmval[2] = Memd[coeff+1]
        parmval[3] = Memd[coeff]
        parmval[4] = Memd[coeff+1]

        # Write output table
        
        call clgstr ("outtable", outtable, SZ_FNAME)

        otp = tbtopn (outtable, NEW_FILE, 0)

        call salloc (out,4,TY_POINTER)        
        call tbcdef (otp, out[1], "coeff1", "", "", TY_DOUBLE, 1, 1)
        call tbcdef (otp, out[2], "coeff2", "", "", TY_DOUBLE, 1, 1)
        call tbcdef (otp, out[3], "deltax", "", "", TY_DOUBLE, 1, 1)
        call tbcdef (otp, out[4], "xmag", "", "", TY_DOUBLE, 1, 1)

        call tbtcre (otp)
        
        # Dump the parameters to the table

        call tbrptd (otp, out, parmval, 4,1)

        call tbtclo (otp)

        call sfree (sp)

 end