自动查找适用于模板匹配的纹理区域

模板匹配是图像处理中常用的精对位手段，但不是所有的区域都适用于模板匹配，本文提出一种启发式的自动搜索纹理区域的方法。

实现目标

1. 初始化配置完成后仅需要输入图像
2. 输出最适用于模板匹配的纹理矩形区域
3. 大小可变
4. 速度不能太慢（几秒内）

核心思想

整体路线

为了在全图范围内找到适合模板匹配的纹理区域，需要在全图范围内的指定尺寸窗口中计算是否适合模板匹配的分值

难点主要在于给定一个窗口如何评价这个窗口内的图像是否适合用于模板匹配，另一个问题是对每个窗口都做复杂的评估需要巨大运算量，如果兼顾不同尺寸的窗口，可能需要很长时间计算，需要在效果和运算量上做取舍和平衡。

基本思想

将需求拆分成三个维度：

1. 目标区域需要和周围图像像素差异大
2. 目标区域纹理丰富
3. 目标区域与周围纹理特征差异大

实现思路

1. 为了兼顾多种尺寸，难以做出细粒度的纹理评估方法，采用固定几种窗口尺寸的策略折中进行评估
2. 在每种尺寸下使用 Hog特征作为基础，生成纹理强度特征和邻域相似度特征，再生成残差特征
3. 将三种特征归一化后合并得到最终适合模板匹配的得分 Map
4. 取得分最高的作为输出结果

特征生成

图像预处理

原图如果太大需要下采样到合适的尺寸

窗口尺寸

64 - 256 * 64 - 256，步长为 64，一共16种尺寸

残差特征

当前图像在 8 个方向上计算梯度，求得 8 组梯度特征，在一个小窗口内对梯度平方求和，取每个位置下最小梯度作为该位置的梯度特征，归一化后得到残差特征。

Hog 特征

45度为步进角度计算 Hog 特征，在所有尺寸下制作 Hog 特征和矩阵，作为生成纹理和邻域相似特征的基础

纹理强度特征

不希望细密的纹理称为精对位的核心纹理，因此计算纹理强度特征时，Hog 特征原始图像经过了中值滤波

在每种尺寸的 Hog 特征和矩阵中，计算每个位置所有方向的 Hog 特征平方和作为纹理强度特征

邻域相似特征

在每种尺寸的 Hog 特征和矩阵中，计算每个位置与周围（上下左右）的 Hog 不同方向特征的最小差异作为该尺寸下该位置的邻域相似特征

结果汇总

将结果做归一化，汇总得到统一的纹理强度 Map，取前几高分值作为结果返回

本质上是局部纹理特征的独特性搜索，并不是全局唯一

问题

• 最主要的问题是没有提出合理的结果评估方法
• 没有评估就无法做消融实验，因此算法设计是启发式的
• 已经尽量做了速度优化，但还是需要几秒，主要是兼容多种尺寸浪费了很多时间

Python 源码

Demo Code

import mtutils as mt
import numpy as np
import cv2
import time
class UniqueZone:
    def __init__(self, crop_step=5, angle_step=45, size_max=256, size_step=64, search_step=4, search_num=5, pattern_dis=6, blur_kernel=7, resize_target=400, max_factor=15) -> None:
        assert size_step < size_max
        self.angle_step = angle_step
        self.crop_step = crop_step
        self.size_max = size_max
        self.size_step = size_step
        self.search_step = search_step
        self.search_num = search_num
        self.pattern_dis = pattern_dis
        self.blur_kernel = blur_kernel
        self.resize_target = resize_target
        self.max_factor = max_factor
    @staticmethod
    def gradient_integral_map(img, angle_step):
        H, W = img.shape[:2]

        #  梯度积分图
        angle_feature_channel_num = 360 // angle_step
        hog_matrix = np.zeros([H+1, W+1, angle_feature_channel_num], dtype='float32')
        hog_gra_matrix = np.zeros([H, W, angle_feature_channel_num], dtype='float32')

        gx = cv2.Sobel(img, cv2.CV_32F, 1, 0, ksize=1)
        gy = cv2.Sobel(img, cv2.CV_32F, 0, 1, ksize=1)

        mag, angle = cv2.cartToPolar(gx, gy, angleInDegrees=True)
        angle_index = (angle / angle_step).round().astype('uint8')

        for channel_index in range(angle_feature_channel_num):
            temp_map = mag * (angle_index == channel_index)
            hog_gra_matrix[:,:,channel_index] = temp_map
            sum_map = cv2.integral(temp_map.astype('float64'))
            hog_matrix[:,:,channel_index] = sum_map

        return hog_matrix, hog_gra_matrix

    @staticmethod
    def make_map_for_size(hog_matrix, crop_width, crop_height):
        H, W, C = hog_matrix.shape[:3]
        if crop_width >= W or crop_height >= H:
            return None
        result_map = np.zeros([H-crop_height, W-crop_width, C])
        result_map = hog_matrix[crop_height:, crop_width:, :] - hog_matrix[crop_height:, :-crop_width] - hog_matrix[:-crop_height, crop_width:] + hog_matrix[:-crop_height, :-crop_width]
        return result_map


    def get_score(self, result_map_normed, x_index, y_index):
        H, W = result_map_normed.shape[:2]
        res_score = 999999

        def cal_score(res_data, res_score):
            score = (1 - (res_data.sum(axis=1) ** 4)).mean()
            if res_score > score:
                res_score = score
            return res_score

        # top
        top_matrix = result_map_normed[max(0, y_index - self.search_num * self.search_step):y_index:self.search_step, x_index , :]
        if min(top_matrix.shape[:2]) > 2:
            res_data = top_matrix * result_map_normed[y_index, x_index, :]
            res_score = cal_score(res_data, res_score)

        # bottom
        bottom_matrix = result_map_normed[y_index:min(H, y_index + self.search_num * self.search_step):self.search_step, x_index , :]
        if min(bottom_matrix.shape[:2]) > 2:
            res_data = bottom_matrix * result_map_normed[y_index, x_index, :]
            res_score = cal_score(res_data, res_score)

        # left
        left_matrix = result_map_normed[y_index, max(0, x_index - self.search_num * self.search_step):x_index:self.search_step, :]
        if min(left_matrix.shape[:2]) > 2:
            res_data = left_matrix * result_map_normed[y_index, x_index, :]
            res_score = cal_score(res_data, res_score)

        # right
        right_matrix = result_map_normed[y_index, x_index:min(W, x_index + self.search_num * self.search_step):self.search_step, :]
        if min(right_matrix.shape[:2]) > 2:
            res_data = right_matrix * result_map_normed[y_index, x_index , :]
            res_score = cal_score(res_data, res_score)

        return res_score

    def edge_map(self, result_map):
        H, W = result_map.shape[:2]
        result_map_normed = result_map / (((result_map ** 2).sum(axis=2)) ** 0.5)[:,:,None]

        # start = time.time()
        win_size = self.search_num * self.search_step
        pad_map = mt.image_border_move(result_map_normed, win_size, value=0.4)

        result_map = np.zeros([H, W, 8], "float32")

        temp_map = np.zeros([H, W, self.search_num-1], dtype='float32')

        # left
        for index in range(1, self.search_num):
            temp_map[:,:,index-1] = np.clip((result_map_normed * pad_map[win_size:H+win_size, win_size - index*self.search_step:win_size+W - index*self.search_step]).sum(axis=2), None, 1)
        result_map[:,:,0] = (1 - temp_map ** 4).mean(axis=2)
        # left top
        for index in range(1, self.search_num):
            temp_map[:,:,index-1] = np.clip((result_map_normed * pad_map[win_size - index*self.search_step:H+win_size - index*self.search_step, win_size - index*self.search_step:win_size+W - index*self.search_step]).sum(axis=2), None, 1)
        result_map[:,:,1] = (1 - temp_map ** 4).mean(axis=2)
        # top
        for index in range(1, self.search_num):
            temp_map[:,:,index-1] = np.clip((result_map_normed * pad_map[win_size - index*self.search_step:H+win_size - index*self.search_step, win_size:win_size+W]).sum(axis=2), None, 1)
        result_map[:,:,2] = (1 - temp_map ** 4).mean(axis=2)
        # right top
        for index in range(1, self.search_num):
            temp_map[:,:,index-1] = np.clip((result_map_normed * pad_map[win_size - index*self.search_step:H+win_size - index*self.search_step, win_size + index*self.search_step:win_size+W + index*self.search_step]).sum(axis=2), None, 1)
        result_map[:,:,3] = (1 - temp_map ** 4).mean(axis=2)
        # right
        for index in range(1, self.search_num):
            temp_map[:,:,index-1] = np.clip((result_map_normed * pad_map[win_size:H+win_size, win_size + index*self.search_step:win_size+W + index*self.search_step]).sum(axis=2), None, 1)
        result_map[:,:,4] = (1 - temp_map ** 4).mean(axis=2)
        # right bottom
        for index in range(1, self.search_num):
            temp_map[:,:,index-1] = np.clip((result_map_normed * pad_map[win_size + index*self.search_step:H+win_size + index*self.search_step, win_size + index*self.search_step:win_size+W + index*self.search_step]).sum(axis=2), None, 1)
        result_map[:,:,5] = (1 - temp_map ** 4).mean(axis=2)
        # bottom
        for index in range(1, self.search_num):
            temp_map[:,:,index-1] = np.clip((result_map_normed * pad_map[win_size + index*self.search_step:H+win_size + index*self.search_step, win_size:win_size+W]).sum(axis=2), None, 1)
        result_map[:,:,6] = (1 - temp_map ** 4).mean(axis=2)
        # left bottom
        for index in range(1, self.search_num):
            temp_map[:,:,index-1] = np.clip((result_map_normed * pad_map[win_size + index*self.search_step:H+win_size + index*self.search_step, win_size - index*self.search_step:win_size+W - index*self.search_step]).sum(axis=2), None, 1)
        result_map[:,:,7] = (1 - temp_map ** 4).mean(axis=2)
        
        score_map = result_map.min(axis=2)
        # print(f"new time {time.time() - start}")
        # start = time.time()
        # score_map2 = np.zeros([H, W], "float32")

        # x_index = 0
        # while x_index < W:
        #     y_index = 0
        #     while y_index < H:
        #         score_map2[y_index, x_index] = self.get_score(result_map_normed, x_index, y_index)
        #         y_index += 1
        #     x_index += 1
        # print(f"old time {time.time() - start}")
        return score_map


    def edge_map2(self, result_map):
        H, W = result_map.shape[:2]

        # start = time.time()
        win_size = self.search_num * self.search_step
        pad_map = mt.image_border_move(result_map, win_size, value=0)

        temp_map = np.zeros([H, W, (self.search_num-1) * 4, result_map.shape[-1]], dtype='float32')

        # left
        g_index = 0
        for index in range(1, self.search_num):
            temp_map[:,:,index-1 + g_index * (self.search_num-1)] = result_map - pad_map[win_size:H+win_size, win_size - index*self.search_step:win_size+W - index*self.search_step]

        # # left top
        # g_index = 1
        # for index in range(1, self.search_num):
        #     temp_map[:,:,index-1 + g_index * (self.search_num-1)] = result_map - pad_map[win_size - index*self.search_step:H+win_size - index*self.search_step, win_size - index*self.search_step:win_size+W - index*self.search_step]

        # top
        g_index = 1
        for index in range(1, self.search_num):
            temp_map[:,:,index-1 + g_index * (self.search_num-1)] = result_map - pad_map[win_size - index*self.search_step:H+win_size - index*self.search_step, win_size:win_size+W]

        # # right top
        # g_index = 3
        # for index in range(1, self.search_num):
        #     temp_map[:,:,index-1 + g_index * (self.search_num-1)] = result_map - pad_map[win_size - index*self.search_step:H+win_size - index*self.search_step, win_size + index*self.search_step:win_size+W + index*self.search_step]

        # right
        g_index = 2
        for index in range(1, self.search_num):
            temp_map[:,:,index-1 + g_index * (self.search_num-1)] = result_map - pad_map[win_size:H+win_size, win_size + index*self.search_step:win_size+W + index*self.search_step]

        # # right bottom
        # g_index = 5
        # for index in range(1, self.search_num):
        #     temp_map[:,:,index-1 + g_index * (self.search_num-1)] = result_map - pad_map[win_size + index*self.search_step:H+win_size + index*self.search_step, win_size + index*self.search_step:win_size+W + index*self.search_step]

        # bottom
        g_index = 3
        for index in range(1, self.search_num):
            temp_map[:,:,index-1 + g_index * (self.search_num-1)] = result_map - pad_map[win_size + index*self.search_step:H+win_size + index*self.search_step, win_size:win_size+W]

        # # left bottom
        # g_index = 7
        # for index in range(1, self.search_num):
        #     temp_map[:,:,index-1 + g_index * (self.search_num-1)] = result_map - pad_map[win_size + index*self.search_step:H+win_size + index*self.search_step, win_size - index*self.search_step:win_size+W - index*self.search_step]

        return temp_map


    def make_edge_map_dict(self, size_map_dict, cur_size_step, cur_size_max):
        edge_map_dict = dict()
        def cal_score(matrix):
            return (matrix ** 2).sum(axis=3).min(axis=2) ** 0.5

        for size_x in size_map_dict.keys():
            if size_x not in edge_map_dict:
                edge_map_dict[size_x] = dict()
            for size_y in size_map_dict[size_x].keys():
                if size_y == cur_size_step:
                    if size_x == cur_size_step:
                        edge_map_dict[size_x][size_y] = self.edge_map2(size_map_dict[size_x][size_y])
                    else:
                        edge_map_dict[size_x][size_y] = edge_map_dict[size_x - cur_size_step][size_y][:, :-cur_size_step, ...] + edge_map_dict[cur_size_step][cur_size_step][:, cur_size_step * ((size_x - cur_size_step) // cur_size_step):, ...]
                        if size_x - cur_size_step > cur_size_step:
                            edge_map_dict[size_x - cur_size_step][size_y] = cal_score(edge_map_dict[size_x - cur_size_step][size_y])
                else:
                    edge_map_dict[size_x][size_y] = edge_map_dict[size_x][size_y - cur_size_step][:-cur_size_step, :, ...] + edge_map_dict[size_x][cur_size_step][cur_size_step * ((size_y - cur_size_step) // cur_size_step):, :, ...]
                    if size_y - cur_size_step > cur_size_step:
                        edge_map_dict[size_x][size_y - cur_size_step] = cal_score(edge_map_dict[size_x][size_y - cur_size_step])
                    if size_y + cur_size_step > cur_size_max:
                        edge_map_dict[size_x][size_y] = cal_score(edge_map_dict[size_x][size_y])

        edge_map_dict[cur_size_step][cur_size_step] = (edge_map_dict[cur_size_step][cur_size_step] ** 2).sum(axis=3).min(axis=2) ** 0.5
        edge_map_dict[size_x][cur_size_step] = cal_score(edge_map_dict[size_x][cur_size_step])

        for size_x in edge_map_dict.keys():
            for size_y in edge_map_dict[size_x].keys():
                edge_map_dict[size_x][size_y] /= (size_x ** 2 + size_y ** 2) ** 0.5

        return edge_map_dict

    @staticmethod
    def diff_map(img, cur_size_step, cur_pattern_dis):
        H, W = img.shape[:2]
        win_size = max(4, cur_size_step // cur_pattern_dis)
        img_resized = mt.image_resize(img, factor=1 / cur_pattern_dis)
        rs_H, rs_W = img_resized.shape[:2]
        pad_img = mt.image_border_move(img_resized, 1).astype('float32')
        feature_map = np.zeros([rs_H, rs_W, 8], 'float32')

        # left
        feature_map[:, :, 0] = pad_img[1:-1, 1:-1] - pad_img[1:-1, 0:-2]
        # left_top
        feature_map[:, :, 1] = pad_img[1:-1, 1:-1] - pad_img[0:-2, 0:-2]
        # top
        feature_map[:, :, 2] = pad_img[1:-1, 1:-1] - pad_img[0:-2, 1:-1]
        # top_right
        feature_map[:, :, 3] = pad_img[1:-1, 1:-1] - pad_img[0:-2, 2:]
        # right
        feature_map[:, :, 4] = pad_img[1:-1, 1:-1] - pad_img[1:-1, 2:]
        # right_bottom
        feature_map[:, :, 5] = pad_img[1:-1, 1:-1] - pad_img[2:, 2:]
        # bottom
        feature_map[:, :, 6] = pad_img[1:-1, 1:-1] - pad_img[2:, 1:-1]
        # left_bottom
        feature_map[:, :, 7] = pad_img[1:-1, 1:-1] - pad_img[2:, 0:-2]

        _, sum_feature = cv2.integral2(feature_map)

        result_map = sum_feature[win_size:, win_size:, :] - sum_feature[win_size:, :-win_size] - sum_feature[:-win_size, win_size:] + sum_feature[:-win_size, :-win_size]
        result_map = result_map ** 0.5
        result_map = result_map.min(axis=2)
        # result_map = (mt.min_max_normalize(result_map + 1e-6) * 255).astype('uint8')
        result_map = mt.image_border_move(result_map, 0, 0, win_size-1, win_size-1)
        result_map = mt.image_resize(result_map, [W, H], uint8=False)
        result_map = mt.min_max_normalize(result_map)
        return result_map

    def result_show(self, color_img, final_map_img, index_map, cur_size_step):
        H, W = color_img.shape[:2]
        map_img = final_map_img.copy()
        score = map_img.max()
        margin = 20
        gap = cur_size_step // 2
        box_num = 0
        while score > 70 and box_num < 3:
            _, value, _, (x, y) = cv2.minMaxLoc(map_img)
            map_index = index_map[y, x]
            size_w = cur_size_step * (map_index // (self.size_max // cur_size_step) + 1)
            size_h = cur_size_step * (map_index % (self.size_max // cur_size_step) + 1)
            bbox = [max(0, x - margin), max(0, y - margin), min(W- 3, x + size_w + margin), min(H-3, y + size_h + margin)]
            color_img = mt.draw_boxes(color_img, [bbox], [0, min(255, value ** 1.13), 0], thickness=2)
            gap_x = max(min((bbox[2] - bbox[0]) * 2, W // 5), W // 10)
            gap_y = max(min((bbox[3] - bbox[1]) * 2, H // 5), H // 10)
            map_img[max(0, y - gap_y):min(H, y+gap_y), max(0, x - gap_x):min(W, x+gap_x)] = 0
            score = map_img.max()
            box_num += 1
        return color_img

    def mask_size_map_dict(self, hog_matrix, cur_size_step, cur_size_max):
        # step_map = self.make_map_for_size(hog_matrix, cur_size_step, cur_size_step)

        # make empty dict
        size_map_dict = dict()
        for size1 in range(cur_size_step, cur_size_max + 1, cur_size_step):
            size_map_dict[size1] = dict()
            for size2 in range(cur_size_step, cur_size_max + 1, cur_size_step):
                size_map_dict[size1][size2] = None

        for size_x in range(cur_size_step, cur_size_max + 1, cur_size_step):
            for size_y in range(cur_size_step, cur_size_max + 1, cur_size_step):
                if size_y == cur_size_step:
                    if size_x == cur_size_step:
                        size_map_dict[size_x][size_y] = self.make_map_for_size(hog_matrix, size_x, size_y)
                    else:
                        size_map_dict[size_x][size_y] = size_map_dict[size_x - cur_size_step][size_y][:, :-cur_size_step, ...] + size_map_dict[cur_size_step][cur_size_step][:, cur_size_step * ((size_x - cur_size_step) // cur_size_step):, ...]
                else:
                    size_map_dict[size_x][size_y] = size_map_dict[size_x][size_y - cur_size_step][:-cur_size_step, :, ...] + size_map_dict[size_x][cur_size_step][cur_size_step * ((size_y - cur_size_step) // cur_size_step):, :, ...]

        return size_map_dict

    @staticmethod
    def make_pattern_map_dict(size_map_dict):
        pattern_map_dict = dict()
        for size_x in size_map_dict.keys():
            pattern_map_dict[size_x] = dict()
            for size_y in size_map_dict[size_x].keys():
                pattern_map_dict[size_x][size_y] = mt.min_max_normalize(size_map_dict[size_x][size_y].sum(axis=2))
        return pattern_map_dict

    def make_final_map(self, ori_img, resized_img, pattern_map_dict, size_map_dict, edge_map_dict, diff_map):
        ori_H, ori_W = ori_img.shape[:2]
        H, W = resized_img.shape[:2]
        
        size_final_map_temp = np.zeros([H, W, (len(pattern_map_dict.keys())) ** 2])
        index = 0
        for key1 in size_map_dict.keys():
            for key2 in size_map_dict[key1].keys():
                edge_map = edge_map_dict[key1][key2]
                pattern_map = pattern_map_dict[key1][key2]
                cur_h, cur_w = pattern_map.shape[:2]
                size_final_map_temp[:cur_h, :cur_w, index] = (edge_map ** 2) * (pattern_map ** 0.5)
                index += 1
        ori_size_final_map_temp = mt.image_resize(size_final_map_temp, shape=[ori_W, ori_H], uint8=False)
        index_map = np.argmax(ori_size_final_map_temp, axis=2)
        final_map = ori_size_final_map_temp.max(axis=2)

        zero_edge = 20
        edge_mask = np.ones([ori_H - zero_edge * 2 - self.size_step, ori_W - zero_edge * 2 - self.size_step])
        edge_mask = mt.image_border_move(edge_mask, zero_edge, zero_edge, self.size_step + zero_edge, self.size_step + zero_edge)

        final_map *= diff_map ** 0.5
        final_map *= edge_mask

        # size_final_map_map = (pattern_map ** 0.5) * (edge_map ** 2.5) * (diff_map[:-height+1, :-width+1] ** 1)
        final_map_img = (mt.min_max_normalize(final_map) * 255).astype('uint8')
        return final_map_img, index_map


    def infer(self, img):
        begin = time.time()
        color_img = mt.to_colorful_image(img)
        img = mt.to_gray_image(img)
        ori_H, ori_W = img.shape[:2]
        print(f"original img size H {ori_H} W {ori_W}")

        cur_factor = max(1, min(self.max_factor, max([ori_H, ori_W]) / max(1, self.resize_target)))

        resized_img = cv2.resize(img, [mt.round(ori_W/cur_factor), mt.round(ori_H/cur_factor)], interpolation=cv2.INTER_LINEAR)
        H, W = resized_img.shape[:2]

        cur_size_max = int(min(max(7, self.size_max // cur_factor), W * 0.8))
        cur_size_step = int(min(max(7, self.size_step // cur_factor), W * 0.4))
        
        cur_pattern_dis = int(max(2, self.pattern_dis // cur_factor))

        cur_blur_kernel = max(mt.round(self.blur_kernel / (cur_factor ** 0.5)) // 2 * 2 + 1, 1)
        mb_img = cv2.medianBlur(resized_img, cur_blur_kernel)

        start = time.time()
        diff_map = self.diff_map(resized_img, cur_size_step, cur_pattern_dis)
        diff_map = mt.image_resize(diff_map, shape=[ori_W, ori_H], uint8=False)
        print(f"diff_map time {time.time() - start}s.")
        
        ## 积分图
        start = time.time()
        hog_matrix, _ = self.gradient_integral_map(resized_img, self.angle_step)
        hog_mb_matrix, _ = self.gradient_integral_map(mb_img, self.angle_step)
        print(f"hog fea extract time {time.time() - start}s.")

        start = time.time()
        size_map_dict = self.mask_size_map_dict(hog_matrix, cur_size_step, cur_size_max)
        print(f"size {W} {H} size_map * 16 time {time.time() - start}s.")

        start = time.time()
        size_mb_map_dict = self.mask_size_map_dict(hog_mb_matrix, cur_size_step, cur_size_max)
        pattern_map_dict = self.make_pattern_map_dict(size_mb_map_dict)
        print(f"size {W} {H} pattern_map time {time.time() - start}s.")

        start = time.time()
        edge_map_dict = self.make_edge_map_dict(size_map_dict, cur_size_step, cur_size_max)
        print(f"size {W} {H} edge_map time {time.time() - start}s.")

        start = time.time()
        final_map_img, index_map = self.make_final_map(color_img, resized_img, pattern_map_dict, size_map_dict, edge_map_dict, diff_map)
        print(f"size {W} {H} final_map time {time.time() - start}s.")
        print(f"Shape {W} {H} total time {time.time() - begin}s.")

        color_img = self.result_show(color_img, final_map_img, index_map, self.size_step)

        return final_map_img, color_img


if __name__ == '__main__':

    img = mt.cv_rgb_imread('test.png')
    obj = UniqueZone()
    final_map_img, color_img = obj.infer(img)
    mt.cv_rgb_imwrite(final_map_img, 'result.png')

文章链接：
https://www.zywvvd.com/notes/study/image-processing/auto-pattern/auto-pattern/