pansharpen

cmd/main.go

@@ -27,18 +27,18 @@ func main() {
 	}
 
 	godal.RegisterAll()
-	if err := reg.CalcDownPhaseCorrelation(); err != nil {
+
+	// reg.SaveOriginalPanToGDALGTiff("data/original_pan.tiff")
+	// reg.SaveFilteredPanToGDALGTiff("data/filtered_pan.tiff")
+
+	if err := reg.DoPhaseCorrelation(); err != nil {
 		panic(err)
 	}
-	// reg.DoCoRegestration()
+	reg.DoCoRegestration()
 
-	// if err := reg.DoUpPhaseCorrelation(); err != nil {
+	// // reg.SaveRegisteredMssToRaw("data/registered_mss.RAW")
-	// 	panic(err)
+	reg.SaveRegisteredMssToGDALGTiff("data/registered_mss.tiff")
-	// }
+	reg.SavePansharpenedToGDALGTiff("data/pansharpened.tiff")
-
-	// reg.SaveRegisteredMssToRaw("data/registered_mss.RAW")
-	// reg.SavePanToGDALGTiff("data/original_pan.tiff")
-	// reg.SaveRegisteredMssToGDALGTiff("data/registered_mss.tiff")
 
 	reg.Clean()
 }

filter.go

@@ -0,0 +1,78 @@
+package imageproc
+
+import (
+	"image"
+	"math"
+
+	log "github.com/sirupsen/logrus"
+
+	"gocv.io/x/gocv"
+)
+
+func PANFilter(panImage16UC1 gocv.Mat) gocv.Mat {
+	log.Println("Applying PAN filter...")
+	// 转换为浮点数类型进行傅里叶变换
+	imgFloat := gocv.NewMat()
+	panImage16UC1.ConvertTo(&imgFloat, gocv.MatTypeCV32F)
+
+	// 使用傅里叶变换
+	planes := gocv.NewMat()
+	gocv.Merge([]gocv.Mat{imgFloat, gocv.NewMatWithSize(imgFloat.Rows(), imgFloat.Cols(), gocv.MatTypeCV32F)}, &planes)
+	dft := gocv.NewMat()
+	gocv.DFT(planes, &dft, gocv.DftComplexOutput)
+
+	// 转换DFT图像，使低频分量位于中心
+	dftShifted := shiftDFT(dft)
+
+	// 创建掩膜
+	rows, cols := panImage16UC1.Rows(), panImage16UC1.Cols()
+	crow, ccol := rows/2, cols/2
+	mask := gocv.NewMatWithSize(rows, cols, gocv.MatTypeCV32FC2)
+	mask.SetTo(gocv.NewScalar(1.0, 1.0, 0.0, 0.0))
+
+	r := 30 // 调整这个参数来改变滤波效果
+	for i := 0; i < rows; i++ {
+		for j := 0; j < cols; j++ {
+			if math.Abs(float64(i-crow)) <= float64(r) && math.Abs(float64(j-ccol)) <= float64(r) {
+				mask.SetFloatAt(i, j*2, 0)
+				mask.SetFloatAt(i, j*2+1, 0)
+			}
+		}
+	}
+
+	// 应用掩膜并进行反向DFT
+	filtered := gocv.NewMat()
+	gocv.MulSpectrums(dftShifted, mask, &filtered, 0)
+	filteredShifted := shiftDFT(filtered)
+	idft := gocv.NewMat()
+	gocv.DFT(filteredShifted, &idft, gocv.DftInverse|gocv.DftRealOutput)
+
+	// 标准化图像到0-65535范围
+	gocv.Normalize(idft, &idft, 0, 65535, gocv.NormMinMax)
+	idft.ConvertTo(&idft, gocv.MatTypeCV16U)
+
+	return idft
+}
+
+// shiftDFT 将DFT结果进行频谱平移
+func shiftDFT(src gocv.Mat) gocv.Mat {
+	rows, cols := src.Rows(), src.Cols()
+	crow, ccol := rows/2, cols/2
+
+	q0 := src.Region(image.Rect(0, 0, ccol, crow))
+	q1 := src.Region(image.Rect(ccol, 0, cols, crow))
+	q2 := src.Region(image.Rect(0, crow, ccol, rows))
+	q3 := src.Region(image.Rect(ccol, crow, cols, rows))
+
+	tmp := gocv.NewMatWithSize(crow, ccol, src.Type())
+
+	q0.CopyTo(&tmp)
+	q3.CopyTo(&q0)
+	tmp.CopyTo(&q3)
+
+	q1.CopyTo(&tmp)
+	q2.CopyTo(&q1)
+	tmp.CopyTo(&q2)
+
+	return src
+}

gdal_pansharpen.go

@@ -1,5 +1,55 @@
 package imageproc
 
+import (
+	"fmt"
+	"image"
+
+	"gocv.io/x/gocv"
+)
+
 const (
 	GDALPansharpenCommand = "gdal_pansharpen.py {{.PANGTiff}} {{.MSSGTiff}} {{.FUSGTiff}} -r cubic -of GTiff"
 )
+
+// IHS 变换是一种将 RGB 空间转换为亮度、色调和饱和度的变换方法
+func PansharpenIHS(panImage, mssImage gocv.Mat) gocv.Mat {
+	pan32 := gocv.NewMat()
+	defer pan32.Close()
+	panImage.ConvertTo(&pan32, gocv.MatTypeCV32FC1)
+
+	mss32 := gocv.NewMat()
+	mssImage.ConvertTo(&mss32, gocv.MatTypeCV32FC4)
+
+	// 调整 MSS 图像的大小到与 PAN 图像相同的大小
+	mssResized := gocv.NewMat()
+	defer mssResized.Close()
+	gocv.Resize(mss32, &mssResized,
+		image.Point{X: panImage.Cols(), Y: panImage.Rows()}, 0, 0, gocv.InterpolationCubic)
+
+	// 将 MSS 图像从 BGR 转换为 HLS
+	mssHLS := gocv.NewMat()
+	defer mssHLS.Close()
+	gocv.CvtColor(mssResized, &mssHLS, gocv.ColorBGRToHLS)
+
+	// 分离 HLS 通道 只有3个通道
+	hlsChannels := gocv.Split(mssHLS)
+
+	fmt.Println("size of channels:", len(hlsChannels), mssImage.Channels(), mssResized.Channels(), mssHLS.Channels())
+	fmt.Println("size of pan32:", pan32.Size(), pan32.Channels())
+	fmt.Println("size of resized mss", mss32.Size(), mss32.Rows(), mss32.Cols())
+
+	// 使用 PAN 图像替换亮度分量（Intensity）- 通常替换的是 G 通道
+	hlsChannels[1] = pan32
+
+	// 合并通道
+	gocv.Merge(hlsChannels, &mssHLS)
+
+	// 将图像从 HLS 转换回 BGR
+	fusedImage := gocv.NewMat()
+	defer fusedImage.Close()
+	gocv.CvtColor(mssHLS, &fusedImage, gocv.ColorHLSToBGR)
+	fused16 := gocv.NewMat()
+	fusedImage.ConvertTo(&fused16, gocv.MatTypeCV16UC3)
+
+	return fused16
+}

image_registration.go

@@ -21,9 +21,9 @@ const (
 	PixelBytes                          = 2
 	PanWidth                            = 9344 // 像素宽度
 	MssWidth                            = 2336
-	BlockNH                             = 8
+	BlockNH                             = 4
-	BlockNW                             = 4
+	BlockNW                             = 16
-	OverlappedBlockLines                = 2000 // 重叠块的行数
+	OverlappedBlockLines                = 3000 // 重叠块的行数
 	DownSampled          ResampleMethod = "down_sample_pan"
 	UpSampled            ResampleMethod = "up_sample_mss"
 )
@@ -41,8 +41,8 @@ type Registrator struct {
 
 	shiftMutex   sync.Mutex
 	phaseShifts  [4][]PhaseShiftM
-	deltaXCoeffs [4][]float64 // Polynomial fitting coefficients: 图像内畸变（非线性变换），捕捉图像在水平方向上引起的垂直方向的变形
+	deltaXCoeffs [4][]float64 // 图像内畸变（非线性变换），捕捉图像在水平方向上引起的垂直方向的变形
-	deltaYCoeffs [4][]float64 // Polynomial fitting coefficients: 图像内畸变（非线性变换），捕捉图像在垂直方向上引起的水平方向的变形
+	deltaYCoeffs [4][]float64 // 图像内畸变（非线性变换），捕捉图像在垂直方向上引起的水平方向的变形
 
 	registeredMssImages [4]gocv.Mat // 配准后的MSS图像
 	rgbirImage          gocv.Mat
@@ -52,6 +52,7 @@ type Registrator struct {
 
 func NewRegistrator() *Registrator {
 	var r Registrator
+	r.resampleMethod = DownSampled
 
 	return &r
 }
@@ -121,6 +122,15 @@ func (r *Registrator) LoadMssRaw(raw string) error {
 	return nil
 }
 
+func (r *Registrator) DoPhaseCorrelation() error {
+	switch r.resampleMethod {
+	case UpSampled:
+		return r.CalcUpPhaseCorrelation()
+	default:
+		return r.CalcDownPhaseCorrelation()
+	}
+}
+
 // 将PAN降采样后计算相位相关的偏移量
 func (r *Registrator) CalcDownPhaseCorrelation() error {
 	// 确保 MSS 高度是 PAN 高度的 1/4
@@ -134,63 +144,14 @@ func (r *Registrator) CalcDownPhaseCorrelation() error {
 	downsampledPanImage := gocv.NewMat()
 	gocv.Resize(r.PanImage, &downsampledPanImage,
 		image.Point{X: r.MssWidth, Y: r.MssHeight}, 0, 0, gocv.InterpolationCubic)
-	fmt.Println("down sampled PAN images size:", downsampledPanImage.Size())
+	log.Println("down sampled PAN images size:", downsampledPanImage.Size())
 
 	// 分块高度
 	blockHeight := r.MssHeight / BlockNH
-	for band := 0; band < MssBands; band++ {
+	blockWidth := r.MssWidth / BlockNW
-		for bh := 0; bh < BlockNH; bh++ {
-			// 读取 PAN 和 MSS 分块数据
-			y1 := (bh+1)*blockHeight + 800
-			if y1 > r.MssHeight {
-				y1 = r.MssHeight
-			}
 
-			var shiftM PhaseShiftM
+	return r.calcPhaseCorrelation(downsampledPanImage, r.MssImages, r.MssHeight, r.MssWidth, blockHeight, blockWidth)
-			shiftM.Block.width = r.MssWidth           // 块宽度
-			shiftM.Block.height = y1 - bh*blockHeight // 块高度
-			shiftM.Block.coord.X = 0                  // 块左上角x坐标
-			shiftM.Block.coord.Y = bh * blockHeight   // 块左上角y坐标
 
-			rect := image.Rect(
-				shiftM.Block.coord.X, shiftM.Block.coord.Y,
-				shiftM.Block.coord.X+shiftM.Block.width, shiftM.Block.coord.Y+shiftM.Block.height,
-			)
-			log.Println("Band", band+1, ", processing block", bh, rect)
-
-			panBlock := downsampledPanImage.Region(rect)
-			mssBlock := r.MssImages[band].Region(rect)
-
-			// 处理每个分块
-			phaseShift, response := r.calculateBlockPhaseShift(panBlock, mssBlock)
-			shiftM.dx = phaseShift.X
-			shiftM.dy = phaseShift.Y
-			shiftM.response = response
-
-			r.phaseShifts[band] = append(r.phaseShifts[band], shiftM)
-
-			panBlock.Close()
-			mssBlock.Close()
-		}
-	}
-
-	// if err := r.DoMssPhaseShift(); err != nil {
-	// 	log.Error("Error calculating MSS phase shift: ", err)
-	// 	return err
-	// }
-
-	for i := 0; i < MssBands; i++ {
-		for j, shift := range r.phaseShifts[i] {
-			if shift.response > 0.4 || shift.dy > 8 {
-				fmt.Printf("Band %d, block %d, dx=%f, dy=%f, response=%f\n",
-					i, j, shift.dx, shift.dy, shift.response)
-			}
-		}
-	}
-
-	r.calcDeltaCoeffs()
-
-	return nil
 }
 
 // 将MSS升采样采样后计算相位相关的偏移量
@@ -218,6 +179,13 @@ func (r *Registrator) CalcUpPhaseCorrelation() error {
 
 	log.Infof("blockHeight=%d, blockWidth=%d", blockHeight, blockWidth)
 
+	return r.calcPhaseCorrelation(r.PanImage, upsampledMssImages, r.PanHeight, r.PanWidth, blockHeight, blockWidth)
+}
+
+func (r *Registrator) calcPhaseCorrelation(panImage gocv.Mat,
+	mssImages [4]gocv.Mat,
+	height, width,
+	blockHeight, blockWidth int) error {
 	var wg sync.WaitGroup
 
 	for bh := 0; bh < BlockNH; bh++ {
@@ -231,12 +199,12 @@ func (r *Registrator) CalcUpPhaseCorrelation() error {
 				y1 := (bh + 1) * blockHeight
 				y1 += OverlappedBlockLines // Y偏移量过大，需要将重叠块的行数加上，以避免边界影响
 
-				if x1 > r.PanWidth || y1 > r.PanHeight {
+				if x1 > width || y1 > height {
 					log.Warnf("Block out of range. x0=%d, y0=%d, x1=%d, y1=%d", x0, y0, x1, y1)
 				}
 
-				if y1 > r.PanHeight {
+				if y1 > height {
-					y1 = r.PanHeight
+					y1 = height
 				}
 
 				var shiftM PhaseShiftM
@@ -250,10 +218,10 @@ func (r *Registrator) CalcUpPhaseCorrelation() error {
 					x1, y1,
 				)
 
-				panBlock := r.PanImage.Region(rect)
+				panBlock := panImage.Region(rect)
 				for band := 0; band < MssBands; band++ {
 					log.Println("Band", band+1, ", processing block", bh, rect)
-					mssBlock := upsampledMssImages[band].Region(rect)
+					mssBlock := mssImages[band].Region(rect)
 
 					// 处理每个分块
 					phaseShift, response := r.calculateBlockPhaseShift(panBlock, mssBlock)
@@ -277,12 +245,14 @@ func (r *Registrator) CalcUpPhaseCorrelation() error {
 	for i := 0; i < MssBands; i++ {
 		for _, shift := range r.phaseShifts[i] {
 			if shift.response > 0.4 || shift.dx > 8 || shift.dy > 8 {
-				fmt.Printf("Band %d, block %d, dx=%f, dy=%f, response=%f\n",
+				log.Printf("Band %d, block %d, dx=%f, dy=%f, response=%f\n",
 					i, shift.Block.coord.X, shift.dx, shift.dy, shift.response)
 			}
 		}
 	}
 
+	r.calcDeltaCoeffs()
+
 	return nil
 }
 
@@ -313,127 +283,6 @@ func (r *Registrator) SaveRegisteredMssToRaw(raw string) error {
 	return nil
 }
 
-func (r *Registrator) bytesToRaw(mssData []byte, filePath string) error {
-	f, err := os.OpenFile(filePath, os.O_CREATE|os.O_TRUNC|os.O_WRONLY, 0777)
-	if err != nil {
-		return err
-	}
-	w := bufio.NewWriter(f)
-	w.Write(mssData)
-
-	return nil
-}
-
-func (r *Registrator) SaveRegisteredMssToGDALGTiff(tiffFile string) error {
-	log.Println("Saving registered MSS to TIFF file:", tiffFile)
-
-	width := r.MssWidth
-	height := r.MssHeight
-
-	// 创建合并后的图像（RGBIR）
-	r.rgbirImage = gocv.NewMatWithSize(height, width, gocv.MatTypeCV16UC4) // 4通道，16位
-
-	for y := 0; y < height; y++ {
-		for x := 0; x < width; x++ {
-			red := r.registeredMssImages[0].GetShortAt(y, x)
-			green := r.registeredMssImages[1].GetShortAt(y, x)
-			blue := r.registeredMssImages[2].GetShortAt(y, x)
-			ir := r.registeredMssImages[3].GetShortAt(y, x)
-			r.rgbirImage.SetShortAt(y, x*4+0, red)
-			r.rgbirImage.SetShortAt(y, x*4+1, green)
-			r.rgbirImage.SetShortAt(y, x*4+2, blue)
-			r.rgbirImage.SetShortAt(y, x*4+3, ir)
-		}
-	}
-
-	// 创建一个二维切片来存储图像数据
-	data := make([][]uint16, MssBands)
-	for i := range data {
-		data[i] = make([]uint16, width*height)
-	}
-
-	// 从gocv.Mat中提取数据
-	for y := 0; y < height; y++ {
-		for x := 0; x < width; x++ {
-			for b := 0; b < MssBands; b++ {
-				data[b][y*width+x] = uint16(r.rgbirImage.GetShortAt(y, x*4+b))
-			}
-		}
-	}
-
-	ds, err := godal.Create(godal.GTiff,
-		tiffFile,
-		MssBands,
-		godal.UInt16,
-		width, height)
-	if err != nil {
-		log.Error("Error creating TIFF file: ", err)
-		return err
-	}
-	defer ds.Close()
-
-	setGeoTransform(ds, 0, 0, float64(width), float64(height), 1.2*4)
-
-	for b := 0; b < MssBands; b++ {
-		band := ds.Bands()[b]
-		err := band.IO(godal.IOWrite,
-			0, 0,
-			data[b],
-			width, height,
-			godal.PixelSpacing(2),
-			godal.LineSpacing(width*2))
-		if err != nil {
-			log.Error("Failed to write data to band:", err)
-			return err
-		}
-	}
-
-	log.Info("Saved registered mss to ", tiffFile)
-
-	return nil
-}
-
-func (r *Registrator) SavePanToGDALGTiff(tiffFile string) error {
-	log.Println("Saving PAN image to TIFF file:", tiffFile)
-
-	width := r.PanWidth
-	height := r.PanHeight
-
-	ds, err := godal.Create(godal.GTiff, tiffFile, 1, godal.UInt16, width, height)
-	if err != nil {
-		log.Error("Error creating TIFF file: ", err)
-		return err
-	}
-	defer ds.Close()
-
-	setGeoTransform(ds, 0, 0, float64(width), float64(height), 1.2)
-	ds.SetMetadata("NBITS", "16")
-
-	// 将通道的数据转换为uint16数组
-	data := make([]uint16, width*height)
-	for y := 0; y < height; y++ {
-		for x := 0; x < width; x++ {
-			data[y*width+x] = uint16(r.PanImage.GetShortAt(y, x))
-		}
-	}
-
-	band := ds.Bands()[0]
-	err = band.IO(godal.IOWrite,
-		0, 0,
-		data,
-		width, height,
-		godal.PixelSpacing(2),
-		godal.LineSpacing(width*2))
-	if err != nil {
-		log.Error("Failed to write data to band:", err)
-		return err
-	}
-
-	log.Info("Saved pan image to ", tiffFile)
-
-	return nil
-}
-
 func (r *Registrator) Clean() {
 	r.PanImage.Close()
 	for i := 0; i < MssBands; i++ {
@@ -454,14 +303,19 @@ func (r *Registrator) calcDeltaCoeffs() error {
 		var dx []float64
 		var dy []float64
 		effectShift := 0
-		for j, shift := range r.phaseShifts[i] {
+		for _, shift := range r.phaseShifts[i] {
 			if math.IsNaN(float64(shift.dx)) || math.IsNaN(float64(shift.dy)) {
 				continue
 			}
 
+			// 经验值过滤
+			if shift.dy < 64.0 {
+				continue
+			}
+
 			effectShift++
-			cx = append(cx, float64(shift.Block.coord.X+j)) // MSS 块在X方向没有分块
+			cx = append(cx, float64(shift.Block.coord.X+shift.Block.width/2)) // MSS 块在X方向没有分块
-			fmt.Println("effectShift:", effectShift, "cx:", shift.Block.coord.X, "dy:", shift.dy)
+			log.Debug("effective shift value:", effectShift, "cx:", shift.Block.coord.X, "dy:", shift.dy)
 			dx = append(dx, float64(shift.dx))
 			dy = append(dy, float64(shift.dy))
 
@@ -487,11 +341,25 @@ func (r *Registrator) DoCoRegestration() error {
 		mapY := gocv.NewMatWithSize(r.MssHeight, r.MssWidth, gocv.MatTypeCV32FC1)
 		for y := 0; y < r.MssHeight; y++ {
 			for x := 0; x < r.MssWidth; x++ {
-				// dx := r.deltaXCoeffs[band][1]*float64(x) + r.deltaXCoeffs[band][0] + float64(x)
+				var dx, dy float64
-				// dy := r.deltaYCoeffs[band][2]*float64(x*x) + r.deltaYCoeffs[band][1]*float64(x) + r.deltaYCoeffs[band][0] + float64(y)
+				if r.resampleMethod == UpSampled {
-				// fmt.Println("x:", x, "dx:", dx, "y:", y, "dy:", dy)
+					xx := float64(x * MssBands)
-				mapX.SetFloatAt(y, x, float32(x)+float32(r.deltaXCoeffs[band][0]))
+					yy := float64(y * MssBands)
-				mapY.SetFloatAt(y, x, float32(y)+float32(r.deltaYCoeffs[band][0]))
+					dx = (r.deltaXCoeffs[band][1]*float64(xx) + r.deltaXCoeffs[band][0] + xx) / MssBands
+					dy = (r.deltaYCoeffs[band][2]*float64(xx*xx) + r.deltaYCoeffs[band][1]*float64(xx) + r.deltaYCoeffs[band][0] + yy) / MssBands
+				} else {
+					dx = r.deltaXCoeffs[band][1]*float64(x) + r.deltaXCoeffs[band][0] + float64(x)
+					dy = r.deltaYCoeffs[band][2]*float64(x*x) + r.deltaYCoeffs[band][1]*float64(x) + r.deltaYCoeffs[band][0] + float64(y)
+				}
+
+				// if band+1 == 4 {
+				// 	fmt.Println("band:", band+1, "x:", x, "map_x:", mx, "y:", y, "map_y:", my)
+				// }
+
+				// mapX.SetFloatAt(y, x, float32(x)+float32(r.deltaXCoeffs[band][0]))
+				// mapY.SetFloatAt(y, x, float32(y)+float32(r.deltaYCoeffs[band][0]))
+				mapX.SetFloatAt(y, x, float32(dx))
+				mapY.SetFloatAt(y, x, float32(dy))
 			}
 		}
 

output.go

@@ -0,0 +1,153 @@
+package imageproc
+
+import (
+	"bufio"
+	"os"
+
+	"github.com/airbusgeo/godal"
+	log "github.com/sirupsen/logrus"
+	"gocv.io/x/gocv"
+)
+
+func (r *Registrator) SaveOriginalPanToGDALGTiff(tiffFile string) error {
+	return savePanToGDALGTiff(r.PanImage, tiffFile)
+}
+
+func (r *Registrator) SaveFilteredPanToGDALGTiff(tiffFile string) error {
+	img := PANFilter(r.PanImage)
+	img.ConvertTo(&img, gocv.MatTypeCV16U)
+	return savePanToGDALGTiff(img, tiffFile)
+}
+
+func savePanToGDALGTiff(pan gocv.Mat, tiffFile string) error {
+	log.Println("Saving PAN image to TIFF file:", tiffFile)
+
+	width := pan.Rows()
+	height := pan.Cols()
+
+	ds, err := godal.Create(godal.GTiff, tiffFile, 1, godal.UInt16, width, height)
+	if err != nil {
+		log.Error("Error creating TIFF file: ", err)
+		return err
+	}
+	defer ds.Close()
+
+	setGeoTransform(ds, 0, 0, float64(width), float64(height), 1.25)
+	ds.SetMetadata("NBITS", "16")
+
+	// 将通道的数据转换为uint16数组
+	data := make([]uint16, width*height)
+	for y := 0; y < height; y++ {
+		for x := 0; x < width; x++ {
+			data[y*width+x] = uint16(pan.GetShortAt(y, x))
+		}
+	}
+
+	band := ds.Bands()[0]
+	err = band.IO(godal.IOWrite,
+		0, 0,
+		data,
+		width, height,
+		godal.PixelSpacing(2),
+		godal.LineSpacing(width*2))
+	if err != nil {
+		log.Error("Failed to write data to band:", err)
+		return err
+	}
+
+	log.Info("Saved pan image to ", tiffFile)
+
+	return nil
+}
+
+func (r *Registrator) SaveRegisteredMssToGDALGTiff(tiffFile string) error {
+	width := r.MssWidth
+	height := r.MssHeight
+
+	// 创建合并后的图像（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, tiffFile)
+}
+
+func (r *Registrator) SavePansharpenedToGDALGTiff(tiffFile string) error {
+	ihsImage := PansharpenIHS(r.PanImage, r.rgbirImage)
+	return SaveBGRToGDALGTiff(ihsImage, 3, 1.25, tiffFile)
+}
+
+func SaveBGRToGDALGTiff(bgr gocv.Mat, bands int, resolution float64, tiffFile string) error {
+	log.Println("Saving BGR to TIFF file:", tiffFile)
+
+	width := bgr.Cols()
+	height := bgr.Rows()
+
+	// 创建一个二维切片来存储图像数据
+	data := make([][]uint16, bands)
+	for i := range data {
+		data[i] = make([]uint16, width*height)
+	}
+
+	// 从gocv.Mat中提取数据
+	for y := 0; y < height; y++ {
+		for x := 0; x < width; x++ {
+			for b := 0; b < bands; b++ {
+				data[b][y*width+x] = uint16(bgr.GetShortAt(y, x*bands+b))
+			}
+		}
+	}
+
+	ds, err := godal.Create(godal.GTiff,
+		tiffFile,
+		bands,
+		godal.UInt16,
+		width, height)
+	if err != nil {
+		log.Error("Error creating TIFF file: ", err)
+		return err
+	}
+	defer ds.Close()
+
+	setGeoTransform(ds, 0, 0, float64(width), float64(height), resolution)
+
+	for b := 0; b < bands; b++ {
+		band := ds.Bands()[b]
+		err := band.IO(godal.IOWrite,
+			0, 0,
+			data[b],
+			width, height,
+			godal.PixelSpacing(2),
+			godal.LineSpacing(width*2))
+		if err != nil {
+			log.Error("Failed to write data to band:", err)
+			return err
+		}
+	}
+
+	log.Info("Saved registered mss to ", tiffFile)
+
+	return nil
+}
+
+func (r *Registrator) BytesToRaw(mssData []byte, filePath string) error {
+	f, err := os.OpenFile(filePath, os.O_CREATE|os.O_TRUNC|os.O_WRONLY, 0777)
+	if err != nil {
+		return err
+	}
+	w := bufio.NewWriter(f)
+	w.Write(mssData)
+
+	return nil
+}

phase_correlation.go

@@ -37,7 +37,7 @@ func (r *Registrator) calculateBlockPhaseShift(panBlock, mssBlock gocv.Mat) (goc
 
 	dx := shift.X
 	dy := shift.Y
-	log.Printf("Block shift: dx = %f, dy = %f. response = %f \n", dx, dy, response)
+	log.Debugf("Block shift: dx = %f, dy = %f. response = %f \n", dx, dy, response)
 
 	return shift, response
 }