default channel order of OpenCV is BGR

cmd/log.go

@@ -1,4 +1,4 @@
-package imageproc
+package main
 
 import (
 	"os"

cmd/proc.go

@@ -56,8 +56,8 @@ var procCmd = &cobra.Command{
 		}
 
 		reg.SaveScenesToTiff(panScenes, mssScenes)
-		reg.SaveOriginalPanToGDALGTiff(reg.Params.PanTiffFile)
+		// reg.SaveOriginalPanToGDALGTiff(reg.Params.PanTiffFile)
-		reg.SaveRegisteredMssToGDALGTiff(reg.Params.MssTiffFile)
+		// reg.SaveRegisteredMssToGDALGTiff(reg.Params.MssTiffFile)
 
 		if params.PansharpenIHS {
 			reg.DoScenePansharpen(panScenes, mssScenes)

producer/gdal_pansharpen.go

@@ -104,7 +104,6 @@ func GDALTranslate(srcTiff, dstJPG string) error {
 	}
 	// gdal_translate -of JPEG -outsize 50% 50% -r bilinear input.tif output.jpg
 	cmd := exec.Command("gdal_translate", "-of", "JPEG",
-		"-scale",
 		srcTiff, dstJPG)
 	log.Println("Run Command: ", cmd.String())
 	cmd.Stdout = os.Stdout

producer/image_proc_test.go

@@ -18,8 +18,23 @@ func TestGDALPansharpen(t *testing.T) {
 
 func TestGTiffToJPG(t *testing.T) {
 	godal.RegisterAll()
-	err := GTiffToJPG("data/052022/010/SJY01_MSS_20240520_115428_052022_103_010.tiff",
+	tiff := "../data/051622/006/SJY01_MSS_20240516_101236_051622_096_006.tiff"
-		"data/052022/SJY01_MSS_20240520_115428_052022_103_010.jpg")
+	jpg := "../data/051622/006/SJY01_MSS_20240516_101236_051622_096_006.jpg"
+	err := GTiffToJPG(tiff, jpg, false)
+	if err != nil {
+		t.Error(err)
+	}
+
+	tiff = "../data/051622/006/SJY01_PAN_20240516_101236_051622_096_006.tiff"
+	jpg = "../data/051622/006/SJY01_PAN_20240516_101236_051622_096_006.jpg"
+	err = GTiffToJPG(tiff, jpg, false)
+	if err != nil {
+		t.Error(err)
+	}
+
+	tiff = "../data/051622/006/SJY01_FUS_20240516_101236_051622_096_006.tiff"
+	jpg = "../data/051622/006/SJY01_FUS_20240516_101236_051622_096_006.jpg"
+	err = GTiffToJPG(tiff, jpg, true)
 	if err != nil {
 		t.Error(err)
 	}

producer/image_registration.go

@@ -201,7 +201,7 @@ func (r *Registrator) calcPhaseCorrelation(panImage gocv.Mat,
 				y1 += OverlappedBlockLines // Y偏移量过大，需要将重叠块的行数加上，以避免边界影响
 
 				if x1 > width || y1 > height {
-					log.Warnf("Block out of range. x0=%d, y0=%d, x1=%d, y1=%d", x0, y0, x1, y1)
+					log.Debugf("Block out of range. x0=%d, y0=%d, x1=%d, y1=%d", x0, y0, x1, y1)
 				}
 
 				if y1 > height {

producer/jpg.go

@@ -9,7 +9,7 @@ import (
 	"gocv.io/x/gocv"
 )
 
-func GTiffToJPG(ftiff, fjpg string) error {
+func GTiffToJPG(ftiff, fjpg string, reversed bool) error {
 	// 打开 TIFF 文件
 	ds, err := godal.Open(ftiff)
 	if err != nil {
@@ -20,62 +20,61 @@ func GTiffToJPG(ftiff, fjpg string) error {
 
 	bands := ds.Bands()
 	log.Infof("TIFF file %s has %d bands", ftiff, len(bands))
-	if len(bands) < 3 {
-		err = fmt.Errorf("TIFF file %s has less than 3 bands", ftiff)
-		log.Error(err)
-		return err
-	}
 
 	// 获取图像大小
 	width := bands[0].Structure().SizeX
 	height := bands[0].Structure().SizeY
+	bandsCnt := len(bands)
 
-	// 读取每个波段并转换为 8 位
+	// 创建 16 位图像矩阵
-	img := gocv.NewMatWithSize(height, width, gocv.MatTypeCV16UC3)
+	var img gocv.Mat
-	for i := 1; i <= 3; i++ {
+	defer img.Close()
-		band := ds.Bands()[i-1]
+	if bandsCnt == 1 {
+		img = gocv.NewMatWithSize(height, width, gocv.MatTypeCV16UC1)
+	} else {
+		bandsCnt = 3 // 只取前三个波段
+		img = gocv.NewMatWithSize(height, width, gocv.MatTypeCV16UC3)
+	}
+
+	// 读取每个波段并存储到图像矩阵中
+	for i := 0; i < bandsCnt; i++ {
+		band := bands[i]
 		data := make([]uint16, width*height)
 		err = band.Read(0, 0, data, width, height)
 		if err != nil {
-			log.Printf("Error reading band %d: %v", i, err)
+			fmt.Printf("Error reading band %d: %v\n", i, err)
 			return err
 		}
 
-		// 计算最小值和最大值
+		// 将 16 位数据存储到图像的相应通道
-		var minVal, maxVal uint16 = 65535, 0
-		for _, value := range data {
-			if value < minVal {
-				minVal = value
-			}
-			if value > maxVal {
-				maxVal = value
-			}
-		}
-
-		// 将 16 位数据缩放到 8 位并存储到图像的相应通道
 		for y := 0; y < height; y++ {
 			for x := 0; x < width; x++ {
 				value := data[y*width+x]
-				scaledValue := uint8(float64(value-minVal) / float64(maxVal-minVal) * 255)
+				img.SetShortAt(y, x*bandsCnt+i, int16(value))
-				img.SetUCharAt(y, x*3+(i-1), scaledValue)
 			}
 		}
-
-		// 将 16 位数据存储到图像的相应通道
-		// for y := 0; y < height; y++ {
-		// 	for x := 0; x < width; x++ {
-		// 		value := data[y*width+x]
-		// 		img.SetShortAt(y, x*3+(i-1), int16(value))
-		// 	}
-		// }
 	}
 
-	// 调整图像大小
+	// 重采样图像
-	resizedImg := gocv.NewMat()
+	gocv.Resize(img, &img, image.Point{X: img.Cols() / 2, Y: img.Rows() / 2}, 0, 0, gocv.InterpolationCubic)
-	gocv.Resize(img, &resizedImg, image.Point{X: 2336, Y: 2336}, 0, 0, gocv.InterpolationCubic)
 
-	// 保存为 JPG 文件
+	// 创建 8 位图像矩阵
-	ok := gocv.IMWrite(fjpg, resizedImg)
+	img8 := gocv.NewMat()
+	defer img8.Close()
+	if bandsCnt == 1 {
+		img.ConvertToWithParams(&img8, gocv.MatTypeCV8UC1, 1.0/256.0*2, 0)
+	} else {
+		img.ConvertToWithParams(&img8, gocv.MatTypeCV8UC3, 1.0/256.0*2, 0)
+		if reversed {
+			channels := gocv.Split(img8)
+			gocv.Merge([]gocv.Mat{channels[2], channels[1], channels[0]}, &img8)
+			for _, ch := range channels {
+				ch.Close()
+			}
+		}
+	}
+
+	ok := gocv.IMWrite(fjpg, img8)
 	if !ok {
 		err = fmt.Errorf("error saving %s", fjpg)
 		return err

producer/meta.go

@@ -2,7 +2,7 @@ package imageproc
 
 import (
 	"encoding/xml"
-	"io/ioutil"
+	"os"
 )
 
 // 定义与XML结构对应的Go结构体
@@ -61,7 +61,7 @@ func WriteProductMeta(productMeta *ProductMeta, filename string) error {
 		return err
 	}
 
-	err = ioutil.WriteFile(filename, output, 0644)
+	err = os.WriteFile(filename, output, 0644)
 	if err != nil {
 		return err
 	}

producer/output.go

@@ -23,7 +23,7 @@ func (r *Registrator) SaveFilteredPanToGDALGTiff(tiffFile string) error {
 }
 
 func savePanToGDALGTiff(pan gocv.Mat, tiffFile string) error {
-	log.Println("Saving PAN image to TIFF file:", tiffFile)
+	// log.Println("Saving PAN image to TIFF file:", tiffFile)
 
 	width := pan.Cols()
 	height := pan.Rows()
@@ -65,28 +65,12 @@ func savePanToGDALGTiff(pan gocv.Mat, tiffFile string) error {
 }
 
 func (r *Registrator) SaveRegisteredMssToGDALGTiff(tiffFile string) error {
-	width := r.MssWidth
+	r.rgbirImage = gocv.NewMat()
-	height := r.MssHeight
+	gocv.Merge(r.registeredMssImages[:], &r.rgbirImage)
-
-	// 创建合并后的图像（RGBIR）
-	r.rgbirImage = gocv.NewMatWithSize(height, width, gocv.MatTypeCV16UC4) // 4通道，16位
-
-	for y := 0; y < height; y++ {
-		for x := 0; x < width; x++ {
-			blue := r.registeredMssImages[0].GetShortAt(y, x)
-			green := r.registeredMssImages[1].GetShortAt(y, x)
-			red := r.registeredMssImages[2].GetShortAt(y, x)
-			ir := r.registeredMssImages[3].GetShortAt(y, x)
-			r.rgbirImage.SetShortAt(y, x*4+0, blue)
-			r.rgbirImage.SetShortAt(y, x*4+1, green)
-			r.rgbirImage.SetShortAt(y, x*4+2, red)
-			r.rgbirImage.SetShortAt(y, x*4+3, ir)
-		}
-	}
 
 	return SaveBGRToGDALGTiff(r.rgbirImage,
 		4, 5,
-		[]godal.ColorInterp{godal.CIBlue, godal.CIGreen, godal.CIRed, godal.CIAlpha},
+		[]godal.ColorInterp{godal.CIBlue, godal.CIGreen, godal.CIRed, godal.CIUndefined},
 		tiffFile)
 }
 
@@ -99,8 +83,6 @@ func (r *Registrator) SavePansharpenedToGDALGTiff(tiffFile string) error {
 }
 
 func SaveBGRToGDALGTiff(bgr gocv.Mat, bands int, resolution float64, colorInterps []godal.ColorInterp, tiffFile string) error {
-	log.Println("Saving BGR to TIFF file:", tiffFile)
-
 	width := bgr.Cols()
 	height := bgr.Rows()
 
@@ -130,7 +112,7 @@ func SaveBGRToGDALGTiff(bgr gocv.Mat, bands int, resolution float64, colorInterp
 	}
 	defer ds.Close()
 
-	ds.SetMetadata("NBITS", "16")
+	// ds.SetMetadata("NBITS", "16")
 
 	setGeoTransform(ds, 0, 0, float64(width), float64(height), resolution)
 

producer/proj.go

@@ -3,13 +3,14 @@ package imageproc
 import "fmt"
 
 func calculateProj(topLeftX, topLeftY, width, height, resolution float64) string {
+	var projections []string
 	// 计算图像的地理范围
 	bottomRightX := topLeftX + float64(width)*resolution
 	bottomRightY := topLeftY + float64(height)*resolution
 
 	// 根据地理范围和分辨率选择适当的投影信息
 	// 这里只是一个示例，你需要根据实际情况选择合适的投影系统和参数
-	projection := fmt.Sprintf(`PROJCS["Custom Projection",
+	projections = append(projections, fmt.Sprintf(`PROJCS["Custom Projection",
 	GEOGCS["WGS 84",
 		DATUM["WGS_1984",
 			SPHEROID["WGS 84",6378137,298.257223563,
@@ -30,9 +31,9 @@ func calculateProj(topLeftX, topLeftY, width, height, resolution float64) string
 		AUTHORITY["EPSG","9001"]],
 	AXIS["Easting",EAST],
 	AXIS["Northing",NORTH]]`,
-		(topLeftX+bottomRightX)/2, (topLeftY+bottomRightY)/2)
+		(topLeftX+bottomRightX)/2, (topLeftY+bottomRightY)/2))
 
-	projection = `PROJCS["WGS 84 / UTM zone 51N",
+	projections = append(projections, `PROJCS["WGS 84 / UTM zone 51N",
 	GEOGCS["WGS 84",
 		DATUM["WGS_1984",
 			SPHEROID["WGS 84",6378137,298.257223563,
@@ -53,9 +54,9 @@ func calculateProj(topLeftX, topLeftY, width, height, resolution float64) string
 		AUTHORITY["EPSG","9001"]],
 	AXIS["Easting",EAST],
 	AXIS["Northing",NORTH],
-	AUTHORITY["EPSG","32651"]]`
+	AUTHORITY["EPSG","32651"]]`)
 
-	projection = `PROJCS["WGS 84 / Pseudo-Mercator",
+	projections = append(projections, `PROJCS["WGS 84 / Pseudo-Mercator",
 	GEOGCS["WGS 84",DATUM["WGS_1984",
 		SPHEROID["WGS 84",6378137,298.257223563,
 			AUTHORITY["EPSG","7030"]],
@@ -72,8 +73,8 @@ func calculateProj(topLeftX, topLeftY, width, height, resolution float64) string
 	AXIS["X",EAST],
 	AXIS["Y",NORTH],
 	EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],
-	AUTHORITY["EPSG","3857"]]`
+	AUTHORITY["EPSG","3857"]]`)
 
 	// 输出投影信息
-	return projection
+	return projections[0]
 }

producer/report.go

@@ -1,6 +1,9 @@
 package imageproc
 
-import "encoding/xml"
+import (
+	"encoding/xml"
+	"os"
+)
 
 // Report represents the root XML element
 type Report struct {
@@ -26,3 +29,17 @@ type ReportScene struct {
 	MetaData       string `xml:"metaData"`
 	QualityData    string `xml:"qualityData"`
 }
+
+func WriteReport(report *Report, filename string) error {
+	output, err := xml.MarshalIndent(report, "", "    ")
+	if err != nil {
+		return err
+	}
+
+	err = os.WriteFile(filename, output, 0644)
+	if err != nil {
+		return err
+	}
+
+	return nil
+}

producer/scenes.go

@@ -114,7 +114,7 @@ func (r *Registrator) SaveScenesToTiff(panScenes []*Scene, mssScenes []*Scene) e
 		rgbirImage, _ := r.MergeMSSToBGRNIR(scene.Mat)
 		err := SaveBGRToGDALGTiff(rgbirImage,
 			4, 5,
-			[]godal.ColorInterp{godal.CIBlue, godal.CIGreen, godal.CIRed, godal.CIAlpha},
+			[]godal.ColorInterp{godal.CIBlue, godal.CIGreen, godal.CIRed, godal.CIUndefined},
 			scene.Tiff)
 		if err != nil {
 			log.Errorf("save mss scene %d to tiff failed: %v", i+1, err)
@@ -143,24 +143,11 @@ func (r *Registrator) MergeMSSToBGRNIR(channels []gocv.Mat) (gocv.Mat, error) {
 	if len(channels) != 4 {
 		return rgbirImage, fmt.Errorf("mss channels count not match")
 	}
-	width := channels[0].Cols()
-	height := channels[0].Rows()
 
-	// 创建合并后的图像（RGBIR）
+	rgbirImage = gocv.NewMat()
-	rgbirImage = gocv.NewMatWithSize(height, width, gocv.MatTypeCV16UC4) // 4通道，16位
+	gocv.Merge(channels, &rgbirImage)
-
+	log.Printf("merge mss to bgr nir image, channels: %d, height: %d, width: %d",
-	for y := 0; y < height; y++ {
+		rgbirImage.Channels(), rgbirImage.Rows(), rgbirImage.Cols())
-		for x := 0; x < width; x++ {
-			blue := channels[0].GetShortAt(y, x)
-			green := channels[1].GetShortAt(y, x)
-			red := channels[2].GetShortAt(y, x)
-			ir := channels[3].GetShortAt(y, x)
-			rgbirImage.SetShortAt(y, x*4+0, blue)
-			rgbirImage.SetShortAt(y, x*4+1, green)
-			rgbirImage.SetShortAt(y, x*4+2, red)
-			rgbirImage.SetShortAt(y, x*4+3, ir)
-		}
-	}
 
 	return rgbirImage, nil