4*4*5

config/config.yaml

@@ -19,5 +19,10 @@ radiation:
   hf_band_stop_width: 24 # 带阻滤波器宽度
   scene_moment_matching: false
 
+rpc:
+  grid_size: 3
+  altitude_layer: 5
+  altitude_range: 1000
+
 dem:
   dem_1km: "dem/gdlebm.tif"

pkg/auxilary/aux.go

@@ -48,9 +48,6 @@ func ExtractAux(auxfile string) ([]*AuxFrameHead, []*AuxFocalBox, []*AuxPlatform
 // 长光卫星姿态和GPS数据更新频率为 4 次/秒
 func transfromGPSandAttMicrosec(microsec uint32) uint32 {
 	unit := uint32(250000)
-	// return microsec
-
 	microsec = (microsec / unit) * unit
-
 	return microsec
 }

pkg/auxilary/image_time.go

@@ -37,7 +37,7 @@ func (it *ImageTime) Extract(aps []*AuxPlatform) {
 // Interp 内插出成像时刻
 // 参数: imgrow 图像行号 cross 图像跨度 aps 辅助平台列表
 // 返回值: 成像时刻(UTC seconds)
-func (it *ImageTime) Interp(imgrow int, cross int) (float64, float64) {
+func (it *ImageTime) Interp0(imgrow int, cross int) (float64, float64) {
 	// 内插出成像时刻
 	u := int(imgrow / cross)
 	var u1, u2 *auxT
@@ -53,6 +53,14 @@ func (it *ImageTime) Interp(imgrow int, cross int) (float64, float64) {
 	return t, dt
 }
 
+func (it *ImageTime) Interp(imgrow int, cross int) (float64, float64) {
+	n := len(it.auxT)
+	rCnt := (it.auxT[n-1].Row-it.auxT[0].Row)*cross + 1
+	dt := (it.auxT[n-1].Tutc - it.auxT[0].Tutc) / float64(rCnt)
+	t := it.auxT[0].Tutc + dt*float64(imgrow)
+	return t, dt
+}
+
 func (it *ImageTime) Store(imgtimeFile string) error {
 	f, err := os.Create(imgtimeFile)
 	if err != nil {

pkg/config/config.go

@@ -11,6 +11,7 @@ type Config struct {
 	BrowserImg     BrowserImgConfig     `yaml:"browser_img" mapstructure:"browser_img"`
 	Log            LogConfig            `yaml:"log" mapstructure:"log"`
 	Dem            DemConfig            `yaml:"dem" mapstructure:"dem"`
+	RPC            RPCConfig            `yaml:"rpc" mapstructure:"rpc"`
 }
 
 type LogConfig struct {
@@ -52,6 +53,12 @@ type BrowserImgConfig struct {
 	CumulativeCutUpper float64           `yaml:"cumulative_cut_upper" mapstructure:"cumulative_cut_upper"`
 }
 
+type RPCConfig struct {
+	GridSize      int `yaml:"grid_size" mapstructure:"grid_size"`
+	AltitudeRange int `yaml:"altitude_range" mapstructure:"altitude_range"`
+	AltitudeLayer int `yaml:"altitude_layer" mapstructure:"altitude_layer"`
+}
+
 var GCONFIG Config
 
 func init() {
@@ -94,5 +101,10 @@ func init() {
 			DemDir:  "dem/SRTM",
 			Dem1Km:  "dem/gdlebm.tif",
 		},
+		RPC: RPCConfig{
+			GridSize:      3,
+			AltitudeRange: 1000,
+			AltitudeLayer: 5,
+		},
 	}
 }

pkg/producer/grid_img.go

@@ -1,25 +1,27 @@
 package producer
 
+import "starwiz.cn/sjy01/image-proc/pkg/config"
+
 type GridPoint struct {
 	Row, Col, H int
 }
 
 // 网格点要覆盖边界，甚至大于边界
 func gridImage2(m, n, height, width, k, hmin, hmax int) []*GridPoint {
-	a := int((height) / (m + 1))
+	a := int((height) / (m))
 	var lines []int
-	for i := 1; i <= m; i++ {
+	for i := 0; i <= m; i++ {
 		lines = append(lines, a*i)
 	}
 
-	b := int((width) / (n + 1))
+	b := int((width) / (n))
 	var samples []int
-	for i := 1; i <= n; i++ {
+	for i := 0; i <= n; i++ {
 		samples = append(samples, b*i)
 	}
 
-	hmax = hmax + 500
+	hmax = (hmax+hmin)/2.0 + config.GCONFIG.RPC.AltitudeRange/2
-	hmin = hmin - 500
+	hmin = (hmax+hmin)/2.0 - config.GCONFIG.RPC.AltitudeRange/2
 	dh := (hmax - hmin) / (k)
 
 	var h []int

pkg/producer/rpc.go

@@ -11,6 +11,7 @@ import (
 
 	log "github.com/sirupsen/logrus"
 	"gonum.org/v1/gonum/mat"
+	"starwiz.cn/sjy01/image-proc/pkg/config"
 	"starwiz.cn/sjy01/image-proc/pkg/dem"
 )
 
@@ -56,8 +57,8 @@ type RPCModel struct {
 // rational polynomial coeffients
 func NewRPC(r *Registrator, scene *Scene, rpb string) *RPC {
 	rpc := RPC{
-		elevationLayer: 5,
+		elevationLayer: config.GCONFIG.RPC.AltitudeLayer,
-		gridsize:       20,
+		gridsize:       config.GCONFIG.RPC.GridSize,
 		registrator:    r,
 		scene:          scene,
 		rpb:            rpb,
@@ -140,22 +141,28 @@ func (rpc *RPC) RPC() error {
 	rpc.saveVec(strings.Replace(rpc.scene.Tiff, ".tiff", ".vec.txt", -1),
 		rowVec, colVec, latVec, lonVec, heightVec)
 
-	rpc.lineOffset = float64(rpc.scene.Height) / 2.0
+	// rpc.lineOffset = float64(rpc.scene.Height) / 2.0
-	rpc.lineScale = float64(rpc.scene.Height) / 2.0
+	// rpc.lineScale = float64(rpc.scene.Height)
-	rpc.sampOffset = float64(rpc.scene.Width) / 2.0
+	// rowVec = normalize2(rowVec, rpc.lineOffset, rpc.lineScale)
-	rpc.sampScale = float64(rpc.scene.Width) / 2.0
+	// rpc.sampOffset = float64(rpc.scene.Width) / 2.0
-	rowVec = normalize2(rowVec, rpc.lineOffset, rpc.lineScale)
+	// rpc.sampScale = float64(rpc.scene.Width)
-	colVec = normalize2(colVec, rpc.sampOffset, rpc.sampScale)
+	// colVec = normalize2(colVec, rpc.sampOffset, rpc.sampScale)
+	// rpc.heightOffset = (rpc.minH + rpc.maxH) / 2.0
+	// rpc.heightScale = 1000.0
+	// heightVec = normalize2(heightVec, rpc.heightOffset, rpc.heightScale)
 
 	latVec, rpc.latOffset, rpc.latScale = normalize(latVec)
 	lonVec, rpc.longOffset, rpc.longScale = normalize(lonVec)
 	heightVec, rpc.heightOffset, rpc.heightScale = normalize(heightVec)
+	rowVec, rpc.lineOffset, rpc.lineScale = normalize(rowVec)
+	colVec, rpc.sampOffset, rpc.sampScale = normalize(colVec)
 
 	rpc.saveVec(strings.Replace(rpc.scene.Tiff, ".tiff", ".vec_norm.txt", -1),
 		rowVec, colVec, latVec, lonVec, heightVec)
 
 	// 设计矩阵 B = [ 20个分子系数 19个分母系数 ] W = 权矩阵 R = 观测结果
 	method := SolveMethodNelderMead
+
 	// x = (B^T * W^2 * B)^-1 * (B^T * W^2 * R), 其中 x = [a1..a20 b2..b20]^T
 	// 行参数
 	// B := setupSystemOfEquations(rowVec, latVec, lonVec, heightVec)

pkg/producer/rpc_helper.go

@@ -58,6 +58,7 @@ func solveCoefficients(f, latVec, lonVec, heightVec *mat.VecDense, method SolveM
 	var MtF mat.VecDense
 	MtF.MulVec(M.T(), f)
 	x0.MulVec(invMtM, &MtF)
+	// return mat.Col(nil, 0, &x0), nil
 
 	if method == SolveMethodNelderMead {
 		numerator := mat.NewVecDense(20, nil)
@@ -83,12 +84,10 @@ func solveCoefficients(f, latVec, lonVec, heightVec *mat.VecDense, method SolveM
 	}
 
 	// 迭代
-	var wm mat.Dense
+
-	var wmx, wf mat.VecDense
 	var x1 mat.VecDense
 
 	var vx []*VX
-	v0 := 0.0
 	iterations := 0
 	maxIterations := 10
 	denominator := mat.NewVecDense(20, nil)
@@ -113,22 +112,36 @@ func solveCoefficients(f, latVec, lonVec, heightVec *mat.VecDense, method SolveM
 		MtW2F.MulVec(&MtW2, f)
 		x1.MulVec(invMtW2M, &MtW2F)
 
-		wm.Mul(weights, M)
+		numerator1 := mat.NewVecDense(20, nil)
-		wmx.MulVec(&wm, &x1)
+		denominator1 := mat.NewVecDense(20, nil)
-		wf.MulVec(weights, f)
+		denominator1.SetVec(0, 1.0)
-		wmx.SubVec(&wmx, &wf)
+		numerator1.SetVec(0, x0.AtVec(0))
-		wmx.MulElemVec(&wmx, &wmx)
+		for i := 1; i < 20; i++ {
-		v := math.Sqrt(mat.Max(&wmx))
+			numerator1.SetVec(i, x1.AtVec(i))
-		log.Println("iteration:", iterations, "v-err:", v)
+			denominator1.SetVec(i, x1.AtVec(i+19))
+		}
 
-		vx = append(vx, &VX{v: v, x: x1})
+		errorSquared := 0.0
-		if math.Abs(v-v0) < epsilon {
+		lambda := 1e-4
-			break
+		for i := 0; i < n; i++ {
+			predictedV := project(numerator1, denominator1, latVec.AtVec(i), lonVec.AtVec(i), heightVec.AtVec(i))
+			errorV := predictedV - f.AtVec(i)
+			errorSquared += errorV * errorV
+		}
+		fmt.Printf("squared error: %.8f\n", errorSquared)
+		var coeffsSquared float64
+		for i := 0; i < 20; i++ {
+			coeffsSquared += lambda * (numerator1.AtVec(i)*numerator1.AtVec(i) + denominator1.AtVec(i)*denominator1.AtVec(i))
 		}
 
 		x0 = x1
-		v0 = v
 		iterations++
+
+		fmt.Printf("squared error+lambda*coeffs: %.8f\n", coeffsSquared)
+		vx = append(vx, &VX{v: errorSquared, x: x1})
+		if errorSquared < 0.001 {
+			break
+		}
 	}
 
 	log.Println("iterations:", iterations)

pkg/producer/rpc_optimize.go

@@ -19,6 +19,12 @@ func objectiveFunc(numerator, denominator, f, latVec, lonVec, heightVec *mat.Vec
 		errorV := predictedV - f.AtVec(i)
 		errorSquared += errorV * errorV
 	}
+	// lambda := 1e-2
+	// coeffsSquaredError := 0.0
+	// for i := 0; i < 20; i++ {
+	// 	coeffsSquaredError += lambda * (numerator.AtVec(i)*numerator.AtVec(i) + denominator.AtVec(i)*denominator.AtVec(i))
+	// }
+
 	return errorSquared
 }