histogram RRC

cmd/proc.go

@@ -44,6 +44,8 @@ var procCmd = &cobra.Command{
 		godal.RegisterAll()
 		os.MkdirAll(params.OutputDir, 0755)
 
+		reg.DoRRC()
+
 		if err := reg.DoPhaseCorrelation(); err != nil {
 			logrus.Fatal(err)
 		}

cmd/rrc.go

@@ -0,0 +1,28 @@
+package main
+
+import (
+	"github.com/spf13/cobra"
+	"starwiz.cn/sjy01/image-proc/pkg/rrc"
+)
+
+var (
+	panDS string
+	mssDS string
+)
+
+var rrcCmd = &cobra.Command{
+	Use:   "rrc",
+	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.StatisticalPAN(panDS)
+	},
+}
+
+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.txt", "path to mss dataset")
+
+	rootCmd.AddCommand(rrcCmd)
+}

pkg/producer/rrc.go

@@ -0,0 +1,41 @@
+package producer
+
+import (
+	"fmt"
+
+	"github.com/sirupsen/logrus"
+	"starwiz.cn/sjy01/image-proc/pkg/rrc"
+)
+
+func (r *Registrator) DoRRC() error {
+	logrus.Println("try to do RRC...")
+	tablePAN, err := rrc.LoadGrayTableMatrix("data/rrc/pan_gray_table.dat")
+	if err != nil {
+		logrus.Error("load pan gray table failed")
+		return err
+	}
+
+	for y := 0; y < r.PanImage.Rows(); y++ {
+		for x := 0; x < r.PanImage.Cols(); x++ {
+			newGray := tablePAN.At(x, int(uint16(r.PanImage.GetShortAt(y, x))))
+			r.PanImage.SetShortAt(y, x, int16(newGray))
+		}
+	}
+
+	for i := 0; i < 4; i++ {
+		tableMSS, err := rrc.LoadGrayTableMatrix(fmt.Sprintf("data/rrc/mss%d_gray_table.dat", i+1))
+		if err != nil {
+			logrus.Error("load mss gray table failed")
+			return err
+		}
+
+		for y := 0; y < r.MssImages[i].Rows(); y++ {
+			for x := 0; x < r.MssImages[i].Cols(); x++ {
+				newGray := tableMSS.At(x, int(uint16(r.MssImages[i].GetShortAt(y, x))))
+				r.MssImages[i].SetShortAt(y, x, int16(newGray))
+			}
+		}
+	}
+
+	return nil
+}

pkg/rrc/histogram.go

@@ -0,0 +1,171 @@
+package rrc
+
+import (
+	"fmt"
+	"math"
+	"os"
+
+	"github.com/sirupsen/logrus"
+	log "github.com/sirupsen/logrus"
+	"gocv.io/x/gocv"
+	"gonum.org/v1/gonum/mat"
+)
+
+type ProbeHistogram struct {
+	probes int         // 探元数
+	N_i    []int64     // N_i 探元像素总数 PAN 0-9343 MSS 0-2335
+	n_ik   [][]int64   // n_ik 第i探元灰度等级为k的像素数统计 PAN 9343 x 65536  MSS 2335 x 65536
+	p_ik   [][]float64 // p_ik = n_ik / N_i 探元灰度概率密度 PAN 9343 x 65536  MSS 2335 x 65536
+	m_l    []int64     // 具有灰度等级l的像素总数 l = 0-65535
+	M      int64       // 参与直方图统计的总像素数
+	P_l    []float64   // P_l = m_l/M 所有探元的期望直方图灰度等级为l的概率密度
+	S_ik   [][]float64 // S_ik = sum(p_ij),j=0..k 第i个探元直方图灰度等级k的累积概率密度
+	V_l    []float64   // V_l = sum(P_j),j=0..l // 期望直方图灰度级l对应的累积概率密度
+	Tmat   *mat.Dense  // 第i个像元的j灰度等级对应的新的灰度值，用于修正图像
+}
+
+func (hist *ProbeHistogram) init(width int) {
+	hist.probes = width
+	hist.M = 0
+	hist.N_i = make([]int64, width)
+	hist.n_ik = make([][]int64, width)
+	hist.p_ik = make([][]float64, width)
+	hist.m_l = make([]int64, MaxGrayLevel)
+	hist.P_l = make([]float64, MaxGrayLevel)
+	hist.S_ik = make([][]float64, width)
+	hist.V_l = make([]float64, MaxGrayLevel)
+
+	for i := 0; i < width; i++ {
+		hist.n_ik[i] = make([]int64, MaxGrayLevel)
+		hist.p_ik[i] = make([]float64, MaxGrayLevel)
+		hist.S_ik[i] = make([]float64, MaxGrayLevel)
+	}
+
+	hist.Tmat = mat.NewDense(width, MaxGrayLevel, nil)
+}
+
+func (hist *ProbeHistogram) statistical(img gocv.Mat) error {
+	hist.M += int64(img.Rows() * img.Cols())
+	fmt.Println("Hist.M:", hist.M)
+
+	// 探元i像素总数
+	for i := 0; i < hist.probes; i++ {
+		hist.N_i[i] += int64(img.Rows())
+	}
+
+	// 探元i灰度等级k的像素数统计
+	for y := 0; y < img.Rows(); y++ {
+		for x := 0; x < img.Cols(); x++ {
+			gray := uint16(img.GetShortAt(x, y))
+			hist.n_ik[x][gray]++
+		}
+	}
+
+	// 灰度等级l的像素总数
+	for gray := 0; gray < MaxGrayLevel; gray++ {
+		for i := 0; i < hist.probes; i++ {
+			hist.m_l[gray] += int64(hist.n_ik[i][gray])
+		}
+	}
+
+	return nil
+}
+
+func (hist *ProbeHistogram) compute() {
+	// 探元i灰度概率密度
+	for i := 0; i < hist.probes; i++ {
+		for k := 0; k < MaxGrayLevel; k++ {
+			hist.p_ik[i][k] = float64(hist.n_ik[i][k]) / float64(hist.N_i[i])
+		}
+	}
+
+	// 所有探元的期望直方图灰度等级为l的概率密度
+	for gray := 0; gray < MaxGrayLevel; gray++ {
+		hist.P_l[gray] = float64(hist.m_l[gray]) / float64(hist.M)
+	}
+
+	// 第i个探元直方图灰度等级k的累积概率密度
+	for i := 0; i < hist.probes; i++ {
+		for k := 0; k < MaxGrayLevel; k++ {
+			hist.S_ik[i][k] = 0
+			for j := 0; j <= k; j++ {
+				hist.S_ik[i][k] += hist.p_ik[i][j]
+			}
+		}
+	}
+
+	// 期望直方图灰度级l对应的累积概率密度
+	for gray := 0; gray < MaxGrayLevel; gray++ {
+		hist.V_l[gray] = 0
+		for j := 0; j <= gray; j++ {
+			hist.V_l[gray] += hist.P_l[j]
+		}
+	}
+
+	// 生成查找表
+	nT := 0
+	for i := 0; i < hist.probes; i++ {
+		for k := 0; k < MaxGrayLevel; k++ {
+			for l := 0; l < MaxGrayLevel-1; l++ {
+				if hist.S_ik[i][k] >= hist.V_l[l] && hist.S_ik[i][k] <= hist.V_l[l+1] {
+					nT += 1
+					if math.Abs(hist.S_ik[i][k]-hist.V_l[l]) <= math.Abs(hist.S_ik[i][k]-hist.V_l[l+1]) {
+						hist.Tmat.Set(i, k, float64(l))
+					} else {
+						hist.Tmat.Set(i, k, float64(l+1))
+					}
+				}
+			}
+		}
+	}
+
+	if nT != hist.probes*MaxGrayLevel {
+		logrus.Error("error in computing Tij table, some values are not satisfied")
+	}
+}
+
+func (hist *ProbeHistogram) save(matrixFile string) error {
+	log.Printf("total pixels: %d", hist.M)
+	f, err := os.OpenFile(matrixFile, os.O_TRUNC|os.O_WRONLY|os.O_CREATE, 0644)
+	if err != nil {
+		return err
+	}
+	defer f.Close()
+
+	_, err = hist.Tmat.MarshalBinaryTo(f)
+	if err != nil {
+		return err
+	}
+
+	return nil
+}
+
+func LoadGrayTableMatrix(matrixFile string) (*mat.Dense, error) {
+	f, err := os.OpenFile(matrixFile, os.O_RDONLY, 0644)
+	if err != nil {
+		return nil, err
+	}
+	defer f.Close()
+
+	matrix := mat.Dense{}
+	if _, err := matrix.UnmarshalBinaryFrom(f); err != nil {
+		return nil, err
+	}
+
+	return &matrix, nil
+}
+
+func (hist *ProbeHistogram) sum(hists []*ProbeHistogram) {
+	for _, h := range hists {
+		hist.M += h.M
+		for i := 0; i < hist.probes; i++ {
+			hist.N_i[i] += h.N_i[i]
+			for k := 0; k < MaxGrayLevel; k++ {
+				hist.n_ik[i][k] += h.n_ik[i][k]
+			}
+		}
+		for gray := 0; gray < MaxGrayLevel; gray++ {
+			hist.m_l[gray] += h.m_l[gray]
+		}
+	}
+}

pkg/rrc/rrc.go

@@ -1,6 +1,8 @@
 package rrc
 
 import (
+	"bufio"
+	"fmt"
 	"os"
 
 	log "github.com/sirupsen/logrus"
@@ -13,55 +15,51 @@ import (
 const (
 	PANCameraProbeNum = 9344  // 全色探元数
 	MSSCameraProbeNum = 2336  // 多光谱探元数
-	MaxGrayLevel      = 65536
+	MaxGrayLevel      = 65536 // 16bit像素值
 )
 
 type RRC struct {
-	PANDataSet []string
+	Histograms [5]ProbeHistogram
-	MSSDataSet []string
-
-	Histograms [5]BandHistogram
 }
 
-type BandHistogram struct {
+func NewRRC() *RRC {
-	width int         // 探元数
+	r := RRC{}
-	N_i   []int       // N_i 探元像素总数 PAN 0-9343 MSS 0-2335
+	r.Histograms[0].init(PANCameraProbeNum)
-	n_ik  [][]int     // n_ik 第i探元灰度等级为k的像素数统计 PAN 9343 x 65536  MSS 2335 x 65536
+	for i := 1; i < 5; i++ {
-	p_ik  [][]float64 // p_ik = n_ik / N_i 探元灰度概率密度 PAN 9343 x 65536  MSS 2335 x 65536
+		r.Histograms[i].init(MSSCameraProbeNum)
-	m_l   []int64     // 具有灰度等级l的像素总数 l = 0-65535
-	M     int64       // 参与直方图统计的总像素数
-	P_l   []float64   // P_l = m_l/M 所有探元的期望直方图灰度等级为l的概率密度
-	S_ik  [][]float64 // S_ik = sum(p_ij),j=0..k 第i个探元直方图灰度等级k的累积概率密度
-	V_l   []float64   // V_l = sum(P_j),j=0..l // 期望直方图灰度级l对应的累积概率密度
-	Tij   [][]float64 // 第i个像元的j灰度等级对应的新的灰度值，用于修正图像
-}
-
-func (hist *BandHistogram) init(width int) {
-	hist.width = width
-	hist.M = 0
-	hist.N_i = make([]int, width)
-	hist.n_ik = make([][]int, width)
-	hist.p_ik = make([][]float64, width)
-	hist.m_l = make([]int64, MaxGrayLevel)
-	hist.P_l = make([]float64, MaxGrayLevel)
-	hist.S_ik = make([][]float64, width)
-	hist.V_l = make([]float64, MaxGrayLevel)
-	hist.Tij = make([][]float64, width)
-
-	for i := 0; i < width; i++ {
-		hist.n_ik[i] = make([]int, MaxGrayLevel)
-		hist.p_ik[i] = make([]float64, MaxGrayLevel)
-		hist.S_ik[i] = make([]float64, MaxGrayLevel)
-		hist.Tij[i] = make([]float64, MaxGrayLevel)
 	}
+
+	os.MkdirAll("data/rrc", 0755)
+
+	return &r
 }
 
+func (rrc *RRC) Statistical(dsPAN, dsMSS string) error {
+	rrc.StatisticalPAN(dsPAN)
+	rrc.StatisticalMSS(dsMSS)
+
+	return nil
+}
+
+func (rrc *RRC) Close() {}
+
 // 统计探元灰度的累积概率密度
-func (rrc *RRC) StatisticalPAN(l0 string) error {
+func (rrc *RRC) StatisticalPAN(dsfile string) error {
-	data, err := os.ReadFile(l0)
+	f, err := os.Open(dsfile)
 	if err != nil {
 		return err
 	}
+	defer f.Close()
+
+	scanner := bufio.NewScanner(f)
+	for scanner.Scan() {
+		l0 := scanner.Text()
+		log.Println("statistical PAN RAW: ", l0)
+		data, err := os.ReadFile(l0)
+		if err != nil {
+			log.Error("Error reading file: ", err)
+			continue
+		}
 
 		height := len(data) / (PANCameraProbeNum * 2)
 		img, err := gocv.NewMatFromBytes(height, PANCameraProbeNum, gocv.MatTypeCV16U, data)
@@ -70,71 +68,64 @@ func (rrc *RRC) StatisticalPAN(l0 string) error {
 			return err
 		}
 
-	hist := BandHistogram{}
+		rrc.Histograms[0].statistical(img)
-	hist.init(PANCameraProbeNum)
+		img.Close()
+	}
 
-	rrc.statistical(img, &hist)
+	log.Println("compute PAN histogram...")
-	rrc.compute(&hist)
+	rrc.Histograms[0].compute()
-	rrc.output(&hist, "pan_gray_table.tif")
+	log.Println("save PAN gray table matrix.")
+	rrc.Histograms[0].save("data/rrc/pan_gray_table.dat")
 
 	return nil
 }
 
-func (rrc *RRC) statistical(img gocv.Mat, hist *BandHistogram) error {
+func (rrc *RRC) StatisticalMSS(dsfile string) error {
-	hist.M += int64(img.Rows() * img.Cols())
+	f, err := os.Open(dsfile)
+	if err != nil {
+		return err
+	}
+	defer f.Close()
 
-	for i := 0; i < hist.width; i++ {
+	scanner := bufio.NewScanner(f)
-		hist.N_i[i] += img.Rows()
+	for scanner.Scan() {
+		l0 := scanner.Text()
+		log.Println("statistical MSS RAW: ", l0)
+		data, err := os.ReadFile(l0)
+		if err != nil {
+			log.Error("Error reading file: ", err)
+			continue
 		}
 
-	for y := 0; y < img.Rows(); y++ {
+		width := MSSCameraProbeNum
-		for x := 0; x < img.Cols(); x++ {
+		height := len(data) / (width * 2)
-			gray := uint16(img.GetShortAt(x, y))
+		mssData := make([][]byte, 4)
-			hist.n_ik[x][gray]++
+		for h := 0; h < height; h++ {
+			row := data[h*width*4*2 : (h+1)*width*4*2]
+			for i := 0; i < 4; i++ {
+				mssData[i] = append(mssData[i], row[i*width*2:(i+1)*width*2]...)
 			}
 		}
 
-	for gray := 0; gray < MaxGrayLevel; gray++ {
+		var mssImages [4]gocv.Mat
-		for i := 0; i < hist.width; i++ {
+		for i := 0; i < 4; i++ {
-			hist.m_l[gray] += int64(hist.n_ik[i][gray])
+			mssImages[i], err = gocv.NewMatFromBytes(height, width, gocv.MatTypeCV16U, mssData[i])
+			if err != nil {
+				log.Error("Error creating Mat from bytes: ", err)
+				return err
 			}
+
+			rrc.Histograms[i+1].statistical(mssImages[i])
+			mssImages[i].Close()
+		}
+	}
+
+	for i := 1; i < 5; i++ {
+		log.Println("compute MSS histogram...")
+		rrc.Histograms[i].compute()
+		log.Println("save MSS gray table matrix.")
+		rrc.Histograms[i].save(fmt.Sprintf("data/rrc/mss%d_gray_table.dat", i))
 	}
 
 	return nil
 }
-
-func (rrc *RRC) compute(hist *BandHistogram) error {
-	width := len(hist.N_i)
-
-	for i := 0; i < width; i++ {
-		for k := 0; k < MaxGrayLevel; k++ {
-			hist.p_ik[i][k] = float64(hist.n_ik[i][k]) / float64(hist.N_i[i])
-		}
-	}
-
-	for gray := 0; gray < MaxGrayLevel; gray++ {
-		hist.P_l[gray] = float64(hist.m_l[gray]) / float64(hist.M)
-	}
-
-	for i := 0; i < width; i++ {
-		for k := 0; k < MaxGrayLevel; k++ {
-			hist.S_ik[i][k] = 0
-			for j := 0; j <= k; j++ {
-				hist.S_ik[i][k] += hist.p_ik[i][j]
-			}
-		}
-	}
-
-	for gray := 0; gray < MaxGrayLevel; gray++ {
-		hist.V_l[gray] = 0
-		for j := 0; j <= gray; j++ {
-			hist.V_l[gray] += hist.P_l[j]
-		}
-	}
-
-	return nil
-}
-
-func (rrc *RRC) output(hist *BandHistogram, referenceTIF string) error {
-	return nil
-}