diff --git a/cmd/proc.go b/cmd/proc.go index 3136cec..d5a7321 100644 --- a/cmd/proc.go +++ b/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) } diff --git a/cmd/rrc.go b/cmd/rrc.go new file mode 100644 index 0000000..d68575f --- /dev/null +++ b/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) +} diff --git a/pkg/producer/rrc.go b/pkg/producer/rrc.go new file mode 100644 index 0000000..60eb10f --- /dev/null +++ b/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 +} diff --git a/pkg/rrc/histogram.go b/pkg/rrc/histogram.go new file mode 100644 index 0000000..c0d1651 --- /dev/null +++ b/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] + } + } +} diff --git a/pkg/rrc/rrc.go b/pkg/rrc/rrc.go index d54b42d..5ae4d47 100644 --- a/pkg/rrc/rrc.go +++ b/pkg/rrc/rrc.go @@ -1,6 +1,8 @@ package rrc import ( + "bufio" + "fmt" "os" log "github.com/sirupsen/logrus" @@ -11,130 +13,119 @@ import ( // 采用在轨统计定标方法。利用卫星在轨后获取的常规影像数据，统计每个探元出现的灰度频次。 const ( - PANCameraProbeNum = 9344 // 全色探元数 - MSSCameraProbeNum = 2336 // 多光谱探元数 - MaxGrayLevel = 65536 + PANCameraProbeNum = 9344 // 全色探元数 + MSSCameraProbeNum = 2336 // 多光谱探元数 + MaxGrayLevel = 65536 // 16bit像素值 ) type RRC struct { - PANDataSet []string - MSSDataSet []string - - Histograms [5]BandHistogram + Histograms [5]ProbeHistogram } -type BandHistogram struct { - width int // 探元数 - N_i []int // N_i 探元像素总数 PAN 0-9343 MSS 0-2335 - n_ik [][]int // 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对应的累积概率密度 - 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) +func NewRRC() *RRC { + r := RRC{} + r.Histograms[0].init(PANCameraProbeNum) + for i := 1; i < 5; i++ { + r.Histograms[i].init(MSSCameraProbeNum) } + + 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 { - data, err := os.ReadFile(l0) +func (rrc *RRC) StatisticalPAN(dsfile string) error { + f, err := os.Open(dsfile) if err != nil { return err } + defer f.Close() - 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 + 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) + if err != nil { + log.Error("Error creating Mat from bytes: ", err) + return err + } + + rrc.Histograms[0].statistical(img) + img.Close() } - hist := BandHistogram{} - hist.init(PANCameraProbeNum) - - rrc.statistical(img, &hist) - rrc.compute(&hist) - rrc.output(&hist, "pan_gray_table.tif") + 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") return nil } -func (rrc *RRC) statistical(img gocv.Mat, hist *BandHistogram) error { - hist.M += int64(img.Rows() * img.Cols()) - - for i := 0; i < hist.width; i++ { - hist.N_i[i] += img.Rows() +func (rrc *RRC) StatisticalMSS(dsfile string) error { + f, err := os.Open(dsfile) + if err != nil { + return err } + defer f.Close() - 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]++ + scanner := bufio.NewScanner(f) + 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 gray := 0; gray < MaxGrayLevel; gray++ { - for i := 0; i < hist.width; i++ { - hist.m_l[gray] += int64(hist.n_ik[i][gray]) - } - } - - 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] + 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]...) } } - } - 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] + 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() } } - return nil -} + 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)) + } -func (rrc *RRC) output(hist *BandHistogram, referenceTIF string) error { return nil }