lankun/sjy01-image-proc
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

adjust

Browse Source
This commit is contained in:
nuknal
2024-05-29 10:20:21 +08:00
parent 5eb5ae6a72
commit e15ae9247b
17 changed files with 276 additions and 45 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
cmd/proc.go
View File

@@ -8,11 +8,11 @@ import (
"github.com/airbusgeo/godal"
log "github.com/sirupsen/logrus"
"github.com/spf13/cobra"
imageproc "starwiz.cn/sjy01/image-proc"
producer "starwiz.cn/sjy01/image-proc/producer"
)
var (
params imageproc.Params
params producer.Params
)
var procCmd = &cobra.Command{
@@ -20,7 +20,7 @@ var procCmd = &cobra.Command{
Short: "process images",
Long: `process images`,
Run: func(cmd *cobra.Command, args []string) {
reg := imageproc.NewRegistrator(imageproc.DownSampled)
reg := producer.NewRegistrator(producer.DownSampled)
reg.Params = params
reg.Params.MssTiffFile = filepath.Join(params.OutputDir, strings.TrimSuffix(filepath.Base(params.MssRawFile), filepath.Ext(params.MssRawFile))+".tiff")
reg.Params.PanTiffFile = filepath.Join(params.OutputDir, strings.TrimSuffix(filepath.Base(params.PanRawFile), filepath.Ext(params.PanRawFile))+".tiff")

0
filter.go → producer/filter.go
View File

0
fit.go → producer/fit.go
View File

0
gdal_pansharpen.go → producer/gdal_pansharpen.go
View File

0
hdr.go → producer/hdr.go
View File

15
gdal_pansharpen_test.go → producer/image_proc_test.go
View File

@@ -1,6 +1,10 @@
package imageproc
import "testing"
import (
"testing"
"github.com/airbusgeo/godal"
)
func TestGDALPansharpen(t *testing.T) {
// err := GDALPansharpen("data/052022/SJY01_PAN_20240520_115428_052022_103_010.tiff",
@@ -11,3 +15,12 @@ func TestGDALPansharpen(t *testing.T) {
GDALTranslate("data/052022/SJY01_MSS_20240520_115428_052022_103_010_FUS.tiff", "")
}
func TestGTiffToJPG(t *testing.T) {
godal.RegisterAll()
err := GTiffToJPG("data/052022/010/SJY01_MSS_20240520_115428_052022_103_010.tiff",
"data/052022/SJY01_MSS_20240520_115428_052022_103_010.jpg")
if err != nil {
t.Error(err)
}
}

8
image_registration.go → producer/image_registration.go
View File

@@ -42,8 +42,8 @@ type Registrator struct {
shiftMutex sync.Mutex
phaseShifts [4][]PhaseShiftM
deltaXCoeffs [4][]float64 // 图像内畸变（非线性变换），捕捉图像在水平方向上引起的垂直方向的变形
deltaYCoeffs [4][]float64 // 图像内畸变（非线性变换），捕捉图像在垂直方向上引起的水平方向的变形
deltaXCoeffs [4][]float64 // 图像内畸变（线性变换），捕捉图像在水平方向上引起的X方向的变形
deltaYCoeffs [4][]float64 // 图像内畸变（非线性变换），捕捉图像在水平方向上引起的Y方向的变形
registeredMssImages [4]gocv.Mat // 配准后的MSS图像
rgbirImage gocv.Mat
@@ -221,7 +221,7 @@ func (r *Registrator) calcPhaseCorrelation(panImage gocv.Mat,
panBlock := panImage.Region(rect)
for band := 0; band < MssBands; band++ {
log.Println("Band", band+1, ", processing block", bh, rect)
log.Debug("Band", band+1, ", processing block", bh, rect)
mssBlock := mssImages[band].Region(rect)
// 处理每个分块
@@ -246,7 +246,7 @@ func (r *Registrator) calcPhaseCorrelation(panImage gocv.Mat,
for i := 0; i < MssBands; i++ {
for _, shift := range r.phaseShifts[i] {
if shift.response > 0.4 || shift.dx > 8 || shift.dy > 8 {
log.Printf("Band %d, block %d, dx=%f, dy=%f, response=%f\n",
log.Debugf("Band %d, block %d, dx=%f, dy=%f, response=%f",
i, shift.Block.coord.X, shift.dx, shift.dy, shift.response)
}
}

85
producer/jpg.go Normal file
View File

@@ -0,0 +1,85 @@
package imageproc
import (
"fmt"
"image"
"github.com/airbusgeo/godal"
log "github.com/sirupsen/logrus"
"gocv.io/x/gocv"
)
func GTiffToJPG(ftiff, fjpg string) error {
// 打开 TIFF 文件
ds, err := godal.Open(ftiff)
if err != nil {
log.Printf("Error opening TIFF file %s: %v", ftiff, err)
return err
}
defer ds.Close()
bands := ds.Bands()
log.Infof("TIFF file %s has %d bands", ftiff, len(bands))
if len(bands) < 3 {
err = fmt.Errorf("TIFF file %s has less than 3 bands", ftiff)
log.Error(err)
return err
}
// 获取图像大小
width := bands[0].Structure().SizeX
height := bands[0].Structure().SizeY
// 读取每个波段并转换为 8 位
img := gocv.NewMatWithSize(height, width, gocv.MatTypeCV16UC3)
for i := 1; i <= 3; i++ {
band := ds.Bands()[i-1]
data := make([]uint16, width*height)
err = band.Read(0, 0, data, width, height)
if err != nil {
log.Printf("Error reading band %d: %v", i, err)
return err
}
// 计算最小值和最大值
var minVal, maxVal uint16 = 65535, 0
for _, value := range data {
if value < minVal {
minVal = value
}
if value > maxVal {
maxVal = value
}
}
// 将 16 位数据缩放到 8 位并存储到图像的相应通道
for y := 0; y < height; y++ {
for x := 0; x < width; x++ {
value := data[y*width+x]
scaledValue := uint8(float64(value-minVal) / float64(maxVal-minVal) * 255)
img.SetUCharAt(y, x*3+(i-1), scaledValue)
}
}
// 将 16 位数据存储到图像的相应通道
// for y := 0; y < height; y++ {
// for x := 0; x < width; x++ {
// value := data[y*width+x]
// img.SetShortAt(y, x*3+(i-1), int16(value))
// }
// }
}
// 调整图像大小
resizedImg := gocv.NewMat()
gocv.Resize(img, &resizedImg, image.Point{X: 2336, Y: 2336}, 0, 0, gocv.InterpolationCubic)
// 保存为 JPG 文件
ok := gocv.IMWrite(fjpg, resizedImg)
if !ok {
err = fmt.Errorf("error saving %s", fjpg)
return err
}
return nil
}

0
log.go → producer/log.go
View File

70
producer/meta.go Normal file
View File

@@ -0,0 +1,70 @@
package imageproc
import (
"encoding/xml"
"io/ioutil"
)
// 定义与XML结构对应的Go结构体
type ProductMeta struct {
XMLName xml.Name `xml:"ProductMeta"`
Satellite string `xml:"Satellite"`
Sensor string `xml:"Sensor"`
ProductID string `xml:"ProductID"`
ProductLevel string `xml:"ProductLevel"`
OutputFormat string `xml:"OutputFormat"`
ProductGenTime string `xml:"ProductGenTime"`
StartTime string `xml:"StartTime"`
EndTime string `xml:"EndTime"`
CentreTime string `xml:"CentreTime"`
Bands int `xml:"Bands"`
CentreLocation Location `xml:"CentreLocation"`
Corners Corners `xml:"Corners"`
SunAzimuth float64 `xml:"SunAzimuth"`
SunElevation float64 `xml:"SunElevation"`
SatAzimuth float64 `xml:"SatAzimuth"`
SatElevation float64 `xml:"SatElevation"`
Gsd float64 `xml:"Gsd"`
CloudRate float64 `xml:"CloudRate"`
Roll float64 `xml:"Roll"`
Pitch float64 `xml:"Pitch"`
Yaw float64 `xml:"Yaw"`
SatPosX float64 `xml:"SatPosX"`
SatPosY float64 `xml:"SatPosY"`
SatPosZ float64 `xml:"SatPosZ"`
Width int `xml:"Width"`
Height int `xml:"Height"`
DataSize int64 `xml:"DataSize"`
MapProjection string `xml:"MapProjection"`
EarthEllipsoid string `xml:"EarthEllipsoid"`
UtmZone int `xml:"UtmZone"`
GainLevel int `xml:"GainLevel"`
IntegratedLevel int `xml:"IntegratedLevel"`
IntegratedTime float64 `xml:"IntegratedTime"`
}
type Location struct {
Longitude float64 `xml:"Longitude"`
Latitude float64 `xml:"Latitude"`
}
type Corners struct {
UpperLeft Location `xml:"UpperLeft"`
UpperRight Location `xml:"UpperRight"`
LowerRight Location `xml:"LowerRight"`
LowerLeft Location `xml:"LowerLeft"`
}
func WriteProductMeta(productMeta *ProductMeta, filename string) error {
output, err := xml.MarshalIndent(productMeta, "", " ")
if err != nil {
return err
}
err = ioutil.WriteFile(filename, output, 0644)
if err != nil {
return err
}
return nil
}

16
output.go → producer/output.go
View File

@@ -47,6 +47,7 @@ func savePanToGDALGTiff(pan gocv.Mat, tiffFile string) error {
}
band := ds.Bands()[0]
band.SetColorInterp(godal.CIGray)
err = band.IO(godal.IOWrite,
0, 0,
data,
@@ -83,15 +84,21 @@ func (r *Registrator) SaveRegisteredMssToGDALGTiff(tiffFile string) error {
}
}
return SaveBGRToGDALGTiff(r.rgbirImage, 4, 5, tiffFile)
return SaveBGRToGDALGTiff(r.rgbirImage,
4, 5,
[]godal.ColorInterp{godal.CIBlue, godal.CIGreen, godal.CIRed, godal.CIAlpha},
tiffFile)
}
func (r *Registrator) SavePansharpenedToGDALGTiff(tiffFile string) error {
ihsImage := PansharpenIHS(r.PanImage, r.rgbirImage)
return SaveBGRToGDALGTiff(ihsImage, 3, 1.25, tiffFile)
return SaveBGRToGDALGTiff(ihsImage,
3, 1.25,
[]godal.ColorInterp{godal.CIRed, godal.CIGreen, godal.CIBlue},
tiffFile)
}
func SaveBGRToGDALGTiff(bgr gocv.Mat, bands int, resolution float64, tiffFile string) error {
func SaveBGRToGDALGTiff(bgr gocv.Mat, bands int, resolution float64, colorInterps []godal.ColorInterp, tiffFile string) error {
log.Println("Saving BGR to TIFF file:", tiffFile)
width := bgr.Cols()
@@ -123,10 +130,13 @@ func SaveBGRToGDALGTiff(bgr gocv.Mat, bands int, resolution float64, tiffFile st
}
defer ds.Close()
ds.SetMetadata("NBITS", "16")
setGeoTransform(ds, 0, 0, float64(width), float64(height), resolution)
for b := 0; b < bands; b++ {
band := ds.Bands()[b]
band.SetColorInterp(colorInterps[b])
err := band.IO(godal.IOWrite,
0, 0,
data[b],

0
params.go → producer/params.go
View File

0
phase_correlation.go → producer/phase_correlation.go
View File

53
util.go → producer/proj.go
View File

@@ -1,12 +1,6 @@
package imageproc
import (
"fmt"
"time"
"github.com/airbusgeo/godal"
log "github.com/sirupsen/logrus"
)
import "fmt"
func calculateProj(topLeftX, topLeftY, width, height, resolution float64) string {
// 计算图像的地理范围
@@ -61,32 +55,25 @@ func calculateProj(topLeftX, topLeftY, width, height, resolution float64) string
AXIS["Northing",NORTH],
AUTHORITY["EPSG","32651"]]`
projection = `PROJCS["WGS 84 / Pseudo-Mercator",
GEOGCS["WGS 84",DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"]],
PROJECTION["Mercator_1SP"],
PARAMETER["central_meridian",0],
PARAMETER["scale_factor",1],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["metre",1,AUTHORITY["EPSG","9001"]],
AXIS["X",EAST],
AXIS["Y",NORTH],
EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],
AUTHORITY["EPSG","3857"]]`
// 输出投影信息
return projection
}
func setGeoTransform(ds *godal.Dataset, topLeftX, topLeftY, width, height, resolution float64) {
// 设置地理变换（假设左上角坐标为(0,0)，PAN每个像素分辨率为1.2米）
geotransform := [6]float64{
0, // top left x
resolution, // w-e pixel resolution
0, // rotation, 0 if image is "north up"
0, // top left y
0, // rotation, 0 if image is "north up"
-resolution, // n-s pixel resolution (negative value)
}
err := ds.SetGeoTransform(geotransform)
if err != nil {
log.Errorf("Failed to set GeoTransform: %v", err)
}
// 设置投影信息（WGS84）
err = ds.SetProjection(calculateProj(topLeftX, topLeftY, width, height, resolution))
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")
}

28
producer/report.go Normal file
View File

@@ -0,0 +1,28 @@
package imageproc
import "encoding/xml"
// Report represents the root XML element
type Report struct {
XMLName xml.Name `xml:"report"`
Satellite string `xml:"satellite,attr"`
Sensor string `xml:"sensor,attr"`
ProductLevel string `xml:"productLevel,attr"`
Scenes []ReportScene `xml:"scene"`
}
// Scene represents each scene in the report
type ReportScene struct {
Name string `xml:"name,attr"`
TiffData string `xml:"tiffData"`
TfwData string `xml:"tfwData"`
RpbData string `xml:"rpbData"`
BrowserData string `xml:"browserData"`
BrowserRpbData string `xml:"browserRpbData"`
ThumbData string `xml:"thumbData"`
ShpData string `xml:"shpData"`
ShxData string `xml:"shxData"`
DbfData string `xml:"dbfData"`
MetaData string `xml:"metaData"`
QualityData string `xml:"qualityData"`
}

6
scenes.go → producer/scenes.go
View File

@@ -7,6 +7,7 @@ import (
"path/filepath"
"strings"
"github.com/airbusgeo/godal"
log "github.com/sirupsen/logrus"
"gocv.io/x/gocv"
)
@@ -111,7 +112,10 @@ func (r *Registrator) SaveScenesToTiff(panScenes []*Scene, mssScenes []*Scene) e
scene.Tiff = filepath.Join(dir, filename)
rgbirImage, _ := r.MergeMSSToBGRNIR(scene.Mat)
err := SaveBGRToGDALGTiff(rgbirImage, 4, 5, scene.Tiff)
err := SaveBGRToGDALGTiff(rgbirImage,
4, 5,
[]godal.ColorInterp{godal.CIBlue, godal.CIGreen, godal.CIRed, godal.CIAlpha},
scene.Tiff)
if err != nil {
log.Errorf("save mss scene %d to tiff failed: %v", i+1, err)
return err

34
producer/util.go Normal file
View File

@@ -0,0 +1,34 @@
package imageproc
import (
"time"
"github.com/airbusgeo/godal"
log "github.com/sirupsen/logrus"
)
func setGeoTransform(ds *godal.Dataset, topLeftX, topLeftY, width, height, resolution float64) {
// 设置地理变换（假设左上角坐标为(0,0)，PAN每个像素分辨率为1.2米）
geotransform := [6]float64{
0, // top left x
resolution, // w-e pixel resolution
0, // rotation, 0 if image is "north up"
0, // top left y
0, // rotation, 0 if image is "north up"
-resolution, // n-s pixel resolution (negative value)
}
err := ds.SetGeoTransform(geotransform)
if err != nil {
log.Errorf("Failed to set GeoTransform: %v", err)
}
// 设置投影信息（WGS84）
err = ds.SetProjection(calculateProj(topLeftX, topLeftY, width, height, resolution))
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")
}