长条带生产需要分段滤波

.gitignore

@@ -12,3 +12,4 @@ build/go
 *.tif
 *.tiff
 dem/*
+license/keys

.vscode/settings.json

@@ -0,0 +1,5 @@
+{
+  "gopls": {
+    "build.buildFlags": ["-tags=local"]
+  }
+}

Makefile

@@ -12,6 +12,8 @@ endif
 export COMMIT_CNT := $(shell git rev-list HEAD | wc -l | sed 's/ //g' )
 export BUILD_NUMBER := ${BRANCH}-${COMMIT_CNT}
 
+export PriKeyVAL := $(shell cat license/keys/v1/prikey.pem)
+
 export COMPILE_LDFLAGS='-s -w \
                           -X "main.BuildDate=${DATE}" \
                           -X "main.LatestCommit=${LATEST_COMMIT}" \
@@ -20,19 +22,26 @@ export COMPILE_LDFLAGS='-s -w \
                           -X "main.BuiltOnOs=${BUILT_ON_OS}" \
                           -X "main.Branch=${BRANCH}" \
                           -X "main.CommitCnt=${COMMIT_CNT}" \
-                          -X "main.RuntimeVer=${RUNTIME_VER}" '
+                          -X "main.RuntimeVer=${RUNTIME_VER}" \
+						  -X "starwiz.cn/sjy01/image-proc/license/public.ECDSA_PRIVATE=${PriKeyVAL}" '
 
 out:
-	GOOS=darwin GOARCH=arm64 go build -o bin/sjy01-imgproc-darwin-arm64 -ldflags=${COMPILE_LDFLAGS} cmd/*.go 
+	GOOS=darwin GOARCH=arm64 go build -tags=local -o bin/sjy01-imgproc-darwin-arm64 -ldflags=${COMPILE_LDFLAGS} cmd/*.go 
 
 clean:
 	rm -rf ./bin/*
 
 linux:
-	GOOS=linux GOARCH=amd64 go build -o bin/sjy01-imgproc-v2 -ldflags=${COMPILE_LDFLAGS} cmd/*.go
+	GOOS=linux GOARCH=amd64 go build -tags=local -o bin/sjy01-imgproc-v2 -ldflags=${COMPILE_LDFLAGS} cmd/*.go
+
+linux-license:
+	GOOS=linux GOARCH=amd64 go build -tags=license -o bin/sjy01-imgproc-v2 -ldflags=${COMPILE_LDFLAGS} cmd/*.go
 
 release:
 	docker run --rm -v .:/src -v /Users/lan/workspace/sjy01/build/go:/build/go nuknal/gdal38-cv49-builder sh -c "cd /src && make linux"
 
+release-license:
+	docker run --rm -v .:/src -v /Users/lan/workspace/sjy01/build/go:/build/go nuknal/gdal38-cv49-builder sh -c "cd /src && make linux-license"
+
 iau:
 	wget -O cmd/finals2000A.all https://maia.usno.navy.mil/ser7/finals2000A.all

cmd/main.go

@@ -2,6 +2,8 @@ package main
 
 import (
 	_ "embed"
+
+	"starwiz.cn/sjy01/image-proc/license"
 )
 
 //go:embed finals2000A.all
@@ -11,5 +13,6 @@ var eopData []byte
 var eopp5Line []byte
 
 func main() {
+	license.VerifyLicense()
 	rootCmd.Execute()
 }

cmd/proc.go

@@ -45,20 +45,20 @@ var procCmd = &cobra.Command{
 		calculator.EOP = calculator.NewEOPTable()
 		calculator.EOP.Load(eopData, eopp5Line)
 
-		reg := producer.NewRegistrator(producer.DownSampled)
-		reg.Params = initParams()
+		processor := producer.NewImgProc(producer.DownSampled)
+		processor.Params = initParams()
 
-		if err := reg.LoadAuxData(); err != nil {
+		if err := processor.LoadAuxData(); err != nil {
 			logrus.Fatal(err)
 		}
 
 		// reg.AuxPrint()
 
-		if err := reg.LoadMssRaw(); err != nil {
+		if err := processor.LoadMssRaw(); err != nil {
 			logrus.Fatal(err)
 		}
 
-		if err := reg.LoadPanRaw(); err != nil {
+		if err := processor.LoadPanRaw(); err != nil {
 			logrus.Fatal(err)
 		}
 
@@ -71,14 +71,14 @@ var procCmd = &cobra.Command{
 		os.MkdirAll(params.OutputDir, 0755)
 
 		if doLUTRRC {
-			reg.DoRRCbyLUT(lutDir)
+			processor.DoRRCbyLUT(lutDir)
 		}
 
 		if doMomentMatching {
-			reg.DoMomentMatching()
+			processor.DoMomentMatching()
 		}
 
-		if err := reg.DoPhaseCorrelation(); err != nil {
+		if err := processor.DoPhaseCorrelation(false); err != nil {
 			logrus.Fatal(err)
 		}
 
@@ -87,33 +87,38 @@ var procCmd = &cobra.Command{
 				params.OutputDir,
 				strings.TrimSuffix(filepath.Base(params.MssRawFile), filepath.Ext(params.MssRawFile))+"_registered.RAW",
 			)
-			reg.SaveRegisteredMssToRaw(registerdMSSRAW)
+			processor.SaveRegisteredMssToRaw(registerdMSSRAW)
 		}
 
 		runtime.GC()
-		panScenes, mssScenes, err := reg.SubScenes()
+		panScenes, mssScenes, err := processor.SubScenes()
 		if err != nil {
 			logrus.Error(err)
 		}
 
-		reg.OutputL1A(panScenes, mssScenes)
+		// scenesCnt := mathutil.Min(len(panScenes), len(mssScenes))
+		// for i := 0; i < scenesCnt; i++ {
+		// 	processor.RegisterScenes(panScenes[i], mssScenes[i])
+		// }
+
+		processor.OutputL1A(panScenes, mssScenes)
 		producer.CleanScenes(panScenes)
 		producer.CleanScenes(mssScenes)
 
 		runtime.GC()
 
 		if saveStrip {
-			reg.SaveOriginalPanToGDALGTiff(reg.Params.PanTiffFile)
-			reg.SaveRegisteredMssToGDALGTiff(reg.Params.MssTiffFile)
+			processor.SaveOriginalPanToGDALGTiff(processor.Params.PanTiffFile)
+			processor.SaveRegisteredMssToGDALGTiff(processor.Params.MssTiffFile)
 		}
 
-		if reg.Params.DoPansharpen {
-			reg.DoScenePansharpen(panScenes, mssScenes)
+		if processor.Params.DoPansharpen {
+			processor.DoScenePansharpen(panScenes, mssScenes)
 		}
 
-		reg.Report()
+		processor.Report()
 
-		reg.Clean()
+		processor.Clean()
 	},
 }
 

cmd/test.go

@@ -1,8 +1,6 @@
 package main
 
 import (
-	"strconv"
-
 	"github.com/airbusgeo/godal"
 	"github.com/sirupsen/logrus"
 	"github.com/spf13/cobra"
@@ -29,10 +27,21 @@ var testCmd = &cobra.Command{
 		}
 
 		mv := gocv.Split(*input0_mat)
-		for i := 0; i < len(mv); i++ {
-			out := producer.CV_Sobel(mv[i])
-			utils.SavePanToGDALGTiff(out, 0, 0, "data/temp/out_"+strconv.Itoa(i)+".tif", 1.3)
+		// for i := 0; i < len(mv); i++ {
+		// 	out := producer.CV_Sobel(mv[i])
+		// 	utils.SavePanToGDALGTiff(out, 0, 0, "data/temp/out_"+strconv.Itoa(i)+".tif", 1.3)
+		// }
+		index := len(mv) - 1
+		for i := index - 1; i >= 0; i-- {
+			mv[i] = producer.CV_ECCAlign(mv[index], mv[i])
 		}
+
+		out := gocv.NewMat()
+		mv[index].ConvertTo(&mv[index], gocv.MatTypeCV16U)
+		gocv.Merge(mv, &out)
+		utils.SaveBGRToGDALGTiff(out, 4, 0, 0, 5.2,
+			[]godal.ColorInterp{godal.CIBlue, godal.CIGreen, godal.CIRed, godal.CIUndefined},
+			"data/temp/ecc.tif")
 	},
 }
 

config/config.yaml

@@ -12,8 +12,8 @@ coregistration:
   fus_band_order: "RGB"
 
 radiation:
-  pan_remove_hf_noise: false
-  mss_remove_hf_noise: false
+  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

go.sum

@@ -766,6 +766,8 @@ cloud.google.com/go/workflows v1.9.0/go.mod h1:ZGkj1aFIOd9c8Gerkjjq7OW7I5+l6cSvT
 cloud.google.com/go/workflows v1.10.0/go.mod h1:fZ8LmRmZQWacon9UCX1r/g/DfAXx5VcPALq2CxzdePw=
 cloud.google.com/go/workflows v1.11.1/go.mod h1:Z+t10G1wF7h8LgdY/EmRcQY8ptBD/nvofaL6FqlET6g=
 dmitri.shuralyov.com/gpu/mtl v0.0.0-20190408044501-666a987793e9/go.mod h1:H6x//7gZCb22OMCxBHrMx7a5I7Hp++hsVxbQ4BYO7hU=
+e.coding.net/zkxrsz/starwiz/starwiz-license.git v0.1.3 h1:jLW8d/W7t1i2GZshjDTRtIB3CsB8hMQZPPCRS1rR2Zo=
+e.coding.net/zkxrsz/starwiz/starwiz-license.git v0.1.3/go.mod h1:6lZBmLJggC/wplylnH7Zfe4TAvwby+uzdRGVorIktKo=
 gioui.org v0.0.0-20210308172011-57750fc8a0a6/go.mod h1:RSH6KIUZ0p2xy5zHDxgAM4zumjgTw83q2ge/PI+yyw8=
 git.sr.ht/~sbinet/cmpimg v0.1.0 h1:E0zPRk2muWuCqSKSVZIWsgtU9pjsw3eKHi8VmQeScxo=
 git.sr.ht/~sbinet/gg v0.3.1/go.mod h1:KGYtlADtqsqANL9ueOFkWymvzUvLMQllU5Ixo+8v3pc=

license/no_license.go

@@ -0,0 +1,10 @@
+//go:build local
+// +build local
+
+package license
+
+import "fmt"
+
+func VerifyLicense() {
+	fmt.Println("license is not required")
+}

license/public/algorithm.go

@@ -0,0 +1,49 @@
+package public
+
+import (
+	"crypto/ecdsa"
+	"encoding/pem"
+)
+
+type NonEquAlgorthm struct {
+	Algorithm  *SigningMethodECDSA
+	PrivateKey *ecdsa.PrivateKey
+	PublicKey  *ecdsa.PublicKey
+}
+
+func GetNonEquAlgorthm(prikey []byte, pubkey []byte) (*NonEquAlgorthm, error) {
+	var err error
+	alg := &NonEquAlgorthm{
+		Algorithm: GetSignVerifyMgr(),
+	}
+
+	if prikey != nil {
+		if block, _ := pem.Decode(prikey); block == nil {
+			prikey, err = FormatPem(prikey, true)
+			if err != nil {
+				return nil, err
+			}
+		}
+
+		alg.PrivateKey, err = ParseECPrivateKeyFromPEM(prikey)
+		if err != nil {
+			return nil, err
+		}
+	}
+
+	if pubkey != nil {
+		if block, _ := pem.Decode(pubkey); block == nil {
+			pubkey, err = FormatPem(pubkey, false)
+			if err != nil {
+				return nil, err
+			}
+		}
+
+		alg.PublicKey, err = ParseECPublicKeyFromPEM(pubkey)
+		if err != nil {
+			return nil, err
+		}
+	}
+
+	return alg, nil
+}

license/public/ecdsa.go

@@ -0,0 +1,149 @@
+package public
+
+import (
+	"crypto"
+	"crypto/ecdsa"
+	"crypto/rand"
+	"encoding/base64"
+	"errors"
+	"fmt"
+	"math/big"
+	"strings"
+)
+
+var (
+	ErrInvalidKey      = errors.New("key is invalid")
+	ErrInvalidKeyType  = errors.New("key is of invalid type")
+	ErrHashUnavailable = errors.New("the requested hash function is unavailable")
+)
+
+var (
+	// Sadly this is missing from crypto/ecdsa compared to crypto/rsa
+	ErrECDSAVerification = errors.New("crypto/ecdsa: verification error")
+)
+
+func GetSignVerifyMgr() *SigningMethodECDSA {
+	return &SigningMethodECDSA{"ES512", crypto.SHA512, 66, 521}
+}
+
+// Implements the ECDSA family of signing methods signing methods
+// Expects *ecdsa.PrivateKey for signing and *ecdsa.PublicKey for verification
+type SigningMethodECDSA struct {
+	Name      string
+	Hash      crypto.Hash
+	KeySize   int
+	CurveBits int
+}
+
+func (m *SigningMethodECDSA) Alg() string {
+	return m.Name
+}
+
+// Implements the Verify method from SigningMethod
+// For this verify method, key must be an ecdsa.PublicKey struct
+func (m *SigningMethodECDSA) Verify(signingString, signature string, key interface{}) error {
+	var err error
+
+	// Decode the signature
+	var sig []byte
+	if sig, err = DecodeSegment(signature); err != nil {
+		return err
+	}
+
+	// Get the key
+	var ecdsaKey *ecdsa.PublicKey
+	switch k := key.(type) {
+	case *ecdsa.PublicKey:
+		ecdsaKey = k
+	default:
+		return ErrInvalidKeyType
+	}
+
+	if len(sig) != 2*m.KeySize {
+		return ErrECDSAVerification
+	}
+
+	r := big.NewInt(0).SetBytes(sig[:m.KeySize])
+	s := big.NewInt(0).SetBytes(sig[m.KeySize:])
+
+	// Create hasher
+	if !m.Hash.Available() {
+		return ErrHashUnavailable
+	}
+	hasher := m.Hash.New()
+	hasher.Write([]byte(signingString))
+
+	// Verify the signature
+	if verifystatus := ecdsa.Verify(ecdsaKey, hasher.Sum(nil), r, s); verifystatus == true {
+		return nil
+	} else {
+		return ErrECDSAVerification
+	}
+}
+
+// Implements the Sign method from SigningMethod
+// For this signing method, key must be an ecdsa.PrivateKey struct
+func (m *SigningMethodECDSA) Sign(signingString string, key interface{}) (string, error) {
+	// Get the key
+	var ecdsaKey *ecdsa.PrivateKey
+	switch k := key.(type) {
+	case *ecdsa.PrivateKey:
+		ecdsaKey = k
+	default:
+		return "", ErrInvalidKeyType
+	}
+
+	// Create the hasher
+	if !m.Hash.Available() {
+		return "", ErrHashUnavailable
+	}
+
+	hasher := m.Hash.New()
+	hasher.Write([]byte(signingString))
+
+	// Sign the string and return r, s
+	if r, s, err := ecdsa.Sign(rand.Reader, ecdsaKey, hasher.Sum(nil)); err == nil {
+		curveBits := ecdsaKey.Curve.Params().BitSize
+
+		if m.CurveBits != curveBits {
+			fmt.Println(m.CurveBits, curveBits)
+			return "", ErrInvalidKey
+		}
+
+		keyBytes := curveBits / 8
+		if curveBits%8 > 0 {
+			keyBytes += 1
+		}
+
+		// We serialize the outpus (r and s) into big-endian byte arrays and pad
+		// them with zeros on the left to make sure the sizes work out. Both arrays
+		// must be keyBytes long, and the output must be 2*keyBytes long.
+		rBytes := r.Bytes()
+		rBytesPadded := make([]byte, keyBytes)
+		copy(rBytesPadded[keyBytes-len(rBytes):], rBytes)
+
+		sBytes := s.Bytes()
+		sBytesPadded := make([]byte, keyBytes)
+		copy(sBytesPadded[keyBytes-len(sBytes):], sBytes)
+
+		out := append(rBytesPadded, sBytesPadded...)
+
+		return EncodeSegment(out), nil
+	} else {
+		return "", err
+	}
+}
+
+// Encode JWT specific base64url encoding with padding stripped
+func EncodeSegment(seg []byte) string {
+	return strings.TrimRight(base64.URLEncoding.EncodeToString(seg), "=")
+}
+
+// Decode JWT specific base64url encoding with padding stripped
+func DecodeSegment(seg string) ([]byte, error) {
+	if l := len(seg) % 4; l > 0 {
+		seg += strings.Repeat("=", 4-l)
+	}
+
+	return base64.URLEncoding.DecodeString(seg)
+}

license/public/errors.go

@@ -0,0 +1,62 @@
+package public
+
+import "fmt"
+
+import "errors"
+
+func New(code int, msg string) *ErrorMsg {
+	return &ErrorMsg{code: code, msg: msg}
+}
+
+type ErrorMsg struct {
+	code int
+	msg  string
+	err  error
+}
+
+func (e *ErrorMsg) GetCode() int {
+	return e.code
+}
+
+func (e *ErrorMsg) GetMsg() string {
+	return e.msg
+}
+
+func (e *ErrorMsg) GetErr() error {
+	return e.err
+}
+
+func (e *ErrorMsg) SetErr(err error) *ErrorMsg {
+	e.err = err
+	return e
+}
+
+func (e *ErrorMsg) SetErrText(text string) *ErrorMsg {
+	e.err = errors.New(text)
+	return e
+}
+
+func (e *ErrorMsg) Error() string {
+	if e.err != nil {
+		return fmt.Sprintf("code:%d, msg:%s, err:%s", e.code, e.msg, e.err.Error())
+	}
+	return fmt.Sprintf("code:%d, msg:%s", e.code, e.msg)
+}
+
+var (
+	ErrNoCreateObj      = New(0, "uninitialized object")
+	ErrUnKnown          = New(-1, "unknown error")
+	ErrDirNoExist       = New(-2, "dir does not exist")
+	ErrNewWatcher       = New(-3, "new watcher object failed")
+	ErrWatcherAdd       = New(-4, "watcher add dir failed")
+	ErrLoadPubKey       = New(-5, "failed to load public key")
+	ErrReadAuthFile     = New(-6, "reading authorization file failed")
+	ErrDecodeAuthFile   = New(-7, "decode authorization file failed")
+	ErrVerifySign       = New(-8, "failed to verify signature")
+	ErrUnmarshalLiObj   = New(-9, "unmarshal license object failed")
+	ErrGetMachineCode   = New(-10, "failed to get machine code")
+	ErrNoMatchProName   = New(-11, "product name does not match")
+	ErrLicenseExpired   = New(-12, "license is expired")
+	ErrBeforeIssued     = New(-13, "license used before issued")
+	ErrNoMatchMachineID = New(-14, "machine id does not match")
+)

license/public/file.go

@@ -0,0 +1,120 @@
+package public
+
+import (
+	"bufio"
+	"encoding/pem"
+	"fmt"
+	"io/ioutil"
+	"os"
+	"path/filepath"
+)
+
+func Exists(path string) bool {
+	_, err := os.Stat(path) //os.Stat获取文件信息
+	if err != nil {
+		if os.IsExist(err) {
+			return true
+		}
+		return false
+	}
+	return true
+}
+
+func SaveLicensePem(dir string, filename string, licenseString string, headers map[string]string) error {
+	isExist := Exists(dir)
+	if isExist == false && dir != "." {
+		err := os.MkdirAll(dir, os.ModePerm)
+		if err != nil {
+			return err
+		}
+	}
+
+	savePath := filepath.Join(dir, filename)
+
+	fd, err := os.OpenFile(savePath, os.O_RDWR|os.O_CREATE|os.O_TRUNC, 0644)
+	if err != nil {
+		return err
+	}
+	defer fd.Close()
+
+	w := bufio.NewWriter(fd)
+	block := &pem.Block{
+		Type:    "LICENSE",
+		Headers: headers,
+		Bytes:   []byte(licenseString),
+	}
+
+	if err := pem.Encode(w, block); err != nil {
+		return err
+	}
+
+	if err := w.Flush(); err != nil {
+		return err
+	}
+	return nil
+}
+
+func CheckPemAndSave(filePath string, licenseBytes []byte) (*pem.Block, error) {
+	var (
+		block *pem.Block
+	)
+
+	saveDir := filepath.Dir(filePath)
+	fmt.Println(saveDir)
+
+	isExist := Exists(saveDir)
+	if isExist == false {
+		err := os.MkdirAll(saveDir, os.ModePerm)
+		if err != nil {
+			return nil, err
+		}
+	}
+
+	if block, _ = pem.Decode(licenseBytes); block == nil {
+		return nil, fmt.Errorf("%s", "license must be PEM")
+	}
+
+	fd, err := os.OpenFile(filePath, os.O_RDWR|os.O_CREATE|os.O_TRUNC, 0644)
+	if err != nil {
+		return nil, err
+	}
+	defer fd.Close()
+
+	_, err = fd.Write(licenseBytes)
+	if err != nil {
+		return nil, err
+	}
+
+	return block, nil
+}
+
+func ReadLicensePem(filePath string) ([]byte, error) {
+	isExist := Exists(filePath)
+	if isExist == false {
+		return nil, fmt.Errorf("%s", "license file does not exist")
+	}
+
+	fd, err := os.OpenFile(filePath, os.O_RDONLY, 0644)
+	if err != nil {
+		return nil, err
+	}
+	defer fd.Close()
+
+	licenseBytes, err := ioutil.ReadAll(fd)
+	if err != nil {
+		return nil, err
+	}
+
+	// block, _ := pem.Decode(licenseBytes)
+
+	return licenseBytes, nil
+}
+
+func LoadKey(pemFile string) ([]byte, error) {
+	context, err := ioutil.ReadFile(pemFile)
+	if err != nil {
+		return nil, fmt.Errorf("failed to read key file: %s, error: %s", pemFile, err)
+	}
+
+	return context, nil
+}

license/public/key.go

@@ -0,0 +1,68 @@
+package public
+
+import (
+	"crypto/ecdsa"
+	"crypto/x509"
+	"encoding/pem"
+	"errors"
+)
+
+var (
+	ErrKeyMustBePEMEncoded = errors.New("Invalid Key: Key must be PEM encoded PKCS1 or PKCS8 private key")
+	ErrNotECPublicKey      = errors.New("Key is not a valid ECDSA public key")
+	ErrNotECPrivateKey     = errors.New("Key is not a valid ECDSA private key")
+)
+
+// Parse PEM encoded Elliptic Curve Private Key Structure
+func ParseECPrivateKeyFromPEM(key []byte) (*ecdsa.PrivateKey, error) {
+	var err error
+
+	// Parse PEM block
+	var block *pem.Block
+	if block, _ = pem.Decode(key); block == nil {
+		return nil, ErrKeyMustBePEMEncoded
+	}
+
+	// Parse the key
+	var parsedKey interface{}
+	if parsedKey, err = x509.ParseECPrivateKey(block.Bytes); err != nil {
+		return nil, err
+	}
+
+	var pkey *ecdsa.PrivateKey
+	var ok bool
+	if pkey, ok = parsedKey.(*ecdsa.PrivateKey); !ok {
+		return nil, ErrNotECPrivateKey
+	}
+
+	return pkey, nil
+}
+
+// Parse PEM encoded PKCS1 or PKCS8 public key
+func ParseECPublicKeyFromPEM(key []byte) (*ecdsa.PublicKey, error) {
+	var err error
+
+	// Parse PEM block
+	var block *pem.Block
+	if block, _ = pem.Decode(key); block == nil {
+		return nil, ErrKeyMustBePEMEncoded
+	}
+
+	// Parse the key
+	var parsedKey interface{}
+	if parsedKey, err = x509.ParsePKIXPublicKey(block.Bytes); err != nil {
+		if cert, err := x509.ParseCertificate(block.Bytes); err == nil {
+			parsedKey = cert.PublicKey
+		} else {
+			return nil, err
+		}
+	}
+
+	var pkey *ecdsa.PublicKey
+	var ok bool
+	if pkey, ok = parsedKey.(*ecdsa.PublicKey); !ok {
+		return nil, ErrNotECPublicKey
+	}
+
+	return pkey, nil
+}

license/public/license.go

@@ -0,0 +1,167 @@
+package public
+
+import (
+	"encoding/json"
+	"fmt"
+	"reflect"
+	"strings"
+	"time"
+)
+
+type License struct {
+	LicenseUUID string            `json:"licensever,omitempty"`  //license 唯一编号
+	ProductName string            `json:"productname,omitempty"` //产品名称
+	MachineID   string            `json:"machineid,omitempty"`   //机器ID
+	ExpiresAt   int64             `json:"expiresat,omitempty"`   //过期时间
+	IssuedAt    int64             `json:"issuedat,omitempty"`    //签发时间
+	CustomKV    map[string]string `json:"customkv,omitempty"`
+}
+
+func GenerateLicense(uuid string, productName string, machineID string, expires int64, kv map[string]string) *License {
+	return &License{
+		LicenseUUID: uuid,
+		ProductName: productName,
+		MachineID:   machineID,
+		ExpiresAt:   expires,
+		IssuedAt:    time.Now().Unix(),
+		CustomKV:    kv,
+	}
+}
+
+func VerifyLicense(productName string, machine string, licenseBytes []byte, isVerify bool) (*License, error) {
+	l := new(License)
+	if err := json.Unmarshal(licenseBytes, l); err != nil {
+		return nil, err
+	}
+
+	if isVerify {
+		err := l.Valid(productName, machine)
+		if err != nil {
+			return nil, err
+		}
+	}
+
+	return l, nil
+}
+
+func (c *License) Valid(productName string, machine string) error {
+	var vErr error
+	now := time.Now().Unix()
+
+	//比较产品名称
+	if c.CompareProductName(productName) == false {
+		vErr = ErrNoMatchProName
+	}
+
+	//比较过期时间
+	if c.VerifyExpiresAt(now, false) == false {
+		// delta := time.Unix(now, 0).Sub(time.Unix(c.ExpiresAt, 0))
+		expStr := time.Unix(c.ExpiresAt, 0).Format("2006-01-02 15:04:05")
+		vErr = ErrLicenseExpired.SetErrText(fmt.Sprintf("license is expired after %s", expStr))
+	}
+
+	//比较签发时间
+	if c.VerifyIssuedAt(now, false) == false {
+		vErr = ErrBeforeIssued
+	}
+
+	//比较机器是否与license匹配
+	if c.CompareMachine(machine) == false {
+		vErr = ErrNoMatchMachineID
+	}
+	return vErr
+}
+
+//已经过期返回0,未过期返回剩余的秒数
+func (c *License) GetExpiresAt() int64 {
+	now := time.Now().Unix()
+	delta := time.Unix(c.ExpiresAt, 0).Sub(time.Unix(now, 0))
+	if delta <= 0 { //
+		return 0
+	}
+	return int64(delta.Seconds())
+}
+
+func (c *License) GetEndTime() string {
+	return time.Unix(c.ExpiresAt, 0).String()
+}
+
+//比较产品名
+func (c *License) CompareProductName(productName string) bool {
+	return strings.Compare(c.ProductName, productName) == 0
+}
+
+//比较过期时间
+func (c *License) VerifyExpiresAt(now int64, req bool) bool {
+	if c.ExpiresAt == 0 {
+		return !req
+	}
+	return now <= c.ExpiresAt
+}
+
+//比较签发时间
+func (c *License) VerifyIssuedAt(now int64, req bool) bool {
+	if c.IssuedAt == 0 {
+		return !req
+	}
+	return now >= c.IssuedAt
+}
+
+//比较机器ID
+func (c *License) CompareMachine(machineID string) bool {
+	return strings.Compare(c.MachineID, machineID) == 0
+}
+
+//编码
+func (c *License) ToBytes() ([]byte, error) {
+	var (
+		jsonValue []byte
+		err       error
+	)
+	if jsonValue, err = json.Marshal(c); err != nil {
+		return nil, err
+	}
+
+	return jsonValue, nil
+}
+
+func ToLicense(lb []byte) (*License, error) {
+	l := new(License)
+	err := json.Unmarshal(lb, l)
+	if err != nil {
+		return nil, err
+	}
+
+	return l, nil
+}
+
+func BytesToLicense(license string) (*License, error) {
+	parts := strings.Split(license, ".")
+	if len(parts) != 2 {
+		return nil, fmt.Errorf("%s", "license contains an invalid number of segments")
+	}
+
+	plainBytes, err := DecodeSegment(parts[0])
+	if err != nil {
+		return nil, err
+	}
+
+	l := new(License)
+	err = json.Unmarshal(plainBytes, l)
+	if err != nil {
+		return nil, err
+	}
+
+	return l, nil
+}
+
+func Struct2Map(obj interface{}) map[string]interface{} {
+	t := reflect.TypeOf(obj)
+	v := reflect.ValueOf(obj)
+
+	var data = make(map[string]interface{})
+	for i := 0; i < t.NumField(); i++ {
+		data[t.Field(i).Name] = v.Field(i).Interface()
+	}
+	return data
+}

license/public/machine.go

@@ -0,0 +1,203 @@
+package public
+
+import (
+	"crypto/md5"
+	"encoding/hex"
+	"encoding/json"
+	"fmt"
+	"io/ioutil"
+	"os"
+	"os/exec"
+	"runtime"
+	"strconv"
+	"strings"
+)
+
+const (
+	DMIDECODE = "dmidecode"
+)
+
+func CheckCmdExists(command string) (string, error) {
+	path, err := exec.LookPath(command)
+	if err != nil {
+		if runtime.GOOS != "darwin" || //macos,windows测试使用,不获取硬盘分区UUID
+			runtime.GOOS != "windows" {
+			fmt.Printf("didn't find 'blkid' executable,err %s\n", err.Error())
+		}
+		return "", err
+	}
+	return path, nil
+}
+
+func IsExistFile(path string) bool {
+	_, err := os.Stat(path) //os.Stat获取文件信息
+	if err != nil {
+		if os.IsExist(err) {
+			return true
+		}
+		return false
+	}
+	return true
+}
+
+func ReadSysFile(filePath string) ([]byte, error) {
+	isExist := IsExistFile(filePath)
+	if isExist == false {
+		return nil, fmt.Errorf("%s file does not exist", filePath)
+	}
+
+	fd, err := os.OpenFile(filePath, os.O_RDONLY, 0644)
+	if err != nil {
+		if os.IsPermission(err) {
+			return nil, fmt.Errorf("%s", "permission denied get machine id")
+		}
+		return nil, fmt.Errorf("%s", "get machine id failed")
+	}
+	defer fd.Close()
+
+	fstabText, err := ioutil.ReadAll(fd)
+	if err != nil {
+		return nil, err
+	}
+
+	return fstabText, nil
+}
+
+//降低数字串长度
+func Sum(data []byte) string {
+	var (
+		sum    uint64
+		length int = len(data)
+		index  int
+	)
+
+	//以32位求和
+	for length >= 4 {
+		sum += uint64(data[index])<<24 + uint64(data[index+1])<<16 + uint64(data[index+2])<<8 + uint64(data[index+3])
+		index += 4
+		length -= 4
+	}
+
+	switch length {
+	case 3:
+		sum += uint64(data[index])<<16 + uint64(data[index+1])<<8 + uint64(data[index+2])
+	case 2:
+		sum += uint64(data[index])<<8 + uint64(data[index+1])
+	case 1:
+		sum += uint64(data[index])
+	case 0:
+		break
+	}
+
+	return strconv.FormatUint(sum, 16)
+}
+
+//获取BIOS出厂UUID
+func GetProductUUID() (string, error) {
+	if runtime.GOOS == "darwin" {
+		return GetProductUUIDMac()
+	}
+
+	if runtime.GOOS == "windows" {
+		return GetProductUUIDWin()
+	}
+
+	cmd := exec.Command("dmidecode", "-s", "system-uuid")
+	output, err := cmd.CombinedOutput()
+	if err != nil {
+		return "", err
+	}
+
+	if len(output) == 0 {
+		return "", fmt.Errorf("%s", "uuid is empty string")
+	}
+
+	system_uuid := string(output)
+
+	cmd = exec.Command("dmidecode", "-s", "system-serial-number")
+	output, err = cmd.CombinedOutput()
+	if err != nil {
+		return "", err
+	}
+
+	if len(output) == 0 {
+		return "", fmt.Errorf("%s", "serial-number is empty string")
+	}
+
+	system_serial_number := string(output)
+
+	deviceID := system_uuid + "_starwiz_" + system_serial_number
+	return GetMd5String(deviceID), nil
+}
+
+/*
+GetMachineID 获取机器ID (分区UUID和网卡地址 or Product UUID)
+*/
+func GetMachineID() (string, error) {
+	var (
+		MachineIDList []string
+	)
+
+	productUUID, err := GetProductUUID()
+	if err != nil {
+		return "", err
+	}
+	MachineIDList = append(MachineIDList, productUUID)
+
+	//编码
+	encByte, err := json.Marshal(MachineIDList)
+	if err != nil {
+		return "", err
+	}
+
+	//降低长度
+	return Sum(encByte), nil
+}
+
+func GetMd5String(s string) string {
+	h := md5.New()
+	h.Write([]byte(s))
+	return hex.EncodeToString(h.Sum(nil))
+}
+
+func GetProductUUIDMac() (string, error) {
+	// 从系统调用获取硬件 UUID
+	out, err := exec.Command("system_profiler", "SPHardwareDataType").Output()
+	if err != nil {
+		return "", err
+	}
+
+	// 解析输出以获取 UUID
+	lines := strings.Split(string(out), "\n")
+	for _, line := range lines {
+		if strings.Contains(line, "Hardware UUID:") {
+			parts := strings.Split(line, ":")
+			if len(parts) == 2 {
+				return strings.TrimSpace(parts[1]), nil
+			}
+		}
+	}
+
+	return "", fmt.Errorf("could not find hardware UUID in system_profiler output")
+}
+
+func GetProductUUIDWin() (string, error) {
+	// 使用系统调用获取计算机的注册表信息
+	out, err := exec.Command("reg", "query", "HKEY_LOCAL_MACHINE\\SOFTWARE\\Microsoft\\Cryptography", "/v", "MachineGuid").Output()
+	if err != nil {
+		return "", err
+	}
+
+	// 解析注册表输出以获取 MachineGuid
+	lines := strings.Split(string(out), "\n")
+	for _, line := range lines {
+		if strings.Contains(line, "REG_SZ") {
+			parts := strings.Fields(line)
+			if len(parts) >= 3 {
+				return strings.TrimSpace(parts[2]), nil
+			}
+		}
+	}
+
+	return "", fmt.Errorf("MachineGuid not found in registry")
+}

license/public/privatekey.go

@@ -0,0 +1,29 @@
+package public
+
+import "strings"
+
+var (
+	ECDSA_PRIVATE = ""
+
+// 	ECDSA_PRIVATE = `-----BEGIN EC PRIVATE KEY-----
+// MIHcAgEBBEIB0pE4uFaWRx7t03BsYlYvF1YvKaBGyvoakxnodm9ou0R9wC+sJAjH
+// QZZJikOg4SwNqgQ/hyrOuDK2oAVHhgVGcYmgBwYFK4EEACOhgYkDgYYABAAJXIuw
+// 12MUzpHggia9POBFYXSxaOGKGbMjIyDI+6q7wi7LMw3HgbaOmgIqFG72o8JBQwYN
+// 4IbXHf+f86CRY1AA2wHzbHvt6IhkCXTNxBEffa1yMUgu8n9cKKF2iLgyQKcKqW33
+// 8fGOw/n3Rm2Yd/EB56u2rnD29qS+nOM9eGS+gy39OQ==
+// -----END EC PRIVATE KEY-----`
+)
+
+func (nea *NonEquAlgorthm) SignedBytes(plainText []byte) (string, error) {
+	var (
+		sig string
+		err error
+	)
+
+	encBase64String := EncodeSegment(plainText)
+	if sig, err = nea.Algorithm.Sign(encBase64String, nea.PrivateKey); err != nil {
+		return "", err
+	}
+
+	return strings.Join([]string{encBase64String, sig}, "."), nil
+}

license/public/public.go

@@ -0,0 +1,100 @@
+package public
+
+import (
+	"bytes"
+	"fmt"
+	"os"
+
+	"github.com/google/uuid"
+)
+
+func GetUUID() string {
+	return uuid.New().String()
+}
+
+func CheckFileIsExist(filename string) bool {
+	var exist = true
+	if _, err := os.Stat(filename); os.IsNotExist(err) {
+		exist = false
+	}
+	return exist
+}
+
+func RemoveDuplicate(list []int) []int {
+	var x []int
+	for _, i := range list {
+		if len(x) == 0 {
+			x = append(x, i)
+		} else {
+			for k, v := range x {
+				if i == v {
+					break
+				}
+				if k == len(x)-1 {
+					x = append(x, i)
+				}
+			}
+		}
+	}
+	return x
+}
+
+func FormatPem(content []byte, isPriKey bool) ([]byte, error) {
+	var (
+		pemPubKeyStart = []byte("-----BEGIN PUBLIC KEY-----")
+		pemPubKeyEnd   = []byte("-----END PUBLIC KEY-----")
+
+		pemPriKeyStart = []byte("-----BEGIN EC PRIVATE KEY-----")
+		pemPriKeyEnd   = []byte("-----END EC PRIVATE KEY-----")
+
+		orgContent []byte
+		pemContent []byte
+		slen       = len(content)
+	)
+
+	if slen == 0 {
+		if isPriKey && slen < (len(pemPriKeyStart)+len(pemPriKeyEnd)) {
+			return nil, fmt.Errorf("%s", "prikey centent length not enough")
+		} else if slen < (len(pemPubKeyStart) + len(pemPubKeyEnd)) {
+			return nil, fmt.Errorf("%s", "pubkey centent length not enough")
+		}
+		return nil, fmt.Errorf("%s", "content not null")
+	} else {
+		orgContent = make([]byte, slen)
+		copy(orgContent, content)
+		orgContent = bytes.Join(bytes.Fields(orgContent), []byte(" "))
+		slen = len(orgContent)
+	}
+
+	if isPriKey {
+		//私钥
+		if bytes.HasPrefix(orgContent, pemPriKeyStart) && bytes.HasSuffix(orgContent, pemPriKeyEnd) {
+			pemContent = append(pemContent, pemPriKeyStart...)
+			for _, v := range orgContent[len(pemPriKeyStart) : slen-len(pemPriKeyEnd)] {
+				if v == ' ' {
+					pemContent = append(pemContent, []byte("\n")...)
+				} else {
+					pemContent = append(pemContent, v)
+				}
+			}
+			pemContent = append(pemContent, pemPriKeyEnd...)
+		} else {
+			return nil, fmt.Errorf("%s", "pem file prefix must be 'BEGIN/END EC PRIVATE KEY' format")
+		}
+	} else {
+		if bytes.HasPrefix(orgContent, pemPubKeyStart) && bytes.HasSuffix(orgContent, pemPubKeyEnd) {
+			pemContent = append(pemContent, pemPubKeyStart...)
+			for _, v := range orgContent[len(pemPubKeyStart) : slen-len(pemPubKeyEnd)] {
+				if v == ' ' {
+					pemContent = append(pemContent, []byte("\n")...)
+				} else {
+					pemContent = append(pemContent, v)
+				}
+			}
+			pemContent = append(pemContent, pemPubKeyEnd...)
+		} else {
+			return nil, fmt.Errorf("%s", "pem file prefix must be 'BEGIN/END PUBLIC KEY' format")
+		}
+	}
+	return pemContent, nil
+}

license/public/publickey.go

@@ -0,0 +1,39 @@
+package public
+
+import (
+	"fmt"
+	"strings"
+)
+
+var (
+	ECDSA_PUBLICKEY = ""
+
+// 	ECDSA_PUBLICKEY = `-----BEGIN PUBLIC KEY-----
+// MIGbMBAGByqGSM49AgEGBSuBBAAjA4GGAAQACVyLsNdjFM6R4IImvTzgRWF0sWjh
+// ihmzIyMgyPuqu8IuyzMNx4G2jpoCKhRu9qPCQUMGDeCG1x3/n/OgkWNQANsB82x7
+// 7eiIZAl0zcQRH32tcjFILvJ/XCihdoi4MkCnCqlt9/HxjsP590ZtmHfxAeertq5w
+// 9vakvpzjPXhkvoMt/Tk=
+// -----END PUBLIC KEY-----`
+)
+
+func (nea *NonEquAlgorthm) VerifySign(cipherText string) ([]byte, error) {
+	var (
+		err error
+	)
+
+	parts := strings.Split(cipherText, ".")
+	if len(parts) != 2 {
+		return nil, fmt.Errorf("%s", "license contains an invalid number of segments")
+	}
+
+	if err = nea.Algorithm.Verify(parts[0], parts[1], nea.PublicKey); err != nil {
+		return nil, err
+	}
+
+	plainBytes, err := DecodeSegment(parts[0])
+	if err != nil {
+		return nil, err
+	}
+
+	return plainBytes, nil
+}

license/public/version.go

@@ -0,0 +1,24 @@
+package public
+
+import (
+	"bytes"
+)
+
+var (
+	buildName     = "switch-license"
+	buildVersion  = "development"
+	buildBranch   = "development"
+	buildCommitID = "development"
+	buildTime     = "development"
+)
+
+func GetAppInfo() string {
+	b := bytes.NewBufferString("\nAppInfo:")
+	b.WriteString("\n	Name: 	   " + buildName)
+	b.WriteString("\n	Version:   " + buildVersion)
+	b.WriteString("\n	Branch:    " + buildBranch)
+	b.WriteString("\n	CommitID:  " + buildCommitID)
+	b.WriteString("\n	BuildTime: " + buildTime)
+	b.WriteString("\n")
+	return b.String()
+}

license/with_license.go

@@ -0,0 +1,53 @@
+//go:build license
+// +build license
+
+package license
+
+import (
+	"encoding/pem"
+	"fmt"
+	"log"
+	"os"
+	"runtime"
+
+	"starwiz.cn/sjy01/image-proc/license/public"
+)
+
+var ProductName = "sjy01-data-processing"
+
+func VerifyLicense() {
+	machineID, err := public.GetMachineID()
+	if err != nil {
+		fmt.Println(err)
+	}
+	fmt.Println("machine ID:", machineID)
+	fmt.Println("product name:", ProductName)
+
+	if runtime.GOOS == "darwin" || runtime.GOOS == "windows" {
+		fmt.Println("license is not required on this platform:", runtime.GOOS)
+		return
+	}
+
+	value, err := public.LoadKey("license.lic")
+	if err != nil {
+		log.Println("invalid license", err)
+		os.Exit(130)
+	}
+
+	var block *pem.Block
+	if block, _ = pem.Decode(value); block == nil {
+		fmt.Printf("%s\n", "license must be PEM")
+		os.Exit(130)
+	}
+
+	l, err := public.BytesToLicense(string(block.Bytes))
+	if err != nil {
+		fmt.Println("invalid license", err)
+		os.Exit(131)
+	}
+
+	if err := l.Valid(ProductName, machineID); err != nil {
+		fmt.Println("invalid license", err)
+		os.Exit(132)
+	}
+}

pkg/cloud-cover/k_means.go

@@ -0,0 +1,176 @@
+package cloudcover
+
+import (
+	"image"
+	"image/color"
+
+	"math"
+	"math/rand"
+	"os"
+	"time"
+
+	log "github.com/sirupsen/logrus"
+	"gocv.io/x/gocv"
+)
+
+const (
+	k                   = 2   // 聚类数目，2表示云和非云
+	maxIters            = 100 // 最大迭代次数
+	tolerance           = 1.0 // 聚类中心变化容忍度
+	brightnessThreshold = 220 // 亮度阈值，值越大要求云更亮
+)
+
+type Pixel struct {
+	r, g, b float64
+}
+
+func CloudPercentByKMeans(imagePath string) float64 {
+	file, err := os.Open(imagePath)
+	if err != nil {
+		log.Errorf("无法打开图像: %v", err)
+		return 0.0
+	}
+	defer file.Close()
+
+	img, _, err := image.Decode(file)
+	if err != nil {
+		log.Errorf("无法解码图像: %v", err)
+		return 0.0
+	}
+
+	return computeCloudCoverByKMeans(img)
+}
+
+func CloudPercentByKMeans2(imgMat gocv.Mat) float64 {
+	img, err := imgMat.ToImage()
+	if err != nil {
+		log.Errorf("无法将Mat转换为Image: %v", err)
+		return 0.0
+	}
+	return computeCloudCoverByKMeans(img)
+}
+
+func computeCloudCoverByKMeans(img image.Image) float64 {
+	log.Info("cloud cover computing...")
+	bounds := img.Bounds()
+	width, height := bounds.Max.X, bounds.Max.Y
+	pixels := make([]Pixel, 0, width*height)
+
+	for y := 0; y < height; y++ {
+		for x := 0; x < width; x++ {
+			r, g, b, _ := img.At(x, y).RGBA()
+			pixels = append(pixels, Pixel{
+				r: float64(r >> 8),
+				g: float64(g >> 8),
+				b: float64(b >> 8),
+			})
+		}
+	}
+
+	rand.Seed(time.Now().UnixNano())
+	centroids := initializeCentroids(pixels)
+	assignments := make([]int, len(pixels))
+
+	for iter := 0; iter < maxIters; iter++ {
+		for i, pixel := range pixels {
+			assignments[i] = nearestCentroid(pixel, centroids)
+		}
+
+		newCentroids := updateCentroids(pixels, assignments)
+		if hasConverged(centroids, newCentroids) {
+			log.Debugf("Convergence after the %d iteration", iter+1)
+			break
+		}
+		centroids = newCentroids
+	}
+
+	cloudPixels := 0
+	totalPixels := width * height
+
+	// 新增：检查哪个聚类的平均亮度更高，确保白色表示云
+	cloudCluster := 0
+	if centroids[1].r+centroids[1].g+centroids[1].b > centroids[0].r+centroids[0].g+centroids[0].b {
+		cloudCluster = 1
+	}
+
+	outputImg := image.NewGray(bounds)
+	for y := 0; y < height; y++ {
+		for x := 0; x < width; x++ {
+			i := y*width + x
+			pixel := pixels[i]
+			// 计算亮度和颜色比值
+			brightness := (pixel.r + pixel.g + pixel.b) / 3
+			blueRedRatio := pixel.b / pixel.r
+			greenRedRatio := pixel.g / pixel.r
+
+			// 判断是否为云的条件：亮度大于阈值并且蓝色/红色比率较高
+			if (assignments[i] == cloudCluster && brightness > brightnessThreshold && blueRedRatio > 0.5) ||
+				(assignments[i] == cloudCluster && greenRedRatio > 0.8) {
+				outputImg.SetGray(x, y, color.Gray{Y: 255}) // 云区域
+				cloudPixels++
+			} else {
+				outputImg.SetGray(x, y, color.Gray{Y: 0}) // 非云区域
+			}
+		}
+	}
+
+	cloudRatio := float64(cloudPixels) / float64(totalPixels)
+	return cloudRatio
+}
+
+func initializeCentroids(pixels []Pixel) []Pixel {
+	centroids := make([]Pixel, k)
+	for i := range centroids {
+		centroids[i] = pixels[rand.Intn(len(pixels))]
+	}
+	return centroids
+}
+
+func nearestCentroid(pixel Pixel, centroids []Pixel) int {
+	minDist := math.Inf(1)
+	bestIndex := 0
+	for i, c := range centroids {
+		dist := math.Sqrt(math.Pow(pixel.r-c.r, 2) + math.Pow(pixel.g-c.g, 2) + math.Pow(pixel.b-c.b, 2))
+		if dist < minDist {
+			minDist = dist
+			bestIndex = i
+		}
+	}
+	return bestIndex
+}
+
+func updateCentroids(pixels []Pixel, assignments []int) []Pixel {
+	sums := make([]Pixel, k)
+	counts := make([]int, k)
+
+	for i, assignment := range assignments {
+		sums[assignment].r += pixels[i].r
+		sums[assignment].g += pixels[i].g
+		sums[assignment].b += pixels[i].b
+		counts[assignment]++
+	}
+
+	newCentroids := make([]Pixel, k)
+	for i := range newCentroids {
+		if counts[i] > 0 {
+			newCentroids[i].r = sums[i].r / float64(counts[i])
+			newCentroids[i].g = sums[i].g / float64(counts[i])
+			newCentroids[i].b = sums[i].b / float64(counts[i])
+		} else {
+			newCentroids[i] = sums[i]
+		}
+	}
+
+	return newCentroids
+}
+
+func hasConverged(oldCentroids, newCentroids []Pixel) bool {
+	for i := range oldCentroids {
+		if math.Abs(oldCentroids[i].r-newCentroids[i].r) > tolerance ||
+			math.Abs(oldCentroids[i].g-newCentroids[i].g) > tolerance ||
+			math.Abs(oldCentroids[i].b-newCentroids[i].b) > tolerance {
+			return false
+		}
+	}
+	return true
+}

pkg/fusion/fusion.go

@@ -24,6 +24,10 @@ const (
 	GDALPansharpen
 	IHS
 	GIHS
+	WeightR   = 1
+	WeightG   = 0.75
+	WeightB   = 0.55
+	WeightNIR = 0.85
 )
 
 func Pansharpen(panTif, mssTif, fusTif string, method PansharpenMethod, weight float32) error {

pkg/producer/aux.go

@@ -21,7 +21,7 @@ import (
 	"starwiz.cn/sjy01/image-proc/pkg/config"
 )
 
-func (r *Registrator) LoadAuxData() error {
+func (r *ImgProc) LoadAuxData() error {
 	var err error
 	r.auxHeads, r.auxBoxes, r.AuxPlatforms, err = auxilary.ExtractAux(r.Params.AuxRawFile)
 
@@ -37,7 +37,7 @@ func (r *Registrator) LoadAuxData() error {
 }
 
 // 数据校验和测试
-func (r *Registrator) AuxPrint() {
+func (r *ImgProc) AuxPrint() {
 	var fcPos84 geojson.FeatureCollection
 	var fcPos84Interp geojson.FeatureCollection
 	for _, p := range r.AuxPlatforms {
@@ -63,7 +63,7 @@ func (r *Registrator) AuxPrint() {
 	f.Write(data)
 }
 
-func (r *Registrator) SceneImageTime(scene *Scene) (start, center, end time.Time) {
+func (r *ImgProc) SceneImageTime(scene *Scene) (start, center, end time.Time) {
 	startPosInAux, endPosInAux := r.SceneInAuxIndex(scene)
 	centerPosInAux := (startPosInAux + endPosInAux) / 2
 
@@ -78,7 +78,7 @@ func (r *Registrator) SceneImageTime(scene *Scene) (start, center, end time.Time
 }
 
 // FIXME: This function is not accurate enough. 四元数、成像时刻、GPS 等需要修改为插值获取
-func (r *Registrator) ComputeMetaAndRPC(scene *Scene) (topLeft, bottomRight orb.Point) {
+func (r *ImgProc) ComputeMetaAndRPC(scene *Scene) (topLeft, bottomRight orb.Point) {
 	log.Info("using attitude quaternion to calculate image boundary...")
 	line0Start := r.calculateLatLonH(scene, 0, 0, 0)
 	line0End := r.calculateLatLonH(scene, 0, scene.Width, 0)
@@ -169,13 +169,13 @@ func (r *Registrator) ComputeMetaAndRPC(scene *Scene) (topLeft, bottomRight orb.
 	return
 }
 
-func (r *Registrator) SceneInAuxIndex(scene *Scene) (int, int) {
+func (r *ImgProc) SceneInAuxIndex(scene *Scene) (int, int) {
 	startPosInAux := r.sceneOffsetInAuxIndex(scene, 0)
 	endPosInAux := r.sceneOffsetInAuxIndex(scene, scene.Height)
 	return startPosInAux, endPosInAux
 }
 
-func (r *Registrator) sceneOffsetInAuxIndex(scene *Scene, offset int) int {
+func (r *ImgProc) sceneOffsetInAuxIndex(scene *Scene, offset int) int {
 	var auxForImageRow int
 	switch scene.Type {
 	case "MSS":
@@ -194,7 +194,7 @@ func (r *Registrator) sceneOffsetInAuxIndex(scene *Scene, offset int) int {
 }
 
 // row, col 相对于图像景左上角, H 为地面目标点高度
-func (r *Registrator) calculateLatLonH(scene *Scene, row, col, H int) calculator.IntersectionPoint {
+func (r *ImgProc) calculateLatLonH(scene *Scene, row, col, H int) calculator.IntersectionPoint {
 	// 内插值获取图像行时刻
 	ucam := col
 	cross := 16

pkg/producer/dark_bias.go

@@ -1,11 +1,11 @@
 package producer
 
 const (
-	DarkBiasValue = 15675.0/6
+	DarkBiasValue = 15675.0 / 6
 )
 
 // 暗场偏置校正 固定值15675校正
-func (r *Registrator) DoDarkBiasCorrection() {
+func (r *ImgProc) DoDarkBiasCorrection() {
 	r.PanImage.SubtractFloat(DarkBiasValue)
 	for i := 0; i < len(r.registeredMssImages); i++ {
 		r.registeredMssImages[i].SubtractFloat(DarkBiasValue)

pkg/producer/ecc_align.go

@@ -1,3 +1,55 @@
 package producer
 
+import (
+	"image"
+
+	log "github.com/sirupsen/logrus"
+	"gocv.io/x/gocv"
+)
+
 // align the image via Enhanced Correlation Coefficient (ECC) algorithm
+
+func CV_ECCAlign(templateImage, inputImage gocv.Mat) gocv.Mat {
+	retMatType := templateImage.Type()
+	warpMode := gocv.MotionHomography
+	var warpMatrix gocv.Mat
+	switch warpMode {
+	case gocv.MotionHomography:
+		warpMatrix = gocv.NewMatWithSize(3, 3, gocv.MatTypeCV32FC1)
+	default:
+		warpMatrix = gocv.NewMatWithSize(2, 3, gocv.MatTypeCV32FC1)
+	}
+	for i := 0; i < warpMatrix.Rows(); i++ {
+		warpMatrix.SetFloatAt(i, i, 1.0)
+	}
+
+	criteria := gocv.NewTermCriteria(gocv.Count+gocv.EPS, 1000, 0.0001)
+
+	var mask = gocv.NewMat()
+	templateImage.ConvertTo(&templateImage, gocv.MatTypeCV32FC1)
+	inputImage.ConvertTo(&inputImage, gocv.MatTypeCV32FC1)
+
+	retval := gocv.FindTransformECC(templateImage, inputImage, &warpMatrix, warpMode, criteria, mask, 5)
+	log.Info("cv::findTransformECC ret:", retval)
+
+	var aligned = gocv.NewMat()
+	sz := image.Point{inputImage.Cols(), inputImage.Rows()}
+	if warpMode == gocv.MotionHomography {
+		gocv.WarpPerspective(inputImage, &aligned, warpMatrix, sz)
+	} else {
+		gocv.WarpAffine(inputImage, &aligned, warpMatrix, sz)
+	}
+
+	aligned.ConvertTo(&aligned, retMatType)
+	templateImage.ConvertTo(&templateImage, retMatType)
+	inputImage.ConvertTo(&inputImage, retMatType)
+	return aligned
+}
+
+func (r *ImgProc) DoECCAlign() error {
+	r.registeredMssImages[0] = r.MssImages[0].Clone()
+	for i := 1; i < len(r.MssImages); i++ {
+		r.registeredMssImages[i] = CV_ECCAlign(r.registeredMssImages[0], r.MssImages[i])
+	}
+	return nil
+}

pkg/producer/edge_phase_correlation.go

@@ -1,8 +1,6 @@
 package producer
 
 import (
-	"fmt"
-
 	"gocv.io/x/gocv"
 )
 
@@ -14,30 +12,30 @@ import (
 
 // SKIP: 多光谱各个波段的边缘检测结果不佳
 func CV_Canny(img0 gocv.Mat) gocv.Mat {
-	fmt.Println(img0.Cols(), img0.Rows(), img0.Type().String())
 	dst8 := gocv.NewMatWithSize(img0.Rows(), img0.Cols(), gocv.MatTypeCV8U)
 	defer dst8.Close()
 	gocv.Normalize(img0, &dst8, 0, 255, gocv.NormMinMax)
 	dst8.ConvertTo(&dst8, gocv.MatTypeCV8U)
 	dstEdge := gocv.NewMat()
 	gocv.Canny(dst8, &dstEdge, 10, 100)
-	dstEdge.ConvertTo(&dstEdge, gocv.MatTypeCV16U)
+	dstEdge.ConvertTo(&dstEdge, img0.Type())
 	return dstEdge
 }
 
 func CV_Sobel(img0 gocv.Mat) gocv.Mat {
+	matType := img0.Type()
 	// x 方向
 	sobelX := gocv.NewMat()
-	gocv.Sobel(img0, &sobelX, gocv.MatTypeCV32F, 1, 0, 5, 1.5, 0, gocv.BorderDefault)
-	sobelX.ConvertTo(&sobelX, gocv.MatTypeCV16U)
+	gocv.Sobel(img0, &sobelX, gocv.MatTypeCV32FC1, 1, 0, 5, 1.5, 0, gocv.BorderDefault)
+	sobelX.ConvertTo(&sobelX, matType)
 	// y 方向
 	sobelY := gocv.NewMat()
-	gocv.Sobel(img0, &sobelY, gocv.MatTypeCV32F, 0, 1, 5, 1.5, 0, gocv.BorderIsolated)
-	sobelY.ConvertTo(&sobelY, gocv.MatTypeCV16U)
+	gocv.Sobel(img0, &sobelY, gocv.MatTypeCV32FC1, 0, 1, 5, 1.5, 0, gocv.BorderIsolated)
+	sobelY.ConvertTo(&sobelY, matType)
 	// 合并
 	sobelXY := gocv.NewMat()
-	gocv.Sobel(img0, &sobelXY, gocv.MatTypeCV32F, 1, 1, 5, 1.5, 0, gocv.BorderDefault)
-	sobelXY.ConvertTo(&sobelXY, gocv.MatTypeCV16U)
+	gocv.Sobel(img0, &sobelXY, gocv.MatTypeCV32FC1, 1, 1, 5, 1.5, 0, gocv.BorderDefault)
+	sobelXY.ConvertTo(&sobelXY, matType)
 
 	return sobelY
 }

pkg/producer/image_registration.go

@@ -22,8 +22,6 @@ const (
 	PixelBytes                          = 2
 	PanWidth                            = payload.PAN_PIXEL_WIDTH // 像素宽度
 	MssWidth                            = payload.MSS_PIXEL_WIDTH
-	BlockNH                             = 8
-	BlockNW                             = 4
 	OverlappedBlockLines                = 1000 // 重叠块的行数
 	DownSampled          ResampleMethod = "down_sample_pan"
 	UpSampled            ResampleMethod = "up_sample_mss"
@@ -33,11 +31,12 @@ const (
 	ReferenceShiftYB2                   = 200.0
 	ReferenceShiftYB3                   = 300.0
 	ReferenceShiftYB4                   = 400.0
+	MaxShiftY                           = 480
 )
 
 type ResampleMethod string
 
-type Registrator struct {
+type ImgProc struct {
 	Params Params
 
 	PanImage  gocv.Mat
@@ -47,6 +46,8 @@ type Registrator struct {
 	MssImages [4]gocv.Mat
 	MssHeight int
 	MssWidth  int
+	blockH    int
+	blockW    int
 
 	shiftMutex   sync.Mutex
 	phaseShifts  [4][]PhaseShiftM
@@ -69,14 +70,16 @@ type Registrator struct {
 	report Report
 }
 
-func NewRegistrator(rsmethod ResampleMethod) *Registrator {
-	var r Registrator
+func NewImgProc(rsmethod ResampleMethod) *ImgProc {
+	var r ImgProc
 	r.resampleMethod = rsmethod
+	r.blockH = 8
+	r.blockW = 4
 
 	return &r
 }
 
-func (r *Registrator) LoadPanRaw() error {
+func (r *ImgProc) LoadPanRaw() error {
 	data, err := os.ReadFile(r.Params.PanRawFile)
 	if err != nil {
 		log.Error("Error reading raw file: ", err)
@@ -111,7 +114,7 @@ func (r *Registrator) LoadPanRaw() error {
 	return nil
 }
 
-func (r *Registrator) LoadMssRaw() error {
+func (r *ImgProc) LoadMssRaw() error {
 	data, err := os.ReadFile(r.Params.MssRawFile)
 	if err != nil {
 		log.Error("Error reading raw file: ", err)
@@ -141,7 +144,14 @@ func (r *Registrator) LoadMssRaw() error {
 	return nil
 }
 
-func (r *Registrator) DoPhaseCorrelation() error {
+func (r *ImgProc) DoPhaseCorrelation(skip bool) error {
+	if skip {
+		for i := 0; i < MssBands; i++ {
+			r.registeredMssImages[i] = r.MssImages[i].Clone()
+		}
+		return nil
+	}
+
 	switch r.resampleMethod {
 	case UpSampled:
 		return r.CalcUpPhaseCorrelation()
@@ -151,7 +161,7 @@ func (r *Registrator) DoPhaseCorrelation() error {
 }
 
 // 将PAN降采样后计算相位相关的偏移量
-func (r *Registrator) CalcDownPhaseCorrelation() error {
+func (r *ImgProc) 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")
@@ -167,8 +177,10 @@ func (r *Registrator) CalcDownPhaseCorrelation() error {
 	log.Println("down sampled PAN images size:", downsampledPanImage.Size())
 
 	// 分块高度
-	blockHeight := r.MssHeight / BlockNH
-	blockWidth := r.MssWidth / BlockNW
+	blockHeight := r.MssHeight / r.blockH
+	blockWidth := r.MssWidth / r.blockW
+
+	log.Infof("blockHeight=%d, blockWidth=%d", blockHeight, blockWidth)
 
 	// 在 MSS 4 个波段上进行配准
 	err := r.doPhaseCorrelation(r.MssImages[0],
@@ -177,23 +189,10 @@ func (r *Registrator) CalcDownPhaseCorrelation() error {
 	if err != nil {
 		return err
 	}
-	r.fileterPhaseShift([]float64{64, 64, 64, 64}, true)
+	r.fileterPhaseShift([]float64{0, 100, 200, 300})
 	r.calcMSSDeltaCoeffs(4)
 	r.DoMSSCoRegistration(false)
 
-	// 边缘检测后再做一次配准
-
-	// var mssEdges []gocv.Mat
-	// for band := 0; band < len(r.registeredMssImages); band++ {
-	// 	edge := CV_Sobel(r.registeredMssImages[band])
-	// 	mssEdges = append(mssEdges, edge)
-	// }
-	// r.doPhaseCorrelation(mssEdges[0], mssEdges,
-	// 	r.MssHeight, r.MssWidth, blockHeight, blockWidth)
-	// r.fileterPhaseShift([]float64{5, 5, 5, 5}, false)
-	// r.calcMSSDeltaCoeffs(4)
-	// r.DoMSSCoRegistration(true)
-
 	// 基于 PAN 图像进行配准
 	err = r.doPhaseCorrelation(downsampledPanImage,
 		[]gocv.Mat{r.registeredMssImages[0]},
@@ -201,16 +200,26 @@ func (r *Registrator) CalcDownPhaseCorrelation() error {
 	if err != nil {
 		return err
 	}
-	r.fileterPhaseShift([]float64{30.0}, true)
+	r.fileterPhaseShift([]float64{100.0})
 	r.calcMSSDeltaCoeffs(1)
 	r.DoPANCoRegistration()
 
+	// 裁掉位移部分
+	r.MssHeight -= MaxShiftY
+	r.PanHeight -= MaxShiftY * 4
+	pMat := r.PanImage.Region(image.Rect(0, 0, r.PanWidth, r.PanHeight))
+	r.PanImage = pMat.Clone()
+	pMat.Close()
+	for i := 0; i < MssBands; i++ {
+		r.registeredMssImages[i] = r.registeredMssImages[i].Region(image.Rect(0, 0, r.MssWidth, r.MssHeight))
+	}
+
 	return nil
 }
 
 // 将MSS升采样采样后计算相位相关的偏移量
-func (r *Registrator) CalcUpPhaseCorrelation() error {
-	log.Fatal("unsuppotted up-resample method")
+func (r *ImgProc) CalcUpPhaseCorrelation() error {
+	log.Fatal("not implemented yet")
 	// 确保 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")
@@ -219,37 +228,39 @@ func (r *Registrator) CalcUpPhaseCorrelation() error {
 	}
 
 	// 在 MSS 4 个波段上进行配准
-	err := r.doPhaseCorrelation(r.MssImages[0],
-		[]gocv.Mat{r.MssImages[0], r.MssImages[1], r.MssImages[2], r.MssImages[3]},
-		r.MssHeight, r.MssWidth, r.MssHeight/BlockNH, r.MssWidth/BlockNW)
-	if err != nil {
-		return err
-	}
-	r.DoMSSCoRegistration(false)
+	// err := r.doPhaseCorrelation(r.MssImages[0],
+	// 	[]gocv.Mat{r.MssImages[0], r.MssImages[1], r.MssImages[2], r.MssImages[3]},
+	// 	r.MssHeight, r.MssWidth, r.MssHeight/r.blockH, r.MssWidth/r.blockW)
+	// if err != nil {
+	// 	return err
+	// }
+	// r.DoMSSCoRegistration(false)
 
 	upsampledMssImages := make([]gocv.Mat, MssBands)
 	for i := 0; i < MssBands; i++ {
 		upsampledMssImages[i] = gocv.NewMat()
-		gocv.Resize(r.registeredMssImages[i], &upsampledMssImages[i],
+		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
+	blockHeight := r.PanHeight / r.blockH
+	blockWidth := r.PanWidth / r.blockW
 
 	log.Infof("blockHeight=%d, blockWidth=%d", blockHeight, blockWidth)
 
 	// 基于 PAN 图像进行配准
-	r.doPhaseCorrelation(r.PanImage,
-		[]gocv.Mat{upsampledMssImages[0]},
-		r.MssHeight, r.MssWidth, blockHeight, blockWidth)
-	return r.DoPANCoRegistration()
+	r.doPhaseCorrelation(r.PanImage, upsampledMssImages,
+		r.PanHeight, r.PanWidth, blockHeight, blockWidth)
+	r.fileterPhaseShift([]float64{64, 64, 64, 64})
+	r.calcMSSDeltaCoeffs(4)
+	r.DoMSSCoRegistration(false)
+	return nil
 }
 
-func (r *Registrator) doPhaseCorrelation(base gocv.Mat,
+func (r *ImgProc) doPhaseCorrelation(base gocv.Mat,
 	mssImages []gocv.Mat,
 	height, width,
 	blockHeight, blockWidth int) error {
@@ -261,8 +272,8 @@ func (r *Registrator) doPhaseCorrelation(base gocv.Mat,
 		r.deltaYCoeffs[band] = make([]float64, 0)
 	}
 
-	for bh := 0; bh < BlockNH; bh++ {
-		for bw := 0; bw < BlockNW; bw++ {
+	for bh := 0; bh < r.blockH; bh++ {
+		for bw := 0; bw < r.blockW; bw++ {
 			wg.Add(1)
 			go func(bh, bw int) {
 				defer wg.Done()
@@ -297,7 +308,7 @@ func (r *Registrator) doPhaseCorrelation(base gocv.Mat,
 					mssBlock := mssImages[band].Region(rect)
 
 					// 处理每个分块
-					phaseShift, response := r.calculateBlockPhaseShift(panBlock, mssBlock)
+					phaseShift, response := CV_PhaseCorrelate(panBlock, mssBlock)
 					shiftM.dx = phaseShift.X
 					shiftM.dy = phaseShift.Y
 					shiftM.response = response
@@ -327,7 +338,7 @@ func (r *Registrator) doPhaseCorrelation(base gocv.Mat,
 	return nil
 }
 
-func (r *Registrator) Clean() {
+func (r *ImgProc) Clean() {
 	r.PanImage.Close()
 	for i := 0; i < MssBands; i++ {
 		r.MssImages[i].Close()
@@ -340,7 +351,7 @@ func (r *Registrator) Clean() {
 	r.rgbirImage.Close()
 }
 
-func (r *Registrator) calcMSSDeltaCoeffs(bands int) error {
+func (r *ImgProc) calcMSSDeltaCoeffs(bands int) error {
 	// 计算每个通道的delta多项式拟合系数
 	for i := 0; i < bands; i++ {
 		var cx []float64
@@ -397,7 +408,7 @@ func (r *Registrator) calcMSSDeltaCoeffs(bands int) error {
 	return nil
 }
 
-func (r *Registrator) DoMSSCoRegistration(byEdge bool) error {
+func (r *ImgProc) DoMSSCoRegistration(byEdge bool) error {
 	for band := 0; band < MssBands; band++ {
 		if len(r.deltaXCoeffs[band]) < 2 || len(r.deltaYCoeffs[band]) < 3 {
 			log.Errorf("delta coefficients not calculated, skip co-registration %d", band+1)
@@ -449,14 +460,10 @@ func (r *Registrator) DoMSSCoRegistration(byEdge bool) error {
 		mapY.Close()
 	}
 
-	// 裁掉末尾的的 MSS 480 行 和 PAN 的 480*4 行
-	r.PanHeight -= 360 * 4
-	r.MssHeight -= 360
-
 	return nil
 }
 
-func (r *Registrator) DoPANCoRegistration() error {
+func (r *ImgProc) DoPANCoRegistration() error {
 	if len(r.deltaXCoeffs[0]) < 2 || len(r.deltaYCoeffs[0]) < 3 {
 		log.Error("delta coefficients not calculated, skip co-registration")
 		return nil
@@ -498,9 +505,6 @@ func (r *Registrator) DoPANCoRegistration() error {
 		r.registeredMssImages[i].Close()
 		r.registeredMssImages[i] = registeredMSS
 	}
-	// 裁掉末尾的的 MSS 480 行 和 PAN 的 480*4 行
-	r.PanHeight -= 120 * 4
-	r.MssHeight -= 120
 
 	return nil
 }

pkg/producer/meta.go

@@ -72,7 +72,7 @@ func (corners Corners) Extend() (lng1, lat1, lng2, lat2 float64) {
 	return
 }
 
-func (r *Registrator) makeProductMeta(scene *Scene) *ProductMeta {
+func (r *ImgProc) makeProductMeta(scene *Scene) *ProductMeta {
 	meta := &ProductMeta{
 		Satellite:      "SJY01",
 		Sensor:         "PMS",

pkg/producer/output.go

@@ -13,7 +13,7 @@ import (
 	"starwiz.cn/sjy01/image-proc/pkg/utils"
 )
 
-func (r *Registrator) SaveOriginalPanToGDALGTiff(tiffFile string) error {
+func (r *ImgProc) SaveOriginalPanToGDALGTiff(tiffFile string) error {
 	err := utils.SavePanToGDALGTiff(r.PanImage, 0, 0, tiffFile, PanResolution)
 	if err != nil {
 		return err
@@ -22,7 +22,7 @@ func (r *Registrator) SaveOriginalPanToGDALGTiff(tiffFile string) error {
 	return nil
 }
 
-func (r *Registrator) SaveRegisteredMssToGDALGTiff(tiffFile string) error {
+func (r *ImgProc) SaveRegisteredMssToGDALGTiff(tiffFile string) error {
 	r.rgbirImage = gocv.NewMat()
 	gocv.Merge(r.registeredMssImages[:], &r.rgbirImage)
 
@@ -37,7 +37,7 @@ func (r *Registrator) SaveRegisteredMssToGDALGTiff(tiffFile string) error {
 	return nil
 }
 
-func (r *Registrator) BytesToRaw(mssData []byte, filePath string) error {
+func (r *ImgProc) 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
@@ -48,15 +48,15 @@ func (r *Registrator) BytesToRaw(mssData []byte, filePath string) error {
 	return nil
 }
 
-func (r *Registrator) SaveRegisteredMssToRaw(raw string) error {
+func (r *ImgProc) SaveRegisteredMssToRaw(raw string) error {
 	return saveRegisteredMssToRaw(raw, r.registeredMssImages)
 }
 
-func (r *Registrator) Report() error {
+func (r *ImgProc) Report() error {
 	return WriteReport(&r.report, r.Params.ReportFile)
 }
 
-func (r *Registrator) rpcKeywordInTif() {
+func (r *ImgProc) rpcKeywordInTif() {
 	// GDAL库对应的RPC关键词
 	// keys := []string{
 	// 	"ERR_BIAS", "ERR_RAND",

pkg/producer/phase_correlation.go

@@ -3,8 +3,10 @@ package producer
 import (
 	"errors"
 	"image"
+	"math/cmplx"
 
 	"github.com/duke-git/lancet/v2/slice"
+	"github.com/mjibson/go-dsp/fft"
 	log "github.com/sirupsen/logrus"
 	"gocv.io/x/gocv"
 )
@@ -22,7 +24,7 @@ type Block struct {
 	coord  image.Point // top-left corner of the block in the original image
 }
 
-func (r *Registrator) calculateBlockPhaseShift(panBlock, mssBlock gocv.Mat) (gocv.Point2f, float64) {
+func CV_PhaseCorrelate(panBlock, mssBlock gocv.Mat) (gocv.Point2f, float64) {
 	pan := gocv.NewMat()
 	mss := gocv.NewMat()
 
@@ -31,9 +33,9 @@ func (r *Registrator) calculateBlockPhaseShift(panBlock, mssBlock gocv.Mat) (goc
 	mssBlock.ConvertTo(&mss, gocv.MatTypeCV32FC1)
 	defer mss.Close()
 
-	hann := gocv.NewMat()
+	hann := gocv.NewMatWithSize(pan.Rows(), pan.Cols(), pan.Type())
 	defer hann.Close()
-
+	gocv.CreateHanningWindow(&hann, image.Point{X: pan.Cols(), Y: pan.Rows()}, pan.Type())
 	shift, response := gocv.PhaseCorrelate(pan, mss, hann)
 
 	dx := shift.X
@@ -43,7 +45,7 @@ func (r *Registrator) calculateBlockPhaseShift(panBlock, mssBlock gocv.Mat) (goc
 	return shift, response
 }
 
-func (r *Registrator) fileterPhaseShift(thredholds []float64, greaterThan bool) error {
+func (r *ImgProc) fileterPhaseShift(thredholds []float64) error {
 	if len(thredholds) > 4 {
 		return errors.New("thredholds length should be less than 4")
 	}
@@ -55,12 +57,83 @@ func (r *Registrator) fileterPhaseShift(thredholds []float64, greaterThan bool)
 				return false
 			}
 
-			if greaterThan {
-				return value.dy > float32(th)
-			}
-			return value.dy < float32(th)
+			return value.dy > float32(th-20) && value.dy < float32(th+20)
 		})
 	}
 
 	return nil
 }
+
+func PhaseCorrelate(panBlock, mssBlock gocv.Mat) (gocv.Point2f, float64) {
+	// 计算傅里叶变换
+	panFreqDomain, _, _, err := toFreqDomain(panBlock)
+	if err != nil {
+		log.Error(err)
+	}
+
+	mssFreqDomain, width, height, err := toFreqDomain(mssBlock)
+	if err != nil {
+		log.Error(err)
+	}
+
+	// 计算共轭乘积并积累相位信息
+	crossPowerSpectrum := make([][]complex128, height)
+	for i := range crossPowerSpectrum {
+		crossPowerSpectrum[i] = make([]complex128, width)
+	}
+
+	for y := 0; y < height; y++ {
+		for x := 0; x < width; x++ {
+			if panFreqDomain[y][x] != 0 && mssFreqDomain[y][x] != 0 {
+				crossPowerSpectrum[y][x] = panFreqDomain[y][x] * cmplx.Conj(mssFreqDomain[y][x])
+			}
+		}
+	}
+
+	// 归一化
+	for y := 0; y < height; y++ {
+		for x := 0; x < width; x++ {
+			if crossPowerSpectrum[y][x] != 0 {
+				magnitude := cmplx.Abs(crossPowerSpectrum[y][x])
+				crossPowerSpectrum[y][x] = crossPowerSpectrum[y][x] / complex(magnitude, 0)
+			}
+		}
+	}
+
+	ifftResult := fft.IFFT2(crossPowerSpectrum)
+
+	// 查找最大值及其对应的平移参数
+	maxVal := 0.0
+	var dx, dy float64
+	for y := 0; y < height; y++ {
+		for x := 0; x < width; x++ {
+			val := real(ifftResult[y][x])
+			if val > maxVal {
+				maxVal = val
+				dx = float64(x)
+				dy = float64(y)
+			}
+		}
+	}
+	log.Debugf("Block shift: dx = %f, dy = %f. response = %f", dx, dy, 0.0)
+
+	return gocv.Point2f{X: float32(dx), Y: float32(dy)}, 0.0
+}
+
+func toFreqDomain(input gocv.Mat) ([][]complex128, int, int, error) {
+	height := input.Rows()
+	width := input.Cols()
+
+	data := make([][]complex128, height)
+	for y := 0; y < height; y++ {
+		data[y] = make([]complex128, width)
+		for x := 0; x < width; x++ {
+			grayColor := float64(uint16(input.GetShortAt(y, x)))
+			data[y][x] = complex(grayColor, 0)
+		}
+	}
+
+	freqDomain := fft.FFT2(data)
+
+	return freqDomain, width, height, nil
+}

pkg/producer/rpc.go

@@ -32,7 +32,7 @@ type RPC struct {
 	gridsize       int
 
 	scene       *Scene
-	registrator *Registrator
+	registrator *ImgProc
 }
 
 // GroundPoint 表示地面点的三维坐标
@@ -53,7 +53,7 @@ type RPCModel struct {
 }
 
 // rational polynomial coeffients
-func NewRPC(r *Registrator, scene *Scene) *RPC {
+func NewRPC(r *ImgProc, scene *Scene) *RPC {
 	rpc := RPC{
 		elevationLayer: config.GCONFIG.RPC.AltitudeLayer,
 		gridsize:       config.GCONFIG.RPC.GridSize,

pkg/producer/rrc.go

@@ -9,7 +9,7 @@ import (
 	"starwiz.cn/sjy01/image-proc/pkg/rrc"
 )
 
-func (r *Registrator) DoRRCbyLUT(lutDir string) error {
+func (r *ImgProc) DoRRCbyLUT(lutDir string) error {
 	logrus.Printf("try to do RRC [%s]...", lutDir)
 	lutPAN, err := rrc.LoadLUT(filepath.Join(lutDir, "B0.LUT"), 9344)
 	if err != nil {
@@ -42,7 +42,7 @@ func (r *Registrator) DoRRCbyLUT(lutDir string) error {
 	return nil
 }
 
-func (r *Registrator) DoMomentMatching() error {
+func (r *ImgProc) DoMomentMatching() error {
 	rrc.DoMomentMatching(r.PanImage)
 	for i := 0; i < 4; i++ {
 		rrc.DoMomentMatching(r.MssImages[i])

pkg/producer/scene_register.go

@@ -0,0 +1,36 @@
+package producer
+
+import (
+	"github.com/sirupsen/logrus"
+	"gocv.io/x/gocv"
+)
+
+// 分景后再做配准
+func (r *ImgProc) RegisterScenes(pan, mss *Scene) error {
+	proc := NewImgProc(DownSampled)
+	proc.PanImage = pan.Mat[0]
+	proc.MssImages = [4]gocv.Mat{mss.Mat[0], mss.Mat[1], mss.Mat[2], mss.Mat[3]}
+	proc.MssHeight = mss.Height
+	proc.MssWidth = mss.Width
+	proc.PanHeight = pan.Height
+	proc.PanWidth = pan.Width
+
+	err := proc.DoPhaseCorrelation(false)
+	if err != nil {
+		logrus.Error("register scene failed: ", err)
+		return err
+	}
+
+	pan.Height = proc.PanHeight
+	pan.Width = proc.PanWidth
+	pan.Mat[0] = proc.PanImage
+
+	mss.Height = proc.MssHeight
+	mss.Width = proc.MssWidth
+	for i := 0; i < 4; i++ {
+		mss.Mat[i].Close()
+		mss.Mat[i] = proc.registeredMssImages[i]
+	}
+
+	return nil
+}

pkg/producer/scenes.go

@@ -14,6 +14,7 @@ import (
 	"github.com/paulmach/orb/geojson"
 	log "github.com/sirupsen/logrus"
 	"gocv.io/x/gocv"
+	cloudcover "starwiz.cn/sjy01/image-proc/pkg/cloud-cover"
 	"starwiz.cn/sjy01/image-proc/pkg/config"
 	"starwiz.cn/sjy01/image-proc/pkg/fusion"
 	"starwiz.cn/sjy01/image-proc/pkg/payload"
@@ -50,7 +51,7 @@ func CleanScenes(scenes []*Scene) {
 // 默认分景大小：
 // MSS  2336 * 2336 - 1764
 // PAN  9344 * 9344 - 7056
-func (r *Registrator) SubScenes() (panScenes []*Scene, mssScenes []*Scene, err error) {
+func (r *ImgProc) SubScenes() (panScenes []*Scene, mssScenes []*Scene, err error) {
 	if len(r.Params.Targets.Targets) > 0 && r.Params.Targets.Targets[0].EndLine > 0 {
 		return r.RoiScenes()
 	}
@@ -159,7 +160,7 @@ func (r *Registrator) SubScenes() (panScenes []*Scene, mssScenes []*Scene, err e
 	return panScenes, mssScenes, err
 }
 
-func (r *Registrator) OutputL1A(panScenes []*Scene, mssScenes []*Scene) error {
+func (r *ImgProc) OutputL1A(panScenes []*Scene, mssScenes []*Scene) error {
 	var fc geojson.FeatureCollection
 	for i, scene := range panScenes {
 		dir := filepath.Join(r.Params.OutputDir, scene.SceneIndex, "PAN")
@@ -231,11 +232,13 @@ func (r *Registrator) OutputL1A(panScenes []*Scene, mssScenes []*Scene) error {
 			return err
 		}
 
+		jpg := strings.Replace(scene.Tiff, ".tiff", ".jpg", 1)
+		GTiffToJPG(scene.Tiff, jpg, "MSS", false)
+		scene.Meta.CloudRate = cloudcover.CloudPercentByKMeans(jpg)
 		scene.Meta.DataSize = sizeOfFile(scene.Tiff)
 		metaFile := strings.Replace(scene.Tiff, ".tiff", ".meta.xml", 1)
 		writeProductMeta(scene.Meta, metaFile)
-		jpg := strings.Replace(scene.Tiff, ".tiff", ".jpg", 1)
-		GTiffToJPG(scene.Tiff, jpg, "MSS", false)
+
 		r.report.Scenes = append(r.report.Scenes, ReportScene{
 			Name:        scene.SceneId,
 			TiffData:    scene.Tiff,
@@ -247,7 +250,7 @@ func (r *Registrator) OutputL1A(panScenes []*Scene, mssScenes []*Scene) error {
 	return nil
 }
 
-func (r *Registrator) DoScenePansharpen(panScenes []*Scene, mssScenes []*Scene) error {
+func (r *ImgProc) DoScenePansharpen(panScenes []*Scene, mssScenes []*Scene) error {
 	for i := 0; i < len(panScenes); i++ {
 		fusedTiff := strings.Replace(mssScenes[i].Tiff, "MSS", "FUS", -1)
 		err := fusion.Pansharpen(panScenes[i].Tiff, mssScenes[i].Tiff, fusedTiff, fusion.ESRI, 0.1)
@@ -266,7 +269,7 @@ func (r *Registrator) DoScenePansharpen(panScenes []*Scene, mssScenes []*Scene)
 	return nil
 }
 
-func (r *Registrator) MergeMSSToBGRNIR(channels []gocv.Mat) (gocv.Mat, error) {
+func (r *ImgProc) MergeMSSToBGRNIR(channels []gocv.Mat) (gocv.Mat, error) {
 	var rgbirImage gocv.Mat
 	if len(channels) != 4 {
 		return rgbirImage, fmt.Errorf("mss channels count not match")
@@ -281,7 +284,7 @@ func (r *Registrator) MergeMSSToBGRNIR(channels []gocv.Mat) (gocv.Mat, error) {
 
 }
 
-func (r *Registrator) RoiScenes() (panScenes []*Scene, mssScenes []*Scene, err error) {
+func (r *ImgProc) RoiScenes() (panScenes []*Scene, mssScenes []*Scene, err error) {
 	log.Println("using target scenes")
 	for _, target := range r.Params.Targets.Targets {
 		y0 := 4 * target.StartLine

pkg/rrc/hf_filter.go

@@ -1,6 +1,7 @@
 package rrc
 
 import (
+	"fmt"
 	"image"
 	"math"
 
@@ -166,3 +167,32 @@ func visualizeDFT(DFT gocv.Mat, file string) {
 	// 保存DFT结果
 	gocv.IMWrite(file, out)
 }
+
+func HFNoiseFilterSep(img gocv.Mat, radius float64) gocv.Mat {
+	const SepHeight = 4096
+	IDFT := gocv.NewMat()
+	IDFTs := []gocv.Mat{}
+	for rows := 0; rows < img.Rows(); rows += SepHeight {
+		row0 := rows
+		row1 := rows + SepHeight
+		if row1 > img.Rows() {
+			row1 = img.Rows()
+		}
+
+		fmt.Println("doing HF noise fileter rows range:", row0, row1)
+
+		mat := img.Region(image.Rect(0, row0, img.Cols(), row1))
+		matFiltered := HFNoiseFilter(mat, radius)
+		IDFTs = append(IDFTs, matFiltered)
+	}
+
+	gocv.Vconcat(IDFTs[0], IDFTs[1], &IDFT)
+	IDFTs[0].Close()
+	IDFTs[0].Close()
+	for i := 2; i < len(IDFTs); i++ {
+		gocv.Vconcat(IDFT, IDFTs[i], &IDFT)
+		IDFTs[i].Close()
+	}
+
+	return IDFT
+}