使用Python进行图像处理

 来源：磐创AI

下面是一个关于使用Python在几行代码中分析城市轮廓线的快速教程。

说一句显而易见的话：轮廓线很美。

在本文中，我们将学习如何从图片中获取轮廓线轮廓。

 

让我们开始吧。

0.理念

这个想法很简单。为了检测轮廓线，我们只检测天空并拍摄互补图像。

在你之前看到的示例中，我们真正做的是识别天空。下一步当然是获取蒙版图像。

那么，为什么探测天空比探测摩天大楼更容易呢？

好吧，这个概念是天空的图片是相对平坦的。另一方面，摩天大楼是颜色、形状、窗户、水泥等的混合体。

从数学上讲，天空的方差比摩天大楼的方差小，并且期望该参数在区分天空和摩天大楼时起决定性作用。

我想用一个例子来证明我所说的话。让我们拍下这张照片。

 

太棒了现在，让我们在这个区域修剪这张图片：

 

现在，让我们从0到50之间取一部分并打印标准偏差：

 

该方差变化可以使用二阶导数来检测。

当我们讨论离散二维情况时，我们实际上是在讨论拉普拉斯算子。拉普拉斯算子可以被视为卷积，这只是使用泰勒近似的导数的定义。

进入下一节的主题。

1.算法

执行图（1）中所示操作的算法如下：

我明白这是一团糟。让我一步一步地解释一下。

1.1将图像转换为黑白

我知道你知道这一切。但重要的是要说明我们为什么要这样做。正如你所知，当你将它们应用于矩阵时，所有的模糊步骤和过滤都是有意义的。彩色图像在技术上是一个张量，因为它具有行数X列数X 3个通道值（红、绿和蓝）。B&W图像是由行数X列数组成的矩阵。

将此应用于彩色图像的一个简单方法是重复上述相同的过程三次，但我认为没有必要。最终，即使使用B&W图像，我们也能分辨出轮廓线。

1.2模糊步骤

中值和归一化滤波器步骤都是用于在保持边的同时对信号的噪声进行滤波的步骤。

1.3拉普拉斯滤波器

拉普拉斯滤波器被认为是离散空间的二阶时间导数。

为什么我们首先需要二阶时间导数？

我们说过，天空和摩天大楼之间的标准差是不同的。这种标准差的变化发生在一个特定的点上，即图像（和摩天大楼）的边缘。

所以我们希望看到图像的快速变化。特别是，我们希望变化最大。这意味着我们需要二阶导数为空的点（或点的邻居）。

当我们讨论离散二维情况时，我们实际上是在讨论拉普拉斯算子。拉普拉斯算子可以被视为卷积，这只是使用泰勒近似的导数的定义。

二阶导数是这样的：

这是一个核，我们将在图像上运行，它将为我们提供二阶导数图像。

1.4应用1/0阈值

我们不关心二阶导数是正还是负。我们所关心的只是我们有0的一小部分，因为这是我们认为的边缘。

这就是为什么我们应用这个1/0阈值。

1.5侵蚀滤波器

侵蚀滤波器是我们用来平滑图像的东西。这背后的想法是，我们希望使图像更清晰。用更专业的话来说，有一个核在图像上传递，并用它们的最小值替换值。同样，由于我们现在有一张1/0的图像，它只是让我们的图像更清晰。

1.6将掩码设置为0，直到找到最后一个索引

这一步有点难解释，但很容易理解。完成所有这些操作后，图像的一列中可能有一个0和1的序列。这没有太多意义，因为你不能再拥有“skyscraper-sky-skyscraper again”这样的东西了。因此，我们在列中找到值为0的最大索引，并将所有值设置为0，直到找到该值。那么，其他的都是0。

2.实际实现

这解释起来有点长，但非常容易实现。

让我们循序渐进：

 

import numpy as np 
import pandas as pd 
from os import listdir
from PIL import Image
import matplotlib.pyplot as plt
from os.path import isfile, join
import cv2
from scipy.signal import medfilt
from scipy import ndimage
import numpy as np
from matplotlib import pyplot as plt
import os
from sklearn.metrics import plot_confusion_matrix
from sklearn.metrics import confusion_matrix

 

2.1导入库：

 

mypath = 'data/images'
subfolders = [ f.path for f in os.scandir(mypath) if f.is_dir() ]
images = []
images_baw = []
labels = []
label=0
size = 256
string_labels = []
for folder in subfolders: 
    onlyfiles = [f for f in listdir(folder) if isfile(join(folder, f))]
    for file in onlyfiles:
        image_file = Image.open(folder+'/'+file).resize((size,size))
        images.append(np.array(image_file))
        images_baw.append(np.array(image_file.convert('1')))
        labels.append(label)
        string_labels.append(folder.split('/')[-1])
    label=label+1
labels = np.array(labels)
images = np.array(images)
image_baw = np.array(images_baw)

 

2.1导入数据：

我从Kaggle那里得到了数据。数据集是开源的，没有版权（CC0：公共域）。特别是，我只下载了数据集的一部分，其中包含12个城市的图像，每个城市有10座摩天大楼：

 

mypath = 'data/images'
subfolders = [ f.path for f in os.scandir(mypath) if f.is_dir() ]
images = []
images_baw = []
labels = []
label=0
size = 256
string_labels = []
for folder in subfolders: 
    onlyfiles = [f for f in listdir(folder) if isfile(join(folder, f))]
    for file in onlyfiles:
        image_file = Image.open(folder+'/'+file).resize((size,size))
        images.append(np.array(image_file))
        images_baw.append(np.array(image_file.convert('1')))
        labels.append(label)
        string_labels.append(folder.split('/')[-1])
    label=label+1
labels = np.array(labels)
images = np.array(images)
image_baw = np.array(images_baw)

 

2.2数据可视化

 

plt.figure(figsize=(20,20))
for i in range(1,13):
    plt.subplot(4,3,i)
    identify_label = np.where(labels==i-1)[0]
    identify_label = np.random.choice(identify_label)
    plt.title('City label = %s'%(string_labels[identify_label]))
    plt.imshow(images[identify_label])

 

2.3定义函数：

以上所有理论都以以下方式实现：

 

def cal_skyline(mask):
    h, w = mask.shape
    for i in range(w):
        raw = mask[:, i]
        after_median = medfilt(raw, 19)
        try:
            first_zero_index = np.where(after_median == 0)[0][0]
            first_one_index = np.where(after_median == 1)[0][0]
            if first_zero_index > 20:
                mask[first_one_index:first_zero_index, i] = 1
                mask[first_zero_index:, i] = 0
                mask[:first_one_index, i] = 0
        except:
            continue
    return mask


def get_sky_region_gradient(img):

    h, w, _ = img.shape

    img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

    img_gray = cv2.blur(img_gray, (9, 3))
    cv2.medianBlur(img_gray, 5)
    lap = cv2.Laplacian(img_gray, cv2.CV_8U)
    gradient_mask = (lap < 6).astype(np.uint8)

    kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (9, 3))

    mask = cv2.morphologyEx(gradient_mask, cv2.MORPH_ERODE, kernel)
    # plt.imshow(mask)
    # plt.show()
    mask = cal_skyline(mask)
    after_img = cv2.bitwise_and(img, img, mask=mask)

    return after_img

 

2.4应用算法

这是算法在整个数据集上的应用：

 

zero_images = []
zero_images_plot = []
for K in range(len(images)):
    zero_image_plot = np.zeros((size,size))
    zero_image = zero_image_plot.copy()
    for i in range(size):
        zero_image_plot[i,tot_profiles[K][i]-2:tot_profiles[K][i]+2] = 1
        zero_image[i,tot_profiles[K][i]] = 1 
    zero_images_plot.append(zero_image_plot.T)
    zero_images.append(zero_image.T)

 

示例：

 

K=0
plt.figure(figsize=(10,10))
plt.suptitle('City = %s'%(string_labels[K]),fontsize=20,y=0.72)
plt.subplot(1,3,1)
plt.title('Original Image')
plt.imshow(images[K])
plt.subplot(1,3,3)
plt.title('Masked Image')
plt.imshow(zero_images_plot[0])
plt.subplot(1,3,2)
plt.title('Profiled Skyline')
plt.imshow(get_sky_region_gradient(images[K]))

 

2.5将其转换为信号

为了获得图像（1）的正确图像，我们应用以下函数：

 

def signal_from_profile(K):
    x = np.arange(size)
    y = np.array(size-np.array(tot_profiles[K]))
    return x,y

 

我们可以将其应用于数据集的所有图像：

 

plt.figure(figsize=(13,10))
#plt.suptitle('City = %s'%(string_labels[K]),fontsize=20)
i=0
for q in range(3):
    K = np.random.choice(len(images))
    #print(K)
    skyline_signal_x,skyline_signal_y =signal_from_profile(K)
    plt.subplot(4,4,i+1)
    plt.title('Original Image')
    plt.imshow(images[K])
    plt.subplot(4,4,i+3)
    plt.title('Masked Image')
    plt.imshow(zero_images_plot[K])
    plt.subplot(4,4,i+2)
    plt.title('Profiled Skyline')
    plt.imshow(get_sky_region_gradient(images[K]))
    plt.subplot(4,4,i+4)
    plt.title('Extracted Signal')
    plt.plot(skyline_signal_x,skyline_signal_y)
    plt.tight_layout()
    i=i+4

 

3.结尾

我认为这项研究之所以有趣，有多种原因。首先，这很有趣，因为有两个理论上合理的理由。

它解释了如何使用拉普拉斯滤波器以非深度学习的方式应用边缘检测

它解释了如何使用图像进行从头到脚的实验，以及如何创建一个有效的图像处理管道

当然，这本身很有趣，因为它为你提供了一个分析不同城市轮廓线的工具！

你可以看到，城市A和城市B有不同的概况，特别是使用提取的信号，我们可以通过以下方式深化这项研究：

提取轮廓线的平均值、中值和标准差

使用深度学习对城市轮廓线进行分类

对轮廓线与时间进行统计研究（轮廓线如何随时间演变？）

记住，这个项目背后的整个想法是，天空的标准差比摩天大楼的标准差低。

我们还可以使用这种方法作为更复杂研究的起点，并且可以使用编码器-解码器来改进这些结果。