图像校正python 图像校正处理

　　本节将从原理和代码两个方面讲解遥感图像的几何校正。

原理

　　首先介绍几何校正的概念：在遥感成像过程中，传感器生成的图像像元相对于地面目标物的实际位置发生了挤压、扭曲、拉伸和偏移等问题，这一现象叫做几何畸变。几何畸变会给遥感图像的定量分析、变化检测、图像融合、地图测量或更新等处理带来的很大误差，所以需要针对图像的几何畸变进行校正，即几何校正。

　　几何校正分为几何粗校正和几何精校正。粗校正是利用空间位置变化关系，采用计算公式和辅助参数进行的校正，叫做系统几何校正；精校正是在此基础上，使图像的几何位置符合某种地理坐标系统，与地图配准，调整亮度值，即利用地面控制点（GCP）做的几何精校正。

　　几何校正步骤：1.空间位置的变换（像元坐标）2.像元灰度值的重新计算，即重采样。

　　1. 坐标变换

　　坐标变换分为直接法和间接法。

　　1）直接法：从原始图像阵列出发，依次计算每个像元在输出图像中的坐标。直接法输出的像元值大小不会发生变化，但输出图像中的像元分布不均匀。

　　2）间接法：从输出图像阵列出发，依次计算每个像元在原始图像中的位置，然后计算原始图像在该位置的像元值，再将计算的像元值赋予输出图像像元。此方法保证校正后的图像的像元在空间上均匀分布，但需要进行灰度重采样。该方法是最常用的几何校正方法。

 

　　由上图可见，直接法直接以原始图像的坐标为基准点，坐标偏移到校正后的图像，坐标的位置有很多出现在了像元的中间位置，所以直接输出像元值大小导致像元分布不均匀。

　　而对于间接法。以输出图像的坐标为基准点，已经定义在了格点的位置上，此时反算出该点在原始图像上对应的图像坐标，坐标多数落在像元的中间位置。这里采用最邻近法、双线性内插和三次卷积法来计算该点的灰度值，达成重采样的目的。

 

　　2. 重采样

　　　图像数据经过坐标变换之后，像元中心的位置发生改变，其在原始图像的位置不一定是整数行\列，需要根据输出图像各像元在原始图像中对应的位置，对原始图像重采样，建立新的栅格矩阵。

　　　主要有最近邻法、双线性内插法、三次卷积法。原理这里不再叙述。其中最近邻法没有改变图像光谱信息，能保证定量分析的准确性，但是图像不够平滑。后两者虽然数据连续平滑了，但是光谱信息遭到了破坏，比较适合于制图，不适用于定量分析。

 

 

 代码实现

 　　　　在熟悉了几何校正的原理后，可以动手实现这一环节。这里主要介绍三种方法，具体如下。

　　　　1.首先介绍最基础也是最底层的方法。(这里是分步切片得到偏移后的坐标，然后重采样）

from osgeo import gdal

def get_indices(source_ds, target_width, target_height):
    source_geotransform = source_ds.GetGeoTransform()
    source_width = source_ds.GetGeoTransform[1]
    source_height = source_ds.GetGeoTransform[5]
    dx = target_width/source_width
    dy = target_height/source_height
    target_x = np.arange(dx/2, source_ds.RasterXSize, dx)
    target_y = np.arange(dy/2, source_ds.RasterYSize, dy)
    return np.meshgrid(target_x, target_y)

def bilinear(in_data, x, y):
    x-=0.5
    y-=0.5
    x0=np.floor(x).astype(int)
    x1=x0+1
    y0=np.floor(y).astype(int)
    y1=y0+1
    
    ul=in_data[y0,x0]*(y1-y)*(x1-x)
    ur=in_data[y0,x1]*(y1-y)*(x-x0)
    ll=in_data[y1,x0]*(y-y0)*(x1-x)
    lr=in_data[y1,x1]*(y-y0)*(x-x0)

    return ul+ur+ll+lr

if __name__ =="__main__":
    in_fn=r''
    out_fn=r''
    cell_size=(0.02,-0.02)

    in_ds = gdal.Open(in_fn)
    x, y=get_indices(in_ds,*cell_size)
    outdata=bilinear(in_ds.ReadAsArray(), x, y)
    
    driver=gdal.GetDriverByName('GTiff')
    rows,cols=outdata.shape
    out_ds=driver.Create(out_fn,cols,rows,1,gdal.GDT_ Int32)

    gt=list(in_ds.GetGeoTransform())
    gt[1]=cell_size[0]
    gt[5]=cell_size[1]
    out_ds.SetGeoTransform(gt)
    
    out_band = out_ds.GetRasterBand(1)
    out_band.WriteArray(outdata)
    out_band.FlushCache()
    out_band.ComputeStatistics(False)

　　2.第二种方法是ReadAsArray()函数，设置buf_xsize, buf_ysize这个缓冲区大小，就可以按这个大小重采样（最邻近法）。同时需要设置新的gt[1]\gt[5]以及行列数，这里代码省略，有兴趣查阅之前的gdal笔记。

gdal.ReprojectImage。这里有重采样的部分，具体如下。

from osgeo import gdal, gdalconst

# 1.设置分辨率投影
# inputfile = r'E:\桌面文件保存路径\GIS\test\quickbird_test.dat'
# outputfile = r'E:\桌面文件保存路径\GIS\test\resampling_test.tif'
#
# inputds = gdal.Open(inputfile, gdal.GA_ReadOnly)  # 待重采样图像
# band1 = inputds.GetRasterBand(1)
# xsize = inputds.RasterXSize
# ysize = inputds.RasterYSize
# nbands = inputds.RasterCount
# inputProj = inputds.GetProjection()
# inputTrans = inputds.GetGeoTransform()
#
#
# new_x = int(xsize * (inputTrans[1] / 10))  # 计算新行列数
# new_y = int(ysize * (inputTrans[5] / -10)) # 注意抵消负号
#
# print(inputTrans[1])
# print(inputTrans[5])
#
# inputTrans = list(inputTrans)
# inputTrans[1] = 10  # 分辨率
# inputTrans[5] = -10 # 注意设置为-，否则无法SetGeoTransform
#
# print(inputTrans)
#
# driver = gdal.GetDriverByName('GTiff')
# out_ds = driver.Create(outputfile, new_x, new_y, nbands, band1.DataType)
# out_ds.SetGeoTransform(inputTrans)
# out_ds.SetProjection(inputProj)
#
# gdal.ReprojectImage(inputds, out_ds, inputProj, out_ds.GetProjection(), gdalconst.GRA_Bilinear)
#
# del out_ds

# 2.从图1重采样到图2

inputfile = r'E:\桌面文件保存路径\IDL\imgs\TMSPOT\TM-30m.img' #高光谱
referencefile = r'E:\桌面文件保存路径\IDL\imgs\TMSPOT\bldr_sp.img' #高分辨率
outputfile = r'E:\test.tif'  #重采样输出图像

inputras = gdal.Open(inputfile, gdal.GA_ReadOnly)  #待重采样图像，高光谱低分辨率
inputProj = inputras.GetProjection()
inputTrans = inputras.GetGeoTransform()
nbands = inputras.RasterCount

reference = gdal.Open(referencefile, gdal.GA_ReadOnly)  #目标图像，高分辨率全色波段
referenceProj = reference.GetProjection()
referenceTrans = reference.GetGeoTransform()
bandreference = reference.GetRasterBand(1)
x = reference.RasterXSize
y = reference.RasterYSize


driver = gdal.GetDriverByName('GTiff')
output = driver.Create(outputfile, x, y, nbands, bandreference.DataType)
output.SetGeoTransform(referenceTrans)
output.SetProjection(referenceProj)

gdal.ReprojectImage(inputras, output, inputProj, referenceProj, gdalconst.GRA_Bilinear)

del output

　　4.此外完整的几何校正和重采样还可以用opencv的一些api。

from matplotlib import pyplot as plt
import cv2
import numpy as np

img = cv2.imread('aier.jpg')
rows,cols = img.shape[:2]
pts1 = np.float32([[50,50],[200,50],[50,200]])  #对应校正点
pts2 = np.float32([[10,100],[200,20],[100,250]])
M = cv2.getAffineTransform(pts1,pts2)  //计算仿射变换矩阵
#第三个参数：变换后的图像大小
res = cv2.warpAffine(img,M,(rows,cols)) //根据矩阵和图像重采样
plt.subplot(121)
plt.imshow(img[:,:,::-1])

plt.subplot(122)
plt.imshow(res[:,:,::-1])  #画出图片，但不显示

plt.show() #显示