diff --git a/cmd/proc.go b/cmd/proc.go index 43a2550..3136cec 100644 --- a/cmd/proc.go +++ b/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), diff --git a/pkg/calculator/camera.go b/pkg/calculator/camera.go index 853bbad..ee87d1f 100644 --- a/pkg/calculator/camera.go +++ b/pkg/calculator/camera.go @@ -7,15 +7,15 @@ import ( ) const ( - FocalLength = 1300.0 // 焦距, mm - FOV = 1.7 // 对角线视场角,degree - FOVCALC = 1.86 // 对角线视场角,degree - PANPixels = 9344.0 - PANCellSize = 3.2 // µm - MSSPixels = 2336.0 - MSSCellSize = 12.8 // µm - AngleCamSatX = 0.0 // 相机与卫星本体X轴的安装角度, degree FIXME: 安装矩阵应该由卫星方提供 - AngleCamSatY = 0.5 // 相机与卫星本体Y轴的安装角度, degree + FocalLength = 1300.0 // 焦距, mm + FOV = 1.7 // 对角线视场角,degree + FOVCALC = 1.86 // 对角线视场角,degree + PANPixels = 9344.0 + PANCellSize = 3.2 // µm + MSSPixels = 2336.0 + MSSCellSize = 12.8 // µm + CameraRoll = 0.0 // 相机与卫星本体X轴的安装角度, degree FIXME: 安装矩阵应该由卫星方提供 + CameraPitch = 0.5 // 相机与卫星本体Y轴的安装角度, degree ) // 计算过程使用PAN分辨率 diff --git a/pkg/calculator/intersection.go b/pkg/calculator/intersection.go index b93cd43..d769912 100644 --- a/pkg/calculator/intersection.go +++ b/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)) diff --git a/pkg/calculator/orbit.go b/pkg/calculator/orbit.go index 6b95fd4..06c98da 100644 --- a/pkg/calculator/orbit.go +++ b/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)) diff --git a/pkg/producer/aux.go b/pkg/producer/aux.go index ac40be7..a5a9125 100644 --- a/pkg/producer/aux.go +++ b/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 diff --git a/pkg/producer/params.go b/pkg/producer/params.go index e142cb6..98d2654 100644 --- a/pkg/producer/params.go +++ b/pkg/producer/params.go @@ -17,6 +17,7 @@ type Params struct { FusTIffFile string SubScenes bool ReportFile string + DataId string } type XMLImageTask struct { diff --git a/pkg/rrc/rrc.go b/pkg/rrc/rrc.go new file mode 100644 index 0000000..d54b42d --- /dev/null +++ b/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 +}