save Tmat in binary format

pkg/rrc/helper.go

@@ -0,0 +1,26 @@
+package rrc
+
+func searchVL(V []float64, Sik float64) int {
+	left, right := 0, len(V)-2
+
+	// if Sik < V[0] || Sik > V[len(V)-1] {
+	// 	return -1
+	// }
+
+	for left <= right {
+		mid := left + (right-left)/2
+		// Check if Sik is between V[mid] and V[mid+1]
+		if V[mid]-1e-6 <= Sik && Sik <= V[mid+1]+1e-6 {
+			return mid
+		} else if Sik < V[mid] {
+			right = mid - 1
+		} else if Sik == V[mid] {
+			return mid
+		} else {
+			left = mid + 1
+		}
+	}
+
+	// Return -1 if no such position is found
+	return -1
+}

pkg/rrc/histogram.go

@@ -1,14 +1,18 @@
 package rrc
 
 import (
+	"bytes"
+	"encoding/binary"
 	"fmt"
 	"math"
 	"os"
 
-	"github.com/sirupsen/logrus"
 	log "github.com/sirupsen/logrus"
 	"gocv.io/x/gocv"
-	"gonum.org/v1/gonum/mat"
+	"gonum.org/v1/plot"
+	"gonum.org/v1/plot/plotter"
+	"gonum.org/v1/plot/plotutil"
+	"gonum.org/v1/plot/vg"
 )
 
 type ProbeHistogram struct {
@@ -21,7 +25,7 @@ type ProbeHistogram struct {
 	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灰度等级对应的新的灰度值，用于修正图像
+	Tmat   [][]uint16  // 第i个像元的j灰度等级对应的新的灰度值，用于修正图像
 }
 
 func (hist *ProbeHistogram) init(width int) {
@@ -34,19 +38,30 @@ func (hist *ProbeHistogram) init(width int) {
 	hist.P_l = make([]float64, MaxGrayLevel)
 	hist.S_ik = make([][]float64, width)
 	hist.V_l = make([]float64, MaxGrayLevel)
+	hist.Tmat = make([][]uint16, width)
 
 	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[i] = make([]uint16, MaxGrayLevel)
 	}
+}
 
-	hist.Tmat = mat.NewDense(width, MaxGrayLevel, nil)
+func (hist *ProbeHistogram) init0(width int) {
+	hist.probes = width
+	hist.M = 0
+	hist.N_i = make([]int64, width)
+	hist.n_ik = make([][]int64, width)
+	hist.m_l = make([]int64, MaxGrayLevel)
+
+	for i := 0; i < width; i++ {
+		hist.n_ik[i] = make([]int64, MaxGrayLevel)
+	}
 }
 
 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++ {
@@ -72,87 +87,118 @@ func (hist *ProbeHistogram) statistical(img gocv.Mat) error {
 }
 
 func (hist *ProbeHistogram) compute() {
-	// 探元i灰度概率密度
+	// 探元i灰度k概率密度
+	log.Info("computing p_ik")
 	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])
 		}
 	}
 
+	points := plotter.XYs{}
+	var maxY int64
 	// 所有探元的期望直方图灰度等级为l的概率密度
+	log.Info("computing p_l")
 	for gray := 0; gray < MaxGrayLevel; gray++ {
 		hist.P_l[gray] = float64(hist.m_l[gray]) / float64(hist.M)
+
+		points = append(points, plotter.XY{X: float64(gray), Y: float64(hist.m_l[gray])})
+		if hist.m_l[gray] > maxY {
+			maxY = hist.m_l[gray]
+		}
 	}
 
+	plt := plot.New()
+	plt.Title.Text = "Gray Level Histogram"
+	plt.X.Label.Text = "Gray Level"
+	plt.Y.Label.Text = "Pixels"
+	plt.Y.Min, plt.X.Min, plt.Y.Max, plt.X.Max = 0, 0, float64(maxY), 65536
+	plotutil.AddLines(plt, points)
+	plt.Save(5*vg.Inch, 5*vg.Inch, "data/rrc/01-gray-level-histogram.png")
+
 	// 第i个探元直方图灰度等级k的累积概率密度
+	log.Info("computing s_ik")
 	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]
+			if k == 0 {
+				hist.S_ik[i][k] = hist.p_ik[i][0]
+			} else {
+				hist.S_ik[i][k] = hist.S_ik[i][k-1] + hist.p_ik[i][k]
 			}
 		}
 	}
 
 	// 期望直方图灰度级l对应的累积概率密度
+
+	log.Info("computing v_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]
+		if gray == 0 {
+			hist.V_l[gray] = hist.P_l[0]
+		} else {
+			hist.V_l[gray] = hist.V_l[gray-1] + hist.P_l[gray]
 		}
 	}
 
 	// 生成查找表
-	nT := 0
+	log.Info("computing Tij table")
+	var nT int64
 	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))
-					}
-				}
+			l := searchVL(hist.V_l, hist.S_ik[i][k])
+			if l == -1 {
+				continue
+			}
+
+			nT++
+			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[i][k] = uint16(l)
+			} else {
+				hist.Tmat[i][k] = uint16(l + 1)
 			}
 		}
 	}
 
-	if nT != hist.probes*MaxGrayLevel {
-		logrus.Error("error in computing Tij table, some values are not satisfied")
+	log.Info("total Tij table entries:", nT)
+	if nT != int64(hist.probes*MaxGrayLevel) {
+		log.Warn("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)
+func (hist *ProbeHistogram) saveLUT(fLUT string) error {
+	file, err := os.Create(fLUT)
 	if err != nil {
 		return err
 	}
-	defer f.Close()
-
-	_, err = hist.Tmat.MarshalBinaryTo(f)
-	if err != nil {
-		return err
+	defer file.Close()
+	for i := 0; i < hist.probes; i++ {
+		binary.Write(file, binary.LittleEndian, hist.Tmat[i])
 	}
 
 	return nil
 }
 
-func LoadGrayTableMatrix(matrixFile string) (*mat.Dense, error) {
-	f, err := os.OpenFile(matrixFile, os.O_RDONLY, 0644)
+const (
+	CheckPointProbe = 1000
+	CheckPointGray  = 15000
+)
+
+func LoadLUT(fLUT string, probes int) ([][]uint16, error) {
+	data, err := os.ReadFile(fLUT)
 	if err != nil {
 		return nil, err
 	}
-	defer f.Close()
 
-	matrix := mat.Dense{}
-	if _, err := matrix.UnmarshalBinaryFrom(f); err != nil {
-		return nil, err
+	Tmat := make([][]uint16, probes)
+	for i := 0; i < probes; i++ {
+		Tmat[i] = make([]uint16, MaxGrayLevel)
+		binary.Read(bytes.NewReader(data[i*MaxGrayLevel*2:i*MaxGrayLevel*2+MaxGrayLevel*2]), binary.LittleEndian, &Tmat[i])
+		if i == CheckPointProbe {
+			log.Infof("Probes[%d], LUT check point [%d][%d]: %d", probes, i, CheckPointGray, Tmat[i][CheckPointGray])
+		}
 	}
 
-	return &matrix, nil
+	return Tmat, nil
 }
 
 func (hist *ProbeHistogram) sum(hists []*ProbeHistogram) {
@@ -168,4 +214,7 @@ func (hist *ProbeHistogram) sum(hists []*ProbeHistogram) {
 			hist.m_l[gray] += h.m_l[gray]
 		}
 	}
+
+	fmt.Println("Hist.M:", hist.M)
+	fmt.Println("Hist.probes:", hist.probes)
 }

pkg/rrc/rrc.go

@@ -4,6 +4,7 @@ import (
 	"bufio"
 	"fmt"
 	"os"
+	"sync"
 
 	log "github.com/sirupsen/logrus"
 	"gocv.io/x/gocv"
@@ -52,30 +53,57 @@ func (rrc *RRC) StatisticalPAN(dsfile string) error {
 	defer f.Close()
 
 	scanner := bufio.NewScanner(f)
+	var files []string
 	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)
-		if err != nil {
-			log.Error("Error creating Mat from bytes: ", err)
-			return err
-		}
-
-		rrc.Histograms[0].statistical(img)
-		img.Close()
+		files = append(files, scanner.Text())
 	}
 
+	// 并发处理
+	var wg sync.WaitGroup
+	jobs := make(chan struct{}, 5)
+	var hists []*ProbeHistogram
+	var mutex sync.Mutex
+	for _, file := range files {
+		wg.Add(1)
+		go func(l0 string) {
+			defer wg.Done()
+
+			jobs <- struct{}{}
+			defer func() { <-jobs }()
+
+			log.Println("statistical PAN RAW: ", l0)
+			data, err := os.ReadFile(l0)
+			if err != nil {
+				log.Error("Error reading file: ", err)
+				return
+			}
+
+			height := len(data) / (PANCameraProbeNum * 2)
+			img, err := gocv.NewMatFromBytes(height, PANCameraProbeNum, gocv.MatTypeCV16U, data)
+			if err != nil {
+				log.Error("Error creating Mat from bytes: ", err)
+				return
+			}
+
+			var hist ProbeHistogram
+			hist.init(PANCameraProbeNum)
+			hist.statistical(img)
+			img.Close()
+
+			mutex.Lock()
+			hists = append(hists, &hist)
+			mutex.Unlock()
+
+		}(file)
+	}
+
+	wg.Wait()
+	log.Println("sum PAN histogram...")
+	rrc.Histograms[0].sum(hists)
 	log.Println("compute PAN histogram...")
 	rrc.Histograms[0].compute()
 	log.Println("save PAN gray table matrix.")
-	rrc.Histograms[0].save("data/rrc/pan_gray_table.dat")
+	rrc.Histograms[0].saveLUT("data/rrc/B0.LUT")
 
 	return nil
 }
@@ -88,43 +116,68 @@ func (rrc *RRC) StatisticalMSS(dsfile string) error {
 	defer f.Close()
 
 	scanner := bufio.NewScanner(f)
+
+	var files []string
 	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
-		}
-
-		width := MSSCameraProbeNum
-		height := len(data) / (width * 2)
-		mssData := make([][]byte, 4)
-		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]...)
-			}
-		}
-
-		var mssImages [4]gocv.Mat
-		for i := 0; i < 4; i++ {
-			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()
-		}
+		files = append(files, scanner.Text())
 	}
 
+	var wg sync.WaitGroup
+	jobs := make(chan struct{}, 5)
+	var hists [4][]*ProbeHistogram
+	var mutex sync.Mutex
+	for _, file := range files {
+		wg.Add(1)
+		go func(l0 string) {
+			defer wg.Done()
+			jobs <- struct{}{}
+			defer func() { <-jobs }()
+
+			log.Println("statistical MSS RAW: ", l0)
+			data, err := os.ReadFile(l0)
+			if err != nil {
+				log.Error("Error reading file: ", err)
+				return
+			}
+
+			width := MSSCameraProbeNum
+			height := len(data) / (width * 2)
+			mssData := make([][]byte, 4)
+			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]...)
+				}
+			}
+
+			var mssImages [4]gocv.Mat
+			for i := 0; i < 4; i++ {
+				mssImages[i], err = gocv.NewMatFromBytes(height, width, gocv.MatTypeCV16U, mssData[i])
+				if err != nil {
+					log.Error("Error creating Mat from bytes: ", err)
+					return
+				}
+
+				var hist ProbeHistogram
+				hist.init(PANCameraProbeNum)
+				hist.statistical(mssImages[i])
+				mssImages[i].Close()
+
+				mutex.Lock()
+				hists[i] = append(hists[i], &hist)
+				mutex.Unlock()
+			}
+		}(file)
+	}
+
+	wg.Wait()
 	for i := 1; i < 5; i++ {
+		log.Println("sum MSS histogram...")
+		rrc.Histograms[i].sum(hists[i-1])
 		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))
+		rrc.Histograms[i].saveLUT(fmt.Sprintf("data/rrc/B%d.LUT", i))
 	}
 
 	return nil