控制点过少时导致图像内畸变较大

pkg/producer/rpc.go

@@ -57,7 +57,7 @@ type RPCModel struct {
 func NewRPC(r *Registrator, scene *Scene, rpb string) *RPC {
 	rpc := RPC{
 		elevationLayer: 3,
-		gridsize:       3,
+		gridsize:       19,
 		registrator:    r,
 		scene:          scene,
 		rpb:            rpb,
@@ -150,14 +150,14 @@ func (rpc *RPC) RPC() error {
 	rpc.saveVec(strings.Replace(rpc.scene.Tiff, ".tiff", ".vec_norm.txt", -1),
 		rowVec, colVec, latVec, lonVec, heightVec)
 
-	// 设计矩阵 B = [ 20个分子系数 19个分母系数 ]
+	// 设计矩阵 B = [ 20个分子系数 19个分母系数 ] W = 权矩阵 R = 观测结果
-
+	method := SolveMethodNelderMead
-	// x = (B^T * B)^-1 * B^T * l, 其中 x = [a1..a20 b2..b20]^T
+	// x = (B^T * W^2 * B)^-1 * (B^T * W^2 * R), 其中 x = [a1..a20 b2..b20]^T
 	// 行参数
 	// B := setupSystemOfEquations(rowVec, latVec, lonVec, heightVec)
 	// J, err := SolveNormalEquation(B, rowVec)
 	log.Println("solving row coefficients")
-	J, err := solveCoefficients(rowVec, latVec, lonVec, heightVec, SolveMethodNelderMead)
+	J, err := solveCoefficients(rowVec, latVec, lonVec, heightVec, method)
 	if err != nil {
 		return err
 	}
@@ -166,7 +166,7 @@ func (rpc *RPC) RPC() error {
 	// D := setupSystemOfEquations(colVec, latVec, lonVec, heightVec)
 	// K, err := SolveNormalEquation(D, colVec)
 	log.Println("solving col coefficients")
-	K, err := solveCoefficients(colVec, latVec, lonVec, heightVec, SolveMethodNelderMead)
+	K, err := solveCoefficients(colVec, latVec, lonVec, heightVec, method)
 	if err != nil {
 		return err
 	}

pkg/producer/rpc_helper.go

@@ -46,28 +46,18 @@ func normalize2(v *mat.VecDense, vOff, vScale float64) *mat.VecDense {
 func solveCoefficients(f, latVec, lonVec, heightVec *mat.VecDense, method SolveMethod) ([]float64, error) {
 	M := setupSystemOfEquations(f, latVec, lonVec, heightVec)
 	n := f.Len()
-	weights := mat.NewDiagDense(n, nil)
-	for i := 0; i < n; i++ {
-		weights.SetDiag(i, 1.0)
-	}
 
-	// 计算初始值
+	// 计算初始值 x = (M^T * M)^-1 * (M^T * f)
 	var x0 mat.VecDense
-	w2 := mat.NewDiagDense(n, nil)
+	var MtM mat.Dense
-	for i := 0; i < n; i++ {
+	MtM.Mul(M.T(), M)
-		w2.SetDiag(i, weights.At(i, i)*weights.At(i, i))
+	invMtM, err := invertRPCMatrix(&MtM)
-	}
-	var MtW2 mat.Dense
-	MtW2.Mul(M.T(), w2)
-	var MtW2M mat.Dense
-	MtW2M.Mul(&MtW2, M)
-	invMtW2M, err := invertRPCMatrix(&MtW2M)
 	if err != nil {
 		return nil, err
 	}
-	var MtW2F mat.VecDense
+	var MtF mat.VecDense
-	MtW2F.MulVec(&MtW2, f)
+	MtF.MulVec(M.T(), f)
-	x0.MulVec(invMtW2M, &MtW2F)
+	x0.MulVec(invMtM, &MtF)
 
 	if method == SolveMethodNelderMead {
 		numerator := mat.NewVecDense(20, nil)
@@ -102,12 +92,15 @@ func solveCoefficients(f, latVec, lonVec, heightVec *mat.VecDense, method SolveM
 	iterations := 0
 	maxIterations := 10
 	denominator := mat.NewVecDense(20, nil)
+	w2 := mat.NewDiagDense(n, nil)
+	var MtW2, MtW2M mat.Dense
+	var MtW2F mat.VecDense
 	for iterations < maxIterations {
 		denominator.SetVec(0, 1.0)
 		for idx := 1; idx < 20; idx++ {
 			denominator.SetVec(idx, x0.AtVec(idx+19))
 		}
-		weights = setupWeightMatrix(denominator, latVec, lonVec, heightVec)
+		weights := setupWeightMatrix(denominator, latVec, lonVec, heightVec)
 		for i := 0; i < n; i++ {
 			w2.SetDiag(i, weights.At(i, i)*weights.At(i, i))
 		}

pkg/producer/rpc_optimize.go

@@ -1,8 +1,6 @@
 package producer
 
 import (
-	"fmt"
-
 	log "github.com/sirupsen/logrus"
 	"gonum.org/v1/gonum/mat"
 	"gonum.org/v1/gonum/optimize"
@@ -54,7 +52,6 @@ func solveNelderMead(num, den, f, latVec, lonVec, heightVec *mat.VecDense) (*mat
 	denominator := mat.NewVecDense(20, nil)
 	denominator.SetVec(0, 1.0)
 	numerator.SetVec(0, result.X[0])
-	fmt.Println("len of result.X: ", len(result.X))
 	for i := 1; i < 20; i++ {
 		numerator.SetVec(i, result.X[i])
 		denominator.SetVec(i, result.X[i+19])