From 2fcbc389c699edfd33508973a084a9795d9dc568 Mon Sep 17 00:00:00 2001 From: nuknal Date: Tue, 10 Sep 2024 14:13:30 +0800 Subject: [PATCH] 4*4*5 --- config/config.yaml | 5 +++++ pkg/auxilary/aux.go | 3 --- pkg/auxilary/image_time.go | 10 ++++++++- pkg/config/config.go | 12 +++++++++++ pkg/producer/grid_img.go | 14 ++++++------ pkg/producer/rpc.go | 23 +++++++++++++------- pkg/producer/rpc_helper.go | 41 ++++++++++++++++++++++++------------ pkg/producer/rpc_optimize.go | 6 ++++++ 8 files changed, 82 insertions(+), 32 deletions(-) diff --git a/config/config.yaml b/config/config.yaml index d0089a3..9c757be 100644 --- a/config/config.yaml +++ b/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" \ No newline at end of file diff --git a/pkg/auxilary/aux.go b/pkg/auxilary/aux.go index d578fcd..b6e6e4a 100644 --- a/pkg/auxilary/aux.go +++ b/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 } diff --git a/pkg/auxilary/image_time.go b/pkg/auxilary/image_time.go index 5739993..21725c3 100644 --- a/pkg/auxilary/image_time.go +++ b/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 { diff --git a/pkg/config/config.go b/pkg/config/config.go index 3b417d9..1e190f2 100644 --- a/pkg/config/config.go +++ b/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, + }, } } diff --git a/pkg/producer/grid_img.go b/pkg/producer/grid_img.go index 093e4ab..ce9d624 100644 --- a/pkg/producer/grid_img.go +++ b/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 - hmin = hmin - 500 + hmax = (hmax+hmin)/2.0 + config.GCONFIG.RPC.AltitudeRange/2 + hmin = (hmax+hmin)/2.0 - config.GCONFIG.RPC.AltitudeRange/2 dh := (hmax - hmin) / (k) var h []int diff --git a/pkg/producer/rpc.go b/pkg/producer/rpc.go index 5491b05..7e39267 100644 --- a/pkg/producer/rpc.go +++ b/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, - gridsize: 20, + elevationLayer: config.GCONFIG.RPC.AltitudeLayer, + 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.lineScale = float64(rpc.scene.Height) / 2.0 - rpc.sampOffset = float64(rpc.scene.Width) / 2.0 - rpc.sampScale = float64(rpc.scene.Width) / 2.0 - rowVec = normalize2(rowVec, rpc.lineOffset, rpc.lineScale) - colVec = normalize2(colVec, rpc.sampOffset, rpc.sampScale) + // rpc.lineOffset = float64(rpc.scene.Height) / 2.0 + // rpc.lineScale = float64(rpc.scene.Height) + // rowVec = normalize2(rowVec, rpc.lineOffset, rpc.lineScale) + // rpc.sampOffset = float64(rpc.scene.Width) / 2.0 + // rpc.sampScale = float64(rpc.scene.Width) + // 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) diff --git a/pkg/producer/rpc_helper.go b/pkg/producer/rpc_helper.go index 00b8533..77ab8f8 100644 --- a/pkg/producer/rpc_helper.go +++ b/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) - wmx.MulVec(&wm, &x1) - wf.MulVec(weights, f) - wmx.SubVec(&wmx, &wf) - wmx.MulElemVec(&wmx, &wmx) - v := math.Sqrt(mat.Max(&wmx)) - log.Println("iteration:", iterations, "v-err:", v) + numerator1 := mat.NewVecDense(20, nil) + denominator1 := mat.NewVecDense(20, nil) + denominator1.SetVec(0, 1.0) + numerator1.SetVec(0, x0.AtVec(0)) + for i := 1; i < 20; i++ { + numerator1.SetVec(i, x1.AtVec(i)) + denominator1.SetVec(i, x1.AtVec(i+19)) + } - vx = append(vx, &VX{v: v, x: x1}) - if math.Abs(v-v0) < epsilon { - break + errorSquared := 0.0 + lambda := 1e-4 + 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) diff --git a/pkg/producer/rpc_optimize.go b/pkg/producer/rpc_optimize.go index 7c14385..eab2246 100644 --- a/pkg/producer/rpc_optimize.go +++ b/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 }