diff --git a/cmd/proc.go b/cmd/proc.go index a0b0245..6ac2155 100644 --- a/cmd/proc.go +++ b/cmd/proc.go @@ -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") diff --git a/filter.go b/producer/filter.go similarity index 100% rename from filter.go rename to producer/filter.go diff --git a/fit.go b/producer/fit.go similarity index 100% rename from fit.go rename to producer/fit.go diff --git a/gdal_pansharpen.go b/producer/gdal_pansharpen.go similarity index 100% rename from gdal_pansharpen.go rename to producer/gdal_pansharpen.go diff --git a/hdr.go b/producer/hdr.go similarity index 100% rename from hdr.go rename to producer/hdr.go diff --git a/gdal_pansharpen_test.go b/producer/image_proc_test.go similarity index 53% rename from gdal_pansharpen_test.go rename to producer/image_proc_test.go index 92af82f..a7351b5 100644 --- a/gdal_pansharpen_test.go +++ b/producer/image_proc_test.go @@ -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) + } +} diff --git a/image_registration.go b/producer/image_registration.go similarity index 96% rename from image_registration.go rename to producer/image_registration.go index 3c4a4c1..707e28a 100644 --- a/image_registration.go +++ b/producer/image_registration.go @@ -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) } } diff --git a/producer/jpg.go b/producer/jpg.go new file mode 100644 index 0000000..eed6f4a --- /dev/null +++ b/producer/jpg.go @@ -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 +} diff --git a/log.go b/producer/log.go similarity index 100% rename from log.go rename to producer/log.go diff --git a/producer/meta.go b/producer/meta.go new file mode 100644 index 0000000..4808cca --- /dev/null +++ b/producer/meta.go @@ -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 +} diff --git a/output.go b/producer/output.go similarity index 91% rename from output.go rename to producer/output.go index 3778d6c..a9280c3 100644 --- a/output.go +++ b/producer/output.go @@ -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], diff --git a/params.go b/producer/params.go similarity index 100% rename from params.go rename to producer/params.go diff --git a/phase_correlation.go b/producer/phase_correlation.go similarity index 100% rename from phase_correlation.go rename to producer/phase_correlation.go diff --git a/util.go b/producer/proj.go similarity index 63% rename from util.go rename to producer/proj.go index 4340462..0af926f 100644 --- a/util.go +++ b/producer/proj.go @@ -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") -} diff --git a/producer/report.go b/producer/report.go new file mode 100644 index 0000000..76a7dd1 --- /dev/null +++ b/producer/report.go @@ -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"` +} diff --git a/scenes.go b/producer/scenes.go similarity index 95% rename from scenes.go rename to producer/scenes.go index ae3cdaa..acb73cc 100644 --- a/scenes.go +++ b/producer/scenes.go @@ -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 diff --git a/producer/util.go b/producer/util.go new file mode 100644 index 0000000..bacc338 --- /dev/null +++ b/producer/util.go @@ -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") +}