gocv v0.36.1 opencv 4.9.0

2024-06-09 20:08:32 +08:00
parent c046da2321
commit 031b76be39
8 changed files with 377 additions and 77 deletions
--- a/pkg/calculator/camera.go
+++ b/pkg/calculator/camera.go
@@ -0,0 +1,59 @@
+package calculator
+
+import (
+	"fmt"
+	"math"
+)
+
+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
+)
+
+// 计算过程使用PAN分辨率
+func CameraDirectionVec(u, v float64) []float64 {
+	w := PANPixels * PANCellSize / 1000.0 // 像素宽度, mm
+	h := w
+	d := math.Sqrt(math.Pow(w, 2) + math.Pow(h, 2)) // 对角线长度, mm
+	fov := 2 * math.Atan2(d/2, FocalLength)         // 对角线视场角
+
+	fmt.Println("Camera Parameters:")
+	fmt.Printf("Focal Length: %.6f mm\n", FocalLength)
+	fmt.Printf("FOV (calculated): %.6f degree\n", fov*180/math.Pi)
+	fmt.Printf("Width: %.6f mm\n", w)
+	fmt.Printf("Height: %.6f mm\n", h)
+	fmt.Printf("Diagonal: %.6f mm\n", d)
+
+	// 从给定FOV计算d
+	dCalcOfFOV := 2 * FocalLength * math.Tan(FOV/2*math.Pi/180)
+	fmt.Printf("Diagonal (calculated from FOV): %.6f mm\n", dCalcOfFOV)
+
+	directionVec := []float64{0, 0, 0}
+	directionVec[0] = 0                                        // x方向, mm，线性CCD每次单行成像
+	directionVec[1] = (v - PANPixels/2) * PANCellSize / 1000.0 // y方向, mm
+	directionVec[2] = -FocalLength                             // z方向, mm
+
+	// 归一化
+	fmt.Printf("Direction Vector: (%.6f, %.6f, %.6f) \n", directionVec[0], directionVec[1], directionVec[2])
+	directionVec = normalizeVec(directionVec)
+	fmt.Printf("Direction Vector (normalized): (%.6f, %.6f, %.6f) \n", directionVec[0], directionVec[1], directionVec[2])
+
+	return directionVec
+}
+
+func normalizeVec(vec []float64) []float64 {
+	var norm float64
+	for i := 0; i < len(vec); i++ {
+		norm += vec[i] * vec[i]
+	}
+
+	for i := 0; i < len(vec); i++ {
+		vec[i] /= math.Sqrt(norm)
+	}
+	return vec
+}
--- a/pkg/calculator/intersection.go
+++ b/pkg/calculator/intersection.go
@@ -7,13 +7,6 @@ import (
 	"gonum.org/v1/gonum/mat"
 )

-// 常量
-const (
-	focal   = 1.3  // 焦距, m
-	FOV     = 1.7  // 视场角,degree
-	nPixels = 9344 // 像素数
-)
-
 type IntersectionPoint struct {
 	Lat float64
 	Lon float64
@@ -22,46 +15,56 @@ type IntersectionPoint struct {

 func Intersection(q Quaternion, satPos84 []float64, satTime time.Time, ucam int) IntersectionPoint {
 	alpha := FOV * math.Pi / 180.0
-	alpha = -alpha/2.0 + float64(ucam)*(alpha/float64(nPixels))
+	alpha = -alpha/2.0 + float64(ucam)*(alpha/float64(PANPixels))
 	direction := []float64{0, math.Tan(alpha), -1.3}

+	// -------- 相机坐标系下CCD成像方向向量转到卫星坐标系 --------
+	Rcam := mat.NewDense(3, 3, []float64{
+		1, 0, 0,
+		0, 1, 0,
+		0, 0, 1,
+	})
+	var dCam mat.VecDense
+	dCam.MulVec(Rcam, mat.NewVecDense(3, direction))
+
+	// -------- 转到ECI坐标系 --------
 	Ratt := q.ToRotationMatrix()
-	RattT := &mat.Dense{}
-	RattT.Inverse(Ratt)
-
-	v := mat.NewVecDense(3, direction)
 	var result mat.VecDense
-	result.MulVec(Ratt, v)
-	eciDirection := result.RawVector().Data
+	result.MulVec(Ratt, &dCam)
+	dECI := result.RawVector().Data

-	// intersection := intersectWithEllipsoid(satPos, eciDirection)
-	// lat, lon, _ := J2000ToWGS84(intersection[0], intersection[1], intersection[2], satTime)
+	// -------- 转到ECEF坐标系 --------
+	x, y, z := ECItoECEF(dECI[0], dECI[1], dECI[2], satTime)
+	dECEF := []float64{x, y, z}

-	x, y, z := ECItoECEF(eciDirection[0], eciDirection[1], eciDirection[2], satTime)
-	ecefDirection := []float64{x, y, z}
-	intersection := intersectWithEllipsoid(satPos84, ecefDirection)
+	// -------- 计算与地球表面的交点 --------
+	intersection := intersectWithEllipsoid(satPos84, dECEF)
 	lat, lon, h := ECEFToGeodetic(intersection[0], intersection[1], intersection[2])
 	return IntersectionPoint{Lat: lat, Lon: lon, H: h}
-
 }

 func Intersection2(Qsat2orbit, Qorbit2eci Quaternion, satPos84 []float64, satTime time.Time, ucam int) IntersectionPoint {
 	alpha := FOV * math.Pi / 180.0
-	alpha = -alpha/2.0 + float64(ucam)*(alpha/float64(nPixels))
+	alpha = -alpha/2.0 + float64(ucam)*(alpha/float64(PANPixels))
 	direction := []float64{0, math.Tan(alpha), -1.3} // 卫星（相机）坐标系下CCD成像方向向量

+	// -------- 相机坐标系下CCD成像方向向量转到卫星坐标系 --------
+	Rcam := mat.NewDense(3, 3, []float64{
+		1, 0, 0,
+		0, 1, 0,
+		0, 0, 1,
+	})
+	var dCam mat.VecDense
+	dCam.MulVec(Rcam, mat.NewVecDense(3, direction))
+
 	// -------- 转到轨道坐标系 --------
 	Rsat2orbit := Qsat2orbit.ToRotationMatrix()
-	Rsat2orbitT := &mat.Dense{}
-	Rsat2orbitT.Inverse(Rsat2orbit)
 	var r0 mat.VecDense
-	r0.MulVec(Rsat2orbit, mat.NewVecDense(3, direction))
+	r0.MulVec(Rsat2orbit, &dCam)
 	dOrbit := r0.RawVector().Data

 	// -------- 转到ECI坐标系 --------
 	Rorbit2eci := Qorbit2eci.ToRotationMatrix()
-	Rorbit2eciT := &mat.Dense{}
-	Rorbit2eciT.Inverse(Rorbit2eci)
 	var r1 mat.VecDense
 	r1.MulVec(Rorbit2eci, mat.NewVecDense(3, dOrbit))
 	dECI := r1.RawVector().Data
--- a/pkg/producer/aux.go
+++ b/pkg/producer/aux.go
@@ -39,12 +39,12 @@ func (r *Registrator) SetSceneBoundary(scene *Scene) (topLeft, bottomRight orb.P
 	startPosInAux, endPosInAux := r.SceneInAuxIndex(scene)

 	as := r.AuxPlatforms[startPosInAux]
-	startPos84 := []float64{as.W84PosX, as.W84PosY, as.W84PosZ}
 	startTime := time.Unix(int64(auxilary.ReferenceTime2000)+int64(as.UTCTimeSec), int64(as.Microsecond)*1000).UTC()
+	startPos84 := []float64{as.W84PosX, as.W84PosY, as.W84PosZ}

 	ae := r.AuxPlatforms[endPosInAux]
-	endPos84 := []float64{ae.W84PosX, ae.W84PosY, ae.W84PosZ}
 	endTime := time.Unix(int64(auxilary.ReferenceTime2000)+int64(ae.UTCTimeSec), int64(ae.Microsecond)*1000).UTC()
+	endPos84 := []float64{ae.W84PosX, ae.W84PosY, ae.W84PosZ}

 	// ------------------ 使用定姿态四元数计算图像边界 ------------------
 	log.Info("using attitude quaternion to calculate image boundary...")
@@ -62,15 +62,15 @@ func (r *Registrator) SetSceneBoundary(scene *Scene) (topLeft, bottomRight orb.P
 	// Qsat2orbit.W = math.Sqrt(1 - Qsat2orbit.X*Qsat2orbit.X - Qsat2orbit.Y*Qsat2orbit.Y - Qsat2orbit.Z*Qsat2orbit.Z)
 	// Qorbit2eci := calculator.Quaternion{X: as.QuatOrbJQ1, Y: as.QuatOrbJQ2, Z: as.QuatOrbJQ3}
 	// Qorbit2eci.W = math.Sqrt(1 - Qorbit2eci.X*Qorbit2eci.X - Qorbit2eci.Y*Qorbit2eci.Y - Qorbit2eci.Z*Qorbit2eci.Z)
-	// line0Start = calculator.Intersection2(Qsat2orbit, Qorbit2eci, startPos84, startTime, 0)
-	// line0End = calculator.Intersection2(Qsat2orbit, Qorbit2eci, startPos84, startTime, 9344)
+	// line0Start := calculator.Intersection2(Qsat2orbit, Qorbit2eci, startPos84, startTime, 0)
+	// line0End := calculator.Intersection2(Qsat2orbit, Qorbit2eci, startPos84, startTime, 9344)

 	// Qsat2orbit = calculator.Quaternion{X: ae.QuatOrbitQ1, Y: ae.QuatOrbitQ2, Z: ae.QuatOrbitQ3}
 	// Qsat2orbit.W = math.Sqrt(1 - Qsat2orbit.X*Qsat2orbit.X - Qsat2orbit.Y*Qsat2orbit.Y - Qsat2orbit.Z*Qsat2orbit.Z)
 	// Qorbit2eci = calculator.Quaternion{X: ae.QuatOrbJQ1, Y: ae.QuatOrbJQ2, Z: ae.QuatOrbJQ3}
 	// Qorbit2eci.W = math.Sqrt(1 - Qorbit2eci.X*Qorbit2eci.X - Qorbit2eci.Y*Qorbit2eci.Y - Qorbit2eci.Z*Qorbit2eci.Z)
-	// lineNStart = calculator.Intersection2(Qsat2orbit, Qorbit2eci, endPos84, endTime, 0)
-	// lineNEnd = calculator.Intersection2(Qsat2orbit, Qorbit2eci, endPos84, endTime, 9344)
+	// lineNStart := calculator.Intersection2(Qsat2orbit, Qorbit2eci, endPos84, endTime, 0)
+	// lineNEnd := calculator.Intersection2(Qsat2orbit, Qorbit2eci, endPos84, endTime, 9344)

 	// ------------------ 计算图像边界距离和分辨率 ------------------
 	W0 := geo.Distance(orb.Point{line0Start.Lon, line0Start.Lat}, orb.Point{line0End.Lon, line0End.Lat})
@@ -79,7 +79,7 @@ func (r *Registrator) SetSceneBoundary(scene *Scene) (topLeft, bottomRight orb.P
 	// HN := geo.Distance(orb.Point{line0End.Lon, line0End.Lat}, orb.Point{lineNEnd.Lon, lineNEnd.Lat})
 	xResolution := W0 / float64(scene.Width)
 	yResolution := H0 / float64(scene.Height)
-	scene.Meta.Gsd = math.Max(xResolution, yResolution)
+	scene.Meta.Gsd = math.Min(xResolution, yResolution)

 	// log.Debug("distance 0: ", W0)
 	// log.Debug("distance N: ", WN)
--- a/pkg/producer/scenes.go
+++ b/pkg/producer/scenes.go
@@ -196,9 +196,12 @@ func (r *Registrator) DoScenePansharpen(panScenes []*Scene, mssScenes []*Scene)
 			return err
 		}

-		GTiffToJPG(fusedTiff, strings.Replace(fusedTiff, ".tiff", ".jpg", 1), true)
+		jpg := strings.Replace(fusedTiff, ".tiff", ".jpg", 1)
+		GTiffToJPG(fusedTiff, jpg, true)
 		r.report.Scenes = append(r.report.Scenes, ReportScene{
-			TiffData: fusedTiff,
+			TiffData:    fusedTiff,
+			Name:        strings.TrimSuffix(filepath.Base(fusedTiff), ".tiff"),
+			BrowserData: jpg,
 		})
 	}
 	return nil