DFT 低通滤波器

2024-06-19 22:42:31 +08:00
parent 6e5ffe1ab2
commit e13038474a
8 changed files with 337 additions and 2 deletions
--- 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
--- 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=
--- 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)
--- a/pkg/rrc/filter_test.go
+++ 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)
+}
--- a/pkg/rrc/hf_filter.go
+++ 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)
+
+}
--- a/pkg/rrc/hf_filter_fft.go
+++ 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
+}
--- a/pkg/rrc/mean_filter.go
+++ 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
+}
--- a/pkg/utils/tiff.go
+++ 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")
+}