Do RRC by moment matching method

cmd/proc.go

@@ -13,10 +13,11 @@ import (
 )
 
 var (
-	params    producer.Params
-	saveStrip bool
-	doRRC     bool
-	lutDir    string
+	params           producer.Params
+	saveStrip        bool
+	doRRC            bool
+	doMomentMatching bool
+	lutDir           string
 )
 
 var procCmd = &cobra.Command{
@@ -50,6 +51,10 @@ var procCmd = &cobra.Command{
 			reg.DoRRC(lutDir)
 		}
 
+		if doMomentMatching {
+			reg.DoMomentMatching()
+		}
+
 		if err := reg.DoPhaseCorrelation(); err != nil {
 			logrus.Fatal(err)
 		}
@@ -99,6 +104,7 @@ func init() {
 	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().BoolVarP(&doMomentMatching, "mm", "", false, "do moment matching")
 }
 
 func initParams() producer.Params {

cmd/rrc.go

@@ -6,8 +6,9 @@ import (
 )
 
 var (
-	panDS string
-	mssDS string
+	panDS  string
+	mssDS  string
+	lutOut string
 )
 
 var rrcCmd = &cobra.Command{
@@ -15,7 +16,7 @@ var rrcCmd = &cobra.Command{
 	Short: "Run RRC algorithm on an image",
 	Long:  `Run RRC algorithm on an image`,
 	Run: func(cmd *cobra.Command, args []string) {
-		rrc := rrc.NewRRC()
+		rrc := rrc.NewRRC(lutOut)
 		rrc.StatisticalPAN(panDS)
 		rrc.StatisticalMSS(mssDS)
 	},
@@ -24,6 +25,7 @@ var rrcCmd = &cobra.Command{
 func init() {
 	rrcCmd.Flags().StringVarP(&panDS, "pan-dataset", "p", "data/RAW/pan.txt", "path to pan dataset")
 	rrcCmd.Flags().StringVarP(&mssDS, "mss-dataset", "m", "data/RAW/mss-none.txt", "path to mss dataset")
+	rrcCmd.Flags().StringVarP(&lutOut, "lut-out-dir", "o", "data/rrc", "path to output lut file")
 
 	rootCmd.AddCommand(rrcCmd)
 }

pkg/producer/rrc.go

@@ -40,3 +40,12 @@ func (r *Registrator) DoRRC(lutDir string) error {
 
 	return nil
 }
+
+func (r *Registrator) DoMomentMatching() error {
+	rrc.DoMomentMatching(r.PanImage)
+	for i := 0; i < 4; i++ {
+		rrc.DoMomentMatching(r.MssImages[i])
+	}
+
+	return nil
+}

pkg/rrc/moment_matching.go

@@ -0,0 +1,56 @@
+package rrc
+
+import (
+	"math"
+
+	log "github.com/sirupsen/logrus"
+	"gocv.io/x/gocv"
+)
+
+// Destriping multisensor imagery with moment matching
+
+func DoMomentMatching(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))
+	}
+
+	// 列参考值和列参考标准差
+	var mu float64
+	var sig float64
+	for x := 0; x < probes; x++ {
+		mu += means[x]
+		sig += stds[x]
+	}
+	mu = mu / float64(probes)
+	sig = sig / float64(probes)
+
+	log.Printf("mean reference value: %f, std reference value: %f", mu, sig)
+	// 修正 DN_adjusted[i] = (DN[i] - means[i]) *sig/stds[i]+mu
+	for x := 0; x < originalImg.Cols(); x++ {
+		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))
+		}
+	}
+}

pkg/rrc/rrc.go

@@ -4,6 +4,7 @@ import (
 	"bufio"
 	"fmt"
 	"os"
+	"path/filepath"
 	"runtime"
 	"sync"
 
@@ -21,17 +22,18 @@ const (
 )
 
 type RRC struct {
+	LUTOutDir  string
 	Histograms [5]ProbeHistogram
 }
 
-func NewRRC() *RRC {
-	r := RRC{}
+func NewRRC(dir string) *RRC {
+	r := RRC{LUTOutDir: dir}
 	r.Histograms[0].init(PANCameraProbeNum)
 	for i := 1; i < 5; i++ {
 		r.Histograms[i].init(MSSCameraProbeNum)
 	}
 
-	os.MkdirAll("data/rrc", 0755)
+	os.MkdirAll(r.LUTOutDir, 0755)
 
 	return &r
 }
@@ -110,7 +112,7 @@ func (rrc *RRC) StatisticalPAN(dsfile string) error {
 	log.Println("compute PAN histogram...")
 	rrc.Histograms[0].compute()
 	log.Println("save PAN gray table matrix.")
-	rrc.Histograms[0].saveLUT("data/rrc/B0.LUT")
+	rrc.Histograms[0].saveLUT(filepath.Join(rrc.LUTOutDir, "B0.LUT"))
 
 	return nil
 }
@@ -189,7 +191,7 @@ func (rrc *RRC) StatisticalMSS(dsfile string) error {
 		log.Println("compute MSS histogram...")
 		rrc.Histograms[i].compute()
 		log.Println("save MSS gray table matrix.")
-		rrc.Histograms[i].saveLUT(fmt.Sprintf("data/rrc/B%d.LUT", i))
+		rrc.Histograms[i].saveLUT(filepath.Join(rrc.LUTOutDir, fmt.Sprintf("B%d.LUT", i)))
 	}
 
 	return nil