diff --git a/cmd/proc.go b/cmd/proc.go index 0938f45..e8bd738 100644 --- a/cmd/proc.go +++ b/cmd/proc.go @@ -74,7 +74,6 @@ var procCmd = &cobra.Command{ if err := reg.DoPhaseCorrelation(); err != nil { logrus.Fatal(err) } - reg.DoCoRegistration() if params.SaveRegisteredMssRaw { registerdMSSRAW := filepath.Join( @@ -94,12 +93,6 @@ var procCmd = &cobra.Command{ producer.CleanScenes(panScenes) producer.CleanScenes(mssScenes) - // err := reg.SubScenesLessMem() - // if err != nil { - // logrus.Error(err) - // return - // } - runtime.GC() if saveStrip { diff --git a/config/config.yaml b/config/config.yaml index d220b8a..d640fe6 100644 --- a/config/config.yaml +++ b/config/config.yaml @@ -1,4 +1,5 @@ -log_level: 5 +log: + log_level: 5 coregistration: mss_bands: 4 pixel_bytes: 2 @@ -11,8 +12,8 @@ coregistration: fus_band_order: "RGB" radiation: - pan_remove_hf_noise: true - mss_remove_hf_noise: true + pan_remove_hf_noise: false + mss_remove_hf_noise: false hf_radius_ratio: 0.49 min_hist_level: 0.3 max_hist_level: 0.6 diff --git a/pkg/producer/edge_phase_correlation.go b/pkg/producer/edge_phase_correlation.go index cba08b0..34e4c1a 100644 --- a/pkg/producer/edge_phase_correlation.go +++ b/pkg/producer/edge_phase_correlation.go @@ -12,11 +12,14 @@ import ( // by performing phase correlation on detected // edges instead of the raw image -func FindEdges(img gocv.Mat) gocv.Mat { - img0 := gocv.IMRead("/Users/lan/workspace/sjy01/image-proc/data/052022/010/PAN/SJY01_PAN_20240520_115428_052022_103_010_L1A.tiff", gocv.IMReadUnchanged) +func FindEdges(img0 gocv.Mat) gocv.Mat { fmt.Println(img0.Cols(), img0.Rows(), img0.Type().String()) - img0.ConvertTo(&img0, gocv.MatTypeCV8U) - dst := gocv.NewMat() - gocv.Canny(img0, &dst, 100, 200) - return dst + 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) + return dstEdge } diff --git a/pkg/producer/image_registration.go b/pkg/producer/image_registration.go index e7e2c6c..7bfefe2 100644 --- a/pkg/producer/image_registration.go +++ b/pkg/producer/image_registration.go @@ -22,9 +22,9 @@ const ( PixelBytes = 2 PanWidth = payload.PAN_PIXEL_WIDTH // 像素宽度 MssWidth = payload.MSS_PIXEL_WIDTH - BlockNH = 4 - BlockNW = 8 - OverlappedBlockLines = 3000 // 重叠块的行数 + BlockNH = 8 + BlockNW = 4 + OverlappedBlockLines = 1000 // 重叠块的行数 DownSampled ResampleMethod = "down_sample_pan" UpSampled ResampleMethod = "up_sample_mss" PanResolution = 1.3 // mm/pixel @@ -170,12 +170,31 @@ func (r *Registrator) CalcDownPhaseCorrelation() error { blockHeight := r.MssHeight / BlockNH blockWidth := r.MssWidth / BlockNW - return r.calcPhaseCorrelation(downsampledPanImage, r.MssImages, r.MssHeight, r.MssWidth, blockHeight, blockWidth) + // 在 MSS 4 个波段上进行配准 + err := r.doMSSPhaseCorrelation(r.MssImages[0], + []gocv.Mat{r.MssImages[0], r.MssImages[1], r.MssImages[2], r.MssImages[3]}, + r.MssHeight, r.MssWidth, blockHeight, blockWidth) + if err != nil { + return err + } + r.DoMSSCoRegistration() + // 基于 PAN 图像进行配准 + r.phaseShifts[0] = make([]PhaseShiftM, 0) + r.deltaXCoeffs[0] = make([]float64, 0) + r.deltaYCoeffs[0] = make([]float64, 0) + err = r.doMSSPhaseCorrelation(downsampledPanImage, + []gocv.Mat{r.registeredMssImages[0]}, + r.MssHeight, r.MssWidth, blockHeight, blockWidth) + if err != nil { + return err + } + return r.DoPANCoRegistration() } // 将MSS升采样采样后计算相位相关的偏移量 func (r *Registrator) CalcUpPhaseCorrelation() error { + log.Fatal("unsuppotted up-resample method") // 确保 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") @@ -183,11 +202,19 @@ func (r *Registrator) CalcUpPhaseCorrelation() error { return err } - // 将PAN将采样作为轮廓匹配基准图像 - var upsampledMssImages [MssBands]gocv.Mat + // 在 MSS 4 个波段上进行配准 + err := r.doMSSPhaseCorrelation(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() + + upsampledMssImages := make([]gocv.Mat, MssBands) for i := 0; i < MssBands; i++ { upsampledMssImages[i] = gocv.NewMat() - gocv.Resize(r.MssImages[i], &upsampledMssImages[i], + gocv.Resize(r.registeredMssImages[i], &upsampledMssImages[i], image.Point{X: r.PanWidth, Y: r.PanHeight}, 0, 0, gocv.InterpolationCubic) } @@ -199,15 +226,23 @@ 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) + // 基于 PAN 图像进行配准 + err = r.doMSSPhaseCorrelation(r.PanImage, + []gocv.Mat{upsampledMssImages[0]}, + r.MssHeight, r.MssWidth, blockHeight, blockWidth) + return err } -func (r *Registrator) calcPhaseCorrelation(panImage gocv.Mat, - mssImages [4]gocv.Mat, +func (r *Registrator) doMSSPhaseCorrelation(base gocv.Mat, + mssImages []gocv.Mat, height, width, blockHeight, blockWidth int) error { var wg sync.WaitGroup + for band := 0; band < len(mssImages); band++ { + r.phaseShifts[band] = make([]PhaseShiftM, 0) + } + for bh := 0; bh < BlockNH; bh++ { for bw := 0; bw < BlockNW; bw++ { wg.Add(1) @@ -238,8 +273,8 @@ func (r *Registrator) calcPhaseCorrelation(panImage gocv.Mat, x1, y1, ) - panBlock := panImage.Region(rect) - for band := 0; band < MssBands; band++ { + panBlock := base.Region(rect) + for band := 0; band < len(mssImages); band++ { log.Debug("processing band:", band+1, ",block:", bh, rect) mssBlock := mssImages[band].Region(rect) @@ -262,7 +297,7 @@ func (r *Registrator) calcPhaseCorrelation(panImage gocv.Mat, wg.Wait() - for i := 0; i < MssBands; i++ { + for i := 0; i < len(mssImages); i++ { for _, shift := range r.phaseShifts[i] { if shift.response > 0.4 || shift.dx > 8 || shift.dy > 8 { log.Debugf("Band %d, block %d, dx=%f, dy=%f, response=%f", @@ -271,7 +306,7 @@ func (r *Registrator) calcPhaseCorrelation(panImage gocv.Mat, } } - return r.calcDeltaCoeffs() + return r.calcMSSDeltaCoeffs(len(mssImages)) } func (r *Registrator) Clean() { @@ -287,9 +322,9 @@ func (r *Registrator) Clean() { r.rgbirImage.Close() } -func (r *Registrator) calcDeltaCoeffs() error { +func (r *Registrator) calcMSSDeltaCoeffs(bands int) error { // 计算每个通道的delta多项式拟合系数 - for i := 0; i < MssBands; i++ { + for i := 0; i < bands; i++ { var cx []float64 var dx []float64 var dy []float64 @@ -330,7 +365,7 @@ func (r *Registrator) calcDeltaCoeffs() error { } } - for i := 0; i < MssBands; i++ { + for i := 0; i < bands; i++ { if len(r.deltaXCoeffs[i]) < 2 || len(r.deltaYCoeffs[i]) < 3 { continue } @@ -344,10 +379,10 @@ func (r *Registrator) calcDeltaCoeffs() error { return nil } -func (r *Registrator) DoCoRegistration() error { +func (r *Registrator) DoMSSCoRegistration() error { for band := 0; band < MssBands; band++ { if len(r.deltaXCoeffs[band]) < 2 || len(r.deltaYCoeffs[band]) < 3 { - log.Error("delta coefficients not calculated, skip co-registration") + log.Errorf("delta coefficients not calculated, skip co-registration %d", band) r.registeredMssImages[band] = r.MssImages[band].Clone() continue } @@ -367,12 +402,6 @@ func (r *Registrator) DoCoRegistration() error { 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)) } @@ -392,8 +421,58 @@ func (r *Registrator) DoCoRegistration() error { } // 裁掉末尾的的 MSS 480 行 和 PAN 的 480*4 行 - r.PanHeight -= 480 * 4 - r.MssHeight -= 480 + r.PanHeight -= 360 * 4 + r.MssHeight -= 360 + + return nil +} + +func (r *Registrator) DoPANCoRegistration() error { + if len(r.deltaXCoeffs[0]) < 2 || len(r.deltaYCoeffs[0]) < 3 { + log.Error("delta coefficients not calculated, skip co-registration") + return nil + } + + mapX := gocv.NewMatWithSize(r.MssHeight, r.MssWidth, gocv.MatTypeCV32FC1) + defer mapX.Close() + mapY := gocv.NewMatWithSize(r.MssHeight, r.MssWidth, gocv.MatTypeCV32FC1) + defer mapY.Close() + + for y := 0; y < r.MssHeight; y++ { + for x := 0; x < r.MssWidth; x++ { + var dx, dy float64 + if r.resampleMethod == UpSampled { + xx := float64(x * MssBands) + yy := float64(y * MssBands) + dx = (r.deltaXCoeffs[0][1]*float64(xx) + r.deltaXCoeffs[0][0] + xx) / MssBands + dy = (r.deltaYCoeffs[0][2]*float64(xx*xx) + r.deltaYCoeffs[0][1]*float64(xx) + r.deltaYCoeffs[0][0] + yy) / MssBands + } else { + dx = r.deltaXCoeffs[0][1]*float64(x) + r.deltaXCoeffs[0][0] + float64(x) + dy = r.deltaYCoeffs[0][2]*float64(x*x) + r.deltaYCoeffs[0][1]*float64(x) + r.deltaYCoeffs[0][0] + float64(y) + } + + mapX.SetFloatAt(y, x, float32(dx)) + mapY.SetFloatAt(y, x, float32(dy)) + } + } + + log.Println("co-registration for MSS (Align with PAN)") + for i := 0; i < MssBands; i++ { + log.Infof("r.registeredMssImages[%d]: %v", i, r.registeredMssImages[i].Size()) + registeredMSS := gocv.NewMatWithSize(r.MssHeight, r.MssWidth, gocv.MatTypeCV16UC1) + gocv.Remap(r.registeredMssImages[i], + ®isteredMSS, + &mapX, &mapY, + gocv.InterpolationCubic, + gocv.BorderConstant, + color.RGBA{0, 0, 0, 0}) + + 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/output.go b/pkg/producer/output.go index 73a5b5f..0bf1fb6 100644 --- a/pkg/producer/output.go +++ b/pkg/producer/output.go @@ -49,6 +49,30 @@ func (r *Registrator) BytesToRaw(mssData []byte, filePath string) error { } func (r *Registrator) SaveRegisteredMssToRaw(raw string) error { + return saveRegisteredMssToRaw(raw, r.registeredMssImages) +} + +func (r *Registrator) Report() error { + return WriteReport(&r.report, r.Params.ReportFile) +} + +func (r *Registrator) rpcKeywordInTif() { + // GDAL库对应的RPC关键词 + // keys := []string{ + // "ERR_BIAS", "ERR_RAND", + // "LINE_OFF", "SAMP_OFF", + // "LAT_OFF", "LONG_OFF", "HEIGHT_OFF", + // "LINE_SCALE", "SAMP_SCALE", + // "LAT_SCALE", "LONG_SCALE", "HEIGHT_SCALE", + // "LINE_NUM_COEFF", "LINE_DEN_COEFF", + // "SAMP_NUM_COEFF", "SAMP_DEN_COEFF", + // } + + // values := map[string]string{} + +} + +func saveRegisteredMssToRaw(raw string, mssImages [4]gocv.Mat) error { f, err := os.OpenFile(raw, os.O_CREATE|os.O_TRUNC|os.O_WRONLY, 0777) if err != nil { return err @@ -56,11 +80,11 @@ func (r *Registrator) SaveRegisteredMssToRaw(raw string) error { var mssData [4][]byte for i := 0; i < MssBands; i++ { - mssData[i] = r.registeredMssImages[i].ToBytes() + mssData[i] = mssImages[i].ToBytes() } - width := r.registeredMssImages[0].Cols() * PixelBytes - height := r.registeredMssImages[0].Rows() + width := mssImages[0].Cols() * PixelBytes + height := mssImages[0].Rows() log.Println("Writing registered MSS to RAW file...", len(mssData[0])*4) log.Println("width:", width) log.Println("height:", height) @@ -85,23 +109,3 @@ func (r *Registrator) SaveRegisteredMssToRaw(raw string) error { return nil } - -func (r *Registrator) Report() error { - return WriteReport(&r.report, r.Params.ReportFile) -} - -func (r *Registrator) rpcKeywordInTif() { - // GDAL库对应的RPC关键词 - // keys := []string{ - // "ERR_BIAS", "ERR_RAND", - // "LINE_OFF", "SAMP_OFF", - // "LAT_OFF", "LONG_OFF", "HEIGHT_OFF", - // "LINE_SCALE", "SAMP_SCALE", - // "LAT_SCALE", "LONG_SCALE", "HEIGHT_SCALE", - // "LINE_NUM_COEFF", "LINE_DEN_COEFF", - // "SAMP_NUM_COEFF", "SAMP_DEN_COEFF", - // } - - // values := map[string]string{} - -}