"Fdrawtps" = function(coords1, coords2, ngrid = 15., xl = "", psize = 0.04, draw1 = T, draw2 = T, number = T)
{
	#thin-plate spline from coords1 to coords2
	#
	if(mode(coords1) == "complex") c1 <- coords1 else {
		c1a <- matrix(coords1, 2.)
		c1 <- complex(real = c1a[1.,  ], imag = c1a[2.,  ])
	}
	if(mode(coords2) == "complex")
		c2 <- coords2
	else {
		c2a <- matrix(coords2, 2.)
		c2 <- complex(real = c2a[1.,  ], imag = c2a[2.,  ])
	}
	m1 <- FLmatrix(c1)
	nland <- length(c2)
	if(nland != length(c1))
		stop("landmark counts don't match")
	#tpscoeff <- complex(real = solve(m1, c(Re(c2), 0., 0., 0.)), imag = solve(m1, c(Im(c2), 0., 0., 0.)))
	m1.svd <- svd(m1)
	m1.inv <- m1.svd$v %*% diag(1/m1.svd$d) %*% t(m1.svd$u)
	#tpscoeff <- complex(real = solve(m1, c(Re(c2), 0., 0., 0.)), imag = solve(m1, c(Im(c2), 0., 0., 0.)))
	tpscoeff <- complex(real = m1.inv %*% c(Re(c2), 0., 0., 0.), imag = m1.inv %*% c(Im(c2), 0., 0., 0.))
	spacing <- max(diff(range(Re(c1))), diff(range(Im(c1))))/(ngrid - 4.)
	gridx0 <- min(Re(c1)) - 1.5 * spacing
	gridy0 <- min(Im(c1)) - 1.5 * spacing
	ngcol <- ceiling(diff(range(Re(c1)))/spacing) + 4.
	ngrow <- ceiling(diff(range(Im(c1)))/spacing) + 4.
	gcol <- gridx0 - spacing + spacing * (1.:ngcol)
	grow <- gridy0 - spacing + spacing * (1.:ngrow)
	gmatrix <- outer(gcol, grow, function(x, y)
	complex(real = x, imag = y))
	affine <- matrix(tpscoeff[nland + 1.], ngcol, ngrow) + tpscoeff[nland + 2.] * Re(gmatrix) + tpscoeff[nland + 3.] * Im(gmatrix)
	kernels <- array(gmatrix, c(ngcol, ngrow, nland))
	centers <- aperm(array(c1, c(nland, ngrow, ngcol)), c(3., 2., 1.))
	sepmtx <- kernels - centers
	terms <- Re(sepmtx * Conj(sepmtx))
	terms <- ifelse(terms > 0., 0.5 * terms * log(terms), 0.)
	coeffs <- aperm(array(tpscoeff[1.:nland], c(nland, ngrow, ngcol)), c(3., 2., 1.))
	partials <- terms * coeffs
	nonaff <- apply(partials, c(1., 2.), sum)
	par(mar = c(1, 1, 1, 1))
	# original grid
	if(draw2) {
		if(draw1) {
			par(mfrow = c(1, 2))
			plotrange <- max(diff(range(Re(gmatrix))), diff(range(Im(gmatrix))))
			par(pty = "s")
			plot(c(min(Re(gmatrix)), min(Re(gmatrix)) + plotrange),
			 c(min(Im(gmatrix)), min(Im(gmatrix)) + plotrange), type = "n", axes = F, xlab = xl, ylab = "")
			for(i in 1.:ngcol)
				lines(gmatrix[i,  ])
			for(i in 1.:ngrow)
				lines(gmatrix[, i])
			par(mkh = psize)
			#rudimentary plot
			if(number) {
				text(Re(c1), Im(c1), 1:length(c1))
			}
			else {
				points(c1, pch = 16)
			}
			par(mkh = 0.01)
		}
	}
	# distorted grid
	tpsgrid <- affine + nonaff
	plotrange <- max(diff(range(Re(tpsgrid))), diff(range(Im(tpsgrid))))
	par(pty = "s")
	plot(c(min(Re(tpsgrid)), min(Re(tpsgrid)) + plotrange), 
		c(min(Im(tpsgrid)), min(Im(tpsgrid)) + plotrange), type = "n", axes = F, xlab = xl, ylab = "")
	for(i in 1.:ngcol)
		lines(tpsgrid[i,  ])
	for(i in 1.:ngrow)
		lines(tpsgrid[, i])
	par(mkh = psize)
	#rudimentary plot
	if(number) {
		text(Re(c2), Im(c2), 1:length(c2))
	}
	else {
		points(c2, pch = 16)
	}
	par(mkh = 0.01)
	tpscoeff
}