From 6e5ffe1ab2f3692a6a0c1cd887993b24c00d1b22 Mon Sep 17 00:00:00 2001
From: nuknal <linktolane@gmail.com>
Date: Tue, 18 Jun 2024 16:19:33 +0800
Subject: [PATCH] =?UTF-8?q?=E7=9B=B8=E9=82=BB=E5=88=97=E5=9D=87=E8=A1=A1?=
 =?UTF-8?q?=E7=AE=97=E6=B3=95=EF=BC=88=E5=BC=A0=E5=85=B5=EF=BC=8C2006?=
 =?UTF-8?q?=EF=BC=89?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 cmd/proc.go                | 10 ++++-----
 pkg/producer/rrc.go        |  2 +-
 pkg/rrc/moment_matching.go | 42 ++++++++++++++++++++++++++++++++++++--
 3 files changed, 46 insertions(+), 8 deletions(-)

diff --git a/cmd/proc.go b/cmd/proc.go
index 34b0e4d..fb597ef 100644
--- a/cmd/proc.go
+++ b/cmd/proc.go
@@ -15,7 +15,7 @@ import (
 var (
 	params           producer.Params
 	saveStrip        bool
-	doRRC            bool
+	doLUTRRC         bool
 	doMomentMatching bool
 	lutDir           string
 )
@@ -47,8 +47,8 @@ var procCmd = &cobra.Command{
 		godal.RegisterAll()
 		os.MkdirAll(params.OutputDir, 0755)
 
-		if doRRC {
-			reg.DoRRC(lutDir)
+		if doLUTRRC {
+			reg.DoRRCbyLUT(lutDir)
 		}
 
 		if doMomentMatching {
@@ -102,8 +102,8 @@ func init() {
 	procCmd.Flags().BoolVarP(&params.SubScenes, "sub-scenes", "", false, "process sub-scenes")
 	procCmd.Flags().BoolVarP(&saveStrip, "save-strip", "", false, "save original and registered images as GDAL GTiff")
 	procCmd.Flags().StringVarP(&paramsXML, "params", "x", "", "params xml file path")
-	procCmd.Flags().StringVarP(&lutDir, "lut", "l", "data/lut", "LUT directory")
-	procCmd.Flags().BoolVarP(&doRRC, "rrc", "", false, "do RRC")
+	procCmd.Flags().StringVarP(&lutDir, "lut-dir", "", "data/lut", "LUT directory")
+	procCmd.Flags().BoolVarP(&doLUTRRC, "lut-rrc", "", false, "do RRC with LUT method")
 	procCmd.Flags().BoolVarP(&doMomentMatching, "mm", "", false, "do moment matching")
 }
 
diff --git a/pkg/producer/rrc.go b/pkg/producer/rrc.go
index 36408b4..11d9e1d 100644
--- a/pkg/producer/rrc.go
+++ b/pkg/producer/rrc.go
@@ -8,7 +8,7 @@ import (
 	"starwiz.cn/sjy01/image-proc/pkg/rrc"
 )
 
-func (r *Registrator) DoRRC(lutDir string) error {
+func (r *Registrator) DoRRCbyLUT(lutDir string) error {
 	logrus.Printf("try to do RRC [%s]...", lutDir)
 	lutPAN, err := rrc.LoadLUT(filepath.Join(lutDir, "B0.LUT"), 9344)
 	if err != nil {
diff --git a/pkg/rrc/moment_matching.go b/pkg/rrc/moment_matching.go
index 1955173..7b64d33 100644
--- a/pkg/rrc/moment_matching.go
+++ b/pkg/rrc/moment_matching.go
@@ -7,8 +7,7 @@ import (
 	"gocv.io/x/gocv"
 )
 
-// Destriping multisensor imagery with moment matching
-
+// Destriping multisensor imagery with moment matching [Gadallah, 2000]
 func DoMomentMatching(originalImg gocv.Mat) {
 	probes := originalImg.Cols()
 	log.Printf("do moment matching for %d probes", probes)
@@ -54,3 +53,42 @@ func DoMomentMatching(originalImg gocv.Mat) {
 		}
 	}
 }
+
+// 相邻列均衡算法（张兵，2006）- 实测结果不如上面的矩匹配算法
+func DoMomentMatching2006(originalImg gocv.Mat) {
+	probes := originalImg.Cols()
+	log.Printf("do moment matching for %d probes", probes)
+	// 第i个探元的像元均值和标准差
+	means := make([]float64, probes)
+	stds := make([]float64, probes)
+	// 计算每个探元的均值和标准差
+	for x := 0; x < originalImg.Cols(); x++ {
+		var total int64
+		n := originalImg.Rows()
+		var dn uint16
+		for y := 0; y < n; y++ {
+			dn = uint16(originalImg.GetShortAt(y, x))
+			total += int64(dn)
+		}
+		means[x] = float64(total) / float64(n)
+
+		var a float64
+		for y := 0; y < n; y++ {
+			dn = uint16(originalImg.GetShortAt(y, x))
+			a += math.Pow(float64(dn)-means[x], 2)
+		}
+		stds[x] = math.Sqrt(a / float64(n))
+	}
+
+	// 修正 DN_adjusted[i] = (DN[i] - means[i]) *sig/stds[i]+mu
+	for x := 1; x < originalImg.Cols()-1; x++ {
+		// 列参考值和列参考标准差
+		mu := (means[x-1] + means[x+1]/2 + means[x]) / 2
+		sig := (stds[x-1] + stds[x+1]/2 + stds[x]) / 2
+		for y := 0; y < originalImg.Rows(); y++ {
+			dn := uint16(originalImg.GetShortAt(y, x))
+			dn_adjusted := uint16((float64(dn)-means[x])*sig/stds[x] + mu)
+			originalImg.SetShortAt(y, x, int16(dn_adjusted))
+		}
+	}
+}