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

使用gamma校正提升jpg亮度

Browse Source
This commit is contained in:
nuknal
2024-05-30 11:22:34 +08:00
parent 7d9ec46750
commit 07ee4d88d4
18 changed files with 174 additions and 98 deletions
Download Patch File Download Diff File Expand all files Collapse all files

78
pkg/producer/filter.go Normal file
View File

@@ -0,0 +1,78 @@
package imageproc
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
}

39
pkg/producer/fit.go Normal file
View File

@@ -0,0 +1,39 @@
package imageproc
import (
"gonum.org/v1/gonum/mat"
)
// PolynomialFit fits a polynomial of given degree to the data points (x, y).
func PolynomialFit(x, y []float64, degree int) ([]float64, error) {
n := len(x)
// Create the Vandermonde matrix
vander := mat.NewDense(n, degree+1, nil)
for i := 0; i < n; i++ {
for j := 0; j <= degree; j++ {
vander.Set(i, j, pow(x[i], j))
}
}
// Create the right-hand side vector
yVec := mat.NewVecDense(n, y)
// Solve the least squares problem
var qr mat.QR
qr.Factorize(vander)
coeffs := mat.NewDense(degree+1, 1, nil)
err := qr.SolveTo(coeffs, false, yVec)
return coeffs.RawMatrix().Data, err
}
// pow is a helper function to calculate power of a number.
func pow(x float64, n int) float64 {
if n == 0 {
return 1
}
res := x
for i := 1; i < n; i++ {
res *= x
}
return res
}

114
pkg/producer/gdal_pansharpen.go Normal file
View File

@@ -0,0 +1,114 @@
package imageproc
import (
"image"
"os"
"os/exec"
"path/filepath"
"strings"
log "github.com/sirupsen/logrus"
"gocv.io/x/gocv"
)
const (
GDALPansharpenCommand = "gdal_pansharpen.py {{.PANGTiff}} {{.MSSGTiff}} {{.FUSGTiff}} -r cubic -of GTiff"
)
// IHS 变换是一种将 RGB 空间转换为亮度、色调和饱和度的变换方法
func PansharpenIHS(panImage, mssImage gocv.Mat) gocv.Mat {
log.Println("start pansharpen IHS")
mss32F := gocv.NewMat()
mssImage.ConvertTo(&mss32F, gocv.MatTypeCV32FC4)
// 调整 MSS 图像的大小到与 PAN 图像相同的大小
log.Println("resize mss to pan size")
mssResized := gocv.NewMat()
defer mssResized.Close()
gocv.Resize(mss32F, &mssResized,
image.Point{X: panImage.Cols(), Y: panImage.Rows()}, 0, 0, gocv.InterpolationCubic)
// 将 MSS 图像从 BGR 转换为 HLS
mssHLS := gocv.NewMat()
defer mssHLS.Close()
gocv.CvtColor(mssResized, &mssHLS, gocv.ColorBGRToHLS)
// 分离 HLS 通道 只有3个通道
// Hue
// Lightness
// Saturation
hlsChannels := gocv.Split(mssHLS)
// 使用 PAN 图像替换亮度分量（Intensity）
// hlsChannels[1] = pan32
// 将 PAN 图像转换为与 HLS 亮度分量相同的范围
panFloat := gocv.NewMat()
defer panFloat.Close()
panImage.ConvertTo(&panFloat, gocv.MatTypeCV32F)
log.Println("normalize pan to 0-255")
gocv.Normalize(panFloat, &panFloat, 0, 255, gocv.NormMinMax)
hlsFloat := gocv.NewMat()
defer hlsFloat.Close()
hlsChannels[1].ConvertTo(&hlsFloat, gocv.MatTypeCV32F)
gocv.Normalize(hlsFloat, &hlsFloat, 0, 255, gocv.NormMinMax)
// 替换亮度分量
hlsChannels[1] = panFloat
// 合并通道
gocv.Merge(hlsChannels, &mssHLS)
// 将图像从 HLS 转换回 BGR
fusedImage := gocv.NewMat()
defer fusedImage.Close()
log.Println("cvt color from HLS to BGR")
gocv.CvtColor(mssHLS, &fusedImage, gocv.ColorHLSToBGR)
// 交换 B 和 R 通道
channels := gocv.Split(fusedImage)
defer channels[0].Close() // B
defer channels[1].Close() // G
defer channels[2].Close() // R
// 交换 B 和 R 通道
swappedChannels := []gocv.Mat{channels[2], channels[1], channels[0]}
gocv.Merge(swappedChannels, &fusedImage)
fused16U := gocv.NewMat()
fusedImage.ConvertTo(&fused16U, gocv.MatTypeCV16UC3)
return fused16U
}
func GDALPansharpen(panTiff, mssTiff string) (string, error) {
fusedTIff := strings.Replace(mssTiff, "MSS", "FUS", -1)
os.MkdirAll(filepath.Dir(fusedTIff), 0755)
// pansharpenCmd := exec.Command("gdal_pansharpen.py", "-w 0.6 -w 0.1 -w 0.3 -w 0 -b 3 -b 2 -b 1", panTiff, mssTiff, fusedTIff, "-r cubic", "-of GTiff")
pansharpenCmd := exec.Command("gdal_pansharpen.py",
"-w", "0.4", "-w", "0.2", "-w", "0.4", "-w", "0",
"-b", "3", "-b", "2", "-b", "1",
"-r", "cubic", "-of", "GTiff",
panTiff,
mssTiff,
fusedTIff)
log.Println("Run Command: ", pansharpenCmd.String())
pansharpenCmd.Stdout = os.Stdout
pansharpenCmd.Stderr = os.Stderr
return fusedTIff, pansharpenCmd.Run()
}
func GDALTranslate(srcTiff, dstJPG string) error {
if dstJPG == "" {
dstJPG = strings.TrimSuffix(srcTiff, ".tiff") + ".jpg"
}
// gdal_translate -of JPEG -outsize 50% 50% -r bilinear input.tif output.jpg
cmd := exec.Command("gdal_translate", "-of", "JPEG",
srcTiff, dstJPG)
log.Println("Run Command: ", cmd.String())
cmd.Stdout = os.Stdout
cmd.Stderr = os.Stderr
return cmd.Run()
}

55
pkg/producer/hdr.go Normal file
View File

@@ -0,0 +1,55 @@
package imageproc
import (
"bytes"
"html"
"os"
"text/template"
)
var HDRTemplate = `ENVI
description = {
File Imported into ENVI.}
samples = {{.Samples}}
lines = {{.Lines}}
bands = {{.Bands}}
header offset = 0
file type = ENVI Standard
data type = 12
interleave = bsq
sensor type = Unknown
byte order = 0
wavelength units = Unknown`
// https://www.nv5geospatialsoftware.com/docs/ENVIHeaderFiles.html
type EnviHdr struct {
Samples int // samples (pixels) each line
Lines int
Bands int
}
func (h *EnviHdr) String() string {
var t *template.Template
t = template.New("template")
t, err := t.Parse(HDRTemplate)
if err != nil {
return ""
}
buf := new(bytes.Buffer)
t.Execute(buf, h)
buf = bytes.NewBufferString(html.UnescapeString(buf.String()))
return buf.String()
}
func (h *EnviHdr) Write(filename string) error {
str := h.String()
f, err := os.Create(filename)
if err != nil {
return err
}
defer f.Close()
_, err = f.WriteString(str)
return err
}

41
pkg/producer/image_proc_test.go Normal file
View File

@@ -0,0 +1,41 @@
package imageproc
import (
"testing"
"github.com/airbusgeo/godal"
)
func TestGDALPansharpen(t *testing.T) {
// err := GDALPansharpen("data/052022/SJY01_PAN_20240520_115428_052022_103_010.tiff",
// "data/052022/SJY01_MSS_20240520_115428_052022_103_010.tiff")
// if err != nil {
// t.Error(err)
// }
GDALTranslate("data/052022/SJY01_MSS_20240520_115428_052022_103_010_FUS.tiff", "")
}
func TestGTiffToJPG(t *testing.T) {
godal.RegisterAll()
tiff := "../data/051622/006/MSS/SJY01_MSS_20240516_101236_051622_096_006.tiff"
jpg := "../data/051622/006/MSS/SJY01_MSS_20240516_101236_051622_096_006.jpg"
err := GTiffToJPG(tiff, jpg, false)
if err != nil {
t.Error(err)
}
tiff = "../data/051622/006/PAN/SJY01_PAN_20240516_101236_051622_096_006.tiff"
jpg = "../data/051622/006/PAN/SJY01_PAN_20240516_101236_051622_096_006.jpg"
err = GTiffToJPG(tiff, jpg, false)
if err != nil {
t.Error(err)
}
tiff = "../data/051622/006/FUS/SJY01_FUS_20240516_101236_051622_096_006.tiff"
jpg = "../data/051622/006/FUS/SJY01_FUS_20240516_101236_051622_096_006.jpg"
err = GTiffToJPG(tiff, jpg, true)
if err != nil {
t.Error(err)
}
}

373
pkg/producer/image_registration.go Normal file
View File

@@ -0,0 +1,373 @@
package imageproc
import (
"fmt"
"image"
"image/color"
"math"
"os"
"sync"
"github.com/airbusgeo/godal"
log "github.com/sirupsen/logrus"
"gocv.io/x/gocv"
)
type Registrate interface{}
const (
MssBands = 4
PixelBytes = 2
PanWidth = 9344 // 像素宽度
MssWidth = 2336
BlockNH = 4
BlockNW = 8
OverlappedBlockLines = 3000 // 重叠块的行数
DownSampled ResampleMethod = "down_sample_pan"
UpSampled ResampleMethod = "up_sample_mss"
)
type ResampleMethod string
type Registrator struct {
Params Params
PanImage gocv.Mat
PanHeight int
PanWidth int
MssImages [4]gocv.Mat
MssHeight int
MssWidth int
shiftMutex sync.Mutex
phaseShifts [4][]PhaseShiftM
deltaXCoeffs [4][]float64 // 图像内畸变（线性变换），捕捉图像在水平方向上引起的X方向的变形
deltaYCoeffs [4][]float64 // 图像内畸变（非线性变换），捕捉图像在水平方向上引起的Y方向的变形
registeredMssImages [4]gocv.Mat // 配准后的MSS图像
rgbirImage gocv.Mat
resampleMethod ResampleMethod
}
func NewRegistrator(rsmethod ResampleMethod) *Registrator {
var r Registrator
r.resampleMethod = rsmethod
return &r
}
func (r *Registrator) LoadPanRaw() error {
data, err := os.ReadFile(r.Params.PanRawFile)
if err != nil {
log.Error("Error reading raw file: ", err)
return err
}
height := len(data) / (PanWidth * PixelBytes)
r.PanImage, err = gocv.NewMatFromBytes(height, PanWidth, gocv.MatTypeCV16U, data)
if err != nil {
log.Error("Error creating Mat from bytes: ", err)
return err
}
r.PanHeight = height
r.PanWidth = PanWidth
godal.RegisterAll()
hDriver, ok := godal.RasterDriver("Gtiff")
if !ok {
panic("Gtiff not found")
}
md := hDriver.Metadatas()
if md["DCAP_CREATE"] == "YES" {
fmt.Printf("Driver GTiff supports Create() method.\n")
}
if md["DCAP_CREATECOPY"] == "YES" {
fmt.Printf("Driver Gtiff supports CreateCopy() method.\n")
}
fmt.Println("Gtiff driver name:", hDriver.LongName(), hDriver.ShortName())
return nil
}
func (r *Registrator) LoadMssRaw() error {
data, err := os.ReadFile(r.Params.MssRawFile)
if err != nil {
log.Error("Error reading raw file: ", err)
return err
}
height := len(data) / (MssWidth * PixelBytes * MssBands)
mssData := make([][]byte, MssBands)
for h := 0; h < height; h++ {
row := data[h*MssWidth*MssBands*PixelBytes : (h+1)*MssWidth*MssBands*PixelBytes]
for i := 0; i < MssBands; i++ {
mssData[i] = append(mssData[i], row[i*MssWidth*PixelBytes:(i+1)*MssWidth*PixelBytes]...)
}
}
for i := 0; i < MssBands; i++ {
r.MssImages[i], err = gocv.NewMatFromBytes(height, MssWidth, gocv.MatTypeCV16U, mssData[i])
if err != nil {
log.Error("Error creating Mat from bytes: ", err)
return err
}
}
r.MssHeight = height
r.MssWidth = MssWidth
return nil
}
func (r *Registrator) DoPhaseCorrelation() error {
switch r.resampleMethod {
case UpSampled:
return r.CalcUpPhaseCorrelation()
default:
return r.CalcDownPhaseCorrelation()
}
}
// 将PAN降采样后计算相位相关的偏移量
func (r *Registrator) CalcDownPhaseCorrelation() error {
// 确保 MSS 高度是 PAN 高度的 1/4
if r.MssHeight*4 != r.PanHeight {
err := fmt.Errorf("MSS height is not 1/4 of PAN height, invalid raw file")
log.Error(err)
return err
}
// 将PAN将采样作为轮廓匹配基准图像
downsampledPanImage := gocv.NewMat()
gocv.Resize(r.PanImage, &downsampledPanImage,
image.Point{X: r.MssWidth, Y: r.MssHeight}, 0, 0, gocv.InterpolationCubic)
log.Println("down sampled PAN images size:", downsampledPanImage.Size())
// 分块高度
blockHeight := r.MssHeight / BlockNH
blockWidth := r.MssWidth / BlockNW
return r.calcPhaseCorrelation(downsampledPanImage, r.MssImages, r.MssHeight, r.MssWidth, blockHeight, blockWidth)
}
// 将MSS升采样采样后计算相位相关的偏移量
func (r *Registrator) CalcUpPhaseCorrelation() error {
// 确保 MSS 高度是 PAN 高度的 1/4
if r.MssHeight*4 != r.PanHeight {
err := fmt.Errorf("MSS height is not 1/4 of PAN height, invalid raw file")
log.Error(err)
return err
}
// 将PAN将采样作为轮廓匹配基准图像
var upsampledMssImages [MssBands]gocv.Mat
for i := 0; i < MssBands; i++ {
upsampledMssImages[i] = gocv.NewMat()
gocv.Resize(r.MssImages[i], &upsampledMssImages[i],
image.Point{X: r.PanWidth, Y: r.PanHeight}, 0, 0, gocv.InterpolationCubic)
}
fmt.Println("up sampled MSS images size:", upsampledMssImages[0].Size())
// 分块高度 - BlockNH, BlockNW % 4 == 0
blockHeight := r.PanHeight / BlockNH
blockWidth := r.PanWidth / BlockNW
log.Infof("blockHeight=%d, blockWidth=%d", blockHeight, blockWidth)
return r.calcPhaseCorrelation(r.PanImage, upsampledMssImages, r.PanHeight, r.PanWidth, blockHeight, blockWidth)
}
func (r *Registrator) calcPhaseCorrelation(panImage gocv.Mat,
mssImages [4]gocv.Mat,
height, width,
blockHeight, blockWidth int) error {
var wg sync.WaitGroup
for bh := 0; bh < BlockNH; bh++ {
for bw := 0; bw < BlockNW; bw++ {
wg.Add(1)
go func(bh, bw int) {
defer wg.Done()
x0 := bw * blockWidth
y0 := bh * blockHeight
x1 := (bw + 1) * blockWidth
y1 := (bh + 1) * blockHeight
y1 += OverlappedBlockLines // Y偏移量过大，需要将重叠块的行数加上，以避免边界影响
if x1 > width || y1 > height {
log.Debugf("Block out of range. x0=%d, y0=%d, x1=%d, y1=%d", x0, y0, x1, y1)
}
if y1 > height {
y1 = height
}
var shiftM PhaseShiftM
shiftM.Block.width = x1 - x0
shiftM.Block.height = y1 - y0
shiftM.Block.coord.X = x0 // 块左上角x坐标
shiftM.Block.coord.Y = y0 // 块左上角y坐标
rect := image.Rect(
x0, y0,
x1, y1,
)
panBlock := panImage.Region(rect)
for band := 0; band < MssBands; band++ {
log.Debug("processing band:", band+1, ",block:", bh, rect)
mssBlock := mssImages[band].Region(rect)
// 处理每个分块
phaseShift, response := r.calculateBlockPhaseShift(panBlock, mssBlock)
shiftM.dx = phaseShift.X
shiftM.dy = phaseShift.Y
shiftM.response = response
r.shiftMutex.Lock()
r.phaseShifts[band] = append(r.phaseShifts[band], shiftM)
r.shiftMutex.Unlock()
mssBlock.Close()
}
panBlock.Close()
}(bh, bw)
}
}
wg.Wait()
for i := 0; i < MssBands; i++ {
for _, shift := range r.phaseShifts[i] {
if shift.response > 0.4 || shift.dx > 8 || shift.dy > 8 {
log.Debugf("Band %d, block %d, dx=%f, dy=%f, response=%f",
i, shift.Block.coord.X, shift.dx, shift.dy, shift.response)
}
}
}
return r.calcDeltaCoeffs()
}
func (r *Registrator) Clean() {
r.PanImage.Close()
for i := 0; i < MssBands; i++ {
r.MssImages[i].Close()
}
for i := 0; i < MssBands; i++ {
r.registeredMssImages[i].Close()
}
r.rgbirImage.Close()
}
func (r *Registrator) calcDeltaCoeffs() error {
// 计算每个通道的delta多项式拟合系数
for i := 0; i < MssBands; i++ {
var cx []float64
var dx []float64
var dy []float64
effectShift := 0
for _, shift := range r.phaseShifts[i] {
if math.IsNaN(float64(shift.dx)) || math.IsNaN(float64(shift.dy)) {
continue
}
// 经验值过滤
if shift.dy < 64.0 {
continue
}
effectShift++
cx = append(cx, float64(shift.Block.coord.X+shift.Block.width/2)) // MSS 块在X方向没有分块
log.Debugf("effective shift value: %v, cx: %v, dx: %v, dy: %v",
effectShift, shift.Block.coord.X, shift.dx, shift.dy)
dx = append(dx, float64(shift.dx))
dy = append(dy, float64(shift.dy))
}
if len(cx) < 3 {
log.Errorf("No effective shift value found for band %d, skip delta coefficients calculation", i+1)
continue
} else {
var err error
if r.deltaXCoeffs[i], err = PolynomialFit(cx, dx, 1); err != nil {
log.Error("Error fitting deltaX coeffs: ", err)
return err
}
if r.deltaYCoeffs[i], err = PolynomialFit(cx, dy, 2); err != nil {
log.Error("Error fitting deltaY coeffs: ", err)
return err
}
}
}
for i := 0; i < MssBands; i++ {
if len(r.deltaXCoeffs[i]) < 2 || len(r.deltaYCoeffs[i]) < 3 {
continue
}
log.Printf("Band %d:\n delta_x = %.6f*x + %.6f, \n delta_y = %.6f*x^2 + %.6f*x + %.6f\n",
i+1,
r.deltaXCoeffs[i][1], r.deltaXCoeffs[i][0],
r.deltaYCoeffs[i][2], r.deltaYCoeffs[i][1], r.deltaYCoeffs[i][0])
}
return nil
}
func (r *Registrator) DoCoRegestration() error {
for band := 0; band < MssBands; band++ {
if len(r.deltaXCoeffs[band]) < 2 || len(r.deltaYCoeffs[band]) < 3 {
log.Error("delta coefficients not calculated, skip co-registration")
r.registeredMssImages[band] = r.MssImages[band].Clone()
continue
}
mapX := gocv.NewMatWithSize(r.MssHeight, r.MssWidth, gocv.MatTypeCV32FC1)
mapY := gocv.NewMatWithSize(r.MssHeight, r.MssWidth, gocv.MatTypeCV32FC1)
for y := 0; y < r.MssHeight; y++ {
for x := 0; x < r.MssWidth; x++ {
var dx, dy float64
if r.resampleMethod == UpSampled {
xx := float64(x * MssBands)
yy := float64(y * MssBands)
dx = (r.deltaXCoeffs[band][1]*float64(xx) + r.deltaXCoeffs[band][0] + xx) / MssBands
dy = (r.deltaYCoeffs[band][2]*float64(xx*xx) + r.deltaYCoeffs[band][1]*float64(xx) + r.deltaYCoeffs[band][0] + yy) / MssBands
} else {
dx = r.deltaXCoeffs[band][1]*float64(x) + r.deltaXCoeffs[band][0] + float64(x)
dy = r.deltaYCoeffs[band][2]*float64(x*x) + r.deltaYCoeffs[band][1]*float64(x) + r.deltaYCoeffs[band][0] + float64(y)
}
// if band+1 == 4 {
// fmt.Println("band:", band+1, "x:", x, "map_x:", mx, "y:", y, "map_y:", my)
// }
// mapX.SetFloatAt(y, x, float32(x)+float32(r.deltaXCoeffs[band][0]))
// mapY.SetFloatAt(y, x, float32(y)+float32(r.deltaYCoeffs[band][0]))
mapX.SetFloatAt(y, x, float32(dx))
mapY.SetFloatAt(y, x, float32(dy))
}
}
log.Println("co-registration for band", band+1)
r.registeredMssImages[band] = gocv.NewMatWithSize(r.MssHeight, r.MssWidth, gocv.MatTypeCV16UC1)
gocv.Remap(r.MssImages[band],
&r.registeredMssImages[band],
&mapX, &mapY,
gocv.InterpolationCubic,
gocv.BorderConstant,
color.RGBA{0, 0, 0, 0})
}
return nil
}

150
pkg/producer/jpg.go Normal file
View File

@@ -0,0 +1,150 @@
package imageproc
import (
"fmt"
"image"
"math"
"github.com/airbusgeo/godal"
log "github.com/sirupsen/logrus"
"gocv.io/x/gocv"
)
func GTiffToJPG(ftiff, fjpg string, reversed bool) 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("creating JPG for TIFF file %s.", ftiff)
// 获取图像大小
width := bands[0].Structure().SizeX
height := bands[0].Structure().SizeY
bandsCnt := len(bands)
// 创建 16 位图像矩阵
var img gocv.Mat
defer img.Close()
if bandsCnt == 1 {
img = gocv.NewMatWithSize(height, width, gocv.MatTypeCV16UC1)
} else {
bandsCnt = 3 // 只取前三个波段
img = gocv.NewMatWithSize(height, width, gocv.MatTypeCV16UC3)
}
// 读取每个波段并存储到图像矩阵中
for i := 0; i < bandsCnt; i++ {
band := bands[i]
data := make([]uint16, width*height)
err = band.Read(0, 0, data, width, height)
if err != nil {
fmt.Printf("Error reading band %d: %v\n", i, err)
return err
}
// 将 16 位数据存储到图像的相应通道
for y := 0; y < height; y++ {
for x := 0; x < width; x++ {
value := data[y*width+x]
img.SetShortAt(y, x*bandsCnt+i, int16(value))
}
}
}
channels := gocv.Split(img)
for i, ch := range channels {
// 2. 计算图像的最小值和最大值
minVal, maxVal, minLoc, maxLoc := gocv.MinMaxLoc(ch)
log.Printf("Min value: %f, Max value: %f, min location: %v, max location: %v", minVal, maxVal, minLoc, maxLoc)
// 3. 将16位图像线性拉伸到8位图像
scale := 255.0 / (maxVal - minVal)
shift := -minVal * scale
ch.ConvertToWithParams(&channels[i], gocv.MatTypeCV8U, scale, shift)
}
img8bit := gocv.NewMat()
defer img8bit.Close()
if reversed {
gocv.Merge([]gocv.Mat{channels[2], channels[1], channels[0]}, &img8bit)
} else {
gocv.Merge(channels, &img8bit)
}
for i := range channels {
channels[i].Close()
}
gocv.Resize(img8bit, &img8bit,
image.Point{X: img8bit.Cols() / 2, Y: img8bit.Rows() / 2},
0, 0, gocv.InterpolationCubic)
// 7. 应用伽玛校正提升亮度
// avgBrightness := calculateAverageBrightness(img8bit)
// gamma := determineGammaValue(avgBrightness)
// log.Printf("Average Brightness: %f, Determined Gamma: %f", avgBrightness, gamma)
gammaCorrected := applyGammaCorrection(img8bit, 1.6) // 伽玛值
defer gammaCorrected.Close()
ok := gocv.IMWriteWithParams(fjpg, gammaCorrected, []int{gocv.IMWriteJpegOptimize, 1})
if !ok {
err = fmt.Errorf("error saving %s", fjpg)
return err
}
// gocv.IMWrite(strings.Replace(fjpg, ".jpg", "_8.jpg", 1), img8bit)
return nil
}
// applyGammaCorrection applies gamma correction to the image
func applyGammaCorrection(src gocv.Mat, gamma float64) gocv.Mat {
lookupTable := make([]uint8, 256)
invGamma := 1.0 / gamma
for i := 0; i < 256; i++ {
lookupTable[i] = uint8(math.Pow(float64(i)/255.0, invGamma) * 255.0)
}
dst := gocv.NewMat()
gammaMat, err := gocv.NewMatFromBytes(1, 256, gocv.MatTypeCV8UC1, lookupTable)
if err != nil {
log.Printf("Error creating gamma correction matrix: %v", err)
return src
}
defer gammaMat.Close()
gocv.LUT(src, gammaMat, &dst)
return dst
}
// calculateAverageBrightness calculates the average brightness of the image
func calculateAverageBrightness(img gocv.Mat) float64 {
if img.Channels() == 1 {
return img.Mean().Val1
}
gray := gocv.NewMat()
defer gray.Close()
gocv.CvtColor(img, &gray, gocv.ColorBGRToGray)
return gray.Mean().Val1
// mean := gocv.NewMat()
// dstStdDev := gocv.NewMat()
// defer mean.Close()
// defer dstStdDev.Close()
// gocv.MeanStdDev(gray, &mean, &dstStdDev)
// return float64(mean.GetFloatAt(0, 0))
}
// determineGammaValue determines an appropriate gamma value based on the average brightness
func determineGammaValue(avgBrightness float64) float64 {
targetBrightness := 180.0
if avgBrightness < 0.1 {
return 1.0
}
gamma := math.Log(targetBrightness/255.0) / math.Log(avgBrightness/255.0)
return gamma
}

70
pkg/producer/meta.go Normal file
View File

@@ -0,0 +1,70 @@
package imageproc
import (
"encoding/xml"
"os"
)
// 定义与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 = os.WriteFile(filename, output, 0644)
if err != nil {
return err
}
return nil
}

196
pkg/producer/output.go Normal file
View File

@@ -0,0 +1,196 @@
package imageproc
import (
"bufio"
"fmt"
"os"
"path/filepath"
"strings"
"github.com/airbusgeo/godal"
log "github.com/sirupsen/logrus"
"gocv.io/x/gocv"
)
func (r *Registrator) SaveOriginalPanToGDALGTiff(tiffFile string) error {
err := savePanToGDALGTiff(r.PanImage, tiffFile)
if err != nil {
return err
}
return nil
}
func (r *Registrator) SaveFilteredPanToGDALGTiff(tiffFile string) error {
img := PANFilter(r.PanImage)
img.ConvertTo(&img, gocv.MatTypeCV16U)
return savePanToGDALGTiff(img, tiffFile)
}
func savePanToGDALGTiff(pan gocv.Mat, tiffFile string) 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, 0, 0, float64(width), float64(height), 1.25)
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,
4, 5,
[]godal.ColorInterp{godal.CIBlue, godal.CIGreen, godal.CIRed, godal.CIUndefined},
tiffFile)
if err != nil {
return err
}
return nil
}
func (r *Registrator) SavePansharpenedToGDALGTiff(tiffFile string) error {
ihsImage := PansharpenIHS(r.PanImage, r.rgbirImage)
return SaveBGRToGDALGTiff(ihsImage,
3, 1.25,
[]godal.ColorInterp{godal.CIRed, godal.CIGreen, godal.CIBlue},
tiffFile)
}
func SaveBGRToGDALGTiff(bgr gocv.Mat, bands int, 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, 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],
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 {
return err
}
w := bufio.NewWriter(f)
w.Write(mssData)
return nil
}
func (r *Registrator) SaveRegisteredMssToRaw(raw string) error {
f, err := os.OpenFile(raw, os.O_CREATE|os.O_TRUNC|os.O_WRONLY, 0777)
if err != nil {
return err
}
var mssData [4][]byte
for i := 0; i < MssBands; i++ {
mssData[i] = r.registeredMssImages[i].ToBytes()
}
width := r.registeredMssImages[0].Cols() * PixelBytes
height := r.registeredMssImages[0].Rows()
log.Println("Writing registered MSS to RAW file...", len(mssData[0])*4)
log.Println("width:", width)
log.Println("height:", height)
w := bufio.NewWriter(f)
for row := 0; row < height; row++ {
w.Write(mssData[0][row*width : (row+1)*width])
w.Write(mssData[1][row*width : (row+1)*width])
w.Write(mssData[2][row*width : (row+1)*width])
w.Write(mssData[3][row*width : (row+1)*width])
}
hdr := EnviHdr{
Samples: 4 * width / PixelBytes,
Lines: height,
Bands: 1,
}
hdrFile := filepath.Join(filepath.Dir(raw),
fmt.Sprintf("%s.HDR", strings.TrimSuffix(filepath.Base(raw), filepath.Ext(raw))),
)
hdr.Write(hdrFile)
return nil
}

14
pkg/producer/params.go Normal file
View File

@@ -0,0 +1,14 @@
package imageproc
type Params struct {
PanRawFile string
MssRawFile string
AuxRawFile string
SaveRegisteredMssRaw bool
PansharpenIHS bool
OutputDir string
PanTiffFile string
MssTiffFile string
FusTIffFile string
SubScenes bool
}

109
pkg/producer/phase_correlation.go Normal file
View File

@@ -0,0 +1,109 @@
package imageproc
import (
"image"
"math"
log "github.com/sirupsen/logrus"
"gocv.io/x/gocv"
)
type PhaseShiftM struct {
dx float32
dy float32
response float64
Block Block
}
type Block struct {
width int
height int
coord image.Point // top-left corner of the block in the original image
}
func (r *Registrator) calculateBlockPhaseShift(panBlock, mssBlock gocv.Mat) (gocv.Point2f, float64) {
pan := gocv.NewMat()
mss := gocv.NewMat()
panBlock.ConvertTo(&pan, gocv.MatTypeCV32FC1)
defer pan.Close()
mssBlock.ConvertTo(&mss, gocv.MatTypeCV32FC1)
defer mss.Close()
hann := gocv.NewMat()
defer hann.Close()
shift, response := gocv.PhaseCorrelate(pan, mss, hann)
dx := shift.X
dy := shift.Y
log.Debugf("Block shift: dx = %f, dy = %f. response = %f", dx, dy, response)
return shift, response
}
func (r *Registrator) DoMssPhaseShift() error {
// 使用平均偏移量来做平移变换
for band := 0; band < MssBands; band++ {
var efficientShiftM int
var xTotal, yTotal float32
for _, shift := range r.phaseShifts[band] {
if math.IsNaN(float64(shift.dx)) || math.IsNaN(float64(shift.dy)) {
continue
}
efficientShiftM += 1
xTotal += shift.dx
yTotal += shift.dy
}
dx := xTotal / float32(efficientShiftM)
dy := yTotal / float32(efficientShiftM)
log.Println("Band", band+1, "average shift:", dx, dy, "efficientShiftM:", efficientShiftM)
translationMat := gocv.NewMatWithSize(2, 3, gocv.MatTypeCV32F)
defer translationMat.Close()
translationMat.SetFloatAt(0, 0, 1)
translationMat.SetFloatAt(0, 1, 0)
translationMat.SetFloatAt(0, 2, dx)
translationMat.SetFloatAt(1, 0, 0)
translationMat.SetFloatAt(1, 1, 1)
translationMat.SetFloatAt(1, 2, dy)
alignedMss := gocv.NewMatWithSize(r.MssHeight, r.MssWidth, gocv.MatTypeCV32FC1)
defer alignedMss.Close()
cvtMss := gocv.NewMat()
defer cvtMss.Close()
r.MssImages[band].ConvertTo(&cvtMss, gocv.MatTypeCV32FC1)
// 手动平移像素
outY := math.MaxInt
for y := 0; y < r.MssHeight; y++ {
for x := 0; x < r.MssWidth; x++ {
// 计算新的坐标
newX := x + int(dx)
newY := y + int(dy)
// 如果新坐标在图像范围内，进行像素值赋值
if newX >= 0 && newX < r.MssWidth && newY >= 0 && newY < r.MssHeight {
alignedMss.SetFloatAt(y, x, cvtMss.GetFloatAt(newY, newX))
} else {
// 如果新坐标不在图像范围内，设置为黑色
alignedMss.SetFloatAt(y, x, 0)
if outY > y {
outY = y
log.Println("Warning: pixel out of range", x, y)
}
}
}
}
// gocv.WarpAffine(cvtMss, &alignedMss, translationMat, image.Pt(cvtMss.Size()[1], cvtMss.Size()[0]))
r.registeredMssImages[band] = gocv.NewMat()
alignedMss.ConvertTo(&r.registeredMssImages[band], gocv.MatTypeCV16U)
log.Println("Band", band+1, "registeredMssImages size:", r.registeredMssImages[band].Size())
}
return nil
}

80
pkg/producer/proj.go Normal file
View File

@@ -0,0 +1,80 @@
package imageproc
import "fmt"
func calculateProj(topLeftX, topLeftY, width, height, resolution float64) string {
var projections []string
// 计算图像的地理范围
bottomRightX := topLeftX + float64(width)*resolution
bottomRightY := topLeftY + float64(height)*resolution
// 根据地理范围和分辨率选择适当的投影信息
// 这里只是一个示例，你需要根据实际情况选择合适的投影系统和参数
projections = append(projections, fmt.Sprintf(`PROJCS["Custom Projection",
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["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",%f],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",%f],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AXIS["Easting",EAST],
AXIS["Northing",NORTH]]`,
(topLeftX+bottomRightX)/2, (topLeftY+bottomRightY)/2))
projections = append(projections, `PROJCS["WGS 84 / UTM zone 51N",
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["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",123],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AXIS["Easting",EAST],
AXIS["Northing",NORTH],
AUTHORITY["EPSG","32651"]]`)
projections = append(projections, `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 projections[0]
}

45
pkg/producer/report.go Normal file
View File

@@ -0,0 +1,45 @@
package imageproc
import (
"encoding/xml"
"os"
)
// 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"`
}
func WriteReport(report *Report, filename string) error {
output, err := xml.MarshalIndent(report, "", " ")
if err != nil {
return err
}
err = os.WriteFile(filename, output, 0644)
if err != nil {
return err
}
return nil
}

163
pkg/producer/scenes.go Normal file
View File

@@ -0,0 +1,163 @@
package imageproc
import (
"fmt"
"image"
"os"
"path/filepath"
"strings"
"github.com/airbusgeo/godal"
log "github.com/sirupsen/logrus"
"gocv.io/x/gocv"
)
type Scene struct {
Width int
Height int
X int // coordinate of the left top corner of the scene
Y int
Mat []gocv.Mat
Tiff string
}
func (s *Scene) Cleanup() {
for i := range s.Mat {
s.Mat[i].Close()
}
}
// 对 PAN 和 配准后的MSS 在 Y 方向进行分景，景之间有25%的重叠
// 默认分景大小：
// MSS 2336 * 2336
// PAN 9344 * 9344
func (r *Registrator) SubScenes() (panScenes []*Scene, mssScenes []*Scene, err error) {
hPAN := (9344 - 2336)
hPANOverlap := 2336
panScenesCnt := r.PanHeight / hPAN
for i := 0; i < panScenesCnt; i++ {
y1 := (i+1)*hPAN + hPANOverlap
if y1 > r.PanHeight {
y1 = r.PanHeight
}
scene := &Scene{
Width: 9344,
Height: y1 - i*hPAN,
X: 0,
Y: i * hPAN,
}
mat := r.PanImage.Region(image.Rect(0, i*hPAN, 9344, y1))
scene.Mat = append(scene.Mat, mat)
panScenes = append(panScenes, scene)
}
hMSS := hPAN / 4
hMSSOverlap := hPANOverlap / 4
mssScenesCnt := r.MssHeight / hMSS
for i := 0; i < mssScenesCnt; i++ {
y1 := (i+1)*hMSS + hMSSOverlap
if y1 > r.MssHeight {
y1 = r.MssHeight
}
scene := &Scene{
Width: 2336,
Height: y1 - i*hMSS,
X: 0,
Y: i * hMSS,
}
for band := 0; band < 4; band++ {
mat := r.registeredMssImages[band].Region(image.Rect(0, i*hMSS, 2336, y1))
scene.Mat = append(scene.Mat, mat)
}
mssScenes = append(mssScenes, scene)
}
if len(panScenes) != len(mssScenes) {
log.Error("pan and mss scenes count not match")
err = fmt.Errorf("pan and mss scenes count not match")
}
return panScenes, mssScenes, err
}
func (r *Registrator) SaveScenesToTiff(panScenes []*Scene, mssScenes []*Scene) error {
for i, scene := range panScenes {
dir := filepath.Join(r.Params.OutputDir, fmt.Sprintf("%03d", i+1), "PAN")
os.MkdirAll(dir, 0755)
tiff := filepath.Base(r.Params.PanTiffFile)
tiff = strings.TrimSuffix(tiff, ".tiff")
filename := fmt.Sprintf("%s_%03d.tiff", tiff, i+1)
scene.Tiff = filepath.Join(dir, filename)
err := savePanToGDALGTiff(scene.Mat[0], scene.Tiff)
if err != nil {
log.Errorf("save pan scene %d to tiff failed: %v", i+1, err)
return err
}
GTiffToJPG(scene.Tiff, strings.Replace(scene.Tiff, ".tiff", ".jpg", 1), false)
}
for i, scene := range mssScenes {
dir := filepath.Join(r.Params.OutputDir, fmt.Sprintf("%03d", i+1), "MSS")
os.MkdirAll(dir, 0755)
tiff := filepath.Base(r.Params.MssTiffFile)
tiff = strings.TrimSuffix(tiff, ".tiff")
filename := fmt.Sprintf("%s_%03d.tiff", tiff, i+1)
scene.Tiff = filepath.Join(dir, filename)
rgbirImage, _ := r.MergeMSSToBGRNIR(scene.Mat)
err := SaveBGRToGDALGTiff(rgbirImage,
4, 5,
[]godal.ColorInterp{godal.CIBlue, godal.CIGreen, godal.CIRed, godal.CIUndefined},
scene.Tiff)
if err != nil {
log.Errorf("save mss scene %d to tiff failed: %v", i+1, err)
return err
}
GTiffToJPG(scene.Tiff, strings.Replace(scene.Tiff, ".tiff", ".jpg", 1), false)
}
return nil
}
func (r *Registrator) DoScenePansharpen(panScenes []*Scene, mssScenes []*Scene) error {
for i := 0; i < len(panScenes); i++ {
// rgbirImg, _ := r.MergeMSSToBGRNIR(mssScenes[i].Mat)
// ihsImage := PansharpenIHS(panScenes[i].Mat[0], rgbirImg)
// path := strings.TrimSuffix(r.Params.MssTiffFile, ".tiff")
// filename := fmt.Sprintf("%s_%03d_FUS.tiff", path, i+1)
// SaveBGRToGDALGTiff(ihsImage, 3, 1.25, filename)
fusedTiff, err := GDALPansharpen(panScenes[i].Tiff, mssScenes[i].Tiff)
if err != nil {
return err
}
GTiffToJPG(fusedTiff, strings.Replace(fusedTiff, ".tiff", ".jpg", 1), true)
}
return nil
}
func (r *Registrator) MergeMSSToBGRNIR(channels []gocv.Mat) (gocv.Mat, error) {
var rgbirImage gocv.Mat
if len(channels) != 4 {
return rgbirImage, fmt.Errorf("mss channels count not match")
}
rgbirImage = gocv.NewMat()
gocv.Merge(channels, &rgbirImage)
log.Printf("merge mss to bgr nir image, channels: %d, height: %d, width: %d",
rgbirImage.Channels(), rgbirImage.Rows(), rgbirImage.Cols())
return rgbirImage, nil
}

34
pkg/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")
}