gpsman 6.4.4.2-2 / usr / lib / gpsman / maptransf.tcl

#
# This file is part of:
#
#  gpsman --- GPS Manager: a manager for GPS receiver data
#
# Copyright (c) 1998-2013 Miguel Filgueiras migfilg@t-online.de
#
#    This program is free software; you can redistribute it and/or modify
#      it under the terms of the GNU General Public License as published by
#      the Free Software Foundation; either version 3 of the License, or
#      (at your option) any later version.
#
#      This program is distributed in the hope that it will be useful,
#      but WITHOUT ANY WARRANTY; without even the implied warranty of
#      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#      GNU General Public License for more details.
#
#      You should have received a copy of the GNU General Public License
#      along with this program.
#
#  File: maptransf.tcl
#  Last change:  6 October 2013
#

## Uses information kindly supplied by
#     Jose Alberto Goncalves, Universidade do Porto
#     J. B. Mehl
##

## transformations for geo-referencing of images
#
# (xt,yt) terrain coordinates (metre)
# (xm,-ym) map coordinates (pixel, ym grows downwards!)
#
# affine transformation, 6 parameters
#  (xt,yt) = [ aij ] x (xm,ym) + (e,f)
#   where  a11 = a   a12 = b   a21 = c   a22 = d
#
# affine conformal transformation, 4 parameters
#  (xt,yt) = lambda [ aij ] x (xm,ym) + (e,f)
#   where  a11 = cos a   a12 = -sin a   a21 = sin a   a22 = cos a
#     rotation angle: a    scaling factor: lambda
#
# affine conformal no rotation transformation, 3 parameters
#  (xt,yt) = lambda (xm,ym) + (e,f)
##

## when adding new transformation procedures the variables MAPTRANSFNPTS and
#   MAPKNOWNTRANSFS must be changed in metadata.tcl
# see also relevant procedures in map.tcl, command_parse.tcl, command.tcl

     # indices of MTData array used for each transformation
array set MAPTRANSFDATA {
    NoRot      {lambda e f}
    Affine     {a b c d e f det k2m1 k4m3}
    AffineConf {a b e f det k2m1 k4m3}
}

proc MapTransfIs {transf} {
    # set global variables so that map transformation is $transf
    #  $transf in $MAPKNOWNTRANSFS
    global MapTransf MapTransfTitle TXT

    set MapTransf $transf ; set MapTransfTitle $TXT(TRNSF$transf)
    return
}

## transformation procs

proc MapInitNoRotTransf {scale xt0 yt0 xm0 ym0} {
    # compute parameters of affine conformal transformation with no rotation
    #  $scale is in metre/pixel
    #  $xt0,$yt0: terrain coordinates of the map point at $xm0,$ym0 pixel
    global MTData

    MapTransfIs NoRot
    # to avoid integer divisions
    set MTData(lambda) [expr $scale*1.0]
    set MTData(e) [expr $xt0-$scale*$xm0]
    set MTData(f) [expr $yt0+$scale*$ym0]
    return 0
}

proc InitNoRotTransf {} {
    # compute representation of affine conformal transformation with no
    #  rotation transformation from 2 WPs
    #  $MapLoadWPs is list of indices of relevant WPs
    #  $MapLoadPos($n,x), $MapLoadPos($n,y)  $n in 0..2 give pixel coordinates
    # set MapScale to appropriate value in metre/pixel
    # return 0 if failed to solve equations
    global MapLoadPos MapScale MTData

    MapTransfIs NoRot
    foreach p "0 1" tcs [MapGeoRefPoints 2] {
	foreach "xt$p yt$p" $tcs {}
    }
    # distance between points in the terrain (metre)
    set dxt [expr $xt0-$xt1] ; set dyt [expr $yt0-$yt1]
    # must have both multiplications by 1.0 !!!
    set lt [expr sqrt(1.0*$dxt*$dxt+1.0*$dyt*$dyt)]
    # distance between points in the map (pixel)
    set dxm [expr $MapLoadPos(0,x)-$MapLoadPos(1,x)]
    set dym [expr $MapLoadPos(0,y)-$MapLoadPos(1,y)]
    # must have both multiplications by 1.0 !!!
    set lm [expr sqrt(1.0*$dxm*$dxm+1.0*$dym*$dym)]
    # scale
    if { [catch {set scale [expr 1.0*$lt/$lm]}] } { return 0 }
    set MapScale $scale
    set MTData(lambda) [expr $scale*1.0]
    set MTData(e) [expr $xt0-$scale*$MapLoadPos(0,x)]
    set MTData(f) [expr $yt0+$scale*$MapLoadPos(0,y)]
    return 1
}

proc MapApplyNoRotTransf {xt yt} {
    # apply affine conformal transformation with no rotation
    global MTData

    set s $MTData(lambda)
    return [list [expr ($xt-$MTData(e))/$s] [expr ($MTData(f)-$yt)/$s]]
}

proc MapInvertNoRotTransf {xm ym} {
    # invert affine conformal transformation with no rotation
    global MTData

    set s $MTData(lambda)
    return [list [expr $s*$xm+$MTData(e)] [expr $MTData(f)-$s*$ym]]
}

proc MapNewScaleNoRotTransf {scale} {
    # set transformation parameters after change in map scale
    # return 1 if possible
    global MTData

    set MTData(lambda) [expr $scale*1.0]
    return 1
}

proc MapInitAffineTransf {args} {
    # compute representation of affine transformation from 3 points
    #  $args==""   points are WPs, with $MapLoadWPs the list of
    #          indices of relevant WPs, used by proc MapGeoRefPoints
    #          that sets MapLoadPos(_,_) to the pixel coordinates
    #       ==list of projected coordinates, MapLoadPos(_,_) was already set
    #  $MapLoadPos($n,x), $MapLoadPos($n,y)  $n in 0..2 give pixel coordinates
    # set MapScale to appropriate value in metre/pixel
    # return 0 if failed to solve equations
    global MapLoadPos MapScale MTData Mat Rx

    MapTransfIs Affine
    if { $args == "" } {
	if { [set tcs [MapGeoRefPoints 3]] == -1 } { return 0 }
    } else { set tcs [lindex $args 0] }
    foreach ps [list "a b e" "c d f"] d "0 1" {
	foreach i "0 1 2" {
	    set Mat($i,0) $MapLoadPos($i,x)
	    set Mat($i,1) $MapLoadPos($i,y)
	    set Mat($i,2) 1
	    set Mat($i,3) [lindex [lindex $tcs $i] $d]
	}
	if { [GaussReduce 3] != 1 } { return 0 }
	foreach p $ps i "0 1 2" {
	    set MTData($p) $Mat($Rx($i),3)
	}
    }
    # the following parameters simplify computations
    set MTData(det) [expr $MTData(a)*$MTData(d)-$MTData(b)*$MTData(c)]
    set MTData(k2m1) [expr $MTData(b)*$MTData(f)-$MTData(d)*$MTData(e)]
    set MTData(k4m3) [expr $MTData(c)*$MTData(e)-$MTData(a)*$MTData(f)]
    # scale along the xm-axis when variation of ym=0
    set MapScale [expr abs($MTData(a))]
    return 1
}

proc MapApplyAffineTransf {xt yt} {
    # apply affine transformation
    global MTData

    set x [expr ($MTData(d)*$xt-$MTData(b)*$yt+$MTData(k2m1))/$MTData(det)]
    set y [expr ($MTData(a)*$yt-$MTData(c)*$xt+$MTData(k4m3))/$MTData(det)]
    return [list $x $y]
}

proc MapInvertAffineTransf {xm ym} {
    # invert affine transformation
    global MTData

    set xt [expr $MTData(a)*$xm+$MTData(b)*$ym+$MTData(e)]
    set yt [expr $MTData(c)*$xm+$MTData(d)*$ym+$MTData(f)]
    return [list $xt $yt]
}

proc MapNewScaleAffineTransf {scale} {
    # set transformation parameters after change in map scale
    # return 1 if possible

    return 0
}

proc MapInitAffineConfTransf {args} {
    # compute representation of affine conformal transformation from 2 points
    #  $args==""   points are WPs, with $MapLoadWPs the list of
    #          indices of relevant WPs, used by proc MapGeoRefPoints
    #          that sets MapLoadPos(_,_) to the pixel coordinates
    #       ==list of projected coordinates, MapLoadPos(_,_) was already set
    #  $MapLoadPos($n,x), $MapLoadPos($n,y)  $n in 0..1 give pixel coordinates
    # set MapScale to appropriate value in metre/pixel
    # return 0 if failed to solve equations
    global MapLoadPos MapScale MTData Mat Rx

    MapTransfIs AffineConf
    catch {unset MTData}
    if { $args == "" } {
	if { [set tcs [MapGeoRefPoints 2]] == -1 } { return 0 }
    } else { set tcs [lindex $args 0] }
    foreach e "0 2" i "0 1" {
	set xyt [lindex $tcs $i]
	set Mat($e,0) $MapLoadPos($i,x)
	set Mat($e,1) [expr -$MapLoadPos($i,y)]
	set Mat($e,2) 1 ; set Mat($e,3) 0
	set Mat($e,4) [lindex $xyt 0]
	incr e
	set Mat($e,0) [expr -$MapLoadPos($i,y)]
	set Mat($e,1) [expr -$MapLoadPos($i,x)]
	set Mat($e,2) 0 ; set Mat($e,3) 1
	set Mat($e,4) [lindex $xyt 1]
    }
    if { [GaussReduce 4] != 1 } { return 0 }
    foreach p "a b e f" i "0 1 2 3" {
	set MTData($p) $Mat($Rx($i),4)
    }
    # the following parameters make calculations easier
    set MTData(det) [expr $MTData(a)*$MTData(a)+$MTData(b)*$MTData(b)]
    set MTData(k2m1) [expr $MTData(b)*$MTData(f)-$MTData(a)*$MTData(e)]
    set MTData(k4m3) [expr $MTData(b)*$MTData(e)+$MTData(a)*$MTData(f)]
    # scale along the xm-axis when variation of ym=0
    set MapScale [expr abs($MTData(a))]
    return 1
}

proc MapApplyAffineConfTransf {xt yt} {
    # apply affine conformal transformation
    global MTData

    set xm [expr ($MTData(a)*$xt-$MTData(b)*$yt+$MTData(k2m1))/$MTData(det)]
    set ym [expr -($MTData(a)*$yt+$MTData(b)*$xt-$MTData(k4m3))/$MTData(det)]
    return [list $xm $ym]
}

proc MapInvertAffineConfTransf {xm ym} {
    # invert affine conformal transformation
    global MTData

    set xt [expr $MTData(a)*$xm-$MTData(b)*$ym+$MTData(e)]
    set yt [expr -$MTData(b)*$xm-$MTData(a)*$ym+$MTData(f)]
    return [list $xt $yt]
}

proc MapNewScaleAffineConfTransf {scale} {
    # set transformation parameters after change in map scale
    # return 1 if possible

    return 0
}

proc MapAffineParams {} {
    # convert parameters of current map transformation to corresponding
    #  parameters of an affine transformation
    # return list of a parameter name followed by its value, where the
    #  parameter names are at least the following ones  a, b, c, d, e, f
    #  as described above for the affine transformation
    global MapTransf MTData

parray MTData
    switch $MapTransf {
	Affine { return [array get MTData] }
	AffineConf {
	    set ps {}
	    foreach p {a b c d e f} cp {a b b a e f} s {1 -1 -1 -1 1 1} {
		lappend ps $p [expr $s*$MTData($cp)]
puts "p=$p, cp=$cp v=$MTData($cp) e=[expr $s*$MTData($cp)]"
	    }
	}
	NoRot {
	    set lambda $MTData(lambda)
	    return [list a $lambda b 0 c 0 d $lambda \
			e $MTData(e) f $MTData(f)]
	}
    }
    return $ps
}


## solving a linear system of equations nxn by Gauss-Jordan elimination
# code adopted from the Slopes Algorithm implementation in C
# by M Filgueiras and A P Tomas / Universidade do Porto / 1996, 1997
#

proc GaussNullFirst {k n} {
    # check first row with non-null element in $k,$k using Rx, Cx and
    #  exchange rows if needs be
    #  $n is dimension of matrix
    # return 0 if non-null element found
    global Rx Cx Mat

    for { set i [expr $k+1] } { $i < $n } { incr i } {
	if { $Mat($Rx($i),$Cx($k)) != 0 } {
	    set l $Rx($k) ; set Rx($k) $Rx($i) ; set Rx($i) $l
	    return 0
	}
    }
    return 1
}

proc GaussElim {i k n p} {
    # eliminate $i,$k element on ?x($n+1) matrix, pivot $p at $k,$k,
    #  using Rx, Cx; the $i,$k element is assumed to be non-null
    global Rx Cx Mat

    set ii $Rx($i) ; set ik $Rx($k)
    set m [expr 1.0*$Mat($ii,$Cx($k))/$p]
    for { set j [expr $k+1] } { $j <= $n } { incr j } {
	set jj $Cx($j)
	set Mat($ii,$jj) [expr $Mat($ii,$jj)-$m*$Mat($ik,$jj)]
    }
    set Mat($ii,$Cx($k)) 0
    return
}

proc GaussSubelim {k n} {
    # eliminate below $k on nx(n+1) matrix, using Rx, Cx[]
    global Rx Cx Mat

    set ck $Cx($k)
    set p $Mat($Rx($k),$ck)
    for { set i [expr $k+1] } { $i < $n } { incr i } {
	if { $Mat($Rx($i),$ck) != 0 } {
	    GaussElim $i $k $n $p
	}
    }
    return
}

proc GaussSupraelim {i n} {
    # eliminate above $i on _x($n1) matrix, using Rx, Cx
    global Rx Cx Mat

    set ci $Cx($i)
    set p $Mat($Rx($i),$ci)
    for { set a 0 } { $a < $i } { incr a } {
	if { $Mat($Rx($a),$ci) != 0 } {
	    GaussElim $a $i $n $p
	}
    }
    return
}

proc GaussReduce {n} {
    # reduction of $nx($n+1) matrix, using Rx, Cx to index rows and columns
    # indices start from 0
    # values in global array $Mat are changed by this procedure
    # return 1 if there are is only one solution
    # solutions to be retrived from $Mat($Rx($i),$n) for each $i from 0 to $n-1
    global Rx Cx Mat

    for { set i 0 } { $i < $n } { incr i } {
	set Rx($i) $i ; set Cx($i) $i
    }
    set Cx($n) $n
    for { set i 0 } { $i < $n } { incr i } {
	if { $Mat($Rx($i),$Cx($i))==0 && [GaussNullFirst $i $n] } { return 0 }
	GaussSubelim $i $n
    }
    for { set i [expr $n-1] } { $i > -1 } { incr i -1 } {
	set Mat($Rx($i),$n) [expr $Mat($Rx($i),$n)/$Mat($Rx($i),$Cx($i))]
	set Mat($Rx($i),$Cx($i)) 1
	GaussSupraelim $i $n
    }
    return 1
}

### new version of least squares fit yielding one of the affine,
#    affine conformal, or affine conformal without rotation  transformation
## using equations worked out by J. B. Mehl 20070712, corrected 20070808
#
# 	(xt,yt) terrain coordinates
# 	(xm,-ym) map coordinates (pixel, ym grows downwards!)
#  with origin of pixel coordinates in the upper left corner of image
#
# 	(e,f) terrain coordinates of upper left corner, map (0,0)
#
#       with   n = number of points
#              Sxm   = sum_i xm_i
#              Sxtxm = sum_i xt_i xm_i
#              ...
#
#  parameters of an affine conformal with no rotation transformation
#
# 	(xt,yt) = lambda (xm,ym) + (e,f)
#
# 	(e,f) = lambda (xm_0,H-ym_0) + (xt_0,yt_0)
#         where H is map height and (xt_0,yt_0) the terrain coordinates of
#         map lower corner
#
#       for all control points replace vertical coordinate ym_i by
#         H-ym-i
#
# 	F = sum_i (xt_i-xt_0-lambda xm_i)^2+(yt_i-yt_0-lambda ym_i)^2
# 	min
#
# 	d(F,xt_0)=0 =>   Sxt - n xt_0 - lambda Sxm = 0
# 		=>   xt_0 = (Sxt - lambda Sxm) / n
# 	d(F,yt_0)=0 =>   Syt - n yt_0 - lambda Sym = 0
# 		=>   yt_0 = (Syt - lambda Sym) / n
# 	d(F,lambda)=0 =>   Sxtxm - xt_0 Sxm - lambda Sxmxm +
# 				Sytm - y_t0 Sym - lambda Symym
#
#       lambda = (n (Sxtxm + Sytym) - Sxt Sxm - Syt Sym) /
# 		   (n (Sxmxm + Symym) - Sxm^2 - Sym^2)
#
#       e = (Sxt - lambda Sxm) / n
#
#       f = (Syt + lambda Sym) / n
#
#  parameters of an affine conformal transformation
#
#                 | a  b |
#       (xt,yt) = | c  d | (xm,ym) + (e,f)
#
#       with   a = -d = s cos r    b = c = s sin r
#       NB: as GPSMan procs use b = c = -s sin r, the equation for b below
#           must be preceded by a minus sign
#
#       minimize sum_i d(xt_i)^2+d(yt_i)^2
#       where
#            d(xt_i) = xt_i - e - a xm_i - b ym_i
#            d(yt_i) = yt_i - f - b xm_i + a ym_i
#
#            |  n    0      Sxm         Sym     | | e |   |     Sxt     |
#            |  0    n     -Sym         Sxm     | | f |   |     Syt     |
#        M = | Sxm -Sym Sxmxm+Symym      0      | | a | = | Sxtxm-Sytym |
#            | Sym  Sxm      0      Sxmxm+Symym | | b |   | Sxtym+Sytxm |
#
#        e = (Sxm (Sxtxm-Sytym)+Sym (Sxtym+Sytxm)-Sxt (Sxmxm+Symym)) /
#               (Sxm^2 + Sym^2 - n (Sxmxm+Symym))
#
#        f = (Sym (Sytym-Sxtxm)+Sxm (Sxtym+Sytxm)-Syt (Sxmxm+Symym)) /
#               (Sxm^2 + Sym^2 - n (Sxmxm+Symym))
#
#        a = (Sxt Sxm - Syt Sym - n (Sxtxm - Sytym)) /
#               (Sxm^2 + Sym^2 - n (Sxmxm+Symym))
#
#        b = (Sxt Sym + Syt Sxm - n (Sxtym + Sytxm)) /
#               (Sxm^2 + Sym^2 - n (Sxmxm+Symym))
#
#  parameters of an affine transformation
#
#                 | a  b |
#       (xt,yt) = | c  d | (xm,ym) + (e,f)
#
#       minimize sum_i d(xt_i)^2+d(yt_i)^2
#       where
#            d(xt_i) = xt_i - xt(xm_i,ym_i,a,...,f)
#            similarly for d(yt_i)
#
#       => minimize sum_i (xt_i-e-a xm_i-ab ym_i)^2 +
#                         (yt_i-f-ac xm_i-ad ym_i)^2
#
#           |  n    Sxm   Sym  |        | e |         |  Sxt  |
#       M = | Sxm  Sxmxm Sxmym |   Vx = | a |    Cx = | Sxtxm |
#           | Sym  Sxmym Symym |        | b |         | Sxtym |
#
#
#                                       | f |         |  Syt  |
#                                  Vy = | c |    Cy = | Sytxm |
#                                       | d |         | Sytym |
#
#       equations:  M Vx = Cx, M Vy = Cy  yield [aij], e, f
#

proc MapInitLeastSquaresTransf {args} {
    # use datum and geodetic and pixel coordinates of points either
    #  from the graphical interface or from a GPSMan least-squares
    #  info file to compute parameters of transformation using a
    #  least-squares fit
    #  $LSqsTransf  in {Affine, AffineConf, NoRot} is the transformation type
    #  $args==""  points are WPs, with $MapLoadWPs the list of
    #          indices of relevant WPs, used by proc MapGeoRefPoints
    #          that sets MapLoadPos(_,_) to the pixel coordinates
    #       ==pair whose 2nd element is the projection and 1st element is
    #          a list with list of latd, longd and datum, a list of pairs
    #          with pixel coordinates (x, y, in Tcl sign convention with
    #          origin at upper left corner), a list with waypoints names (empty
    #          for positions given explictly), and the input file path
    # return list of WP names to be displayed or 0 on error or operation being
    #  cancelled
    global MPData MTData MapScale MPData MapLoadWPs MapLoadPos MapImageHeight \
	 LSqsTransf Mat Rx SUBDTUNIT DSCALE SUBDSCALE TXT GFShowInfo GFVisible

    set dispWPs ""
    if { $args == "" } {
	# from interface: always show fit information, no WPs to display
	set showinfo 1
	set names ""
	set fpath ""
	set n [llength $MapLoadWPs]
	if { [set tcs [MapGeoRefPoints $n]] == -1 } { return 0 }
	set xys ""
	for { set k 0 } { $k < $n } { incr k } {
	    lappend xys [list $MapLoadPos($k,x) $MapLoadPos($k,y)]
	}
	foreach ix $MapLoadWPs {
	    lappend names [NameOf WP $ix]
	}
    } else {
	# from file: show fit information and display WPs if required
	set showinfo $GFShowInfo
	foreach "p proj" $args { break }
	foreach "ldds xys names fpath" $p { break }
	set tcs ""
	foreach ldd $ldds {
	    foreach "latd longd datum" $ldd {}
	    lappend tcs [Proj${proj}Point MPData $latd $longd $datum]
	}
	set n [llength $names]
	if { $GFVisible } { set dispWPs $names }
    }
    # least-squares on $tcs and $xys
    MapTransfIs $LSqsTransf
    if { $LSqsTransf == "NoRot" } {
	set ymoff $MapImageHeight ; set ymsign -1
    } else { set ymoff 0 ; set ymsign 1 }
    foreach i "xt xtxm xtym yt ytxm ytym xm xmxm xmym ym ymym" {
	set s($i) 0.0
    }
    foreach tc $tcs pc $xys {
	foreach "xt yt" $tc {}
	foreach "xm ym" $pc {}
	set ym [expr $ymoff+$ymsign*$ym]
	foreach i "xt yt xm ym" {
	    set s($i) [expr $s($i)+[set $i]]
	}
	set s(xmym) [expr $s(xmym)+1.0*$xm*$ym]
	foreach j "x y" {
	    set z [set ${j}m]
	    set w "${j}m${j}m"
	    set s($w) [expr $s($w)+1.0*$z*$z]
	    foreach i "x y" {
		set w "${i}t${j}m"
		set s($w) [expr $s($w)+1.0*[set ${i}t]*$z]
	    }
	}
    }
    set dfreedom [expr $n+$n]
    set n [expr 1.0*$n]
    switch $LSqsTransf {
	Affine {
	    incr dfreedom -6
	    foreach d "xt yt" ps [list "e a b" "f c d"] {
		catch {unset Mat}
		set Mat(0,0) $n
		set Mat(0,1) [set Mat(1,0) $s(xm)]
		set Mat(0,2) [set Mat(2,0) $s(ym)]
		set Mat(1,1) $s(xmxm)
		set Mat(1,2) [set Mat(2,1) $s(xmym)]
		set Mat(2,2) $s(ymym)
		set Mat(0,3) $s($d)
		set Mat(1,3) $s(${d}xm)
		set Mat(2,3) $s(${d}ym)
		if { [GaussReduce 3] != 1 } { return 0 }
		foreach p $ps i "0 1 2" {
		    set MTData($p) $Mat($Rx($i),3)
		}
	    }
	    # the following parameters simplify computations
	    set MTData(det) [expr $MTData(a)*$MTData(d)-$MTData(b)*$MTData(c)]
	    set MTData(k2m1) [expr $MTData(b)*$MTData(f)-$MTData(d)*$MTData(e)]
	    set MTData(k4m3) [expr $MTData(c)*$MTData(e)-$MTData(a)*$MTData(f)]
	    # scale along the xm-axis when variation of ym=0
	    set MapScale [expr abs($MTData(a))]
	}
	AffineConf {
	    incr dfreedom -4
	    set zxym [expr $s(xmxm)+$s(ymym)]
	    set zuyvx [expr $s(xtym)+$s(ytxm)]
	    set zuxMvy [expr $s(xtxm)-$s(ytym)]
	    if { abs([set den [expr $s(xm)*$s(xm)+$s(ym)*$s(ym)-$n*$zxym]]) \
		     < 1e-50 } { return 0 }
 	    set MTData(a) [expr ($s(xt)*$s(xm)-$s(yt)*$s(ym)-$n*$zuxMvy) \
 			       / $den]
	    # b sign is changed in GPSMan affine conformal procs
 	    set MTData(b) [expr -($s(xt)*$s(ym)+$s(yt)*$s(xm)-$n*$zuyvx) /$den]
 	    set MTData(e) [expr ($s(xm)*$zuxMvy+$s(ym)*$zuyvx-$s(xt)*$zxym) \
 			       / $den]
 	    set MTData(f) [expr ($s(xm)*$zuyvx-$s(ym)*$zuxMvy-$s(yt)*$zxym) \
 			       / $den]
	    # the following parameters make calculations easier
 	    set MTData(det) [expr $MTData(a)*$MTData(a)+$MTData(b)*$MTData(b)]
 	    set MTData(k2m1) [expr $MTData(b)*$MTData(f)-$MTData(a)*$MTData(e)]
 	    set MTData(k4m3) [expr $MTData(b)*$MTData(e)+$MTData(a)*$MTData(f)]
	    # scale along the xm-axis when variation of ym=0
	    set MapScale [expr abs($MTData(a))]
	}
	NoRot {
	    incr dfreedom -3
	    if { [catch {set k [expr 1.0*($n*($s(xtxm)+$s(ytym))- \
					      $s(xt)*$s(xm)-$s(yt)*$s(ym)) /  \
				    ($n*($s(xmxm)+$s(ymym)) - $s(xm)*$s(xm) - \
					 $s(ym)*$s(ym))]}] } {
		return 0
	    }
	    # assuming:
	    # set MapLoadPos(origin,x) 0
	    # set MapLoadPos(origin,y) 0
	    set MTData(lambda) [set MapScale $k]
	    set MTData(e) [expr 1.0*($s(xt)-$k*$s(xm))/$n]
	    set MTData(f) [expr 1.0*($s(yt)-$k*$s(ym))/$n+$k*$MapImageHeight]
	}
    }
    if { $showinfo } {
	# deviations of control points in user units (m, ft)
	set sc [expr $MapScale*$DSCALE/$SUBDSCALE/1000.0]
    
	append info $TXT(lstsqs) " / " $TXT(TRNSF$LSqsTransf) "\n"
	if { $fpath != "" } {
	    append info $TXT(file) ": " $fpath "\n"
	}
	set xt $TXT(xtcoord) ; set yt $TXT(ytcoord)
	set delta $TXT(delta)
	append info \
	    "\t\t$xt\t$yt\t$delta $xt\t$delta $yt\t$delta $TXT(residual)\n" \
	    "\t\tm\tm\t$SUBDTUNIT\t$SUBDTUNIT\t$SUBDTUNIT\n"
	set sumdrt2 0
	foreach tc $tcs pc $xys name $names {
	    foreach "xt yt" $tc {}
	    foreach "xm ym" $pc {}
	    foreach "nxm nym" [MapApply${LSqsTransf}Transf $xt $yt] {}
	    set dxt [expr $sc*($nxm-$xm)]
	    set dyt [expr $sc*($nym-$ym)]
	    set drt [expr sqrt($dxt*$dxt+$dyt*$dyt)]
	    set sumdrt2 [expr $sumdrt2+$drt*$drt]
	    foreach v "dxt dyt drt" { set $v [format %7.2f [set $v]] }
	    append info $name "\t\t" [format %.2f $xt] "\t" [format %.2f $yt] \
		"\t" $dxt "\t" $dyt "\t" $drt "\n"
	}
	set rms [expr sqrt($sumdrt2/($n+$n))]
	set rmsdf [expr sqrt(1.0*$sumdrt2/$dfreedom)]
	append info "$TXT(rmsxydev) = " [format %.2f $rms] $SUBDTUNIT "\n" \
	    "$TXT(resstderr) = " [format %.2f $rmsdf] $SUBDTUNIT "\n"
	DisplayInfo $info tabs \
	    [list -10 right -24 right -38 right -46 right -54 right -62 right]
    }
    return $dispWPs
}

##### TFW metadata file

proc MapInitTFWTransf {data proj} {
    # compute parameters of affine transformation from the values in a
    #  TFW metadata file
    # based on the ESRI ArcGIS 2.9 manual (consulted November 2009) at
    #  http://webhelp.esri.com/arcgisdesktop/9.2/index.cfm?id=2676&pid=2664&topicname=World_files_for_raster_datasets
    # this is an affine transformation with
    # (xt,yt) = [ Aij ] x (xm,ym) + (C,F)
    #  where  A11 = A = tfw1
    #         A12 = B = tfw3
    #         A21 = D = tfw2
    #         A22 = E = tfw4 (correct sign for y map coordinates growing down)
    #         C = tfw5, F = tfw6
    #
    #  therefore in the GPSMan notation
    #     a = A,  b = B,  e = C
    #     c = D,  d = E,  f = F
    #
    #  $data is a pair with an empty list and a list with the 6 parameters
    #        (tfw_ in the formulae above)
    #  $proj is the projection
    # return 0 if determinant is 0
    global MTData MapScale

    MapTransfIs Affine
    foreach m {a c b d e f} v [lindex $data 1] {
	# to avoid integer divisions
	set MTData($m) [expr 0.0+$v]
    }
    set MTData(det) [expr $MTData(a)*$MTData(d)-$MTData(b)*$MTData(c)]
    if { abs($MTData(det)) < 1e-30 } {
	return 0
    }
    # the following parameters simplify computations
    set MTData(k2m1) [expr $MTData(b)*$MTData(f)-$MTData(d)*$MTData(e)]
    set MTData(k4m3) [expr $MTData(c)*$MTData(e)-$MTData(a)*$MTData(f)]
    # scale along the xm-axis when variation of ym=0
    set MapScale [expr abs($MTData(a))]
    return 1
}

##### partial support for OziExplorer .map files

proc MapInitOziMapTransf {data proj} {
    # use datum and geodetic and pixel coordinates of points in Ozi .map
    #  file to compute parameters of affine or affine conformal transformation
    #  $data is pair with list of latd,longd,datum and
    #        list of pairs with pixel coordinates (in Tcl sign convention)
    #        aligned with the previous one
    #  $proj is the projection
    # return 0 on error
    global MTData MapScale MPData MapLoadPos

    foreach "opds opxs" $data {}
    if { [set sel [MapReduceCPoints $opxs]] == "" } { return 0 }
    set n 0
    set projcs ""
    foreach ix $sel {
	foreach "x y" [lindex $opxs $ix] {}
	set MapLoadPos($n,x) $x ; set MapLoadPos($n,y) $y
	set lld [lindex $opds $ix]
	lappend projcs [eval Proj${proj}Point MPData $lld]
	incr n
    }
    if { $n == 3 } {
	set r [MapInitAffineTransf $projcs]
    } else { set r [MapInitAffineConfTransf $projcs] }
    if { $r != 0 } {
	FixMapScale $proj
    }
    return $r
}

proc MapReduceCPoints {xys} {
    # reduce a set of points given by a list of their planar
    #  coordinates to 3 or 2
    # return list of indices in $xys of selected points or "" if less than 2
    #
    # reduction principle (suggested by Luis Damas): maximize the
    #  minimum distance between each 2 points in a triangle
    # algorithm: find the leftmost 3 edges forming a triangle
    #  in a list of triples with distance between two points and
    #  the 2 points indices, sorted by decreasing distance

    if { [set n [llength $xys]] < 2 } { return "" }
    if { $n == 2 } { return "0 1" }
    # form list of triples
    set ds ""
    set ix 0
    foreach p [lrange $xys 0 end-1] {
	foreach "x y" $p {}
	set ixi [expr $ix+1]
	foreach pi [lrange $xys 1 end] {
	    foreach "xi yi" $pi {}
	    set dx [expr $x-$xi] ; set dy [expr $y-$yi]
	    lappend ds [list [expr sqrt($dx*$dx+$dy*$dy)] $ix $ixi]
	    incr ixi
	}
	incr ix
    }
    set ds [lsort -decreasing -real -index 0 $ds]
    foreach t $ds {
	foreach "d ix0 ix1" $t {}
	# $ix0 < $ix1 by construction
	if { [catch {set s0 $seen($ix0)}] } {
	    set seen($ix0) $ix1
	    if { [catch {set seen($ix1)}] } {
		set seen($ix1) $ix0
	    } else { lappend seen($ix1) $ix0 }
	} elseif { [catch {set s1 $seen($ix1)}] } {
	    set seen($ix1) $ix0
	    lappend seen($ix0) $ix1
	} elseif { [set ix2 [Intersect1 $s0 $s1]] == "" } {
	    lappend seen($ix0) $ix1
	    lappend seen($ix1) $ix0
	} else {
	    # distance between $ix0 and $ix1 is the minimum of the
	    #  3 distances, therefore if the points are colinear
	    #  $ix2 is an extreme
	    break
	}	
    }
    foreach n "0 1 2" {
	foreach "x$n y$n" [lindex $xys [set ix$n]] {}
    }
    # colinear?
    if { $x0 == $x2 || $x1 == $x2 } {
	if { abs ($x0-$x1) < 10 } {
	    if { abs($y0-$y2) > abs($y1-$y2) } {
		return [list $ix2 $ix0]
	    }
	    return [list $ix2 $ix1]
	}
    } else {
	set dd02 [expr 1.0*($y0-$y2)/($x0-$x2)]
	set dd12 [expr 1.0*($y1-$y2)/($x1-$x2)]
	set ddm [expr { abs($dd02) > abs($dd12) } ? $dd02 : $dd12]
	if { abs(($dd02-$dd12)/$ddm) < 0.5 } {
	    if { abs($x0-$x2) > abs($x1-$x2) } {
		return [list $ix2 $ix0]
	    }
	    return [list $ix2 $ix1]
	}
    }
    return [list $ix2 $ix1 $ix0]
}