RRC

cmd/proc.go

@@ -31,6 +31,7 @@ var procCmd = &cobra.Command{
 		if err := reg.LoadAuxData(); err != nil {
 			logrus.Fatal(err)
 		}
+		reg.AuxPrint()
 
 		if err := reg.LoadMssRaw(); err != nil {
 			logrus.Fatal(err)
@@ -108,6 +109,7 @@ func initParams() producer.Params {
 		taskParams.DoPansharpen = task.Params.DoPansharpen
 		taskParams.OutputDir = task.Params.OutputPath
 		taskParams.ReportFile = task.Params.ReportFile
+		taskParams.DataId = task.Params.DataID
 	}
 
 	taskParams.MssTiffFile = filepath.Join(taskParams.OutputDir, strings.TrimSuffix(filepath.Base(taskParams.MssRawFile),

pkg/calculator/camera.go

@@ -14,8 +14,8 @@ const (
 	PANCellSize = 3.2 // µm
 	MSSPixels   = 2336.0
 	MSSCellSize = 12.8 // µm
-	AngleCamSatX = 0.0  // 相机与卫星本体X轴的安装角度, degree FIXME: 安装矩阵应该由卫星方提供
+	CameraRoll  = 0.0  // 相机与卫星本体X轴的安装角度, degree FIXME: 安装矩阵应该由卫星方提供
-	AngleCamSatY = 0.5  // 相机与卫星本体Y轴的安装角度, degree
+	CameraPitch = 0.5  // 相机与卫星本体Y轴的安装角度, degree
 )
 
 // 计算过程使用PAN分辨率

pkg/calculator/intersection.go

@@ -22,7 +22,7 @@ func IntersectionAttitude(q Quaternion, satPos84 []float64, satTime time.Time, u
 	direction := CameraDirectionVec(0, float64(ucam))
 
 	// -------- 相机坐标系下CCD成像方向向量转到卫星坐标系 --------
-	Rcam := CameraRotMatrix(AngleCamSatX*math.Pi/180.0, AngleCamSatY*math.Pi/180.0, 0)
+	Rcam := CameraRotMatrix(CameraRoll*math.Pi/180.0, CameraPitch*math.Pi/180.0, 0)
 	var dCam mat.VecDense
 	dCam.MulVec(Rcam, mat.NewVecDense(3, direction))
 
@@ -52,7 +52,7 @@ func IntersectionECI(Qsat2orbit, Qorbit2eci Quaternion, satPos84 []float64, satT
 	direction := []float64{0, math.Tan(alpha), -1.3} // 卫星（相机）坐标系下CCD成像方向向量
 
 	// -------- 相机坐标系下CCD成像方向向量转到卫星坐标系 --------
-	Rcam := CameraRotMatrix(AngleCamSatX*math.Pi/180.0, AngleCamSatY*math.Pi/180.0, 0)
+	Rcam := CameraRotMatrix(CameraRoll*math.Pi/180.0, CameraPitch*math.Pi/180.0, 0)
 	var dCam mat.VecDense
 	dCam.MulVec(Rcam, mat.NewVecDense(3, direction))
 

pkg/calculator/orbit.go

@@ -37,7 +37,7 @@ func IntersectionECEF(Qsat2orbit Quaternion, satPos84, vec84 []float64, ucam int
 	direction := []float64{0, math.Tan(alpha), -1.3} // 卫星（相机）坐标系下CCD成像方向向量
 
 	// -------- 相机坐标系下CCD成像方向向量转到卫星坐标系 --------
-	Rcam := CameraRotMatrix(AngleCamSatX*math.Pi/180.0, AngleCamSatY*math.Pi/180.0, 0)
+	Rcam := CameraRotMatrix(CameraRoll*math.Pi/180.0, CameraPitch*math.Pi/180.0, 0)
 	var dCam mat.VecDense
 	dCam.MulVec(Rcam, mat.NewVecDense(3, direction))
 

pkg/producer/aux.go

@@ -1,7 +1,10 @@
 package producer
 
 import (
+	"encoding/json"
+	"fmt"
 	"math"
+	"os"
 	"time"
 
 	log "github.com/sirupsen/logrus"
@@ -9,6 +12,7 @@ import (
 	"github.com/duke-git/lancet/v2/mathutil"
 	"github.com/paulmach/orb"
 	"github.com/paulmach/orb/geo"
+	"github.com/paulmach/orb/geojson"
 	"github.com/paulmach/orb/planar"
 	"starwiz.cn/sjy01/image-proc/pkg/auxilary"
 	"starwiz.cn/sjy01/image-proc/pkg/calculator"
@@ -20,6 +24,20 @@ func (r *Registrator) LoadAuxData() error {
 	return err
 }
 
+// 数据校验和测试
+func (r *Registrator) AuxPrint() {
+	var fcPos84 geojson.FeatureCollection
+	for _, p := range r.AuxPlatforms {
+		lat, lon, _ := calculator.WGS84XYZtoLatLngH(p.W84PosX, p.W84PosY, p.W84PosZ)
+		point := orb.Point{lon, lat}
+		fcPos84.Features = append(fcPos84.Features, geojson.NewFeature(point))
+	}
+	data, _ := json.Marshal(fcPos84)
+	f, _ := os.Create(fmt.Sprintf("log/%s_aux_pos_84.geojson", r.Params.DataId))
+	defer f.Close()
+	f.Write(data)
+}
+
 func (r *Registrator) SceneImageTime(scene *Scene) (start, center, end time.Time) {
 	startPosInAux, endPosInAux := r.SceneInAuxIndex(scene)
 	centerPosInAux := (startPosInAux + endPosInAux) / 2

pkg/producer/params.go

@@ -17,6 +17,7 @@ type Params struct {
 	FusTIffFile          string
 	SubScenes            bool
 	ReportFile           string
+	DataId               string
 }
 
 type XMLImageTask struct {

pkg/rrc/rrc.go

@@ -0,0 +1,140 @@
+package rrc
+
+import (
+	"os"
+
+	log "github.com/sirupsen/logrus"
+	"gocv.io/x/gocv"
+)
+
+// Relative Radiation Correction
+// 采用在轨统计定标方法。利用卫星在轨后获取的常规影像数据，统计每个探元出现的灰度频次。
+
+const (
+	PANCameraProbeNum = 9344 // 全色探元数
+	MSSCameraProbeNum = 2336 // 多光谱探元数
+	MaxGrayLevel      = 65536
+)
+
+type RRC struct {
+	PANDataSet []string
+	MSSDataSet []string
+
+	Histograms [5]BandHistogram
+}
+
+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 (rrc *RRC) StatisticalPAN(l0 string) error {
+	data, err := os.ReadFile(l0)
+	if err != nil {
+		return err
+	}
+
+	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
+	}
+
+	hist := BandHistogram{}
+	hist.init(PANCameraProbeNum)
+
+	rrc.statistical(img, &hist)
+	rrc.compute(&hist)
+	rrc.output(&hist, "pan_gray_table.tif")
+
+	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()
+	}
+
+	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]++
+		}
+	}
+
+	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]
+			}
+		}
+	}
+
+	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
+}