Python+Numpy+CV2/GDAL实现对图像的Wallis匀色

Wallis匀色原理：

# f(x,y)：Wallis匀色后结果
# g(x,y):输入的待匀色影像
# mg:待处理影像的灰度均值
# mf:参考影像的灰度均值
# sg:待处理影像和的标准偏差        
# sf:参考影像的标准偏差
f(x,y)=(g(x,y)−mg)⋅(sf/sg)+mf
1
2
3
4
5
6
7

匀色代码逻辑解释：

（1）使用变异系数计算影像的分块数目；
（2）分块计算各块的均值、标准差；
（3）均值、标准差图重采样（双线性）成与输入影像相同行列数；
（4）代入Wallis匀色计算公式计算匀色后的图像数组并保存结果。

代码使用注意：

（1）输入影像与参考影像一定要行列数一致，后面采用GDAL的算法做了重采样，但是GDAL重采样要求输入的影像一定要有坐标；
（2）代码里给Wallis匀色后的值取了绝对值，因为保存成8bit的时候一些负值变成255了；
（3）处理的数据必须是8位的，输出也被固定成8位的了（band_i.astype(np.uint)），如果输入别的位深的数据需要修改一下输出时的数值转换。

算法脚本：

cv2进行Wallis匀色处理的代码

cv2适合对没有坐标、数据量小的图片进行处理，带坐标且数据量极大的卫星影像等往下看GDAL的算法。

"""Wallis匀光——cv"""
import cv2
import numpy as np
from osgeo import gdal
import matplotlib.pyplot as plt
org_file = r"输入影像.tif"
ref_file = r"参考影像.tif"
img_org = cv2.imread(org_file)
infer_img = cv2.imread(ref_file)
width,height,bands = img_org.shape
# 将影像分块进行处理
# 计算变异系数
cv_org = np.std(img_org)/np.mean(img_org)
cv_ref = np.std(infer_img)/np.mean(infer_img)
r = cv_org/cv_ref
num = int(np.ceil(8*r))
w = int(np.ceil(width/num))
h = int(np.ceil(height/num))
# mg:待处理影像的灰度均值
# mf:参考影像的灰度均值
# sg:待处理影像和的标准偏差        
# sf:参考影像的标准偏差 
mg = np.zeros((num,num,bands),dtype=np.float)
mf = np.zeros_like(mg)
sg = np.zeros_like(mg)
sf = np.zeros_like(mg)
for b in range(bands):
    for i in range(num):
        for j in range(num):
            orgin_x = min(i*w,width - w)
            orgin_y = min(j*h, height - h)
            end_x = orgin_x + w
            end_y = orgin_y + h
            img = img_org[orgin_x:end_x,orgin_y:end_y,b]
            ref = infer_img[orgin_x:end_x,orgin_y:end_y,b]
            mg[i,j,b] = np.mean(img)
            sg[i,j,b] = np.std(img)
            mf[i,j,b] = np.mean(ref)
            sf[i,j,b] = np.std(ref)

"""Wallis公式:f(x,y)=(g(x,y)−mg)⋅(sf/sg)+mf"""
eps = 1e-8
waillisImg = np.zeros_like(img_org)
for i in range(bands):
    mg_res = cv2.resize(mg[:,:,i],(height,width),interpolation=cv2.INTER_LINEAR)
    mf_res = cv2.resize(mf[:,:,i],(height,width),interpolation=cv2.INTER_LINEAR)
    sf_res = cv2.resize(sf[:,:,i],(height,width),interpolation=cv2.INTER_LINEAR)
    sg_res = cv2.resize(sg[:,:,i],(height,width),interpolation=cv2.INTER_LINEAR)
    band_i = np.abs((img_org[:,:,i] - mg_res) * (sf_res / (sg_res+ eps)) + mf_res)
    waillisImg[:,:,i] = band_i.astype(np.uint)
cv2.imwrite(r"waillis匀色结果.tif"，waillisImg)
plt.subplot(1,3,1)
plt.imshow(img_org)
plt.subplot(1,3,2)
plt.imshow(infer_img)
plt.subplot(1,3,3)
plt.imshow(waillisImg)
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58

贴一下匀色结果：

以上代码是直接对波段灰度值进行Wallis匀色的方法，另有将RGB转为HSV后单独对亮度V进行匀色，得到的新亮度与H/S组合并反算新的RGB。
亮度V的计算公式为（R,G,B为8位无符号整型数据）：

V = max(R/255.,G/255.,B/255.)
1

计算亮度V的均值、标准差—>计算匀色后的V—>反算成RGB的代码如下：

"""计算亮度V的均值、标准差"""
# mg:待处理影像的灰度均值
# mf:参考影像的灰度均值
# sg:待处理影像和的标准偏差        
# sf:参考影像的标准偏差 
res_out = np.zeros((4,num,num),dtype=np.float)
for i in range(num):
    for j in range(num_w):
        orgin_x = min(i*w,width - w)
        orgin_y = min(j*h, height - h)
        end_x = orgin_x + w
        end_y = orgin_y + h
        img = img_org[orgin_x:end_x,orgin_y:end_y,:]
        ref = infer_img[orgin_x:end_x,orgin_y:end_y,:]
        v1 = np.max((img/255.).astype(np.float),axis=-1)
        v2 = np.max((ref/255.).astype(np.float),axis=-1)
        res_out[0,i,j] = np.mean(v1)#mg
        res_out[1,i,j] = np.std(v1)#sg
        res_out[2,i,j] = np.mean(v2)#mf
        res_out[3,i,j] = np.std(v2)#sf
        del img,ref,v1,v2

import matplotlib
eps = 1e-8
gx = (img_org/255.).astype(np.float)
# rgb_to_hsv
hsv = matplotlib.colors.rgb_to_hsv(gx)
h = hsv[:,:,0]
s = hsv[:,:,1]
v = hsv[:,:,2]
# 双线性采样
mg = cv2.resize(res_out[0,:,:],(height,width),interpolation=cv2.INTER_LINEAR)
mf = cv2.resize(res_out[1,:,:],(height,width),interpolation=cv2.INTER_LINEAR)
sf = cv2.resize(res_out[2,:,:],(height,width),interpolation=cv2.INTER_LINEAR)
sg = cv2.resize(res_out[3,:,:],(height,width),interpolation=cv2.INTER_LINEAR)
# 计算匀色亮度
out = np.abs((v - mg) * (sf / (sg+ eps)) + mf)
# hsv_to_rgb
new_hsv = (np.stack((h,s,out))).transpose(1,2,0)
rgb = matplotlib.colors.hsv_to_rgb(new_hsv)
waillisImg = (rgb * 255).astype(np.uint)
cv2.imwrite(r"waillis匀色结果.tif"，waillisImg)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42

GDAL的Wallis匀色算法代码

GDAL的算法就没有办法像上面cv2一样把全图读完计算变异系数了（计算量太大了），采用的是经典的分块处理，将图像分成固定大小的方形切片，计算均值和标准差，并使用gdal.warp进行重采样，后面就是简单的分波段计算、保存与输出了。

"""Wallis匀光——GDAL"""

from osgeo import gdal,gdalconst
import numpy as np

org_file = r"输入影像.tif"
ref_file = r"参考影像.tif"

raster = gdal.Open(org_file)
rows = raster.RasterYSize
cols = raster.RasterXSize
bands = raster.RasterCount
print(cols,rows,bands)
OriginX,psX,_,OriginY,_,psY = raster.GetGeoTransform()
EndX = OriginX + cols * psX
EndY = OriginY + rows * psY
extent = [OriginX,EndY,EndX,OriginY]
# 分块大小定义为512
bk_size = 512
num_w = int(np.ceil(cols / bk_size))
num_h = int(np.ceil(rows/ bk_size))
print(num_w,num_h)
ref_raster = gdal.Open(ref_file)
# 输入影像对齐（对参考影像重采样成相同行列数）
if ref_raster.RasterXSize != cols and ref_raster.RasterYSize != rows:
    new_ref = ref_file[0:-4]+"_resample.tif"
    warp_ds = gdal.Warp(new_ref,ref_file,width = cols,height = rows) 
    warp_ds = None
    ref_raster = gdal.Open(new_ref)
# 计算Wallis需要的参数
# mg:待处理影像的灰度均值
# mf:参考影像的灰度均值
# sg:待处理影像和的标准偏差        
# sf:参考影像的标准偏差 
res_out = np.zeros((4,num_h,num_w,bands),dtype=np.float)
for b in range(bands):
    img_band = raster.GetRasterBand(b+1)
    ref_img = ref_raster.GetRasterBand(b+1)
    for i in range(num_h):
        for j in range(num_w):
            orgin_x = min(j*bk_size,cols -  bk_size)
            orgin_y = min(i*bk_size,rows - bk_size)
            img = img_band.ReadAsArray(orgin_x,orgin_y, bk_size, bk_size)
            ref = ref_img.ReadAsArray(orgin_x,orgin_y, bk_size, bk_size)
            res_out[0,i,j,b] = np.mean(img)#mg
            res_out[1,i,j,b] = np.std(img)#sg
            res_out[2,i,j,b] = np.mean(ref)#mf
            res_out[3,i,j,b] = np.std(ref)#sf

# 重采样
# 对输入影像重采样是为了获得采样后的Projection和GeoTransform,用来赋给mg/sg/mf/sf进行上采样
# 测试中发现gdal.warp无法对没有坐标的图像重采样
outimg = r"输入影像重采样.tif"
warp_ds = gdal.Warp(outimg,org_file,resampleAlg=gdalconst.GRA_Average,width = num_w,height = num_h) 
del warp_ds
temp_ref = gdal.Open(outimg)
in_list = []
for i in range(4):
    driver = gdal.GetDriverByName("GTiff")
    temp_out = r"temp%d.tif" % i # 过程数据 完成后可以删除
    temp_ds = driver.Create(temp_out,num_w,num_h,bands,gdal.GDT_Float32)
    temp_ds.SetGeoTransform(temp_ref.GetGeoTransform())
    temp_ds.SetProjection(temp_ref.GetProjection())
    for tb in range(bands):
        temp_ds.GetRasterBand(tb+1).WriteArray(res_out[i,:,:,tb])
    temp_ds.FlushCache()
    del temp_ds
    temp_res = r"temp_res%d.tif" % i # 过程数据 完成后可以删除
    warp_ds = gdal.Warp(temp_res,temp_out,resampleAlg=gdalconst.GRA_Bilinear,outputBounds = extent,xRes = psX,yRes =psY,targetAlignedPixels=True)                         
    del warp_ds
    in_raster = gdal.Open(temp_res)
    in_list.append(in_raster)

# 输出匀色结果
eps = 1e-8
[mean_raster,std_raster,mean_ref_raster,std_ref_raster] = in_list
driver = gdal.GetDriverByName("GTiff")
out_ds= driver.Create(r"Wallis匀色结果.tif",cols,rows,bands,gdal.GDT_Byte)
out_ds.SetGeoTransform(raster.GetGeoTransform())
out_ds.SetProjection(raster.GetProjection())
for b in range(bands):
    # 输入影像
    in_band = raster.GetRasterBand(b+1)
    mean_band = mean_raster.GetRasterBand(b+1)
    std_band = std_raster.GetRasterBand(b+1)
    # 参考影像
    mean_ref_band = mean_ref_raster.GetRasterBand(b+1)
    std_ref_band = std_ref_raster.GetRasterBand(b+1)
    # 输出影像
    out_band = out_ds.GetRasterBand(b+1)
    # 分块处理
    for i in range(num_h):
        for j in range(num_w):
            orgin_x = min(j*bk_size,cols -  bk_size)
            orgin_y = min(i*bk_size,rows - bk_size)
            # 读取输入参数
            gx = in_band.ReadAsArray(orgin_x, orgin_y, bk_size, bk_size)
            mg = mean_band.ReadAsArray(orgin_x, orgin_y, bk_size, bk_size)
            sg = std_band.ReadAsArray(orgin_x, orgin_y, bk_size, bk_size)
            mf = mean_ref_band.ReadAsArray(orgin_x, orgin_y, bk_size, bk_size)
            sf = std_ref_band.ReadAsArray(orgin_x, orgin_y, bk_size, bk_size)
            # 计算匀色影像
            wallis = np.abs((gx - mg) * (sf / (sg+ eps)) + mf)
            # 保存匀色结果
            out_band.WriteArray(wallis.astype(np.uint),orgin_x,orgin_y)
    out_band.FlushCache()
del out_band
out_ds.FlushCache()
del out_ds
print("Done")
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110

GDAL对亮度V匀色的代码如下（从上片代码中计算mg/sg/mf/sf的部分替换即可）：

# mg:待处理影像的灰度均值
# mf:参考影像的灰度均值
# sg:待处理影像和的标准偏差        
# sf:参考影像的标准偏差 
res_out = np.zeros((4,num_h,num_w),dtype=np.float)
for i in range(num_h):
    for j in range(num_w):
        orgin_x = min(j*bk_size,cols -  bk_size)
        orgin_y = min(i*bk_size,rows - bk_size)
        img = raster.ReadAsArray(orgin_x,orgin_y, bk_size, bk_size)
        img = (img.transpose(1,2,0))[:,:,0:3]
        ref = ref_raster.ReadAsArray(orgin_x,orgin_y, bk_size, bk_size)
        ref = (ref.transpose(1,2,0))[:,:,0:3]
        v1 = np.max((img/255.).astype(np.float),axis=-1)
        v2 = np.max((ref/255.).astype(np.float),axis=-1)
        if img.max() > 0:
            mask = ((v1 > 0) * (v2 > 0)).astype(np.uint)
            res_out[0,i,j] = np.mean(v1*mask)#mg
            res_out[1,i,j] = np.std(v1*mask)#sg
            res_out[2,i,j] = np.mean(v2*mask)#mf
            res_out[3,i,j] = np.std(v2*mask)#sf
        del img,ref,v1,v2
del ref_raster

# 重采样
outimg = r"img_resample.tif"
warp_ds = gdal.Warp(outimg,org_file,resampleAlg=gdalconst.GRA_Average,width = num_w,height = num_h) 
del warp_ds
print(outimg)
temp_ref = gdal.Open(outimg)
temp_geotrans = temp_ref.GetGeoTransform()
temp_proj = temp_ref.GetProjection()
del temp_ref

driver = gdal.GetDriverByName("GTiff")
temp_out = r"temp00.tif"
temp_ds = driver.Create(temp_out,num_w,num_h,4,gdal.GDT_Float32)
temp_ds.SetGeoTransform(temp_geotrans)
temp_ds.SetProjection(temp_proj)
for tb in range(4):
    temp_ds.GetRasterBand(tb+1).WriteArray(res_out[tb,:,:])
temp_ds.FlushCache()
del temp_ds
print(temp_out)
temp_res = r"temp_res00.tif"
warp_ds = gdal.Warp(temp_res,temp_out,resampleAlg=gdalconst.GRA_Bilinear,outputBounds = extent,
                    xRes = psX,yRes =psY,targetAlignedPixels=True) 
                    #width = cols,height = rows,                    
del warp_ds
print(temp_res)
in_raster = gdal.Open(temp_res)
# 输出匀色结果
eps = 1e-8
driver = gdal.GetDriverByName("GTiff")
out_ds= driver.Create(out_wallis_file,cols,rows,3,gdal.GDT_Byte)
out_ds.SetGeoTransform(raster.GetGeoTransform())
out_ds.SetProjection(raster.GetProjection())
# 分块处理
for i in range(num_h):
    for j in range(num_w):
        orgin_x = min(j*bk_size,cols -  bk_size)
        orgin_y = min(i*bk_size,rows - bk_size)
        # 读取输入影像
        gx = raster.ReadAsArray(orgin_x, orgin_y, bk_size, bk_size)
        gx = (gx.transpose(1,2,0))[:,:,0:3]
        hsv = matplotlib.colors.rgb_to_hsv((gx/255.).astype(np.float))
        h = hsv[:,:,0]
        s = hsv[:,:,1]
        v = hsv[:,:,2]
        mg = in_raster.GetRasterBand(1).ReadAsArray(orgin_x, orgin_y, bk_size, bk_size)
        sg = in_raster.GetRasterBand(2).ReadAsArray(orgin_x, orgin_y, bk_size, bk_size)
        # 读取参考影像
        mf = in_raster.GetRasterBand(3).ReadAsArray(orgin_x, orgin_y, bk_size, bk_size)
        sf = in_raster.GetRasterBand(4).ReadAsArray(orgin_x, orgin_y, bk_size, bk_size)
        # 计算匀色影像
        out = np.abs((v - mg) * (sf / (sg+ eps)) + mf)
        # hsv_to_rgb
        new_hsv = (np.stack((h,s,out))).transpose(1,2,0)
        rgb = matplotlib.colors.hsv_to_rgb(new_hsv)
        rgb = (rgb  * 255).astype(np.uint)
        for b in range(3):
            # 输出影像
            out_band = out_ds.GetRasterBand(b+1)
            # 保存匀色结果
            out_band.WriteArray(rgb[:,:,b],orgin_x,orgin_y)
            out_band.FlushCache()
        del mask,hsv,h,s,v,mg,sg,mf,sf,out,new_hsv,rgb
        del gx
del out_band
out_ds.FlushCache()
del out_ds
print("Done")
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92