From ce0a4fc370c54e258b1507143e401d4e622f1e55 Mon Sep 17 00:00:00 2001 From: nuknal Date: Tue, 27 Aug 2024 17:24:07 +0800 Subject: [PATCH] rpc --- cmd/proc.go | 4 +- pkg/dem/dem.go | 2 + pkg/producer/aux.go | 1 - pkg/producer/grid_img.go | 24 ++--- pkg/producer/output.go | 16 +++ pkg/producer/rpc.go | 203 +++++++++++++++++++++++++++------------ pkg/utils/tiff.go | 2 +- 7 files changed, 173 insertions(+), 79 deletions(-) diff --git a/cmd/proc.go b/cmd/proc.go index d5b72ef..1b8cefb 100644 --- a/cmd/proc.go +++ b/cmd/proc.go @@ -36,6 +36,8 @@ var procCmd = &cobra.Command{ logrus.SetLevel(config.GCONFIG.Log.LogLevel) + godal.RegisterAll() + // dem dem.Dem1KmLT = dem.NewDem1Km(config.GCONFIG.Dem.Dem1Km) @@ -51,7 +53,6 @@ var procCmd = &cobra.Command{ // reg.AuxPrint() - if err := reg.LoadMssRaw(); err != nil { logrus.Fatal(err) } @@ -60,7 +61,6 @@ var procCmd = &cobra.Command{ logrus.Fatal(err) } - godal.RegisterAll() os.MkdirAll(params.OutputDir, 0755) if doLUTRRC { diff --git a/pkg/dem/dem.go b/pkg/dem/dem.go index 30646be..a50bfd5 100644 --- a/pkg/dem/dem.go +++ b/pkg/dem/dem.go @@ -51,6 +51,7 @@ func (d *Dem1Km) Load() error { d.hUnit = 180.0 / float32(structure.SizeY) d.width = structure.SizeX d.height = structure.SizeY + hDataset.Close() return nil } @@ -97,5 +98,6 @@ func (d *Dem1Km) MinMaxElevationInRect(lng1, lat1, lng2, lat2 float64) (float64, } } } + return minElev, maxElev } diff --git a/pkg/producer/aux.go b/pkg/producer/aux.go index bde3f74..ac8b331 100644 --- a/pkg/producer/aux.go +++ b/pkg/producer/aux.go @@ -78,7 +78,6 @@ func (r *Registrator) SceneImageTime(scene *Scene) (start, center, end time.Time // FIXME: This function is not accurate enough. 四元数、成像时刻、GPS 等需要修改为插值获取 func (r *Registrator) SetSceneBoundary(scene *Scene) (topLeft, bottomRight orb.Point) { log.Info("using attitude quaternion to calculate image boundary...") - line0Start := r.calculateLatLonH(scene, 0, 0, 0) line0End := r.calculateLatLonH(scene, 0, scene.Width, 0) lineNStart := r.calculateLatLonH(scene, scene.Height, 0, 0) diff --git a/pkg/producer/grid_img.go b/pkg/producer/grid_img.go index 3beb887..3189866 100644 --- a/pkg/producer/grid_img.go +++ b/pkg/producer/grid_img.go @@ -4,24 +4,26 @@ type GridPoint struct { Row, Col, H int } -func gridImage(m, n, height, width, k, hmin, hmax int) []*GridPoint { - a := int(height / (m + 1)) +// 网格点要覆盖边界，甚至大于边界 +func gridImage2(m, n, height, width, k, hmin, hmax int) []*GridPoint { + a := int((height) / (m)) var lines []int - for i := 0; i < m; i++ { - lines = append(lines, a*(i+1)) + 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 := 0; i < n; i++ { - samples = append(samples, b*(i+1)) + for i := 1; i <= n; i++ { + samples = append(samples, b*i) } - averageH := (hmax - hmin) / 2 - dh := 500 // 高度差500m + hmax = hmax + 500 + hmin = hmin - 500 + dh := (hmax - hmin) / (k) var h []int - for i := 0; i < k; i++ { - h = append(h, averageH+(i-k/2)*dh) + for i := 1; i <= k; i++ { + h = append(h, hmin+dh*i) } var points []*GridPoint diff --git a/pkg/producer/output.go b/pkg/producer/output.go index b0caf8c..73a5b5f 100644 --- a/pkg/producer/output.go +++ b/pkg/producer/output.go @@ -89,3 +89,19 @@ func (r *Registrator) SaveRegisteredMssToRaw(raw string) error { func (r *Registrator) Report() error { return WriteReport(&r.report, r.Params.ReportFile) } + +func (r *Registrator) rpcKeywordInTif() { + // GDAL库对应的RPC关键词 + // keys := []string{ + // "ERR_BIAS", "ERR_RAND", + // "LINE_OFF", "SAMP_OFF", + // "LAT_OFF", "LONG_OFF", "HEIGHT_OFF", + // "LINE_SCALE", "SAMP_SCALE", + // "LAT_SCALE", "LONG_SCALE", "HEIGHT_SCALE", + // "LINE_NUM_COEFF", "LINE_DEN_COEFF", + // "SAMP_NUM_COEFF", "SAMP_DEN_COEFF", + // } + + // values := map[string]string{} + +} diff --git a/pkg/producer/rpc.go b/pkg/producer/rpc.go index 17d3dce..7fe9d0a 100644 --- a/pkg/producer/rpc.go +++ b/pkg/producer/rpc.go @@ -1,12 +1,15 @@ package producer import ( + "encoding/json" "fmt" "math" "os" + "strings" "github.com/duke-git/lancet/v2/mathutil" - "github.com/duke-git/lancet/v2/slice" + "github.com/paulmach/orb" + "github.com/paulmach/orb/geojson" log "github.com/sirupsen/logrus" "gonum.org/v1/gonum/mat" @@ -55,8 +58,8 @@ type RPCModel struct { // rational polynomial coeffients func NewRPC(r *Registrator, scene *Scene, rpb string) *RPC { rpc := RPC{ - elevationLayer: 3, - gridsize: 20, + elevationLayer: 9, + gridsize: 19, registrator: r, scene: scene, rpb: rpb, @@ -98,33 +101,29 @@ func (rpc *RPC) init() { rpc.latOffset = (rpc.minLat + rpc.maxLat) / 2.0 rpc.longOffset = (rpc.minLon + rpc.maxLon) / 2.0 - var H []float64 - H = append(H, float64(dem.Dem1KmLT.Elevation(rpc.scene.Meta.Corners.LowerLeft.Longitude, - rpc.scene.Meta.Corners.LowerLeft.Latitude))) - H = append(H, float64(dem.Dem1KmLT.Elevation(rpc.scene.Meta.Corners.LowerRight.Longitude, - rpc.scene.Meta.Corners.LowerRight.Latitude))) - H = append(H, float64(dem.Dem1KmLT.Elevation(rpc.scene.Meta.Corners.UpperLeft.Longitude, - rpc.scene.Meta.Corners.UpperLeft.Latitude))) - H = append(H, float64(dem.Dem1KmLT.Elevation(rpc.scene.Meta.Corners.UpperRight.Longitude, - rpc.scene.Meta.Corners.UpperRight.Latitude))) - H = append(H, float64(dem.Dem1KmLT.Elevation(rpc.longOffset, rpc.latOffset))) - slice.Sort(H, "asc") - rpc.minH = H[0] - rpc.maxH = H[len(H)-1] - rpc.minH, rpc.maxH = dem.Dem1KmLT.MinMaxElevationInRect( rpc.scene.Meta.Corners.UpperLeft.Longitude, rpc.scene.Meta.Corners.UpperLeft.Latitude, rpc.scene.Meta.Corners.LowerRight.Longitude, rpc.scene.Meta.Corners.LowerRight.Latitude, ) + + if rpc.minH < -10.0 { + rpc.minH = 0.0 + } + + if rpc.maxH < -10.0 { + rpc.maxH = 0.0 + } + rpc.heightOffset = (rpc.minH + rpc.maxH) / 2.0 } // 虚拟控制点 func (rpc *RPC) generateVirtualGCP() { - log.Info("Generating virtual GCPs...") - points := gridImage(rpc.gridsize, rpc.gridsize, + log.Infof("Generating virtual GCPs, %d x %d x %d", + rpc.gridsize+1, rpc.gridsize+1, rpc.elevationLayer+1) + points := gridImage2(rpc.gridsize, rpc.gridsize, rpc.scene.Height, rpc.scene.Width, rpc.elevationLayer, int(rpc.minH), int(rpc.maxH)) @@ -143,8 +142,8 @@ func (rpc *RPC) generateVirtualGCP() { func (rpc *RPC) RPC() error { rpc.generateVirtualGCP() n := len(rpc.GCPs) - log.Info("num of virtual GCPs: ", n) + rpc.saveGCP() rowVec := mat.NewVecDense(n, nil) colVec := mat.NewVecDense(n, nil) @@ -160,33 +159,33 @@ func (rpc *RPC) RPC() error { heightVec.SetVec(i, ip.H) } - rowVec, rowOff, rowScale := normalize(rowVec) - colVec, colOff, colScale := normalize(colVec) - latVec, latOff, latScale := normalize(latVec) - lonVec, lonOff, lonScale := normalize(lonVec) - heightVec, heightOff, heightScale := normalize(heightVec) + rpc.saveVec(strings.Replace(rpc.scene.Tiff, ".tiff", ".vec.txt", -1), + rowVec, colVec, latVec, lonVec, heightVec) - rpc.lineOffset = rowOff - rpc.lineScale = rowScale - rpc.sampOffset = colOff - rpc.sampScale = colScale - rpc.latOffset = latOff - rpc.latScale = latScale - rpc.longOffset = lonOff - rpc.longScale = lonScale - rpc.heightOffset = heightOff - rpc.heightScale = heightScale + rpc.lineOffset = float64(rpc.scene.Height / 2) + rpc.lineScale = float64(rpc.scene.Height) + rpc.sampOffset = float64(rpc.scene.Width / 2) + rpc.sampScale = float64(rpc.scene.Width) + rowVec = normalize2(rowVec, rpc.lineOffset, rpc.lineScale) + colVec = normalize2(colVec, rpc.sampOffset, rpc.sampScale) + // rowVec.ScaleVec(1.0/rpc.lineScale, rowVec) + // colVec.ScaleVec(1.0/rpc.sampScale, colVec) - // fmt.Printf("lineOffset: %f, lineScale: %f\n", rpc.lineOffset, rpc.lineScale) - // fmt.Printf("sampOffset: %f, sampScale: %f\n", rpc.sampOffset, rpc.sampScale) - // fmt.Printf("latOffset: %f, latScale: %f\n", rpc.latOffset, rpc.latScale) - // fmt.Printf("longOffset: %f, longScale: %f\n", rpc.longOffset, rpc.longScale) - // fmt.Printf("heightOffset: %f, heightScale: %f\n", rpc.heightOffset, rpc.heightScale) - // fmt.Printf("X0: %f, Y0: %f\n", colVec.At(111, 0), rowVec.At(111, 0)) - // fmt.Printf("lat0: %f, lon0: %f, height0: %f\n", latVec.At(111, 0), lonVec.At(111, 0), heightVec.At(111, 0)) + // rowVec, rpc.lineOffset, rpc.lineScale = normalize(rowVec) + // colVec, rpc.sampOffset, rpc.sampScale = normalize(colVec) + latVec, rpc.latOffset, rpc.latScale = normalize(latVec) + lonVec, rpc.longOffset, rpc.longScale = normalize(lonVec) + heightVec, rpc.heightOffset, rpc.heightScale = normalize(heightVec) + + // rpc.latOffset, rpc.latScale = (rpc.maxLat+rpc.minLat)/2.0, (rpc.maxLat - rpc.minLat) + // rpc.longOffset, rpc.longScale = (rpc.maxLon+rpc.minLon)/2.0, (rpc.maxLon - rpc.minLon) + // rpc.heightOffset, rpc.heightScale = (rpc.maxH+rpc.minH)/2.0+500.0, 500.0 + + rpc.saveVec(strings.Replace(rpc.scene.Tiff, ".tiff", ".vec_norm.txt", -1), + rowVec, colVec, latVec, lonVec, heightVec) // 设计矩阵 B = [ 20个分子系数 19个分母系数 ] - B := buildDesignMatrix(latVec, lonVec, heightVec) + B := buildDesignMatrix(rowVec, latVec, lonVec, heightVec) // x = (B^T * B)^-1 * B^T * l, 其中 x = [a1..a20 b2..b20]^T // 行参数 @@ -204,7 +203,8 @@ func (rpc *RPC) RPC() error { } // 列参数 - K, err := SolveNormalEquation(B, colVec) + D := buildDesignMatrix(colVec, latVec, lonVec, heightVec) + K, err := SolveNormalEquation(D, colVec) if err != nil { return err } @@ -230,7 +230,15 @@ func normalize(v *mat.VecDense) (*mat.VecDense, float64, float64) { return v, vOff, vScale } -func buildDesignMatrix(latVec, lonVec, heightVec *mat.VecDense) *mat.Dense { +func normalize2(v *mat.VecDense, vOff, vScale float64) *mat.VecDense { + for i := 0; i < v.Len(); i++ { + v.SetVec(i, (v.AtVec(i)-vOff)/vScale) + } + + return v +} + +func buildDesignMatrix(vec, latVec, lonVec, heightVec *mat.VecDense) *mat.Dense { n := latVec.Len() // 设计矩阵 B = [ 20个分子系数 19个分母系数 ] B := mat.NewDense(n, 39, nil) @@ -238,6 +246,8 @@ func buildDesignMatrix(latVec, lonVec, heightVec *mat.VecDense) *mat.Dense { P := latVec.AtVec(i) L := lonVec.AtVec(i) H := heightVec.AtVec(i) + r_c := vec.AtVec(i) + B.Set(i, 0, 1) B.Set(i, 1, L) B.Set(i, 2, P) @@ -258,25 +268,25 @@ func buildDesignMatrix(latVec, lonVec, heightVec *mat.VecDense) *mat.Dense { B.Set(i, 17, L*L*H) B.Set(i, 18, P*P*H) B.Set(i, 19, H*H*H) - B.Set(i, 20, -L) - B.Set(i, 21, -P) - B.Set(i, 22, -H) - B.Set(i, 23, -L*P) - B.Set(i, 24, -L*H) - B.Set(i, 25, -P*H) - B.Set(i, 26, -L*L) - B.Set(i, 27, -P*P) - B.Set(i, 28, -H*H) - B.Set(i, 29, -P*L*H) - B.Set(i, 30, -L*L*L) - B.Set(i, 31, -L*P*P) - B.Set(i, 32, -L*H*H) - B.Set(i, 33, -L*L*P) - B.Set(i, 34, -P*P*P) - B.Set(i, 35, -P*H*H) - B.Set(i, 36, -L*L*H) - B.Set(i, 37, -P*P*H) - B.Set(i, 38, -H*H*H) + B.Set(i, 20, -L*r_c) + B.Set(i, 21, -P*r_c) + B.Set(i, 22, -H*r_c) + B.Set(i, 23, -L*P*r_c) + B.Set(i, 24, -L*H*r_c) + B.Set(i, 25, -P*H*r_c) + B.Set(i, 26, -L*L*r_c) + B.Set(i, 27, -P*P*r_c) + B.Set(i, 28, -H*H*r_c) + B.Set(i, 29, -P*L*H*r_c) + B.Set(i, 30, -L*L*L*r_c) + B.Set(i, 31, -L*P*P*r_c) + B.Set(i, 32, -L*H*H*r_c) + B.Set(i, 33, -L*L*P*r_c) + B.Set(i, 34, -P*P*P*r_c) + B.Set(i, 35, -P*H*H*r_c) + B.Set(i, 36, -L*L*H*r_c) + B.Set(i, 37, -P*P*H*r_c) + B.Set(i, 38, -H*H*H*r_c) } return B @@ -303,7 +313,24 @@ func SolveNormalEquation(A *mat.Dense, b *mat.VecDense) ([]float64, error) { // 求解 (A^T * A)^-1 * (A^T * b) var AtInv mat.Dense err := AtInv.Inverse(&At) + if err != nil { + // 岭估计方法调整法方程状态，使得矩阵非奇异，最小二乘平差可以收敛 + r, c := At.Dims() + log.Infof("cannot inverse design matrix(%d*%d): %v", r, c, err) + log.Info("try to adjust design matrix with +kI, k=0.0000001") + k := 0.0000001 // [0.00000005, 0.000005] + I := mat.NewDiagDense(r, nil) + for i := 0; i < r; i++ { + I.SetDiag(i, k) + } + At.Add(&At, I) + + err = AtInv.Inverse(&At) + } + + if err != nil { + log.Infof("cannot inverse design matrix: %v, try SVD method", err) // 计算矩阵的 SVD 分解 var svd mat.SVD ok := svd.Factorize(&At, mat.SVDThin) @@ -403,5 +430,53 @@ func (rpc *RPC) SaveRpb() error { } defer f.Close() _, err = f.WriteString(model) + if err != nil { + log.Errorf("Failed to write RPC model file: %v", err) + } + f.Sync() return err } + +func (rpc *RPC) saveGCP() error { + name := strings.Replace(rpc.scene.Tiff, ".tiff", ".gcp.geojson", -1) + log.Infof("save gcp to %s", name) + f, err := os.Create(name) + if err != nil { + log.Errorf("Failed to create GCP file: %v", err) + return err + } + defer f.Close() + + var gcp geojson.FeatureCollection + for _, p := range rpc.GCPs { + point := orb.Point{p.L, p.P} + feature := geojson.NewFeature(point) + feature.Properties = map[string]interface{}{ + "H": p.H, + "Y": p.Y, + "X": p.X, + } + gcp.Features = append(gcp.Features, feature) + } + data, _ := json.Marshal(gcp) + f.Write(data) + f.Sync() + + return nil +} + +func (rpc *RPC) saveVec(name string, rowVec, colVec, latVec, lonVec, heightVec *mat.VecDense) error { + f, err := os.Create(name) + + if err != nil { + log.Errorf("Failed to create vec file: %v", err) + return err + } + defer f.Close() + + for i := 0; i < rowVec.Len(); i++ { + f.WriteString(fmt.Sprintf("%.8f\t%.8f\t%.8f\t%.8f\t%.8f\n", + rowVec.AtVec(i), colVec.AtVec(i), latVec.AtVec(i), lonVec.AtVec(i), heightVec.AtVec(i))) + } + return nil +} diff --git a/pkg/utils/tiff.go b/pkg/utils/tiff.go index f62ef87..e0c3de0 100644 --- a/pkg/utils/tiff.go +++ b/pkg/utils/tiff.go @@ -23,7 +23,7 @@ func SavePanToGDALGTiff(pan gocv.Mat, topLeftX, topLeftY float64, tiffFile strin } defer ds.Close() - setGeoTransform(ds, topLeftX, topLeftY, resolution) + // setGeoTransform(ds, topLeftX, topLeftY, resolution) ds.SetMetadata("NBITS", "16") // 将通道的数据转换为uint16数组