diff --git a/cmd/fus.go b/cmd/fus.go index b1a6812..327209d 100644 --- a/cmd/fus.go +++ b/cmd/fus.go @@ -7,6 +7,7 @@ import ( "github.com/airbusgeo/godal" "github.com/sirupsen/logrus" "github.com/spf13/cobra" + "starwiz.cn/sjy01/image-proc/pkg/fusion" producer "starwiz.cn/sjy01/image-proc/pkg/producer" ) @@ -30,12 +31,18 @@ var fusCmd = &cobra.Command{ fusReport = filepath.Join(outputDir, id+"_report.xml") } - err := producer.GDALPansharpen(panImage, mssImage, filepath.Join(outputDir, fusedTiff)) + godal.RegisterAll() + + // err := producer.GDALPansharpen(panImage, mssImage, filepath.Join(outputDir, fusedTiff)) + // if err != nil { + // logrus.Fatal(err) + // } + + err := fusion.Pansharpen(panImage, mssImage, filepath.Join(outputDir, fusedTiff), fusion.ESRI) if err != nil { logrus.Fatal(err) } - godal.RegisterAll() producer.GTiffToJPG(filepath.Join(outputDir, fusedTiff), filepath.Join(outputDir, id+".jpg"), "FUS", true) var report producer.Report diff --git a/config/config.yaml b/config/config.yaml index 00572f8..d220b8a 100644 --- a/config/config.yaml +++ b/config/config.yaml @@ -11,8 +11,8 @@ coregistration: fus_band_order: "RGB" radiation: - pan_remove_hf_noise: false - mss_remove_hf_noise: false + pan_remove_hf_noise: true + mss_remove_hf_noise: true hf_radius_ratio: 0.49 min_hist_level: 0.3 max_hist_level: 0.6 diff --git a/pkg/fusion/fusion.go b/pkg/fusion/fusion.go new file mode 100644 index 0000000..b67626d --- /dev/null +++ b/pkg/fusion/fusion.go @@ -0,0 +1,123 @@ +package fusion + +import ( + "image" + "os" + "strings" + + "github.com/airbusgeo/godal" + log "github.com/sirupsen/logrus" + "gocv.io/x/gocv" + "starwiz.cn/sjy01/image-proc/pkg/config" + "starwiz.cn/sjy01/image-proc/pkg/utils" +) + +type PansharpenMethod int + +const ( + SimpleBrovey PansharpenMethod = iota + Brovey + SimpleMean + ESRI + GDALPansharpen +) + +func Pansharpen(panTif, mssTif, fusTif string, method PansharpenMethod) error { + log.Infof("Pansharpening %s and %s", panTif, mssTif) + pan, err := utils.ReadTifftoMat(panTif) + if err != nil { + return err + } + defer pan.Close() + + log.Infof("PAN row %d, col %d", pan.Rows(), pan.Cols()) + + mss, err := utils.ReadTifftoMat(mssTif) + if err != nil { + return err + } + defer mss.Close() + + log.Infof("MSS bands %d row %d, col %d", mss.Channels(), mss.Rows(), mss.Cols()) + + mssResized := gocv.NewMat() + defer mssResized.Close() + gocv.Resize(*mss, &mssResized, image.Point{X: pan.Cols(), Y: pan.Rows()}, 0, 0, gocv.InterpolationCubic) + + BGRI := gocv.Split(mssResized) + + var fus gocv.Mat + switch method { + case ESRI: + fus = PansharpenESRI(*pan, BGRI[0], BGRI[1], BGRI[2], BGRI[3]) + default: + fus = PansharpenESRI(*pan, BGRI[0], BGRI[1], BGRI[2], BGRI[3]) + } + + defer fus.Close() + + var mt gocv.MatType + switch fus.Channels() { + case 1: + mt = gocv.MatTypeCV16UC1 + case 3: + mt = gocv.MatTypeCV16UC3 + case 4: + mt = gocv.MatTypeCV16UC4 + } + fus.ConvertTo(&fus, mt) + + var colorInterps []godal.ColorInterp + if config.GCONFIG.CoRegistration.FUSBandOrder == "RGB" { + colorInterps = []godal.ColorInterp{godal.CIRed, godal.CIGreen, godal.CIBlue, godal.CIUndefined} + } else { + colorInterps = []godal.ColorInterp{godal.CIBlue, godal.CIGreen, godal.CIRed, godal.CIUndefined} + } + err = utils.SaveBGRToGDALGTiff(fus, fus.Channels(), 0, 0, 1.3, colorInterps, fusTif) + + for i := 0; i < len(BGRI); i++ { + BGRI[i].Close() + } + + // use rpc of PAN directly + rpb := strings.Replace(panTif, ".tiff", ".rpb", -1) + if _, err := os.Stat(rpb); err == nil { + fusRPB := strings.Replace(fusTif, ".tiff", ".rpb", -1) + if input, err := os.ReadFile(rpb); err == nil { + os.WriteFile(fusRPB, input, 0644) + } + } + + return err +} + +func PansharpenSimpleBrovey(pan, B, G, R, I gocv.Mat) error { + return nil +} + +func PansharpenBrovey(pan, B, G, R, I gocv.Mat) error { return nil } + +func PansharpenSimpleMean(pan, B, G, R, I gocv.Mat) error { return nil } + +func PansharpenESRI(pan, B, G, R, I gocv.Mat) gocv.Mat { + bgrn_mean := gocv.NewMatWithSize(pan.Rows(), pan.Cols(), pan.Type()) + defer bgrn_mean.Close() + ADJ := gocv.NewMat() + defer ADJ.Close() + + bgrn_mean.MultiplyFloat(-1) + gocv.Add(pan, bgrn_mean, &ADJ) + gocv.Add(B, ADJ, &B) + gocv.Add(G, ADJ, &G) + gocv.Add(R, ADJ, &R) + gocv.Add(I, ADJ, &I) + + mss := gocv.NewMat() + if config.GCONFIG.CoRegistration.FUSBandOrder == "RGB" { + gocv.Merge([]gocv.Mat{R, G, B}, &mss) + } else { + gocv.Merge([]gocv.Mat{B, G, R}, &mss) + } + + return mss +} diff --git a/pkg/producer/rpc_helper.go b/pkg/producer/rpc_helper.go index 693a915..c87dbca 100644 --- a/pkg/producer/rpc_helper.go +++ b/pkg/producer/rpc_helper.go @@ -128,7 +128,6 @@ func solveCoefficients(f, latVec, lonVec, heightVec *mat.VecDense, method SolveM 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)) @@ -137,7 +136,6 @@ func solveCoefficients(f, latVec, lonVec, heightVec *mat.VecDense, method SolveM x0 = x1 iterations++ - fmt.Printf("squared error+lambda*coeffs: %.8f\n", coeffsSquared) vx = append(vx, &VX{v: errorSquared, x: x1}) if errorSquared < 0.001 { break diff --git a/pkg/utils/tiff.go b/pkg/utils/tiff.go index 76a9bb0..e95c4b2 100644 --- a/pkg/utils/tiff.go +++ b/pkg/utils/tiff.go @@ -1,6 +1,8 @@ package utils import ( + "os" + "path/filepath" "time" "github.com/airbusgeo/godal" @@ -12,7 +14,7 @@ import ( func SavePanToGDALGTiff(pan gocv.Mat, topLeftX, topLeftY float64, tiffFile string, resolution float64) error { // log.Println("Saving PAN image to TIFF file:", tiffFile) - + os.MkdirAll(filepath.Dir(tiffFile), 0755) width := pan.Cols() height := pan.Rows() @@ -60,6 +62,8 @@ func SaveBGRToGDALGTiff(bgr gocv.Mat, width := bgr.Cols() height := bgr.Rows() + os.MkdirAll(filepath.Dir(tiffFile), 0755) + // 创建一个二维切片来存储图像数据 data := make([][]uint16, bands) for i := range data { @@ -174,3 +178,46 @@ func ReadTiff(tifFile string) ([][][]float32, error) { } return data, nil } + +func ReadTifftoMat(ftiff string) (*gocv.Mat, error) { + ds, err := godal.Open(ftiff) + if err != nil { + log.Printf("Error opening TIFF file %s: %v", ftiff, err) + return nil, err + } + defer ds.Close() + + bands := ds.Bands() + + // 获取图像大小 + width := bands[0].Structure().SizeX + height := bands[0].Structure().SizeY + bandsCnt := len(bands) + + var img gocv.Mat + if bandsCnt == 1 { + img = gocv.NewMatWithSize(height, width, gocv.MatTypeCV16UC1) + } else { + img = gocv.NewMatWithSize(height, width, gocv.MatTypeCV16UC4) + } + + // 读取每个波段并存储到图像矩阵中 + for i := 0; i < bandsCnt; i++ { + band := bands[i] + data := make([]uint16, width*height) + err = band.Read(0, 0, data, width, height) + if err != nil { + return nil, err + } + + // 将 16 位数据存储到图像的相应通道 + for y := 0; y < height; y++ { + for x := 0; x < width; x++ { + value := data[y*width+x] + img.SetShortAt(y, x*bandsCnt+i, int16(value)) + } + } + } + + return &img, nil +}