来源：
??????? 那些年“南方测绘杯”科技论文比赛做的工作
先看效果：
实现思路：
??????? 由于遮挡物的存在，会使辐射源（太阳）辐射能量不能到达地面某些区域，这些区域就是遥感影像上的阴影区域，它们通常会和遮挡物有一种伴生关系，而建筑物就是遮挡物之一。因此，建筑物区域附近会存在阴影，且建筑物区域与对应阴影区域的位置关系和太阳方位角有关。鉴于此，本文就利用建筑物与阴影之间的位置关系进行消失建筑物检测，同时获取近似建筑物区域。本文阴影提取方法包括以下四个步骤：
??????? 在可见光遥感影像上，辐射能量中的绝大部分仅由太阳光组成，阴影区域色度应与直接照射时相同，因此，阴影区域与高亮度区域都不会受归一化彩色空间的影响 [1]。于是，我们通过求取影像上原彩色空间与归一化彩色空间的差异来表征阴影特征i；
|
| (1) |
|
| (2) |
|
| (3) |
??????? 由于遮挡物的存在，阴影在光学遥感影像上表现出较低亮度，因此可以用亮度作为特征ii。目前计算亮度的方式有很多，如均值法、最大值法、转换到HSV空间等，但这些方法都没有考虑到人眼对RGB三种光的敏感程度，因此根据人眼对RGB三种光的不同敏感度设计了以下亮度计算方式[2]
|
| (4) |
??????? 由于植被区域树叶之间缝隙的存在，会在植被区域形成斑点状阴影区域，因此需要构造植被区域特征将这些斑点去掉以得到较纯净的阴影。在可见光范围内，植被一般呈绿色，因此可以以绿光波段减去红光波段和蓝光波段的最小值作为特征iii[3]
| ? | (5) |
|
| (6) |
??????? 综合决策部分，经过上述分析，对于特征i，阴影区域和高亮区域呈现较大的值；对于特征ii，阴影区域和植被区域呈现较小的值、高亮区域呈现较大的值；对于特征iii，植被呈现出较大的值、其它部分呈现出较小的值。因此，可以构造以下公式作为最终决策项（三个特征在计算前均已归化到[0，1]），再利用大津阈值分割即可得到较好的阴影提取效果。
|
| (7) |
|
| (8) |
其中α、β、λ为三项特征分别对应的权值，可根据实际情况进行调节。
源代码-matlab：
%yangzhen
%2019.5.11
clear;
clc;
close all
%%
pathname = 'data/';
filenames = dir(pathname);
for i=3:length(filenames)
fileID = [pathname, filenames(i).name];
img = imread(fileID);
img = im2double(img);
%% 计算色度归一化
idist = GetColor(img);
%% 根据人眼视觉计算亮度
ilight = GetLight(img);
%% 计算植被维特征
ivege = GetVege(img);
%% 综合决策
final = GetL_D_V(idist, ilight, ivege);
%% 结果后处理
shadow = FinalTrate(final);
% shadow=final;
%% 定义可视化颜色
icolormap = [0,0,0; 76,180,231]/255;
%% 可视化
figure;
colormap(icolormap);
imshow(img);
hold on
h = imagesc(shadow);
set(h,'alphadata',shadow)
title('阴影提取结果', 'Interpreter', 'None', 'FontSize',12,'FontWeight','bold');
colorbar('horiz', 'Ticks',[0.25, 0.75],...
'TickLabels',{'背景','阴影'}, ...
'FontSize',10,'FontWeight','bold');
%% 保存结果
% save('shadow.mat', 'shadow');
end
%% 标准化函数
function result = standard( data )
[irow, icol, idim] = size(data);
data = reshape(data, [irow*icol, idim]);
temp1 = bsxfun(@minus, data, min(data));
result = bsxfun(@rdivide, temp1, max(data) - min(data));
result = reshape(result, [irow, icol,idim]);
end
%% 根据人眼视觉特性计算影像亮度
function result = GetLight(img)
R = standard( img(:, :, 1) );
G = standard( img(:, :, 2) );
B = standard( img(:, :, 3) );
result = 0.04*R+0.5*G+0.46*B;
end
%% 色度空间归一化
function result = GetColor(img)
img1 = img;
misc = (img(:, :, 1) + img(:, :, 2) + img(:, :, 3));
misc(misc == 0) = 0.0000001;
img(:, :, 1) = img(:, :, 1) ./ misc;
img(:, :, 2) = img(:, :, 2) ./ misc;
result = imabsdiff(img, img1);
result = mean(result(:, :, 1:2), 3);
end
%% 获取植被特征
function result = GetVege(img)
R = standard( img(:, :, 1) );
G = standard( img(:, :, 2) );
B = standard( img(:, :, 3) );
result = (G-min(R, B));
result( result<0 ) = 0;
end
%% 总决策
function result = GetL_D_V(idist, ilight, ivege)
idist = standard(idist);
ilight = standard(ilight);
ivege = standard(ivege);
result = (idist-ilight-ivege);
result( result<0 ) = 0;
end
%% 结果后处理,这里有些参数是可以根据实际情况进行调节的
function result = FinalTrate(img)
T = graythresh(img);
result = imbinarize(img, T);
%中值滤波
result=medfilt2(result, [7, 7]);
end源代码-python：
#yangzhen
#2020.4.13
"""get the shadow proportion form images
of remote sensing"""
import numpy as np
import cv2
import os
import glob
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.colors import LinearSegmentedColormap
from pylab import mpl
import random
mpl.rcParams['font.sans-serif'] = ['FangSong']
mpl.rcParams['axes.unicode_minus'] = False
def standard(data):
'''影像文件标准化
输入单通道影像
输出标准化后单通道影像'''
mdata = data.copy()
irow, icol = mdata.shape[0:2]
mdata = np.reshape(mdata, [irow*icol, 1])
temp1 = mdata - np.min(data)
result = temp1/(np.max(data)-np.min(data))
result = np.reshape(result, [irow, icol])
return result
def GetLight(img):
'''计算人眼视觉特性亮度'''
mimg = img.copy()
B = mimg[:,:,0]
G = mimg[:,:,1]
R = mimg[:,:,2]
result = 0.04*R+0.5*G+0.46*B
return result
def GetColor(img):
'''色度空间归一化'''
mimg = img.copy()
misc = mimg[:,:,0]+mimg[:,:,1]+mimg[:,:,2]
misc[misc == 0] = 0.0000001
mimg[:,:,0] = img[:,:,0]/misc
mimg[:,:,1] = img[:,:,1]/misc
result = np.abs(mimg - img)
result = (result[:,:,0]+result[:,:,1])/2
return result
def GetVege(img):
'''获取植被特征'''
mimg = img.copy()
B = mimg[:,:,0]
G = mimg[:,:,1]
R = mimg[:,:,2]
result = G-np.minimum(R, B)
result[result<0] = 0
return result
def GetLDV(idist, ilight, ivege):
'''总决策'''
idist = standard(idist)
ilight = standard(ilight)
ivege = standard(ivege)
result = idist-ilight-ivege
result[result<0]=0
return result
def FinalTrare(img):
'''结果后处理'''
mimg = img.copy()
mimg = np.uint8(standard(mimg)*255)
T, result = cv2.threshold(mimg, 0, 255, cv2.THRESH_OTSU)
result = cv2.medianBlur(result, 7)
return result
if __name__ == "__main__":
#获取输入图片路径
filepath = 'data/'
filenames = glob.glob(filepath + '*')
print(filenames)
for filename in filenames:
img = cv2.imread(filename)
# 获取阴影
img1 = img.copy()
img1 = img1.astype(np.float)
img1[:,:,0] = standard(img1[:,:,0])
img1[:,:,1] = standard(img1[:,:,1])
img1[:,:,2] = standard(img1[:,:,2])
idist = GetColor(img1)
ilight = GetLight(img1)
ivege = GetVege(img1)
final = GetLDV(idist, ilight, ivege)
shadow = FinalTrare(final)
# 可视化
color_list = ['#00000000', '#4cb4e7']
my_cmap = LinearSegmentedColormap.from_list('mcmp', color_list)
cm.register_cmap(cmap=my_cmap)
fig = plt.figure()
plt.title('阴影提取结果', fontsize=12, fontweight='bold')
plt.imshow(img)
plt.imshow(shadow, cmap='mcmp')
plt.show()参考文献：
- Xu L Q, Landabaso, J.L., Pardas, M.Shadow removal with blob-based morphological reconstruction for error correction[C]. Acoustics, Speech, and Signal Processing, 2005. Proceedings. (ICASSP '05). IEEE International Conference on
- 柳稼航.基于视觉特征的高分辨率光学遥感影像目标识别与提取技术研究[D]. 上海交通大学, 2011.
- 沈小乐.视觉注意机制下面向对象高分辨率遥感影像建筑物提取[D]. 武汉大学,2014.
测试数据：
链接：https://pan.baidu.com/s/1Mtew7ESJLFeqeGmvrOwvhQ?
提取码：1234