232 lines
6.0 KiB
Go
232 lines
6.0 KiB
Go
package rrc
|
||
|
||
import (
|
||
"bytes"
|
||
"encoding/binary"
|
||
"math"
|
||
"os"
|
||
|
||
log "github.com/sirupsen/logrus"
|
||
"gocv.io/x/gocv"
|
||
"gonum.org/v1/plot"
|
||
"gonum.org/v1/plot/plotter"
|
||
"gonum.org/v1/plot/plotutil"
|
||
"gonum.org/v1/plot/vg"
|
||
)
|
||
|
||
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 [][]uint16 // 第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)
|
||
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)
|
||
}
|
||
}
|
||
|
||
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())
|
||
|
||
// 探元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灰度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++ {
|
||
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
|
||
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]
|
||
}
|
||
}
|
||
|
||
// 生成查找表
|
||
log.Info("computing Tij table")
|
||
var nT int64
|
||
for i := 0; i < hist.probes; i++ {
|
||
for k := 0; k < MaxGrayLevel; k++ {
|
||
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)
|
||
}
|
||
}
|
||
}
|
||
|
||
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, nT:", nT,
|
||
"probes*MaxGrayLevel:", hist.probes*MaxGrayLevel)
|
||
}
|
||
}
|
||
|
||
const (
|
||
CheckPointProbe = 1000
|
||
CheckPointGray = 15000
|
||
)
|
||
|
||
func (hist *ProbeHistogram) saveLUT(fLUT string) error {
|
||
file, err := os.Create(fLUT)
|
||
if err != nil {
|
||
return err
|
||
}
|
||
defer file.Close()
|
||
for i := 0; i < hist.probes; i++ {
|
||
binary.Write(file, binary.LittleEndian, hist.Tmat[i])
|
||
if i == CheckPointProbe {
|
||
log.Infof("Probes *, LUT check point [%d][%d]: %d", i, CheckPointGray, hist.Tmat[i][CheckPointGray])
|
||
}
|
||
}
|
||
|
||
return nil
|
||
}
|
||
|
||
func LoadLUT(fLUT string, probes int) ([][]uint16, error) {
|
||
data, err := os.ReadFile(fLUT)
|
||
if 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 Tmat, 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]
|
||
}
|
||
}
|
||
}
|
||
|
||
func (hist *ProbeHistogram) free() {
|
||
hist.N_i = nil
|
||
hist.n_ik = nil
|
||
hist.p_ik = nil
|
||
hist.m_l = nil
|
||
hist.P_l = nil
|
||
hist.S_ik = nil
|
||
hist.V_l = nil
|
||
hist.Tmat = nil
|
||
}
|