diff --git a/cmd/proc.go b/cmd/proc.go index 2dba148..4259ace 100644 --- a/cmd/proc.go +++ b/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) diff --git a/cmd/test.go b/cmd/test.go index b3d448c..0eba978 100644 --- a/cmd/test.go +++ b/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++ { - mv[i] = producer.CV_ECCAlign(mv[0], mv[i]) + index := len(mv) - 1 + 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}, diff --git a/pkg/producer/ecc_align.go b/pkg/producer/ecc_align.go index bba9e9e..287ef47 100644 --- a/pkg/producer/ecc_align.go +++ b/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) - templateImage.ConvertTo(&templateImage, gocv.MatTypeCV16U) - inputImage.ConvertTo(&inputImage, gocv.MatTypeCV16U) + aligned.ConvertTo(&aligned, retMatType) + templateImage.ConvertTo(&templateImage, retMatType) + inputImage.ConvertTo(&inputImage, retMatType) return aligned } diff --git a/pkg/producer/edge_phase_correlation.go b/pkg/producer/edge_phase_correlation.go index b5cbcb2..fd80acc 100644 --- a/pkg/producer/edge_phase_correlation.go +++ b/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) - sobelX.ConvertTo(&sobelX, gocv.MatTypeCV16U) + gocv.Sobel(img0, &sobelX, gocv.MatTypeCV32FC1, 1, 0, 5, 1.5, 0, gocv.BorderDefault) + sobelX.ConvertTo(&sobelX, matType) // y 方向 sobelY := gocv.NewMat() - gocv.Sobel(img0, &sobelY, gocv.MatTypeCV32F, 0, 1, 5, 1.5, 0, gocv.BorderIsolated) - sobelY.ConvertTo(&sobelY, gocv.MatTypeCV16U) + gocv.Sobel(img0, &sobelY, gocv.MatTypeCV32FC1, 0, 1, 5, 1.5, 0, gocv.BorderIsolated) + sobelY.ConvertTo(&sobelY, matType) // 合并 sobelXY := gocv.NewMat() - gocv.Sobel(img0, &sobelXY, gocv.MatTypeCV32F, 1, 1, 5, 1.5, 0, gocv.BorderDefault) - sobelXY.ConvertTo(&sobelXY, gocv.MatTypeCV16U) + gocv.Sobel(img0, &sobelXY, gocv.MatTypeCV32FC1, 1, 1, 5, 1.5, 0, gocv.BorderDefault) + sobelXY.ConvertTo(&sobelXY, matType) return sobelY } diff --git a/pkg/producer/image_registration.go b/pkg/producer/image_registration.go index 5693038..5cb0802 100644 --- a/pkg/producer/image_registration.go +++ b/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 - blockWidth := r.MssWidth / BlockNW + blockHeight := r.MssHeight / r.blockH + 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, - []gocv.Mat{r.registeredMssImages[0]}, - r.MssHeight, r.MssWidth, blockHeight, blockWidth) - if err != nil { - return err + // err = r.doPhaseCorrelation(downsampledPanImage, + // []gocv.Mat{r.registeredMssImages[0]}, + // r.MssHeight, r.MssWidth, blockHeight, blockWidth) + // if err != nil { + // 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], - []gocv.Mat{r.MssImages[0], r.MssImages[1], r.MssImages[2], r.MssImages[3]}, - r.MssHeight, r.MssWidth, r.MssHeight/BlockNH, r.MssWidth/BlockNW) - if err != nil { - return err - } - r.DoMSSCoRegistration(false) + // err := r.doPhaseCorrelation(r.MssImages[0], + // []gocv.Mat{r.MssImages[0], r.MssImages[1], r.MssImages[2], r.MssImages[3]}, + // r.MssHeight, r.MssWidth, r.MssHeight/r.blockH, r.MssWidth/r.blockW) + // if err != nil { + // return err + // } + // 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 - blockWidth := r.PanWidth / BlockNW + blockHeight := r.PanHeight / r.blockH + blockWidth := r.PanWidth / r.blockW log.Infof("blockHeight=%d, blockWidth=%d", blockHeight, blockWidth) // 基于 PAN 图像进行配准 - r.doPhaseCorrelation(r.PanImage, - []gocv.Mat{upsampledMssImages[0]}, - r.MssHeight, r.MssWidth, blockHeight, blockWidth) - return r.DoPANCoRegistration() + r.doPhaseCorrelation(r.PanImage, upsampledMssImages, + r.PanHeight, r.PanWidth, blockHeight, blockWidth) + r.fileterPhaseShift([]float64{64, 64, 64, 64}, true) + 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 bw := 0; bw < BlockNW; bw++ { + for bh := 0; bh < r.blockH; bh++ { + 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 } diff --git a/pkg/producer/phase_correlation.go b/pkg/producer/phase_correlation.go index a078351..13a0a0e 100644 --- a/pkg/producer/phase_correlation.go +++ b/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 +} diff --git a/pkg/producer/scene_register.go b/pkg/producer/scene_register.go new file mode 100644 index 0000000..38c198d --- /dev/null +++ b/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 +}