From b9f61cd26f463313a35d180121cb51cb658a3d2f Mon Sep 17 00:00:00 2001 From: nuknal Date: Sat, 25 May 2024 09:24:19 +0800 Subject: [PATCH] phase shift --- go.mod | 6 +- go.sum | 4 + image_registration.go | 320 ++++++++++++++++++++++++++++++++++++++++++ main.go | 2 +- phase_correlation.go | 127 +++++++++++++++++ registrate.go | 269 ----------------------------------- 6 files changed, 454 insertions(+), 274 deletions(-) create mode 100644 image_registration.go create mode 100644 phase_correlation.go delete mode 100644 registrate.go diff --git a/go.mod b/go.mod index 44cd87a..fe201be 100644 --- a/go.mod +++ b/go.mod @@ -4,10 +4,8 @@ go 1.20 require ( github.com/chai2010/tiff v0.0.0-20211005095045-4ec2aa243943 + github.com/sirupsen/logrus v1.9.3 gocv.io/x/gocv v0.36.1 ) -require ( - github.com/sirupsen/logrus v1.9.3 // indirect - golang.org/x/sys v0.0.0-20220715151400-c0bba94af5f8 // indirect -) +require golang.org/x/sys v0.0.0-20220715151400-c0bba94af5f8 // indirect diff --git a/go.sum b/go.sum index 84170fb..b0adea7 100644 --- a/go.sum +++ b/go.sum @@ -1,15 +1,19 @@ github.com/chai2010/tiff v0.0.0-20211005095045-4ec2aa243943 h1:CjuhVIUiyWQZVY4rmcvm/9R+60e/Wi6LkXyHU38MqXI= github.com/chai2010/tiff v0.0.0-20211005095045-4ec2aa243943/go.mod h1:FhMMqekobM33oGdTfbi65oQ9P7bnQ5/0EDfmleW35RE= github.com/davecgh/go-spew v1.1.0/go.mod h1:J7Y8YcW2NihsgmVo/mv3lAwl/skON4iLHjSsI+c5H38= +github.com/davecgh/go-spew v1.1.1 h1:vj9j/u1bqnvCEfJOwUhtlOARqs3+rkHYY13jYWTU97c= github.com/davecgh/go-spew v1.1.1/go.mod h1:J7Y8YcW2NihsgmVo/mv3lAwl/skON4iLHjSsI+c5H38= +github.com/pmezard/go-difflib v1.0.0 h1:4DBwDE0NGyQoBHbLQYPwSUPoCMWR5BEzIk/f1lZbAQM= github.com/pmezard/go-difflib v1.0.0/go.mod h1:iKH77koFhYxTK1pcRnkKkqfTogsbg7gZNVY4sRDYZ/4= github.com/sirupsen/logrus v1.9.3 h1:dueUQJ1C2q9oE3F7wvmSGAaVtTmUizReu6fjN8uqzbQ= github.com/sirupsen/logrus v1.9.3/go.mod h1:naHLuLoDiP4jHNo9R0sCBMtWGeIprob74mVsIT4qYEQ= github.com/stretchr/objx v0.1.0/go.mod h1:HFkY916IF+rwdDfMAkV7OtwuqBVzrE8GR6GFx+wExME= +github.com/stretchr/testify v1.7.0 h1:nwc3DEeHmmLAfoZucVR881uASk0Mfjw8xYJ99tb5CcY= github.com/stretchr/testify v1.7.0/go.mod h1:6Fq8oRcR53rry900zMqJjRRixrwX3KX962/h/Wwjteg= gocv.io/x/gocv v0.36.1 h1:6XkEaPOk7h/umjy+MXgSEtSeCIgcPJhccUjrJFhjdTY= gocv.io/x/gocv v0.36.1/go.mod h1:lmS802zoQmnNvXETpmGriBqWrENPei2GxYx5KUxJsMA= golang.org/x/sys v0.0.0-20220715151400-c0bba94af5f8 h1:0A+M6Uqn+Eje4kHMK80dtF3JCXC4ykBgQG4Fe06QRhQ= golang.org/x/sys v0.0.0-20220715151400-c0bba94af5f8/go.mod h1:oPkhp1MJrh7nUepCBck5+mAzfO9JrbApNNgaTdGDITg= gopkg.in/check.v1 v0.0.0-20161208181325-20d25e280405/go.mod h1:Co6ibVJAznAaIkqp8huTwlJQCZ016jof/cbN4VW5Yz0= +gopkg.in/yaml.v3 v3.0.0-20200313102051-9f266ea9e77c h1:dUUwHk2QECo/6vqA44rthZ8ie2QXMNeKRTHCNY2nXvo= gopkg.in/yaml.v3 v3.0.0-20200313102051-9f266ea9e77c/go.mod h1:K4uyk7z7BCEPqu6E+C64Yfv1cQ7kz7rIZviUmN+EgEM= diff --git a/image_registration.go b/image_registration.go new file mode 100644 index 0000000..147a306 --- /dev/null +++ b/image_registration.go @@ -0,0 +1,320 @@ +package main + +import ( + "bufio" + "fmt" + "image" + "math" + "os" + + "github.com/chai2010/tiff" + log "github.com/sirupsen/logrus" + "gocv.io/x/gocv" +) + +type Registrate interface{} + +const ( + MssBands = 4 + PixelBytes = 2 + PanWidth = 9344 // 像素宽度 + MssWidth = 2336 + BlockNH = 5 + BlockNW = 10 + DownSampled ResampleMethod = "down_sample_pan" + UpSampled ResampleMethod = "up_sample_mss" +) + +type ResampleMethod string + +type Registrator struct { + PanImage gocv.Mat + PanHeight int + PanWidth int + + MssImages [4]gocv.Mat + MssHeight int + MssWidth int + + phaseShifts [4][]*PhaseShiftM + registeredImages [4]gocv.Mat + rgbirImage gocv.Mat + + resampleMethod ResampleMethod +} + +func NewRegistrator() *Registrator { + var r Registrator + + return &r +} + +func (r *Registrator) LoadPanRaw(raw string) error { + data, err := os.ReadFile(raw) + 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 + + return nil +} + +func (r *Registrator) LoadMssRaw(raw string) error { + data, err := os.ReadFile(raw) + 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 +} + +// 将PAN降采样后计算相位相关的偏移量 +func (r *Registrator) DoDownPhaseCorrelation() 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.InterpolationLinear) + + // 对每个 MSS 波段图像进行上采样 + // upsampledMSSImages := make([]gocv.Mat, MssBands) + + 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.InterpolationLinear) + + // r.msToRaw(upsampledMSSImages[i].ToBytes(), fmt.Sprintf("MSS%d.RAW", i+1)) + } + + fmt.Println("down sampled PAN images size:", downsampledPanImage.Size()) + + // 分块高度 + blockHeight := r.MssHeight / BlockNH + + alignedMssData := make([][]byte, MssBands) + + for band := 0; band < MssBands; band++ { + alignedMSSImage := gocv.NewMatWithSize(r.MssHeight, r.MssWidth, gocv.MatTypeCV16U) + for bh := 0; bh < BlockNH; bh++ { + // 读取 PAN 和 MSS 分块数据 + y1 := (bh + 1) * blockHeight + if bh == BlockNH-1 { + y1 = r.MssHeight + } + + var shiftM PhaseShiftM + shiftM.Block.width = r.MssWidth // 块宽度 + shiftM.Block.height = y1 - bh*blockHeight // 块高度 + shiftM.Block.coord.X = 0 // 块左上角x坐标 + shiftM.Block.coord.Y = bh * blockHeight // 块左上角y坐标 + + rect := image.Rect( + shiftM.Block.coord.X, shiftM.Block.coord.Y, + shiftM.Block.coord.X+shiftM.Block.width, shiftM.Block.coord.Y+shiftM.Block.height, + ) + log.Println("Band", band+1, ", processing block", bh, rect) + + panBlock := downsampledPanImage.Region(rect) + mssBlock := r.MssImages[band].Region(rect) + + // 处理每个分块 + alignedBlock, phaseShift := r.processBlock(panBlock, mssBlock) + shiftM.dx = phaseShift.X + shiftM.dy = phaseShift.Y + + r.phaseShifts[band] = append(r.phaseShifts[band], &shiftM) + + // alignedBlockData := alignedBlock.ToBytes() + // alignedMssData[band] = append(alignedMssData[band], alignedBlockData...) + + // if alignedMSSImage.Empty() { + // alignedMSSImage = alignedBlock.Clone() + // } else { + // gocv.Vconcat(alignedMSSImage, alignedBlock, &alignedMSSImage) + // alignedBlock.Close() + // } + + panBlock.Close() + mssBlock.Close() + alignedBlock.Close() + } + + r.registeredImages[band] = alignedMSSImage + } + + // 使用平均偏移量来做平移变换 + 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) + 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) + cvtMss := gocv.NewMat() + 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])) + alignedMss.ConvertTo(&alignedMss, gocv.MatTypeCV16U) + alignedMssData[band] = append(alignedMssData[band], alignedMss.ToBytes()...) + + translationMat.Close() + cvtMss.Close() + alignedMss.Close() + } + + r.mssToRaw(alignedMssData) + + return nil +} + +func (r *Registrator) panToTiff(panImage gocv.Mat, filePath string) error { + return nil +} + +func (r *Registrator) mssToTiff(registeredImages [4]gocv.Mat, filePath string) error { + // 创建合并后的图像（RGBIR） + rgbirImage := gocv.NewMatWithSize(r.PanHeight, r.PanWidth, gocv.MatTypeCV16UC4) // 4通道，16位 + + for y := 0; y < r.PanHeight; y++ { + for x := 0; x < r.PanWidth; x++ { + r := registeredImages[0].GetShortAt(y, x) + g := registeredImages[1].GetShortAt(y, x) + b := registeredImages[2].GetShortAt(y, x) + ir := registeredImages[3].GetShortAt(y, x) + rgbirImage.SetShortAt(y, x*4+0, r) + rgbirImage.SetShortAt(y, x*4+1, g) + rgbirImage.SetShortAt(y, x*4+2, b) + rgbirImage.SetShortAt(y, x*4+3, ir) + } + } + + // 将合并后的图像保存为TIFF文件 + fileName := "data/registered_rgbir.tiff" + tiffFile, err := os.Create(fileName) + if err != nil { + fmt.Println("Error creating TIFF file:", err) + return err + } + defer tiffFile.Close() + + // 使用tiff库保存图像 + img, err := rgbirImage.ToImage() + if err != nil { + fmt.Println("Error converting Mat to image:", err) + return err + } + if err := tiff.Encode(tiffFile, img, nil); err != nil { + fmt.Println("Error encoding TIFF file:", err) + return err + } + + fmt.Println("Saved", fileName) + return nil +} + +func (r *Registrator) mssToRaw(mssData [][]byte) error { + f, err := os.OpenFile("data/downsampled_registered_mss.RAW", os.O_CREATE|os.O_TRUNC|os.O_WRONLY, 0777) + if err != nil { + return err + } + + log.Println("Writing downsampled registered MSS to RAW file...", len(mssData[0])*4) + log.Println("width:", r.MssWidth*PixelBytes*4) + log.Println("height:", r.MssHeight) + w := bufio.NewWriter(f) + for row := 0; row < r.MssHeight; row++ { + w.Write(mssData[0][row*r.MssWidth*PixelBytes : (row+1)*r.MssWidth*PixelBytes]) + w.Write(mssData[1][row*r.MssWidth*PixelBytes : (row+1)*r.MssWidth*PixelBytes]) + w.Write(mssData[2][row*r.MssWidth*PixelBytes : (row+1)*r.MssWidth*PixelBytes]) + w.Write(mssData[3][row*r.MssWidth*PixelBytes : (row+1)*r.MssWidth*PixelBytes]) + } + + 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 +} diff --git a/main.go b/main.go index d378938..9ca759b 100644 --- a/main.go +++ b/main.go @@ -8,7 +8,7 @@ func main() { if err := reg.LoadPanRaw("/Users/lan/workspace/temp/ccv/pan_image.raw"); err != nil { panic(err) } - if err := reg.DoPhaseCorrelation(); err != nil { + if err := reg.DoDownPhaseCorrelation(); err != nil { panic(err) } } diff --git a/phase_correlation.go b/phase_correlation.go new file mode 100644 index 0000000..b666aa4 --- /dev/null +++ b/phase_correlation.go @@ -0,0 +1,127 @@ +package main + +import ( + "image" + "math" + + log "github.com/sirupsen/logrus" + "gocv.io/x/gocv" +) + +type PhaseShiftM struct { + dx float32 + dy float32 + 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) processBlock(panBlock, mssBlock gocv.Mat) (gocv.Mat, gocv.Point2f) { + 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() + shift, response := gocv.PhaseCorrelate(pan, mss, hann) + hann.Close() + + dx := shift.X + dy := shift.Y + log.Printf("Block shift: dx = %f, dy = %f. response = %f \n", dx, dy, response) + + if math.IsNaN(float64(dx)) { + dx = 0.0 + } + + if math.IsNaN(float64(dy)) { + dy = 0.0 + } + + // cv::Mat translationMatrix = (cv::Mat_(2, 3) << 1, 0, shift.x, 0, 1, shift.y); + translationMat := gocv.NewMatWithSize(2, 3, gocv.MatTypeCV32F) + 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) + defer translationMat.Close() + + alignedBlock := gocv.NewMat() + gocv.WarpAffine(mss, &alignedBlock, translationMat, image.Pt(mss.Size()[1], mss.Size()[0])) + + alignedBlock.ConvertTo(&alignedBlock, gocv.MatTypeCV16U) + + return alignedBlock, shift +} + +func (r *Registrator) DoMssPhaseShift() ([][]byte, error) { + alignedMssData := make([][]byte, MssBands) + // 使用平均偏移量来做平移变换 + 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) + 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) + cvtMss := gocv.NewMat() + 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])) + alignedMss.ConvertTo(&alignedMss, gocv.MatTypeCV16U) + alignedMssData[band] = append(alignedMssData[band], alignedMss.ToBytes()...) + + translationMat.Close() + cvtMss.Close() + alignedMss.Close() + } + + return alignedMssData, nil +} diff --git a/registrate.go b/registrate.go deleted file mode 100644 index ad3c1b5..0000000 --- a/registrate.go +++ /dev/null @@ -1,269 +0,0 @@ -package main - -import ( - "bufio" - "fmt" - "image" - "os" - - "github.com/chai2010/tiff" - log "github.com/sirupsen/logrus" - "gocv.io/x/gocv" -) - -type Registrate interface{} - -const ( - MssBands = 4 - PixelBytes = 2 - PanWidth = 9344 // 像素宽度 - MssWidth = 2336 - BlockNum = 5 -) - -type Registrator struct { - PanImage gocv.Mat - PanHeight int - PanWidth int - - MssImages [4]gocv.Mat - MssHeight int - MssWidth int - registeredImages [4]gocv.Mat - rgbirImage gocv.Mat -} - -func NewRegistrator() *Registrator { - var r Registrator - - return &r -} - -func (r *Registrator) LoadPanRaw(raw string) error { - data, err := os.ReadFile(raw) - 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 - - return nil -} - -func (r *Registrator) LoadMssRaw(raw string) error { - data, err := os.ReadFile(raw) - 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 { - // 确保 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.InterpolationLinear) - - // 对每个 MSS 波段图像进行上采样 - upsampledMSSImages := make([]gocv.Mat, MssBands) - 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.InterpolationLinear) - - // r.msToRaw(upsampledMSSImages[i].ToBytes(), fmt.Sprintf("MSS%d.RAW", i+1)) - } - - fmt.Println("PAN images size:", r.PanImage.Size()) - fmt.Println("upsampled MSS1 images size:", upsampledMSSImages[0].Size()) - fmt.Println("upsampled MSS2 images size:", upsampledMSSImages[0].Size()) - fmt.Println("upsampled MSS3 images size:", upsampledMSSImages[0].Size()) - fmt.Println("upsampled MSS4 images size:", upsampledMSSImages[0].Size()) - - // 分块高度 - blockHeight := r.PanHeight / BlockNum - - alignedMssData := make([][]byte, MssBands) - - for band := 0; band < MssBands; band++ { - var alignedMSSImage gocv.Mat - for block := 0; block < BlockNum; block++ { - // 读取 PAN 和 MSS 分块数据 - y1 := (block + 1) * blockHeight - if block == BlockNum-1 { - y1 = r.PanHeight - } - - log.Println("Processing block", block, image.Rect(0, block*blockHeight, r.PanWidth, y1)) - panBlock := r.PanImage.Region(image.Rect(0, block*blockHeight, r.PanWidth, y1)) - mssBlock := upsampledMSSImages[band].Region(image.Rect(0, block*blockHeight, r.PanWidth, y1)) - - // 处理每个分块 - alignedBlock := r.processBlock(panBlock, mssBlock) - - alignedBlockData := alignedBlock.ToBytes() - alignedMssData[band] = append(alignedMssData[band], alignedBlockData...) - - if alignedMSSImage.Empty() { - alignedMSSImage = alignedBlock.Clone() - } else { - gocv.Vconcat(alignedMSSImage, alignedBlock, &alignedMSSImage) - alignedBlock.Close() - } - - panBlock.Close() - mssBlock.Close() - } - - r.registeredImages[band] = alignedMSSImage - } - - r.mssToRaw(alignedMssData, "mss.raw") - - return nil -} - -func (r *Registrator) processBlock(panBlock, mssBlock gocv.Mat) gocv.Mat { - 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() - shift, response := gocv.PhaseCorrelate(pan, mss, hann) - hann.Close() - - dx := shift.X - dy := shift.Y - log.Printf("Block shift: dx = %f, dy = %f. response = %f \n", dx, dy, response) - - // cv::Mat translationMatrix = (cv::Mat_(2, 3) << 1, 0, shift.x, 0, 1, shift.y); - translationMat := gocv.NewMatWithSize(2, 3, gocv.MatTypeCV32F) - translationMat.SetFloatAt(0, 0, 1) - translationMat.SetFloatAt(0, 1, 0) - translationMat.SetFloatAt(0, 2, shift.X) - translationMat.SetFloatAt(1, 0, 0) - translationMat.SetFloatAt(1, 1, 1) - translationMat.SetFloatAt(1, 2, shift.Y) - - alignedBlock := gocv.NewMat() - gocv.WarpAffine(mss, &alignedBlock, translationMat, image.Pt(panBlock.Size()[0], panBlock.Size()[1])) - - translationMat.Close() - alignedBlock.ConvertTo(&alignedBlock, gocv.MatTypeCV16U) - - return alignedBlock -} - -func (r *Registrator) panToTiff(panImage gocv.Mat, filePath string) error { - return nil -} - -func (r *Registrator) mssToTiff(registeredImages [4]gocv.Mat, filePath string) error { - // 创建合并后的图像（RGBIR） - rgbirImage := gocv.NewMatWithSize(r.PanHeight, r.PanWidth, gocv.MatTypeCV16UC4) // 4通道，16位 - - for y := 0; y < r.PanHeight; y++ { - for x := 0; x < r.PanWidth; x++ { - r := registeredImages[0].GetShortAt(y, x) - g := registeredImages[1].GetShortAt(y, x) - b := registeredImages[2].GetShortAt(y, x) - ir := registeredImages[3].GetShortAt(y, x) - rgbirImage.SetShortAt(y, x*4+0, r) - rgbirImage.SetShortAt(y, x*4+1, g) - rgbirImage.SetShortAt(y, x*4+2, b) - rgbirImage.SetShortAt(y, x*4+3, ir) - } - } - - // 将合并后的图像保存为TIFF文件 - fileName := "registered_rgbir.tiff" - tiffFile, err := os.Create(fileName) - if err != nil { - fmt.Println("Error creating TIFF file:", err) - return err - } - defer tiffFile.Close() - - // 使用tiff库保存图像 - img, err := rgbirImage.ToImage() - if err != nil { - fmt.Println("Error converting Mat to image:", err) - return err - } - if err := tiff.Encode(tiffFile, img, nil); err != nil { - fmt.Println("Error encoding TIFF file:", err) - return err - } - - fmt.Println("Saved", fileName) - return nil -} - -func (r *Registrator) mssToRaw(mssData [][]byte, filePath string) error { - f, err := os.OpenFile("mss.RAW", os.O_CREATE|os.O_TRUNC|os.O_WRONLY, 0777) - if err != nil { - return err - } - w := bufio.NewWriter(f) - for row := 0; row < r.PanHeight; row++ { - w.Write(mssData[0][row*r.PanWidth*PixelBytes : (row+1)*r.PanWidth*PixelBytes]) - w.Write(mssData[1][row*r.PanWidth*PixelBytes : (row+1)*r.PanWidth*PixelBytes]) - w.Write(mssData[2][row*r.PanWidth*PixelBytes : (row+1)*r.PanWidth*PixelBytes]) - w.Write(mssData[3][row*r.PanWidth*PixelBytes : (row+1)*r.PanWidth*PixelBytes]) - } - - return nil -} - -func (r *Registrator) msToRaw(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 -}