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）代入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 + wend_y = orgin_y + himg = 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()

贴一下匀色结果：

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

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

计算亮度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 + wend_y = orgin_y + himg = 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)#mgres_out[1,i,j] = np.std(v1)#sgres_out[2,i,j] = np.mean(v2)#mfres_out[3,i,j] = np.std(v2)#sfdel img,ref,v1,v2import 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)

GDAL的Wallis匀色算法代码

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

"""Wallis匀光——GDAL"""from osgeo import gdal,gdalconst
import numpy as nporg_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 = Noneref_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)#mgres_out[1,i,j,b] = np.std(img)#sgres_out[2,i,j,b] = np.mean(ref)#mfres_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_dstemp_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_dsin_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")

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)#mgres_out[1,i,j] = np.std(v1*mask)#sgres_out[2,i,j] = np.mean(v2*mask)#mfres_out[3,i,j] = np.std(v2*mask)#sfdel 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_refdriver = 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_rgbnew_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,rgbdel gx
del out_band
out_ds.FlushCache()
del out_ds
print("Done")