diff --git a/go.mod b/go.mod index a0e306f..75941c0 100644 --- a/go.mod +++ b/go.mod @@ -34,6 +34,7 @@ require ( github.com/mattn/go-isatty v0.0.17 // indirect github.com/mgutz/ansi v0.0.0-20200706080929-d51e80ef957d // indirect github.com/mitchellh/mapstructure v1.5.0 // indirect + github.com/mjibson/go-dsp v0.0.0-20180508042940-11479a337f12 // indirect github.com/onsi/ginkgo v1.16.5 // indirect github.com/onsi/gomega v1.33.1 // indirect github.com/pelletier/go-toml/v2 v2.1.0 // indirect diff --git a/go.sum b/go.sum index cd7075f..665fe0f 100644 --- a/go.sum +++ b/go.sum @@ -1131,6 +1131,8 @@ github.com/mitchellh/mapstructure v1.1.2/go.mod h1:FVVH3fgwuzCH5S8UJGiWEs2h04kUh github.com/mitchellh/mapstructure v1.4.1/go.mod h1:bFUtVrKA4DC2yAKiSyO/QUcy7e+RRV2QTWOzhPopBRo= github.com/mitchellh/mapstructure v1.5.0 h1:jeMsZIYE/09sWLaz43PL7Gy6RuMjD2eJVyuac5Z2hdY= github.com/mitchellh/mapstructure v1.5.0/go.mod h1:bFUtVrKA4DC2yAKiSyO/QUcy7e+RRV2QTWOzhPopBRo= +github.com/mjibson/go-dsp v0.0.0-20180508042940-11479a337f12 h1:dd7vnTDfjtwCETZDrRe+GPYNLA1jBtbZeyfyE8eZCyk= +github.com/mjibson/go-dsp v0.0.0-20180508042940-11479a337f12/go.mod h1:i/KKcxEWEO8Yyl11DYafRPKOPVYTrhxiTRigjtEEXZU= github.com/modern-go/concurrent v0.0.0-20180228061459-e0a39a4cb421/go.mod h1:6dJC0mAP4ikYIbvyc7fijjWJddQyLn8Ig3JB5CqoB9Q= github.com/modern-go/reflect2 v0.0.0-20180701023420-4b7aa43c6742/go.mod h1:bx2lNnkwVCuqBIxFjflWJWanXIb3RllmbCylyMrvgv0= github.com/modern-go/reflect2 v1.0.1/go.mod h1:bx2lNnkwVCuqBIxFjflWJWanXIb3RllmbCylyMrvgv0= diff --git a/pkg/producer/scenes.go b/pkg/producer/scenes.go index 0a52b7c..82e31ae 100644 --- a/pkg/producer/scenes.go +++ b/pkg/producer/scenes.go @@ -13,6 +13,7 @@ import ( "github.com/paulmach/orb/geojson" log "github.com/sirupsen/logrus" "gocv.io/x/gocv" + "starwiz.cn/sjy01/image-proc/pkg/rrc" ) type Scene struct { @@ -66,7 +67,9 @@ func (r *Registrator) SubScenes() (panScenes []*Scene, mssScenes []*Scene, err e scene.SceneId = fmt.Sprintf("%s_%03d", name, i+1) mat := r.PanImage.Region(image.Rect(0, i*hPAN, 9344, y1)) - scene.Mat = append(scene.Mat, mat) + matFiltered := rrc.HFNoiseFilter(mat, float64(mat.Cols())*0.45) + scene.Mat = append(scene.Mat, matFiltered) + mat.Close() panScenes = append(panScenes, scene) } @@ -94,7 +97,9 @@ func (r *Registrator) SubScenes() (panScenes []*Scene, mssScenes []*Scene, err e for band := 0; band < 4; band++ { mat := r.registeredMssImages[band].Region(image.Rect(0, i*hMSS, 2336, y1)) - scene.Mat = append(scene.Mat, mat) + matFiltered := rrc.HFNoiseFilter(mat, float64(mat.Cols())*0.45) + mat.Close() + scene.Mat = append(scene.Mat, matFiltered) } name := filepath.Base(r.Params.MssTiffFile) diff --git a/pkg/rrc/filter_test.go b/pkg/rrc/filter_test.go new file mode 100644 index 0000000..f52fd1a --- /dev/null +++ b/pkg/rrc/filter_test.go @@ -0,0 +1,43 @@ +package rrc + +import ( + "testing" + + "github.com/airbusgeo/godal" + "gocv.io/x/gocv" + "starwiz.cn/sjy01/image-proc/pkg/utils" +) + +func TestFilter(t *testing.T) { + godal.RegisterAll() + tif := "/Users/lan/workspace/sjy01/image-proc/data/052022-no-rrc/010/PAN/SJY01_PAN_20240520_115428_052022_103_010_L1A.tiff" + tif = "/Users/lan/workspace/sjy01/image-proc/data/051622/007/PAN/SJY01_PAN_20240516_101236_051622_096_007_L1A.tiff" + // tif = "/Users/lan/workspace/sjy01/image-proc/data/060622-mm/003/PAN/SJY01_PAN_20240606_012004_060622_065_003_L1A.tiff" + + ds, err := godal.Open(tif) + if err != nil { + t.Error(err) + } + bands := ds.Bands() + width := bands[0].Structure().SizeX + height := bands[0].Structure().SizeY + img := gocv.NewMatWithSize(height, width, gocv.MatTypeCV16UC1) + + data := make([]uint16, width*height) + err = bands[0].Read(0, 0, data, width, height) + if err != nil { + t.Error(err) + } + + // 将 16 位数据存储到图像的相应通道 + for y := 0; y < height; y++ { + for x := 0; x < width; x++ { + value := data[y*width+x] + img.SetShortAt(y, x, int16(value)) + } + } + + m0 := HFNoiseFilter(img, float64(img.Cols())*0.45) + defer m0.Close() + utils.SavePanToGDALGTiff(m0, 0, 0, "/Users/lan/workspace/sjy01/image-proc/data/fft/filtered.tiff", 1.3) +} diff --git a/pkg/rrc/hf_filter.go b/pkg/rrc/hf_filter.go new file mode 100644 index 0000000..0377ebf --- /dev/null +++ b/pkg/rrc/hf_filter.go @@ -0,0 +1,118 @@ +package rrc + +import ( + "image" + "math" + + "github.com/sirupsen/logrus" + "gocv.io/x/gocv" +) + +// filter high frequence noise +func HFNoiseFilter(img gocv.Mat, radius float64) gocv.Mat { + logrus.Println("applying hf noise filter...") + imgFloat := gocv.NewMat() + img.ConvertTo(&imgFloat, gocv.MatTypeCV32F) + + // rows, cols := img.Rows(), img.Cols() + // utils.SavePanToGDALGTiff(img, 0, 0, "/Users/lan/workspace/sjy01/image-proc/data/fft/original.tiff", 1.3) + + // 将图像转换为32位浮点型 + img32f := gocv.NewMat() + img.ConvertTo(&img32f, gocv.MatTypeCV32F) + defer img32f.Close() + + // 获取图像尺寸 + rows, cols := img32f.Rows(), img32f.Cols() + + // zeros := gocv.NewMatFromScalar(gocv.NewScalar(0, 0, 0, 0), gocv.MatTypeCV32F) + + // 执行 DFT + DFT := gocv.NewMat() + defer DFT.Close() + + planes := gocv.NewMatWithSize(rows, cols, gocv.MatTypeCV32F) + img32f.ConvertTo(&planes, gocv.MatTypeCV32F) + gocv.Merge([]gocv.Mat{planes, gocv.Zeros(rows, cols, gocv.MatTypeCV32F)}, &DFT) + gocv.DFT(DFT, &DFT, gocv.DftComplexOutput) + fftshift(&DFT) + // 将复数部分分离并计算振幅光谱 + pss := gocv.Split(DFT) + // SaveDFT(DFT, "/Users/lan/workspace/sjy01/image-proc/data/fft/dft.tiff") + + // 应用滤波器 + for y := 0; y < rows; y++ { + for x := 0; x < cols; x++ { + if math.Sqrt(math.Pow(float64(x-cols/2), 2)+math.Pow(float64(y-rows/2), 2)) > float64(radius) { + // if math.Abs(float64(x-cols/2)) < float64(radius) || math.Abs(float64(y-rows/2)) < float64(radius) { + pss[0].SetFloatAt(y, x, 0.0) + pss[1].SetFloatAt(y, x, 0.0) + } + } + } + gocv.Merge(pss, &DFT) + // SaveDFT(DFT, "/Users/lan/workspace/sjy01/image-proc/data/fft/dft_filtered.tiff") + fftshift(&DFT) + // 执行 IDFT + IDFT := gocv.NewMat() + gocv.DFT(DFT, &IDFT, gocv.DftInverse|gocv.DftRealOutput) + + // 归一化到0-65535 + gocv.Normalize(IDFT, &IDFT, 0, 8192, gocv.NormMinMax) + IDFT.ConvertTo(&IDFT, gocv.MatTypeCV16U) + + // 保存 IDFT 结果 + return IDFT +} + +// 创建低通滤波器 +func createLowPassFilter(filter *gocv.Mat, radius int) { + rows, cols := filter.Rows(), filter.Cols() + center := image.Pt(cols/2, rows/2) + for y := 0; y < rows; y++ { + for x := 0; x < cols; x++ { + dist := math.Sqrt(math.Pow(float64(x-center.X), 2) + math.Pow(float64(y-center.Y), 2)) + if dist <= float64(radius) { + filter.SetFloatAt(y, x, 1.0) + } else { + filter.SetFloatAt(y, x, 0.0) + } + } + } +} + +// fftshift 将频率直流分量移动到图像中心 +func fftshift(input *gocv.Mat) { + cx, cy := input.Cols()/2, input.Rows()/2 + + q0 := input.Region(image.Rect(0, 0, cx, cy)) // Top-Left + q1 := input.Region(image.Rect(cx, 0, input.Cols(), cy)) // Top-Right + q2 := input.Region(image.Rect(0, cy, cx, input.Rows())) // Bottom-Left + q3 := input.Region(image.Rect(cx, cy, input.Cols(), input.Rows())) // Bottom-Right + + tmp := gocv.NewMatWithSize(q0.Rows(), q0.Cols(), input.Type()) + q0.CopyTo(&tmp) + q3.CopyTo(&q0) + tmp.CopyTo(&q3) + + q1.CopyTo(&tmp) + q2.CopyTo(&q1) + tmp.CopyTo(&q2) +} + +func SaveDFT(DFT gocv.Mat, file string) { + pss := gocv.Split(DFT) + planes := gocv.NewMat() + gocv.Magnitude(pss[0], pss[1], &planes) + // 进行对数尺度变换,便于可视化 + ones := gocv.NewMatFromScalar(gocv.NewScalar(1, 0, 0, 0), gocv.MatTypeCV64F) + gocv.Add(planes, ones, &planes) + gocv.Log(planes, &planes) + + // 归一化 + gocv.Normalize(planes, &planes, 0, 1, gocv.NormMinMax) + + // 保存DFT结果 + gocv.IMWrite(file, planes) + +} diff --git a/pkg/rrc/hf_filter_fft.go b/pkg/rrc/hf_filter_fft.go new file mode 100644 index 0000000..723780f --- /dev/null +++ b/pkg/rrc/hf_filter_fft.go @@ -0,0 +1,56 @@ +package rrc + +import ( + "math" + + "github.com/mjibson/go-dsp/fft" + log "github.com/sirupsen/logrus" + "gocv.io/x/gocv" +) + +func HFNoiseFilterFFT(input gocv.Mat, radius int) gocv.Mat { + // 将灰度图像转换为浮点数数组 + 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) + } + } + + // 进行二维DFT + log.Println("Performing 2D DFT...") + dft := fft.FFT2(data) + // 应用低通滤波器 + log.Println("Applying low-frequency-pass filter...") + for y := 0; y < height; y++ { + for x := 0; x < width; x++ { + if math.Abs(float64(x-width/2)) < float64(radius) || math.Abs(float64(y-height/2)) < float64(radius) { + dft[y][x] = 0 + } + } + } + + // 进行二维逆DFT + log.Println("Performing 2D inverse DFT...") + idft := fft.IFFT2(dft) + + // 将浮点数数组转换回灰度图像 + for y := 0; y < height; y++ { + for x := 0; x < width; x++ { + realPart := real(idft[y][x]) + if realPart < 0 { + realPart = 0 + } else if realPart > 65535 { + realPart = 65535 + } + input.SetShortAt(y, x, int16(uint16(realPart))) + } + } + + return input +} diff --git a/pkg/rrc/mean_filter.go b/pkg/rrc/mean_filter.go new file mode 100644 index 0000000..c131e61 --- /dev/null +++ b/pkg/rrc/mean_filter.go @@ -0,0 +1,15 @@ +package rrc + +import ( + "image" + + log "github.com/sirupsen/logrus" + "gocv.io/x/gocv" +) + +func MeanFilter(img gocv.Mat, ksize int) gocv.Mat { + log.Println("Applying mean filter, kernel size:", ksize) + meanFilteredImg := gocv.NewMat() + gocv.Blur(img, &meanFilteredImg, image.Pt(ksize, ksize)) + return meanFilteredImg +} diff --git a/pkg/utils/tiff.go b/pkg/utils/tiff.go new file mode 100644 index 0000000..53b6a19 --- /dev/null +++ b/pkg/utils/tiff.go @@ -0,0 +1,95 @@ +package utils + +import ( + "time" + + "github.com/airbusgeo/godal" + log "github.com/sirupsen/logrus" + "starwiz.cn/sjy01/image-proc/pkg/calculator" + + "gocv.io/x/gocv" +) + +func SavePanToGDALGTiff(pan gocv.Mat, topLeftX, topLeftY float64, tiffFile string, resolution float64) error { + // log.Println("Saving PAN image to TIFF file:", tiffFile) + + width := pan.Cols() + height := pan.Rows() + + 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, topLeftX, topLeftY, resolution) + 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] + band.SetColorInterp(godal.CIGray) + 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 setGeoTransform(ds *godal.Dataset, topLeftLng, topLeftLat, resolution float64) { + // 转换分辨率(米到度) + resLat := resolution / calculator.MetersPerDegreeLatitude + resLng := calculator.MetersToDegreesLongitude(resolution, topLeftLat) + // 设置地理变换(假设左上角坐标为(0,0),PAN每个像素分辨率为1.2米) + geotransform := [6]float64{ + topLeftLng, // top left x + resLng, // w-e pixel resolution + 0, // rotation, 0 if image is "north up" + topLeftLat, // top left y + 0, // rotation, 0 if image is "north up" + -resLat, // n-s pixel resolution (negative value) + } + err := ds.SetGeoTransform(geotransform) + if err != nil { + log.Errorf("Failed to set GeoTransform: %v", err) + return + } + + // 设置投影为 WGS84 + srs, err := godal.NewSpatialRefFromEPSG(4326) // gdal.CreateSpatialReference("") + if err != nil { + log.Errorf("Failed to set spatial reference: %v", err) + return + } + + projWKT, err := srs.WKT() + if err != nil { + log.Errorf("Failed to convert spatial reference to WKT: %v", err) + return + } + + err = ds.SetProjection(projWKT) + if err != nil { + log.Errorf("Failed to set projection: %v", err) + } + + // 设置一些常见的元数据(可选) + ds.SetMetadata("TIFFTAG_DATETIME", time.Now().String()) + ds.SetMetadata("TIFFTAG_SOFTWARE", "StarWiz-SJY01-IMAGE-PROC") +}