From d99b8a974063625206b692f0fa83343ebd6ba9e9 Mon Sep 17 00:00:00 2001 From: nuknal Date: Fri, 21 Jun 2024 12:47:14 +0800 Subject: [PATCH] =?UTF-8?q?DFT=20=E5=B8=A6=E9=98=BB=E6=BB=A4=E6=B3=A2?= =?UTF-8?q?=E5=99=A8?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- cmd/proc.go | 4 ++ config/config.yaml | 20 +++++++ pkg/config/config.go | 28 +++++++-- pkg/producer/aux.go | 7 +++ pkg/producer/filter.go | 78 ------------------------ pkg/producer/output.go | 104 +------------------------------- pkg/producer/scenes.go | 51 ++++++++++++---- pkg/producer/util.go | 45 -------------- pkg/rrc/filter_test.go | 26 +++++++- pkg/rrc/helper.go | 20 ------- pkg/rrc/hf_filter.go | 119 +++++++++++++++++++++++++------------ pkg/rrc/hf_filter_fft.go | 2 + pkg/rrc/moment_matching.go | 8 ++- pkg/rrc/rrc.go | 9 +-- pkg/utils/memory.go | 21 +++++++ pkg/utils/tiff.go | 59 +++++++++++++++++- 16 files changed, 292 insertions(+), 309 deletions(-) create mode 100644 config/config.yaml delete mode 100644 pkg/producer/filter.go create mode 100644 pkg/utils/memory.go diff --git a/cmd/proc.go b/cmd/proc.go index fb597ef..5e613ac 100644 --- a/cmd/proc.go +++ b/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) diff --git a/config/config.yaml b/config/config.yaml new file mode 100644 index 0000000..46b451f --- /dev/null +++ b/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 \ No newline at end of file diff --git a/pkg/config/config.go b/pkg/config/config.go index 6a1406a..7047cfd 100644 --- a/pkg/config/config.go +++ b/pkg/config/config.go @@ -3,8 +3,9 @@ package config import "github.com/sirupsen/logrus" type Config struct { - CRConfig CoRegistrationConfig `yaml:"cr_config" mapstructure:"cr_config"` - LogLevel logrus.Level `yaml:"log_level" mapstructure:"log_level"` + CoRegistration CoRegistrationConfig `yaml:"coregistration" mapstructure:"coregistration"` + Radiation RadiationConfig `yaml:"radiation" mapstructure:"radiation"` + LogLevel logrus.Level `yaml:"log_level" mapstructure:"log_level"` } type CoRegistrationConfig struct { @@ -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, + }, } } diff --git a/pkg/producer/aux.go b/pkg/producer/aux.go index a5a9125..6300c3f 100644 --- a/pkg/producer/aux.go +++ b/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) diff --git a/pkg/producer/filter.go b/pkg/producer/filter.go deleted file mode 100644 index c370d0f..0000000 --- a/pkg/producer/filter.go +++ /dev/null @@ -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 -} diff --git a/pkg/producer/output.go b/pkg/producer/output.go index f5a2d76..b0caf8c 100644 --- a/pkg/producer/output.go +++ b/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 { diff --git a/pkg/producer/scenes.go b/pkg/producer/scenes.go index 82e31ae..d0e9f2b 100644 --- a/pkg/producer/scenes.go +++ b/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) - scene.Mat = append(scene.Mat, matFiltered) - mat.Close() + if config.GCONFIG.Radiation.PANRemoveHfNoise { + 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, diff --git a/pkg/producer/util.go b/pkg/producer/util.go index 4516965..1399982 100644 --- a/pkg/producer/util.go +++ b/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 { diff --git a/pkg/rrc/filter_test.go b/pkg/rrc/filter_test.go index f52fd1a..7b7d801 100644 --- a/pkg/rrc/filter_test.go +++ b/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) +} diff --git a/pkg/rrc/helper.go b/pkg/rrc/helper.go index 4b9d3fa..c6b27ed 100644 --- a/pkg/rrc/helper.go +++ b/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 -} diff --git a/pkg/rrc/hf_filter.go b/pkg/rrc/hf_filter.go index 0377ebf..9cfbfe9 100644 --- a/pkg/rrc/hf_filter.go +++ b/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...") - 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) + minGrayVal, maxGrayVal, _, _ := gocv.MinMaxLoc(img) + log.Println("minGrayVal:", minGrayVal, "maxGrayVal:", maxGrayVal) // 将图像转换为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() - - 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.Merge([]gocv.Mat{img32f, 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) } diff --git a/pkg/rrc/hf_filter_fft.go b/pkg/rrc/hf_filter_fft.go index 723780f..569cd67 100644 --- a/pkg/rrc/hf_filter_fft.go +++ b/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() diff --git a/pkg/rrc/moment_matching.go b/pkg/rrc/moment_matching.go index 7b64d33..f8c9f05 100644 --- a/pkg/rrc/moment_matching.go +++ b/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 sig float64 + var mu, 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++ { diff --git a/pkg/rrc/rrc.go b/pkg/rrc/rrc.go index 14b4d4d..3cf42fc 100644 --- a/pkg/rrc/rrc.go +++ b/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() } diff --git a/pkg/utils/memory.go b/pkg/utils/memory.go new file mode 100644 index 0000000..ec6c60a --- /dev/null +++ b/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 +} diff --git a/pkg/utils/tiff.go b/pkg/utils/tiff.go index 53b6a19..4737199 100644 --- a/pkg/utils/tiff.go +++ b/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 }