人群密度圖生成-Python實現

07-18

人群密度圖生成-Python實現

來自專欄深度學習填坑

最近在做人群密度估計，是一種以密度圖的形式表示人群密集程度的方法，大概如下圖

左圖是某公共場所的人群圖，右圖是人群密度圖，也就是用來訓練生成人群密度圖網路的標籤。首先要將原圖中左右人頭的坐標標定出來，然後生成密度圖像。

原理

如果一個標註點的位置為 $x_i$ ，我們可以將它表示為 $delta(x-x_i)$ ，因此有 $N$ 個人頭的標籤可以表示為：

$H(x)=sum_{i=1}^{N}{delta(x-x_i)}$

我們可以使用高斯核 $G_{sigma}$ 對這個函數卷積得到密度函數 $F(x)=H(x)*G_{sigma}(x)$ 。然而這是假設每個 $x_i$ 在樣本空間中是獨立的。實際上，每個 $x_i$ 是3D場景中的人群密度樣本，由於透視失真，像素與周邊樣本在不同場景區域尺度不一致。

因此為了精確估計群體密度 $F$ 需要考慮透視變換，我們需要考慮真實和圖像的失真。通常我們不能確定場景的幾乎形狀。但是，如果我們每個頭部周圍的人群是均勻分布的，那麼頭部與最近k個鄰居之間的圖像中的平均距離可以給出幾何失真的合理估計。

因此，我們應該根據圖像中每個人頭部大小來確定參數 $sigma$ ，但是在實際情況下，由於遮擋，不可能準確的獲得頭部的尺寸，而且很難找到頭部尺寸和密度圖之間的關係。我們發現頭部大小通常與擁擠場景中兩個頭部中心距離有關，所以我們使用相鄰頭部的平均距離作為參數。

對於每個人頭 $x_i$ ，給出 $k$ 個近鄰頭部距離 ${d^i_1,d^i_2,...,d^i_m}$ 的平均值 $ar{d_i}=frac{1}{m}sum_{j=1}^{m}{d_j^i}$ , $x_i$ 相關的像素對應場景中地面上的一個區域，這個區域的半徑與 $ar{d_i}$ 成正比。因此，為了估計像素 $x_i$ 周圍的人群密度，我們需要用一個自適應的高斯核與 $H(x)$ 卷積，這個高斯核的方差 $sigma_i$ 可變且與 $ar{d_i}$ 成比例。

$F(x)=sum_{i=1}^{N}{delta(x-x_i)}*G_{sigma_i}(x), quad withquad sigma_i=etaar{d_i}$

這裡 $eta=0.3$ 。

代碼

import numpy as npimport osimport matplotlib.image as mpimgimport scipy.io as siofrom scipy.ndimage import filtersimport matplotlib.pyplot as pltfrom sklearn.neighbors import NearestNeighborsfrom PIL import Imageimport mathdef gaussian_filter_density(gt): pts = np.array(list(zip(np.nonzero(gt)[1], np.nonzero(gt)[0]))) neighbors = NearestNeighbors(n_neighbors=4, algorithm=kd_tree, leaf_size=1200) neighbors.fit(pts.copy()) distances, _ = neighbors.kneighbors() density = np.zeros(gt.shape, dtype=np.float32) type(distances) sigmas = distances.sum(axis=1) * 0.075 for i in range(len(pts)): pt = pts[i] pt2d = np.zeros(shape=gt.shape, dtype=np.float32) pt2d[pt[1]][pt[0]] = 1 # starttime = datetime.datetime.now() density += filters.gaussian_filter(pt2d, sigmas[i], mode=constant) # endtime = datetime.datetime.now() # # interval = (endtime - starttime) # print(interval) return densitydef create_density(gts, d_map_h, d_map_w): res = np.zeros(shape=[d_map_h, d_map_w]) bool_res = (gts[:, 0] < d_map_w) & (gts[:, 1] < d_map_h) for k in range(len(gts)): gt = gts[k] if (bool_res[k] == True): res[int(gt[1])][int(gt[0])] = 1 pts = np.array(list(zip(np.nonzero(res)[1], np.nonzero(res)[0]))) neighbors = NearestNeighbors(n_neighbors=4, algorithm=kd_tree, leaf_size=1200) neighbors.fit(pts.copy()) distances, _ = neighbors.kneighbors() map_shape = [d_map_h, d_map_w] density = np.zeros(shape=map_shape, dtype=np.float32) sigmas = distances.sum(axis=1) * 0.075 for i in range(len(pts)): pt = pts[i] pt2d = np.zeros(shape=map_shape, dtype=np.float32) pt2d[pt[1]][pt[0]] = 1 # starttime = datetime.datetime.now() density += filters.gaussian_filter(pt2d, sigmas[i], mode=constant) return densityif __name__ == __main__: train_img = ./ShanghaiTech_Crowd_Counting_Dataset/part_A_final/train_data/images train_gt = ./ShanghaiTech_Crowd_Counting_Dataset/part_A_final/train_data/ground_truth out_path = ./ShanghaiTech_Crowd_Counting_Dataset/true_crowd_counting/ validation_num = 15 img_names = os.listdir(train_img) num = len(img_names) num_list = np.arange(1, num + 1) # random.shuffle(num_list) global_step = 1 for i in num_list: full_img = train_img + /IMG_ + str(i) + .jpg full_gt = train_gt + /GT_IMG_ + str(i) + .mat img = mpimg.imread(full_img) data = sio.loadmat(full_gt) gts = data[image_info][0][0][0][0][0] # shape like (num_count, 2) count = 1 # fig1 = plt.figure(fig1) # plt.imshow(img) shape = img.shape if (len(shape) < 3): img = img.reshape([shape[0], shape[1], 1]) d_map_h = math.floor(math.floor(float(img.shape[0]) / 2.0) / 2.0) d_map_w = math.floor(math.floor(float(img.shape[1]) / 2.0) / 2.0) # starttime = datetime.datetime.now() # den_map = gaussian_filter_density(res) if (global_step == 4): print(1) den_map = create_density(gts / 4, d_map_h, d_map_w) # endtime = datetime.datetime.now() # interval = (endtime - starttime).seconds # print(interval) p_h = math.floor(float(img.shape[0]) / 3.0) p_w = math.floor(float(img.shape[1]) / 3.0) d_map_ph = math.floor(math.floor(p_h / 2.0) / 2.0) d_map_pw = math.floor(math.floor(p_w / 2.0) / 2.0) if (global_step < validation_num): mode = val else: mode = train py = 1 py2 = 1 for j in range(1, 4): px = 1 px2 = 1 for k in range(1, 4): print(global + str(global_step)) # print(j + str(j)) # print(k +str(k)) print(----------) if (global_step == 4 & j == 3 & k == 4): print(global + str(global_step)) final_image = img[py - 1: py + p_h - 1, px - 1: px + p_w - 1, :] final_gt = den_map[py2 - 1: py2 + d_map_ph - 1, px2 - 1: px2 + d_map_pw - 1] px = px + p_w px2 = px2 + d_map_pw if final_image.shape[2] < 3: final_image = np.tile(final_image, [1, 1, 3]) image_final_name = out_path + mode + _img/ IMG_ + str(i) + _ + str(count) + .jpg gt_final_name = out_path + mode + _gt/ + GT_IMG_ + str(i) + _ + str(count) Image.fromarray(final_image).convert(RGB).save(image_final_name) np.save(gt_final_name, final_gt) count = count + 1 py = py + p_h py2 = py2 + d_map_ph global_step = global_step + 1 # fig2 = plt.figure(fig2) # plt.imshow(den_map) # plt.show()

這裡是以ShanghaiTech_Crowd_Counting_Dataset為例。