Wallis匀色原理：
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 = 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()
贴一下匀色结果：
以上代码是直接对波段灰度值进行Wallis匀色的方法，另有将RGB转为HSV后单独对亮度V进行匀色，得到的新亮度与H/S组合并反算新的RGB。
亮度V的计算公式为（R,G,B为8位无符号整型数据）：
V = max(R/255.,G/255.,B/255.)
计算亮度V的均值、标准差—>计算匀色后的V—>反算成RGB的代码如下：
"""计算亮度V的均值、标准差"""
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)
res_out[1,i,j] = np.std(v1)
res_out[2,i,j] = np.mean(v2)
res_out[3,i,j] = np.std(v2)
del img,ref,v1,v2
import matplotlib
eps = 1e-8
gx = (img_org/255.).astype(np.float)
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)
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 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]
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)
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)
res_out[1,i,j,b] = np.std(img)
res_out[2,i,j,b] = np.mean(ref)
res_out[3,i,j,b] = np.std(ref)
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")
GDAL对亮度V匀色的代码如下（从上片代码中计算mg/sg/mf/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)
res_out[1,i,j] = np.std(v1*mask)
res_out[2,i,j] = np.mean(v2*mask)
res_out[3,i,j] = np.std(v2*mask)
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)
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)
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")
Original: https://blog.csdn.net/qq_33339770/article/details/127938687
Author: 一只大笨猪
Title: Python+Numpy+CV2/GDAL实现对图像的Wallis匀色
原创文章受到原创版权保护。转载请注明出处：https://www.johngo689.com/757916/
转载文章受原作者版权保护。转载请注明原作者出处！