分景后配准

cmd/proc.go

@@ -7,6 +7,7 @@ import (
 	"strings"
 
 	"github.com/airbusgeo/godal"
+	"github.com/duke-git/lancet/v2/mathutil"
 	"github.com/sirupsen/logrus"
 	"github.com/spf13/cobra"
 	"starwiz.cn/sjy01/image-proc/pkg/calculator"
@@ -78,7 +79,7 @@ var procCmd = &cobra.Command{
 			processor.DoMomentMatching()
 		}
 
-		if err := processor.DoPhaseCorrelation(); err != nil {
+		if err := processor.DoPhaseCorrelation(true); err != nil {
 			logrus.Fatal(err)
 		}
 
@@ -96,6 +97,11 @@ var procCmd = &cobra.Command{
 			logrus.Error(err)
 		}
 
+		scenesCnt := mathutil.Min(len(panScenes), len(mssScenes))
+		for i := 0; i < scenesCnt; i++ {
+			processor.RegisterScenes(panScenes[i], mssScenes[i])
+		}
+
 		processor.OutputL1A(panScenes, mssScenes)
 		producer.CleanScenes(panScenes)
 		producer.CleanScenes(mssScenes)

cmd/test.go

@@ -31,12 +31,13 @@ var testCmd = &cobra.Command{
 		// 	out := producer.CV_Sobel(mv[i])
 		// 	utils.SavePanToGDALGTiff(out, 0, 0, "data/temp/out_"+strconv.Itoa(i)+".tif", 1.3)
 		// }
-		for i := 1; i < len(mv); i++ {
+		index := len(mv) - 1
-			mv[i] = producer.CV_ECCAlign(mv[0], mv[i])
+		for i := index - 1; i >= 0; i-- {
+			mv[i] = producer.CV_ECCAlign(mv[index], mv[i])
 		}
 
 		out := gocv.NewMat()
-		mv[0].ConvertTo(&mv[0], gocv.MatTypeCV16U)
+		mv[index].ConvertTo(&mv[index], gocv.MatTypeCV16U)
 		gocv.Merge(mv, &out)
 		utils.SaveBGRToGDALGTiff(out, 4, 0, 0, 5.2,
 			[]godal.ColorInterp{godal.CIBlue, godal.CIGreen, godal.CIRed, godal.CIUndefined},

pkg/producer/ecc_align.go

@@ -10,6 +10,7 @@ import (
 // align the image via Enhanced Correlation Coefficient (ECC) algorithm
 
 func CV_ECCAlign(templateImage, inputImage gocv.Mat) gocv.Mat {
+	retMatType := templateImage.Type()
 	warpMode := gocv.MotionHomography
 	var warpMatrix gocv.Mat
 	switch warpMode {
@@ -27,7 +28,7 @@ func CV_ECCAlign(templateImage, inputImage gocv.Mat) gocv.Mat {
 	var mask = gocv.NewMat()
 	templateImage.ConvertTo(&templateImage, gocv.MatTypeCV32FC1)
 	inputImage.ConvertTo(&inputImage, gocv.MatTypeCV32FC1)
-	log.Info("cv::findTransformECC...")
+
 	retval := gocv.FindTransformECC(templateImage, inputImage, &warpMatrix, warpMode, criteria, mask, 5)
 	log.Info("cv::findTransformECC ret:", retval)
 
@@ -39,9 +40,9 @@ func CV_ECCAlign(templateImage, inputImage gocv.Mat) gocv.Mat {
 		gocv.WarpAffine(inputImage, &aligned, warpMatrix, sz)
 	}
 
-	aligned.ConvertTo(&aligned, gocv.MatTypeCV16U)
+	aligned.ConvertTo(&aligned, retMatType)
-	templateImage.ConvertTo(&templateImage, gocv.MatTypeCV16U)
+	templateImage.ConvertTo(&templateImage, retMatType)
-	inputImage.ConvertTo(&inputImage, gocv.MatTypeCV16U)
+	inputImage.ConvertTo(&inputImage, retMatType)
 	return aligned
 }
 

pkg/producer/edge_phase_correlation.go

@@ -1,8 +1,6 @@
 package producer
 
 import (
-	"fmt"
-
 	"gocv.io/x/gocv"
 )
 
@@ -14,30 +12,30 @@ import (
 
 // SKIP: 多光谱各个波段的边缘检测结果不佳
 func CV_Canny(img0 gocv.Mat) gocv.Mat {
-	fmt.Println(img0.Cols(), img0.Rows(), img0.Type().String())
 	dst8 := gocv.NewMatWithSize(img0.Rows(), img0.Cols(), gocv.MatTypeCV8U)
 	defer dst8.Close()
 	gocv.Normalize(img0, &dst8, 0, 255, gocv.NormMinMax)
 	dst8.ConvertTo(&dst8, gocv.MatTypeCV8U)
 	dstEdge := gocv.NewMat()
 	gocv.Canny(dst8, &dstEdge, 10, 100)
-	dstEdge.ConvertTo(&dstEdge, gocv.MatTypeCV16U)
+	dstEdge.ConvertTo(&dstEdge, img0.Type())
 	return dstEdge
 }
 
 func CV_Sobel(img0 gocv.Mat) gocv.Mat {
+	matType := img0.Type()
 	// x 方向
 	sobelX := gocv.NewMat()
-	gocv.Sobel(img0, &sobelX, gocv.MatTypeCV32F, 1, 0, 5, 1.5, 0, gocv.BorderDefault)
+	gocv.Sobel(img0, &sobelX, gocv.MatTypeCV32FC1, 1, 0, 5, 1.5, 0, gocv.BorderDefault)
-	sobelX.ConvertTo(&sobelX, gocv.MatTypeCV16U)
+	sobelX.ConvertTo(&sobelX, matType)
 	// y 方向
 	sobelY := gocv.NewMat()
-	gocv.Sobel(img0, &sobelY, gocv.MatTypeCV32F, 0, 1, 5, 1.5, 0, gocv.BorderIsolated)
+	gocv.Sobel(img0, &sobelY, gocv.MatTypeCV32FC1, 0, 1, 5, 1.5, 0, gocv.BorderIsolated)
-	sobelY.ConvertTo(&sobelY, gocv.MatTypeCV16U)
+	sobelY.ConvertTo(&sobelY, matType)
 	// 合并
 	sobelXY := gocv.NewMat()
-	gocv.Sobel(img0, &sobelXY, gocv.MatTypeCV32F, 1, 1, 5, 1.5, 0, gocv.BorderDefault)
+	gocv.Sobel(img0, &sobelXY, gocv.MatTypeCV32FC1, 1, 1, 5, 1.5, 0, gocv.BorderDefault)
-	sobelXY.ConvertTo(&sobelXY, gocv.MatTypeCV16U)
+	sobelXY.ConvertTo(&sobelXY, matType)
 
 	return sobelY
 }

pkg/producer/image_registration.go

@@ -22,8 +22,6 @@ const (
 	PixelBytes                          = 2
 	PanWidth                            = payload.PAN_PIXEL_WIDTH // 像素宽度
 	MssWidth                            = payload.MSS_PIXEL_WIDTH
-	BlockNH                             = 8
-	BlockNW                             = 4
 	OverlappedBlockLines                = 1000 // 重叠块的行数
 	DownSampled          ResampleMethod = "down_sample_pan"
 	UpSampled            ResampleMethod = "up_sample_mss"
@@ -33,6 +31,7 @@ const (
 	ReferenceShiftYB2                   = 200.0
 	ReferenceShiftYB3                   = 300.0
 	ReferenceShiftYB4                   = 400.0
+	MaxShiftY                           = 480
 )
 
 type ResampleMethod string
@@ -47,6 +46,8 @@ type ImgProc struct {
 	MssImages [4]gocv.Mat
 	MssHeight int
 	MssWidth  int
+	blockH    int
+	blockW    int
 
 	shiftMutex   sync.Mutex
 	phaseShifts  [4][]PhaseShiftM
@@ -72,6 +73,8 @@ type ImgProc struct {
 func NewImgProc(rsmethod ResampleMethod) *ImgProc {
 	var r ImgProc
 	r.resampleMethod = rsmethod
+	r.blockH = 8
+	r.blockW = 8
 
 	return &r
 }
@@ -141,7 +144,14 @@ func (r *ImgProc) LoadMssRaw() error {
 	return nil
 }
 
-func (r *ImgProc) DoPhaseCorrelation() error {
+func (r *ImgProc) DoPhaseCorrelation(skip bool) error {
+	if skip {
+		for i := 0; i < MssBands; i++ {
+			r.registeredMssImages[i] = r.MssImages[i].Clone()
+		}
+		return nil
+	}
+
 	switch r.resampleMethod {
 	case UpSampled:
 		return r.CalcUpPhaseCorrelation()
@@ -167,11 +177,13 @@ func (r *ImgProc) CalcDownPhaseCorrelation() error {
 	log.Println("down sampled PAN images size:", downsampledPanImage.Size())
 
 	// 分块高度
-	blockHeight := r.MssHeight / BlockNH
+	blockHeight := r.MssHeight / r.blockH
-	blockWidth := r.MssWidth / BlockNW
+	blockWidth := r.MssWidth / r.blockW
+
+	log.Infof("blockHeight=%d, blockWidth=%d", blockHeight, blockWidth)
 
 	// 在 MSS 4 个波段上进行配准
-	err := r.doPhaseCorrelation(r.MssImages[0],
+	err := r.doPhaseCorrelation(downsampledPanImage,
 		[]gocv.Mat{r.MssImages[0], r.MssImages[1], r.MssImages[2], r.MssImages[3]},
 		r.MssHeight, r.MssWidth, blockHeight, blockWidth)
 	if err != nil {
@@ -195,22 +207,30 @@ func (r *ImgProc) CalcDownPhaseCorrelation() error {
 	// r.DoMSSCoRegistration(true)
 
 	// 基于 PAN 图像进行配准
-	err = r.doPhaseCorrelation(downsampledPanImage,
+	// err = r.doPhaseCorrelation(downsampledPanImage,
-		[]gocv.Mat{r.registeredMssImages[0]},
+	// 	[]gocv.Mat{r.registeredMssImages[0]},
-		r.MssHeight, r.MssWidth, blockHeight, blockWidth)
+	// 	r.MssHeight, r.MssWidth, blockHeight, blockWidth)
-	if err != nil {
+	// if err != nil {
-		return err
+	// 	return err
+	// }
+	// r.fileterPhaseShift([]float64{30.0}, true)
+	// r.calcMSSDeltaCoeffs(1)
+	// r.DoPANCoRegistration()
+
+	// 裁掉位移部分
+	r.MssHeight -= MaxShiftY
+	r.PanHeight -= MaxShiftY * 4
+	r.PanImage = r.PanImage.Region(image.Rect(0, 0, r.PanWidth, r.PanHeight))
+	for i := 0; i < MssBands; i++ {
+		r.registeredMssImages[i] = r.registeredMssImages[i].Region(image.Rect(0, 0, r.MssWidth, r.MssHeight))
 	}
-	r.fileterPhaseShift([]float64{30.0}, true)
-	r.calcMSSDeltaCoeffs(1)
-	r.DoPANCoRegistration()
 
 	return nil
 }
 
 // 将MSS升采样采样后计算相位相关的偏移量
 func (r *ImgProc) CalcUpPhaseCorrelation() error {
-	log.Fatal("unsuppotted up-resample method")
+	log.Fatal("not implemented yet")
 	// 确保 MSS 高度是 PAN 高度的 1/4
 	if r.MssHeight*4 != r.PanHeight {
 		err := fmt.Errorf("MSS height is not 1/4 of PAN height, invalid raw file")
@@ -219,34 +239,36 @@ func (r *ImgProc) CalcUpPhaseCorrelation() error {
 	}
 
 	// 在 MSS 4 个波段上进行配准
-	err := r.doPhaseCorrelation(r.MssImages[0],
+	// err := r.doPhaseCorrelation(r.MssImages[0],
-		[]gocv.Mat{r.MssImages[0], r.MssImages[1], r.MssImages[2], r.MssImages[3]},
+	// 	[]gocv.Mat{r.MssImages[0], r.MssImages[1], r.MssImages[2], r.MssImages[3]},
-		r.MssHeight, r.MssWidth, r.MssHeight/BlockNH, r.MssWidth/BlockNW)
+	// 	r.MssHeight, r.MssWidth, r.MssHeight/r.blockH, r.MssWidth/r.blockW)
-	if err != nil {
+	// if err != nil {
-		return err
+	// 	return err
-	}
+	// }
-	r.DoMSSCoRegistration(false)
+	// r.DoMSSCoRegistration(false)
 
 	upsampledMssImages := make([]gocv.Mat, MssBands)
 	for i := 0; i < MssBands; i++ {
 		upsampledMssImages[i] = gocv.NewMat()
-		gocv.Resize(r.registeredMssImages[i], &upsampledMssImages[i],
+		gocv.Resize(r.MssImages[i], &upsampledMssImages[i],
 			image.Point{X: r.PanWidth, Y: r.PanHeight}, 0, 0, gocv.InterpolationCubic)
 	}
 
 	fmt.Println("up sampled MSS images size:", upsampledMssImages[0].Size())
 
 	// 分块高度 - BlockNH, BlockNW % 4 == 0
-	blockHeight := r.PanHeight / BlockNH
+	blockHeight := r.PanHeight / r.blockH
-	blockWidth := r.PanWidth / BlockNW
+	blockWidth := r.PanWidth / r.blockW
 
 	log.Infof("blockHeight=%d, blockWidth=%d", blockHeight, blockWidth)
 
 	// 基于 PAN 图像进行配准
-	r.doPhaseCorrelation(r.PanImage,
+	r.doPhaseCorrelation(r.PanImage, upsampledMssImages,
-		[]gocv.Mat{upsampledMssImages[0]},
+		r.PanHeight, r.PanWidth, blockHeight, blockWidth)
-		r.MssHeight, r.MssWidth, blockHeight, blockWidth)
+	r.fileterPhaseShift([]float64{64, 64, 64, 64}, true)
-	return r.DoPANCoRegistration()
+	r.calcMSSDeltaCoeffs(4)
+	r.DoMSSCoRegistration(false)
+	return nil
 }
 
 func (r *ImgProc) doPhaseCorrelation(base gocv.Mat,
@@ -261,8 +283,8 @@ func (r *ImgProc) doPhaseCorrelation(base gocv.Mat,
 		r.deltaYCoeffs[band] = make([]float64, 0)
 	}
 
-	for bh := 0; bh < BlockNH; bh++ {
+	for bh := 0; bh < r.blockH; bh++ {
-		for bw := 0; bw < BlockNW; bw++ {
+		for bw := 0; bw < r.blockW; bw++ {
 			wg.Add(1)
 			go func(bh, bw int) {
 				defer wg.Done()
@@ -449,10 +471,6 @@ func (r *ImgProc) DoMSSCoRegistration(byEdge bool) error {
 		mapY.Close()
 	}
 
-	// 裁掉末尾的的 MSS 480 行 和 PAN 的 480*4 行
-	r.PanHeight -= 360 * 4
-	r.MssHeight -= 360
-
 	return nil
 }
 
@@ -498,9 +516,6 @@ func (r *ImgProc) DoPANCoRegistration() error {
 		r.registeredMssImages[i].Close()
 		r.registeredMssImages[i] = registeredMSS
 	}
-	// 裁掉末尾的的 MSS 480 行 和 PAN 的 480*4 行
-	r.PanHeight -= 120 * 4
-	r.MssHeight -= 120
 
 	return nil
 }

pkg/producer/phase_correlation.go

@@ -3,8 +3,10 @@ package producer
 import (
 	"errors"
 	"image"
+	"math/cmplx"
 
 	"github.com/duke-git/lancet/v2/slice"
+	"github.com/mjibson/go-dsp/fft"
 	log "github.com/sirupsen/logrus"
 	"gocv.io/x/gocv"
 )
@@ -31,9 +33,9 @@ func CV_PhaseCorrelate(panBlock, mssBlock gocv.Mat) (gocv.Point2f, float64) {
 	mssBlock.ConvertTo(&mss, gocv.MatTypeCV32FC1)
 	defer mss.Close()
 
-	hann := gocv.NewMat()
+	hann := gocv.NewMatWithSize(pan.Rows(), pan.Cols(), pan.Type())
 	defer hann.Close()
-
+	gocv.CreateHanningWindow(&hann, image.Point{X: pan.Cols(), Y: pan.Rows()}, pan.Type())
 	shift, response := gocv.PhaseCorrelate(pan, mss, hann)
 
 	dx := shift.X
@@ -64,3 +66,77 @@ func (r *ImgProc) fileterPhaseShift(thredholds []float64, greaterThan bool) erro
 
 	return nil
 }
+
+func PhaseCorrelate(panBlock, mssBlock gocv.Mat) (gocv.Point2f, float64) {
+	// 计算傅里叶变换
+	panFreqDomain, _, _, err := toFreqDomain(panBlock)
+	if err != nil {
+		log.Error(err)
+	}
+
+	mssFreqDomain, width, height, err := toFreqDomain(mssBlock)
+	if err != nil {
+		log.Error(err)
+	}
+
+	// 计算共轭乘积并积累相位信息
+	crossPowerSpectrum := make([][]complex128, height)
+	for i := range crossPowerSpectrum {
+		crossPowerSpectrum[i] = make([]complex128, width)
+	}
+
+	for y := 0; y < height; y++ {
+		for x := 0; x < width; x++ {
+			if panFreqDomain[y][x] != 0 && mssFreqDomain[y][x] != 0 {
+				crossPowerSpectrum[y][x] = panFreqDomain[y][x] * cmplx.Conj(mssFreqDomain[y][x])
+			}
+		}
+	}
+
+	// 归一化
+	for y := 0; y < height; y++ {
+		for x := 0; x < width; x++ {
+			if crossPowerSpectrum[y][x] != 0 {
+				magnitude := cmplx.Abs(crossPowerSpectrum[y][x])
+				crossPowerSpectrum[y][x] = crossPowerSpectrum[y][x] / complex(magnitude, 0)
+			}
+		}
+	}
+
+	ifftResult := fft.IFFT2(crossPowerSpectrum)
+
+	// 查找最大值及其对应的平移参数
+	maxVal := 0.0
+	var dx, dy float64
+	for y := 0; y < height; y++ {
+		for x := 0; x < width; x++ {
+			val := real(ifftResult[y][x])
+			if val > maxVal {
+				maxVal = val
+				dx = float64(x)
+				dy = float64(y)
+			}
+		}
+	}
+	log.Debugf("Block shift: dx = %f, dy = %f. response = %f", dx, dy, 0.0)
+
+	return gocv.Point2f{X: float32(dx), Y: float32(dy)}, 0.0
+}
+
+func toFreqDomain(input gocv.Mat) ([][]complex128, int, int, error) {
+	height := input.Rows()
+	width := input.Cols()
+
+	data := make([][]complex128, height)
+	for y := 0; y < height; y++ {
+		data[y] = make([]complex128, width)
+		for x := 0; x < width; x++ {
+			grayColor := float64(uint16(input.GetShortAt(y, x)))
+			data[y][x] = complex(grayColor, 0)
+		}
+	}
+
+	freqDomain := fft.FFT2(data)
+
+	return freqDomain, width, height, nil
+}

pkg/producer/scene_register.go

@@ -0,0 +1,36 @@
+package producer
+
+import (
+	"github.com/sirupsen/logrus"
+	"gocv.io/x/gocv"
+)
+
+// 分景后再做配准
+func (r *ImgProc) RegisterScenes(pan, mss *Scene) error {
+	proc := NewImgProc(DownSampled)
+	proc.PanImage = pan.Mat[0]
+	proc.MssImages = [4]gocv.Mat{mss.Mat[0], mss.Mat[1], mss.Mat[2], mss.Mat[3]}
+	proc.MssHeight = mss.Height
+	proc.MssWidth = mss.Width
+	proc.PanHeight = pan.Height
+	proc.PanWidth = pan.Width
+
+	err := proc.DoPhaseCorrelation(false)
+	if err != nil {
+		logrus.Error("register scene failed: ", err)
+		return err
+	}
+
+	pan.Height = proc.PanHeight
+	pan.Width = proc.PanWidth
+	pan.Mat[0] = proc.PanImage
+
+	mss.Height = proc.MssHeight
+	mss.Width = proc.MssWidth
+	for i := 0; i < 4; i++ {
+		mss.Mat[i].Close()
+		mss.Mat[i] = proc.registeredMssImages[i]
+	}
+
+	return nil
+}