pansharpen by ESRI method

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

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

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
+}

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

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
+}