IHS pansharpen

cmd/fus.go

@@ -16,6 +16,7 @@ var (
 	mssImage  string
 	outputDir string
 	fusReport string
+	fusMethod int
 )
 
 var fusCmd = &cobra.Command{
@@ -38,7 +39,9 @@ var fusCmd = &cobra.Command{
 		// 	logrus.Fatal(err)
 		// }
 
-		err := fusion.Pansharpen(panImage, mssImage, filepath.Join(outputDir, fusedTiff), fusion.ESRI)
+		err := fusion.Pansharpen(panImage, mssImage,
+			filepath.Join(outputDir, fusedTiff),
+			fusion.PansharpenMethod(fusMethod))
 		if err != nil {
 			logrus.Fatal(err)
 		}
@@ -63,6 +66,7 @@ func init() {
 	fusCmd.Flags().StringVarP(&mssImage, "mss", "m", "", "path to the MSS image")
 	fusCmd.Flags().StringVarP(&outputDir, "output", "o", "", "path to the output directory")
 	fusCmd.Flags().StringVarP(&paramsXML, "params", "x", "", "path to the XML file containing parameters for GDAL_Pansharpen")
+	fusCmd.Flags().IntVarP(&fusMethod, "method", "t", int(fusion.ESRI), "method to use for fusion (5: IHS, 3: ESRI)")
 }
 
 func initFUSParams() {

pkg/fusion/fusion.go

@@ -20,6 +20,8 @@ const (
 	SimpleMean
 	ESRI
 	GDALPansharpen
+	IHS
+	GIHS
 )
 
 func Pansharpen(panTif, mssTif, fusTif string, method PansharpenMethod) error {
@@ -50,6 +52,8 @@ func Pansharpen(panTif, mssTif, fusTif string, method PansharpenMethod) error {
 	switch method {
 	case ESRI:
 		fus = PansharpenESRI(*pan, BGRI[0], BGRI[1], BGRI[2], BGRI[3])
+	case IHS:
+		fus = PansharpenIHS(*pan, BGRI[0], BGRI[1], BGRI[2], BGRI[3])
 	default:
 		fus = PansharpenESRI(*pan, BGRI[0], BGRI[1], BGRI[2], BGRI[3])
 	}
@@ -97,7 +101,11 @@ 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 PansharpenBrovey(pan, B, G, R, I gocv.Mat, w float32) gocv.Mat {
+	// DNF = (PAN - w*I)/(R*w + G*w + B*w)
+	var fus gocv.Mat
+	return fus
+}
 
 func PansharpenSimpleMean(pan, B, G, R, I gocv.Mat) error { return nil }
 
@@ -107,6 +115,13 @@ func PansharpenESRI(pan, B, G, R, I gocv.Mat) gocv.Mat {
 	ADJ := gocv.NewMat()
 	defer ADJ.Close()
 
+	for y := 0; y < pan.Rows(); y++ {
+		for x := 0; x < pan.Cols(); x++ {
+			v := uint16((R.GetFloatAt(y, x) + G.GetFloatAt(y, x) + B.GetFloatAt(y, x) + I.GetFloatAt(y, x)) / 4)
+			bgrn_mean.SetShortAt(y, x, int16(v))
+		}
+	}
+
 	bgrn_mean.MultiplyFloat(-1)
 	gocv.Add(pan, bgrn_mean, &ADJ)
 	gocv.Add(B, ADJ, &B)
@@ -130,3 +145,52 @@ func PansharpenESRI(pan, B, G, R, I gocv.Mat) gocv.Mat {
 
 	return mss
 }
+
+func PansharpenIHS(pan, B, G, R, I gocv.Mat) gocv.Mat {
+	mss := gocv.NewMat()
+	gocv.Merge([]gocv.Mat{B, G, R}, &mss)
+	mssHLS := gocv.NewMat()
+	defer mssHLS.Close()
+
+	mss.ConvertTo(&mss, gocv.MatTypeCV32FC3)
+	gocv.CvtColor(mss, &mssHLS, gocv.ColorBGRToHLS)
+	hlsChannels := gocv.Split(mssHLS)
+
+	// 将 PAN 图像转换为与 HLS 亮度分量相同的范围
+	panFloat := gocv.NewMat()
+	defer panFloat.Close()
+	pan.ConvertTo(&panFloat, gocv.MatTypeCV32F)
+	gocv.Normalize(panFloat, &panFloat, 0, 65535, gocv.NormMinMax)
+
+	hlsFloat := gocv.NewMat()
+	defer hlsFloat.Close()
+	hlsChannels[1].ConvertTo(&hlsFloat, gocv.MatTypeCV32F)
+	gocv.Normalize(hlsFloat, &hlsFloat, 0, 65535, gocv.NormMinMax)
+
+	// 替换亮度分量
+	hlsChannels[1] = panFloat
+
+	// 合并通道
+	gocv.Merge(hlsChannels, &mssHLS)
+
+	// 将图像从 HLS 转换回 BGR
+	fusedImage := gocv.NewMat()
+
+	log.Println("cvt color from HLS to BGR")
+	gocv.CvtColor(mssHLS, &fusedImage, gocv.ColorHLSToBGR)
+
+	channels := gocv.Split(fusedImage)
+	defer channels[0].Close() // B
+	defer channels[1].Close() // G
+	defer channels[2].Close() // R
+
+	if config.GCONFIG.CoRegistration.FUSBandOrder == "RGBI" || config.GCONFIG.CoRegistration.FUSBandOrder == "RGB" {
+		gocv.Merge([]gocv.Mat{channels[2], channels[1], channels[0]}, &fusedImage)
+	} else {
+		gocv.Merge([]gocv.Mat{channels[0], channels[1], channels[2]}, &fusedImage)
+	}
+
+	fusedImage.ConvertTo(&fusedImage, gocv.MatTypeCV16UC3)
+
+	return fusedImage
+}