Density Map

Posted on

This is some implementations of create ground truth density maps for crowd counting.

This one is from Single Image Crowd Counting via MCNN (Unofficial Implementation). The author uses fixed gaussian kernel size. 

1


2


3


4


5


6


7


8


9


10


11


12


13


14


15


16


17


18


19


20


21


22


23


24


25


26


27


28


29


30


31


32


33


34


clc; 


dataset = 'A';


dataset_name = ['shanghaitech_part_' dataset ];


path = ['../data/original/shanghaitech/part_' dataset '_final/test_data/images/'];


gt_path = ['../data/original/shanghaitech/part_' dataset '_final/test_data/ground-truth/'];


gt_path_csv = ['../data/original/shanghaitech/part_' dataset '_final/test_data/ground_truth_csv/'];




mkdir(gt_path_csv);


if (dataset == 'A')


    num_images = 182;


else


    num_images = 316;


end




gt_people_count = [];




for i = 1:num_images  


    if (mod(i,10)==0)


        fprintf(1,'Processing %3d/%d files\n', i, num_images);


    end


    gt_file_name = [gt_path,'GT_IMG_',num2str(i),'.mat'];


    load(gt_file_name);


    input_img_name = strcat(path,'IMG_',num2str(i),'.jpg');


    im = imread(input_img_name);


    [h, w, c] = size(im);


    if (c == 3)


        im = rgb2gray(im);


    end     


    annPoints =  image_info{1}.location;


    gt_people_count = [ gt_people_count;length(annPoints)];


    [h, w, c] = size(im);


    im_density = get_density_map_gaussian(im,annPoints);    


    csvwrite([gt_path_csv ,'IMG_',num2str(i) '.csv'], im_density);       


end


1


2


3


4


5


6


7


8


9


10


11


12


13


14


15


16


17


18


19


20


21


22


23


24


25


26


27


28


29


30


31


32


33


34


35


36


37


38


39


40


41


42


43


44


45


46


47


48


49


50


51


52


53


54


55


56


57


58


59


60


61


function im_density = get_density_map_gaussian(im,points)






im_density = zeros(size(im)); 


[h,w] = size(im_density);




if( isempty(points))


    return;


end


disp('The total points are');


disp(length(points));


if(length(points(:,1))==1)


    x1 = max(1,min(w,round(points(1,1))));


    y1 = max(1,min(h,round(points(1,2))));


    im_density(y1,x1) = 255;


    return;


end




for j = 1:length(points) 	


    f_sz = 15;


    sigma = 4.0;


    H = fspecial('Gaussian',[f_sz, f_sz],sigma);


    x = min(w,max(1,abs(int32(floor(points(j,1))))));


    y = min(h,max(1,abs(int32(floor(points(j,2))))));


    if(x > w || y > h)


        continue;


    end


    


    x1 = x - int32(floor(f_sz/2)); y1 = y - int32(floor(f_sz/2));


    x2 = x + int32(floor(f_sz/2)); y2 = y + int32(floor(f_sz/2));


    dfx1 = 0; dfy1 = 0; dfx2 = 0; dfy2 = 0;


    change_H = false;


    


    if(x1 < 1)


        dfx1 = abs(x1)+1;


        x1 = 1;


        change_H = true;


    end


    if(y1 < 1)


        dfy1 = abs(y1)+1;


        y1 = 1;


        change_H = true;


    end


    if(x2 > w)


        dfx2 = x2 - w;


        x2 = w;


        change_H = true;


    end


    if(y2 > h)


        dfy2 = y2 - h;


        y2 = h;


        change_H = true;


    end


    x1h = 1+dfx1; y1h = 1+dfy1; x2h = f_sz - dfx2; y2h = f_sz - dfy2;


    if (change_H == true)


        H =  fspecial('Gaussian',[double(y2h-y1h+1), double(x2h-x1h+1)],sigma);


    end


    im_density(y1:y2,x1:x2) = im_density(y1:y2,x1:x2) +  H;


     


end


end

The second one is modified by me based on the code above. In this version, it is created for UCSD dataset where the intensities outside ROI is set to 0.

1


2


3


4


5


6


7


8


9


10


11


12


13


14


15


16


17


18


19


20


21


22


23


24


25


26


27


28


29


30


31


32


33


34


35


function im_density = ucsd_density_map(im,points,xi,yi)




im_density = zeros(size(im)); 


[h,w] = size(im_density); %h=160,w=240


if( isempty(points))


    return;


end


if(length(points(:,1))==1)


    x1 = max(1,min(w,round(points(1,1))));


    y1 = max(1,min(h,round(points(1,2))));


    im_density(y1,x1) = 255;


    return;


end




for j = 1:length(points) 


    sigma = 4.0;


    H = fspecial('Gaussian',sigma);


    x = min(w,max(1,abs(int32(floor(points(j,1))))));


    y = min(h,max(1,abs(int32(floor(points(j,2))))));


    if(x > w || y > h)


        continue;


    end 


    


    mask_area = roipoly(im,floor(xi),floor(yi));


    


    imm = zeros(h,w);


    imm(y,x) = 1; % attention: not imm(x,y) = 1 


    H_ROI = roifilt2(H,imm,mask_area);


    im_density = im_density +  H_ROI;


         


end


fprintf("the sum of sum is %f",sum(sum(im_density)));


fprintf("  the length is %f\n",length(points));




end

The third one is from 人群密度图生成-Python实现. In this one, the author uses geometry-adaptive Gaussian kernel for density map creation.

1


2


3


4


5


6


7


8


9


10


11


12


13


14


15


16


17


18


19


20


21


22


23


24


25


26


27


28


29


30


31


32


33


34


35


36


37


38


39


40


41


42


43


44


45


46


47


48


49


50


51


52


53


54


55


56


57


58


59


60


61


62


63


64


65


66


67


68


69


70


71


72


73


74


75


76


77


78


79


80


81


82


83


84


85


86


87


88


89


90


91


92


93


94


95


96


97


98


99


100


101


102


103


104


105


import numpy as np


import os


import matplotlib.image as mpimg


import scipy.io as sio


from scipy.ndimage import filters


from sklearn.neighbors import NearestNeighbors


from PIL import Image


import math




def gaussian_filter_density(gt):


    pts = np.array(list(zip(np.nonzero(gt)[1], np.nonzero(gt)[0]))) #np.nonzero return two arrays, one is x_index of nonzero value and the other is y_index of nonzero value


    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 # 0.075 = 0.3/4


    for i in range(len(pts)):


        pt = pts[i]


        pt2d = np.zeros(shape=gt.shape, dtype=np.float32)


        pt2d[pt[1]][pt[0]] = 1


        density += filters.gaussian_filter(pt2d, sigmas[i], mode='constant')


    return density






def create_density(gts, d_map_h, d_map_w):


    res = np.zeros(shape=[d_map_h, d_map_w]) # res store head positions


    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


        t1 = filters.gaussian_filter(pt2d, sigmas[i])


        density += t1 #bigger the sigma, blurrer the picture


    return density






if __name__ == '__main__':


    train_img = '../data/original/shanghaitech/part_A_final/train_data/images'


    train_gt = '../data/original/shanghaitech/part_A_final/train_data/ground-truth'


    out_path = '../data/formatted_trainval/shanghaitech_adaptive_gaussian/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 is (num_count, 2)


        count = 1




        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)


        den_map = create_density(gts / 4, d_map_h, d_map_w)






        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):


                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]) #tile : repeat "final_image" once in axis=0 and axis = 1 while repeat three times in axis=2


                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

Some other implementation 密度图生成1