DFT 带阻滤波器

cmd/proc.go

@@ -3,6 +3,7 @@ package main
 import (
 	"os"
 	"path/filepath"
+	"runtime"
 	"strings"
 
 	"github.com/airbusgeo/godal"
@@ -74,6 +75,9 @@ var procCmd = &cobra.Command{
 		}
 
 		reg.SaveScenesToTiff(panScenes, mssScenes)
+		producer.CleanScenes(panScenes)
+		producer.CleanScenes(mssScenes)
+		runtime.GC()
 
 		if saveStrip {
 			reg.SaveOriginalPanToGDALGTiff(reg.Params.PanTiffFile)

config/config.yaml

@@ -0,0 +1,20 @@
+log_level: 5
+coregistration:
+  mss_bands: 4
+  pixel_bytes: 2
+  mss_width: 2336
+  pan_width: 9344
+  cr_block_nw: 8
+  cr_block_nh: 4
+  overlapped_block_lines: 3000
+  cr_resample_method: "down_sample_pan"
+  fus_band_order: "RGB"
+
+radiation:
+  pan_remove_hf_noise: true
+  mss_remove_hf_noise: true
+  hf_radius_ratio: 0.49
+  min_hist_level: 0.3
+  max_hist_level: 0.6
+  hf_band_stop_width: 24 # 带阻滤波器宽度
+  scene_moment_matching: false

pkg/config/config.go

@@ -3,7 +3,8 @@ package config
 import "github.com/sirupsen/logrus"
 
 type Config struct {
-	CRConfig CoRegistrationConfig `yaml:"cr_config" mapstructure:"cr_config"`
+	CoRegistration CoRegistrationConfig `yaml:"coregistration" mapstructure:"coregistration"`
+	Radiation      RadiationConfig      `yaml:"radiation" mapstructure:"radiation"`
 	LogLevel       logrus.Level         `yaml:"log_level" mapstructure:"log_level"`
 }
 
@@ -16,7 +17,17 @@ type CoRegistrationConfig struct {
 	CRBlockNH            int    `yaml:"cr_block_nh" mapstructure:"cr_block_nh"`
 	OverlappedBlockLines int    `yaml:"overlapped_block_lines" mapstructure:"overlapped_block_lines"`
 	CRResampleMethod     string `yaml:"cr_resample_method" mapstructure:"cr_resample_method"`
-	FUSBandOrder         string `yaml:"fu_band_order" mapstructure:"fu_band_order"`
+	FUSBandOrder         string `yaml:"fus_band_order" mapstructure:"fus_band_order"`
+}
+
+type RadiationConfig struct {
+	PANRemoveHfNoise    bool    `yaml:"pan_remove_hf_noise" mapstructure:"pan_remove_hf_noise"`
+	MSSRemoveHfNoise    bool    `yaml:"mss_remove_hf_noise" mapstructure:"mss_remove_hf_noise"`
+	SceneMomentMatching bool    `yaml:"scene_moment_matching" mapstructure:"scene_moment_matching"`
+	HfRadiusRatio       float64 `yaml:"hf_radius_ratio" mapstructure:"hf_radius_ratio"`
+	MinHistLevel        float64 `yaml:"min_hist_level" mapstructure:"min_hist_level"`
+	MaxHistLevel        float64 `yaml:"max_hist_level" mapstructure:"max_hist_level"`
+	HFBandStopWidth     int     `yaml:"hf_band_stop_width" mapstructure:"hf_band_stop_width"`
 }
 
 var GCONFIG Config
@@ -24,7 +35,7 @@ var GCONFIG Config
 func init() {
 	GCONFIG = Config{
 		LogLevel: logrus.DebugLevel,
-		CRConfig: CoRegistrationConfig{
+		CoRegistration: CoRegistrationConfig{
 			MssBands:             4,
 			PixelBytes:           2,
 			PanWidth:             9344,
@@ -35,5 +46,14 @@ func init() {
 			CRResampleMethod:     "down_sample_pan",
 			FUSBandOrder:         "RGB",
 		},
+		Radiation: RadiationConfig{
+			PANRemoveHfNoise:    true,
+			MSSRemoveHfNoise:    true,
+			SceneMomentMatching: false,
+			MinHistLevel:        0.3,
+			MaxHistLevel:        0.6,
+			HfRadiusRatio:       0.49,
+			HFBandStopWidth:     24,
+		},
 	}
 }

pkg/producer/aux.go

@@ -113,6 +113,13 @@ func (r *Registrator) SetSceneBoundary(scene *Scene) (topLeft, bottomRight orb.P
 	xResolution := W0 / float64(scene.Width)
 	yResolution := H0 / float64(scene.Height)
 	scene.Meta.Gsd = math.Min(xResolution, yResolution)
+
+	// FIXME: 临时设置分辨率
+	if scene.Meta.Gsd < 2 {
+		scene.Meta.Gsd = 1.3
+	} else {
+		scene.Meta.Gsd = 5.2
+	}
 	log.Debug("resolution x: ", xResolution)
 	log.Debug("resolution y: ", yResolution)
 

pkg/producer/filter.go

@@ -1,78 +0,0 @@
-package producer
-
-import (
-	"image"
-	"math"
-
-	log "github.com/sirupsen/logrus"
-
-	"gocv.io/x/gocv"
-)
-
-func PANFilter(panImage16UC1 gocv.Mat) gocv.Mat {
-	log.Println("Applying PAN filter...")
-	// 转换为浮点数类型进行傅里叶变换
-	imgFloat := gocv.NewMat()
-	panImage16UC1.ConvertTo(&imgFloat, gocv.MatTypeCV32F)
-
-	// 使用傅里叶变换
-	planes := gocv.NewMat()
-	gocv.Merge([]gocv.Mat{imgFloat, gocv.NewMatWithSize(imgFloat.Rows(), imgFloat.Cols(), gocv.MatTypeCV32F)}, &planes)
-	dft := gocv.NewMat()
-	gocv.DFT(planes, &dft, gocv.DftComplexOutput)
-
-	// 转换DFT图像，使低频分量位于中心
-	dftShifted := shiftDFT(dft)
-
-	// 创建掩膜
-	rows, cols := panImage16UC1.Rows(), panImage16UC1.Cols()
-	crow, ccol := rows/2, cols/2
-	mask := gocv.NewMatWithSize(rows, cols, gocv.MatTypeCV32FC2)
-	mask.SetTo(gocv.NewScalar(1.0, 1.0, 0.0, 0.0))
-
-	r := 30 // 调整这个参数来改变滤波效果
-	for i := 0; i < rows; i++ {
-		for j := 0; j < cols; j++ {
-			if math.Abs(float64(i-crow)) <= float64(r) && math.Abs(float64(j-ccol)) <= float64(r) {
-				mask.SetFloatAt(i, j*2, 0)
-				mask.SetFloatAt(i, j*2+1, 0)
-			}
-		}
-	}
-
-	// 应用掩膜并进行反向DFT
-	filtered := gocv.NewMat()
-	gocv.MulSpectrums(dftShifted, mask, &filtered, 0)
-	filteredShifted := shiftDFT(filtered)
-	idft := gocv.NewMat()
-	gocv.DFT(filteredShifted, &idft, gocv.DftInverse|gocv.DftRealOutput)
-
-	// 标准化图像到0-65535范围
-	gocv.Normalize(idft, &idft, 0, 65535, gocv.NormMinMax)
-	idft.ConvertTo(&idft, gocv.MatTypeCV16U)
-
-	return idft
-}
-
-// shiftDFT 将DFT结果进行频谱平移
-func shiftDFT(src gocv.Mat) gocv.Mat {
-	rows, cols := src.Rows(), src.Cols()
-	crow, ccol := rows/2, cols/2
-
-	q0 := src.Region(image.Rect(0, 0, ccol, crow))
-	q1 := src.Region(image.Rect(ccol, 0, cols, crow))
-	q2 := src.Region(image.Rect(0, crow, ccol, rows))
-	q3 := src.Region(image.Rect(ccol, crow, cols, rows))
-
-	tmp := gocv.NewMatWithSize(crow, ccol, src.Type())
-
-	q0.CopyTo(&tmp)
-	q3.CopyTo(&q0)
-	tmp.CopyTo(&q3)
-
-	q1.CopyTo(&tmp)
-	q2.CopyTo(&q1)
-	tmp.CopyTo(&q2)
-
-	return src
-}

pkg/producer/output.go

@@ -10,10 +10,11 @@ import (
 	"github.com/airbusgeo/godal"
 	log "github.com/sirupsen/logrus"
 	"gocv.io/x/gocv"
+	"starwiz.cn/sjy01/image-proc/pkg/utils"
 )
 
 func (r *Registrator) SaveOriginalPanToGDALGTiff(tiffFile string) error {
-	err := savePanToGDALGTiff(r.PanImage, 0, 0, tiffFile, PanResolution)
+	err := utils.SavePanToGDALGTiff(r.PanImage, 0, 0, tiffFile, PanResolution)
 	if err != nil {
 		return err
 	}
@@ -21,53 +22,11 @@ func (r *Registrator) SaveOriginalPanToGDALGTiff(tiffFile string) error {
 	return nil
 }
 
-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 (r *Registrator) SaveRegisteredMssToGDALGTiff(tiffFile string) error {
 	r.rgbirImage = gocv.NewMat()
 	gocv.Merge(r.registeredMssImages[:], &r.rgbirImage)
 
-	err := SaveBGRToGDALGTiff(r.rgbirImage,
+	err := utils.SaveBGRToGDALGTiff(r.rgbirImage,
 		4, 0, 0, MssResolution,
 		[]godal.ColorInterp{godal.CIBlue, godal.CIGreen, godal.CIRed, godal.CIUndefined},
 		tiffFile)
@@ -78,63 +37,6 @@ func (r *Registrator) SaveRegisteredMssToGDALGTiff(tiffFile string) error {
 	return nil
 }
 
-func SaveBGRToGDALGTiff(bgr gocv.Mat,
-	bands int,
-	topLeftX, topLeftY float64,
-	resolution float64,
-	colorInterps []godal.ColorInterp, tiffFile string) error {
-	width := bgr.Cols()
-	height := bgr.Rows()
-
-	// 创建一个二维切片来存储图像数据
-	data := make([][]uint16, bands)
-	for i := range data {
-		data[i] = make([]uint16, width*height)
-	}
-
-	// 从gocv.Mat中提取数据
-	for y := 0; y < height; y++ {
-		for x := 0; x < width; x++ {
-			for b := 0; b < bands; b++ {
-				data[b][y*width+x] = uint16(bgr.GetShortAt(y, x*bands+b))
-			}
-		}
-	}
-
-	ds, err := godal.Create(godal.GTiff,
-		tiffFile,
-		bands,
-		godal.UInt16,
-		width, height)
-	if err != nil {
-		log.Error("Error creating TIFF file: ", err)
-		return err
-	}
-	defer ds.Close()
-
-	// ds.SetMetadata("NBITS", "16")
-	setGeoTransform(ds, topLeftX, topLeftY, resolution)
-
-	for b := 0; b < bands; b++ {
-		band := ds.Bands()[b]
-		band.SetColorInterp(colorInterps[b])
-		err := band.IO(godal.IOWrite,
-			0, 0,
-			data[b],
-			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 BGR MSS to ", tiffFile)
-
-	return nil
-}
-
 func (r *Registrator) BytesToRaw(mssData []byte, filePath string) error {
 	f, err := os.OpenFile(filePath, os.O_CREATE|os.O_TRUNC|os.O_WRONLY, 0777)
 	if err != nil {

pkg/producer/scenes.go

@@ -13,7 +13,9 @@ import (
 	"github.com/paulmach/orb/geojson"
 	log "github.com/sirupsen/logrus"
 	"gocv.io/x/gocv"
+	"starwiz.cn/sjy01/image-proc/pkg/config"
 	"starwiz.cn/sjy01/image-proc/pkg/rrc"
+	"starwiz.cn/sjy01/image-proc/pkg/utils"
 )
 
 type Scene struct {
@@ -34,6 +36,12 @@ func (s *Scene) Cleanup() {
 	}
 }
 
+func CleanScenes(scenes []*Scene) {
+	for _, scene := range scenes {
+		scene.Cleanup()
+	}
+}
+
 // 对 PAN 和 配准后的MSS 在 Y 方向进行分景，景之间有25%的重叠
 // 默认分景大小：
 // MSS  2336 * 2336 - 1764
@@ -67,9 +75,18 @@ 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))
-		matFiltered := rrc.HFNoiseFilter(mat, float64(mat.Cols())*0.45)
+		if config.GCONFIG.Radiation.PANRemoveHfNoise {
-		scene.Mat = append(scene.Mat, matFiltered)
+			log.Println("applying hf noise filter on", scene.SceneId)
+			matFiltered := rrc.HFNoiseFilter(mat, float64(mat.Cols())*config.GCONFIG.Radiation.HfRadiusRatio)
 			mat.Close()
+			mat = matFiltered
+		}
+
+		if config.GCONFIG.Radiation.SceneMomentMatching {
+			rrc.DoMomentMatching(mat)
+		}
+
+		scene.Mat = append(scene.Mat, mat)
 		panScenes = append(panScenes, scene)
 	}
 
@@ -95,16 +112,26 @@ func (r *Registrator) SubScenes() (panScenes []*Scene, mssScenes []*Scene, err e
 			continue
 		}
 
-		for band := 0; band < 4; band++ {
-			mat := r.registeredMssImages[band].Region(image.Rect(0, i*hMSS, 2336, y1))
-			matFiltered := rrc.HFNoiseFilter(mat, float64(mat.Cols())*0.45)
-			mat.Close()
-			scene.Mat = append(scene.Mat, matFiltered)
-		}
-
 		name := filepath.Base(r.Params.MssTiffFile)
 		name = strings.TrimSuffix(name, ".tiff")
 		scene.SceneId = fmt.Sprintf("%s_%03d", name, i+1)
+
+		for band := 0; band < 4; band++ {
+			mat := r.registeredMssImages[band].Region(image.Rect(0, i*hMSS, 2336, y1))
+			if config.GCONFIG.Radiation.MSSRemoveHfNoise {
+				log.Println("applying hf noise filter on", scene.SceneId)
+				matFiltered := rrc.HFNoiseFilter(mat, float64(mat.Cols())*config.GCONFIG.Radiation.HfRadiusRatio/4)
+				mat.Close()
+				mat = matFiltered
+			}
+
+			if config.GCONFIG.Radiation.SceneMomentMatching {
+				rrc.DoMomentMatching(mat)
+			}
+
+			scene.Mat = append(scene.Mat, mat)
+		}
+
 		mssScenes = append(mssScenes, scene)
 	}
 
@@ -127,7 +154,7 @@ func (r *Registrator) SaveScenesToTiff(panScenes []*Scene, mssScenes []*Scene) e
 		scene.Meta = r.makeProductMeta(scene)
 		r.SetSceneBoundary(scene)
 
-		err := savePanToGDALGTiff(scene.Mat[0],
+		err := utils.SavePanToGDALGTiff(scene.Mat[0],
 			scene.Meta.Corners.UpperLeft.Longitude,
 			scene.Meta.Corners.UpperLeft.Latitude,
 			scene.Tiff,
@@ -175,7 +202,7 @@ func (r *Registrator) SaveScenesToTiff(panScenes []*Scene, mssScenes []*Scene) e
 		r.SetSceneBoundary(scene)
 
 		rgbirImage, _ := r.MergeMSSToBGRNIR(scene.Mat)
-		err := SaveBGRToGDALGTiff(rgbirImage,
+		err := utils.SaveBGRToGDALGTiff(rgbirImage,
 			4,
 			scene.Meta.Corners.UpperLeft.Longitude,
 			scene.Meta.Corners.UpperLeft.Latitude,

pkg/producer/util.go

@@ -2,55 +2,10 @@ package producer
 
 import (
 	"os"
-	"time"
 
-	"github.com/airbusgeo/godal"
 	log "github.com/sirupsen/logrus"
-	"starwiz.cn/sjy01/image-proc/pkg/calculator"
 )
 
-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")
-}
-
 func sizeOfFile(file string) int64 {
 	fileInfo, err := os.Stat(file)
 	if err != nil {

pkg/rrc/filter_test.go

@@ -1,6 +1,7 @@
 package rrc
 
 import (
+	"os"
 	"testing"
 
 	"github.com/airbusgeo/godal"
@@ -9,11 +10,14 @@ import (
 )
 
 func TestFilter(t *testing.T) {
+	os.MkdirAll("/Users/lan/workspace/sjy01/image-proc/data/fft", 0755)
 	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/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"
-
+	// tif = "/Users/lan/workspace/sjy01/image-proc/data/052022-no-rrc/010/MSS/SJY01_MSS_20240520_115428_052022_103_010_L1A.tiff"
+	tif = "/Users/lan/workspace/sjy01/image-proc/data/052022-no-rrc/006/PAN/SJY01_PAN_20240520_115428_052022_103_006_L1A.tiff"
+	tif = "/Users/lan/workspace/sjy01/image-proc/data/060622-RRC/002/PAN/SJY01_PAN_20240606_012004_060622_065_002_L1A.tiff"
 	ds, err := godal.Open(tif)
 	if err != nil {
 		t.Error(err)
@@ -37,7 +41,23 @@ func TestFilter(t *testing.T) {
 		}
 	}
 
-	m0 := HFNoiseFilter(img, float64(img.Cols())*0.45)
+	m0 := HFNoiseFilter(img, float64(img.Cols())*0.49)
 	defer m0.Close()
 	utils.SavePanToGDALGTiff(m0, 0, 0, "/Users/lan/workspace/sjy01/image-proc/data/fft/filtered.tiff", 1.3)
 }
+
+func TestMomentMatching(t *testing.T) {
+	raw := "/Users/lan/workspace/sjy01/preprocessing/demo/output/051622/SJY01_PAN_20240516_101236_051622_096.RAW"
+	data, err := os.ReadFile(raw)
+	if err != nil {
+		t.Error(err)
+	}
+
+	height := len(data) / (9344 * 2)
+	img, err := gocv.NewMatFromBytes(height, 9344, gocv.MatTypeCV16U, data)
+	if err != nil {
+		t.Error(err)
+	}
+
+	DoMomentMatching(img)
+}

pkg/rrc/helper.go

@@ -1,12 +1,5 @@
 package rrc
 
-import (
-	"runtime"
-
-	"github.com/dustin/go-humanize"
-	log "github.com/sirupsen/logrus"
-)
-
 func searchVL(V []float64, Sik float64) int {
 	left, right := 0, len(V)-2
 
@@ -31,16 +24,3 @@ func searchVL(V []float64, Sik float64) int {
 	// Return -1 if no such position is found
 	return -1
 }
-
-var lastTotalFreed uint64
-
-func printMemStats() {
-	var m runtime.MemStats
-	runtime.ReadMemStats(&m)
-	log.Printf("[Memory] Alloc = %v TotalAlloc=%v Just Freed = %v Sys = %v NumGc=%v",
-		humanize.Bytes(m.Alloc),
-		humanize.Bytes(m.TotalAlloc),
-		humanize.Bytes(((m.TotalAlloc - m.Alloc) - lastTotalFreed)),
-		humanize.Bytes(m.Sys), m.NumGC)
-	lastTotalFreed = m.TotalAlloc - m.Alloc
-}

pkg/rrc/hf_filter.go

@@ -4,69 +4,74 @@ import (
 	"image"
 	"math"
 
-	"github.com/sirupsen/logrus"
+	log "github.com/sirupsen/logrus"
 	"gocv.io/x/gocv"
+	"starwiz.cn/sjy01/image-proc/pkg/config"
 )
 
 // filter high frequence noise
 func HFNoiseFilter(img gocv.Mat, radius float64) gocv.Mat {
-	logrus.Println("applying hf noise filter...")
+	minGrayVal, maxGrayVal, _, _ := gocv.MinMaxLoc(img)
-	imgFloat := gocv.NewMat()
+	log.Println("minGrayVal:", minGrayVal, "maxGrayVal:", maxGrayVal)
-	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()
 
+	// 扩展图像尺寸，使其宽高为2的幂次方
+	// optRows := gocv.GetOptimalDFTSize(img32f.Rows())
+	// optCols := gocv.GetOptimalDFTSize(img32f.Cols())
+	// padded := gocv.NewMatWithSize(optRows, optCols, gocv.MatTypeCV32F)
+	// defer padded.Close()
+	// gocv.CopyMakeBorder(img32f, &padded,
+	// 	0, optRows-img32f.Rows(), 0, optCols-img32f.Cols(),
+	// 	gocv.BorderConstant, color.RGBA{0, 0, 0, 0})
+
 	// 获取图像尺寸
 	rows, cols := img32f.Rows(), img32f.Cols()
 
-	// zeros := gocv.NewMatFromScalar(gocv.NewScalar(0, 0, 0, 0), gocv.MatTypeCV32F)
-
 	// 执行 DFT
 	DFT := gocv.NewMat()
 	defer DFT.Close()
-
+	gocv.Merge([]gocv.Mat{img32f, gocv.Zeros(rows, cols, gocv.MatTypeCV32F)}, &DFT)
-	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)
+
+	// visualizeDFT(DFT, "/Users/lan/workspace/sjy01/image-proc/data/fft/dft.jpg")
+
+	// 低通滤波器
+	pss := gocv.Split(DFT)
+	defer pss[0].Close()
+	defer pss[1].Close()
+	filter := gocv.NewMatWithSize(rows, cols, gocv.MatTypeCV32F)
+	defer filter.Close()
+	// createLowPassFilter(&filter, radius)
+	// createBandBlockFilter(&filter, radius, radius+0.015*float64(cols))
+	createFixedBlockFilter(&filter, config.GCONFIG.Radiation.HFBandStopWidth)
+	gocv.Multiply(pss[0], filter, &pss[0])
+	gocv.Multiply(pss[1], filter, &pss[1])
+	gocv.Merge(pss, &DFT)
+
+	// visualizeDFT(DFT, "/Users/lan/workspace/sjy01/image-proc/data/fft/dft_filtered.jpg")
+
+	// 频谱位置还原
+	fftshift(&DFT)
+
 	// 执行 IDFT
 	IDFT := gocv.NewMat()
 	gocv.DFT(DFT, &IDFT, gocv.DftInverse|gocv.DftRealOutput)
 
-	// 归一化到0-65535
+	// 归一化到原始图像范围
-	gocv.Normalize(IDFT, &IDFT, 0, 8192, gocv.NormMinMax)
+	gocv.Normalize(IDFT, &IDFT, float64(minGrayVal), float64(maxGrayVal), gocv.NormMinMax)
 	IDFT.ConvertTo(&IDFT, gocv.MatTypeCV16U)
-
-	// 保存 IDFT 结果
 	return IDFT
 }
 
 // 创建低通滤波器
-func createLowPassFilter(filter *gocv.Mat, radius int) {
+func createLowPassFilter(filter *gocv.Mat, radius float64) {
 	rows, cols := filter.Rows(), filter.Cols()
 	center := image.Pt(cols/2, rows/2)
 	for y := 0; y < rows; y++ {
@@ -81,6 +86,38 @@ func createLowPassFilter(filter *gocv.Mat, radius int) {
 	}
 }
 
+// 创建带阻滤波器
+func createBandBlockFilter(filter *gocv.Mat, radius0, radius1 float64) {
+	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(radius0) || dist >= float64(radius1) {
+				filter.SetFloatAt(y, x, 1.0)
+			} else {
+				filter.SetFloatAt(y, x, 0.0)
+			}
+		}
+	}
+}
+
+func createFixedBlockFilter(filter *gocv.Mat, ksize 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++ {
+			if (x >= center.X-ksize/2 && x <= center.X+ksize/2 &&
+				(y < int(ksize) || y >= rows-int(ksize))) ||
+				((x < ksize || x > cols-ksize) && (y >= center.Y-ksize/2 && y <= center.Y+ksize/2)) {
+				filter.SetFloatAt(y, x, 0.0)
+			} else {
+				filter.SetFloatAt(y, x, 1.0)
+			}
+		}
+	}
+}
+
 // fftshift 将频率直流分量移动到图像中心
 func fftshift(input *gocv.Mat) {
 	cx, cy := input.Cols()/2, input.Rows()/2
@@ -100,19 +137,25 @@ func fftshift(input *gocv.Mat) {
 	tmp.CopyTo(&q2)
 }
 
-func SaveDFT(DFT gocv.Mat, file string) {
+func visualizeDFT(DFT gocv.Mat, file string) {
 	pss := gocv.Split(DFT)
 	planes := gocv.NewMat()
+	defer planes.Close()
+
 	gocv.Magnitude(pss[0], pss[1], &planes)
 	// 进行对数尺度变换，便于可视化
 	ones := gocv.NewMatFromScalar(gocv.NewScalar(1, 0, 0, 0), gocv.MatTypeCV64F)
+	defer ones.Close()
 	gocv.Add(planes, ones, &planes)
 	gocv.Log(planes, &planes)
 
 	// 归一化
-	gocv.Normalize(planes, &planes, 0, 1, gocv.NormMinMax)
+	gocv.Normalize(planes, &planes, 0, 255, gocv.NormMinMax)
+	planes.ConvertTo(&planes, gocv.MatTypeCV8U)
+	out := gocv.NewMatWithSize(planes.Rows()/4, planes.Cols()/4, gocv.MatTypeCV8U)
+	defer out.Close()
+	gocv.Resize(planes, &out, image.Point{planes.Cols() / 4, planes.Rows() / 4}, 0, 0, gocv.InterpolationLinear)
 
 	// 保存DFT结果
-	gocv.IMWrite(file, planes)
+	gocv.IMWrite(file, out)
-
 }

pkg/rrc/hf_filter_fft.go

@@ -8,6 +8,8 @@ import (
 	"gocv.io/x/gocv"
 )
 
+// FIXME: still need to implement the HFNoiseFilterFFT function
+
 func HFNoiseFilterFFT(input gocv.Mat, radius int) gocv.Mat {
 	// 将灰度图像转换为浮点数数组
 	height := input.Rows()

pkg/rrc/moment_matching.go

@@ -10,10 +10,11 @@ import (
 // Destriping multisensor imagery with moment matching [Gadallah, 2000]
 func DoMomentMatching(originalImg gocv.Mat) {
 	probes := originalImg.Cols()
-	log.Printf("do moment matching for %d probes", probes)
+	log.Printf("do moment matching for %d probes, %d rows", probes, originalImg.Rows())
 	// 第i个探元的像元均值和标准差
 	means := make([]float64, probes)
 	stds := make([]float64, probes)
+
 	// 计算每个探元的均值和标准差
 	for x := 0; x < originalImg.Cols(); x++ {
 		var total int64
@@ -34,16 +35,17 @@ func DoMomentMatching(originalImg gocv.Mat) {
 	}
 
 	// 列参考值和列参考标准差
-	var mu float64
+	var mu, sig float64
-	var sig float64
 	for x := 0; x < probes; x++ {
 		mu += means[x]
 		sig += stds[x]
 	}
+
 	mu = mu / float64(probes)
 	sig = sig / float64(probes)
 
 	log.Printf("mean reference value: %f, std reference value: %f", mu, sig)
+
 	// 修正 DN_adjusted[i] = (DN[i] - means[i]) *sig/stds[i]+mu
 	for x := 0; x < originalImg.Cols(); x++ {
 		for y := 0; y < originalImg.Rows(); y++ {

pkg/rrc/rrc.go

@@ -10,6 +10,7 @@ import (
 
 	log "github.com/sirupsen/logrus"
 	"gocv.io/x/gocv"
+	"starwiz.cn/sjy01/image-proc/pkg/utils"
 )
 
 // Relative Radiation Correction
@@ -49,7 +50,7 @@ func (rrc *RRC) Close() {}
 
 // 统计探元灰度的累积概率密度
 func (rrc *RRC) StatisticalPAN(dsfile string) error {
-	printMemStats()
+	utils.PrintMemStats()
 	f, err := os.Open(dsfile)
 	if err != nil {
 		return err
@@ -100,7 +101,7 @@ func (rrc *RRC) StatisticalPAN(dsfile string) error {
 			rrc.Histograms[0].sum([]*ProbeHistogram{&hist})
 			hist.free()
 			runtime.GC()
-			printMemStats()
+			utils.PrintMemStats()
 
 			mutex.Unlock()
 
@@ -118,7 +119,7 @@ func (rrc *RRC) StatisticalPAN(dsfile string) error {
 }
 
 func (rrc *RRC) StatisticalMSS(dsfile string) error {
-	printMemStats()
+	utils.PrintMemStats()
 	f, err := os.Open(dsfile)
 	if err != nil {
 		return err
@@ -179,7 +180,7 @@ func (rrc *RRC) StatisticalMSS(dsfile string) error {
 				rrc.Histograms[i+1].sum([]*ProbeHistogram{&hist})
 				hist.free()
 				runtime.GC()
-				printMemStats()
+				utils.PrintMemStats()
 
 				mutex.Unlock()
 			}

pkg/utils/memory.go

@@ -0,0 +1,21 @@
+package utils
+
+import (
+	"log"
+	"runtime"
+
+	"github.com/dustin/go-humanize"
+)
+
+var lastTotalFreed uint64
+
+func PrintMemStats() {
+	var m runtime.MemStats
+	runtime.ReadMemStats(&m)
+	log.Printf("[Memory] Alloc = %v TotalAlloc=%v Just Freed = %v Sys = %v NumGc=%v",
+		humanize.Bytes(m.Alloc),
+		humanize.Bytes(m.TotalAlloc),
+		humanize.Bytes(((m.TotalAlloc - m.Alloc) - lastTotalFreed)),
+		humanize.Bytes(m.Sys), m.NumGC)
+	lastTotalFreed = m.TotalAlloc - m.Alloc
+}

pkg/utils/tiff.go

@@ -47,7 +47,64 @@ func SavePanToGDALGTiff(pan gocv.Mat, topLeftX, topLeftY float64, tiffFile strin
 		return err
 	}
 
-	log.Info("Saved pan image to ", tiffFile)
+	log.Info("Saved image to ", tiffFile)
+
+	return nil
+}
+
+func SaveBGRToGDALGTiff(bgr gocv.Mat,
+	bands int,
+	topLeftX, topLeftY float64,
+	resolution float64,
+	colorInterps []godal.ColorInterp, tiffFile string) error {
+	width := bgr.Cols()
+	height := bgr.Rows()
+
+	// 创建一个二维切片来存储图像数据
+	data := make([][]uint16, bands)
+	for i := range data {
+		data[i] = make([]uint16, width*height)
+	}
+
+	// 从gocv.Mat中提取数据
+	for y := 0; y < height; y++ {
+		for x := 0; x < width; x++ {
+			for b := 0; b < bands; b++ {
+				data[b][y*width+x] = uint16(bgr.GetShortAt(y, x*bands+b))
+			}
+		}
+	}
+
+	ds, err := godal.Create(godal.GTiff,
+		tiffFile,
+		bands,
+		godal.UInt16,
+		width, height)
+	if err != nil {
+		log.Error("Error creating TIFF file: ", err)
+		return err
+	}
+	defer ds.Close()
+
+	// ds.SetMetadata("NBITS", "16")
+	setGeoTransform(ds, topLeftX, topLeftY, resolution)
+
+	for b := 0; b < bands; b++ {
+		band := ds.Bands()[b]
+		band.SetColorInterp(colorInterps[b])
+		err := band.IO(godal.IOWrite,
+			0, 0,
+			data[b],
+			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 BGR MSS to ", tiffFile)
 
 	return nil
 }