|
| 1 | +import numpy as np |
| 2 | +from time import time |
| 3 | +from numba import njit |
| 4 | + |
| 5 | +def neighbors(shape, core): |
| 6 | + shp = [slice(0,i) for i in core] |
| 7 | + idx = np.mgrid[tuple(shp)] |
| 8 | + idx = idx.reshape((len(core),-1)) |
| 9 | + offset = np.array(core)//2 |
| 10 | + offset[0] = 0 |
| 11 | + idx -= offset.reshape((-1,1)) |
| 12 | + acc = np.cumprod((1,)+shape[::-1][:-1]) |
| 13 | + return np.dot(idx.T, acc[::-1]) |
| 14 | + |
| 15 | +@njit |
| 16 | +def fill_col(pdimg, idx, colimg): |
| 17 | + s = 0 |
| 18 | + for i in range(len(pdimg)): |
| 19 | + if pdimg[i]&1==0: continue |
| 20 | + for j in idx: |
| 21 | + colimg[s] = pdimg[i+j] |
| 22 | + s += 1 |
| 23 | + return colimg |
| 24 | + |
| 25 | +def img2col(img, core, stride=(1,1), buf=[np.zeros(1, dtype=np.int32)]): |
| 26 | + # new the col_img, if needed |
| 27 | + strh, strw = stride |
| 28 | + cimg_w = np.cumprod(core.shape[1:])[-1] |
| 29 | + print(img.shape) |
| 30 | + n,c,h,w = img.shape |
| 31 | + cimg_h = n*(h//strh)*(w//strw) |
| 32 | + |
| 33 | + if len(buf[0])<cimg_h*cimg_w: |
| 34 | + buf[0] = col_img = np.zeros(cimg_h*cimg_w, dtype=np.int32) |
| 35 | + else: |
| 36 | + col_img = buf[0][:cimg_h*cimg_w] |
| 37 | + col_img[:] = 0 |
| 38 | + |
| 39 | + # mark where need |
| 40 | + iimg = img.view(dtype=np.int32) |
| 41 | + #iimg.view(dtype=np.uint8).ravel()[::4]&=0xfe |
| 42 | + iimg &= 0xfffffffe |
| 43 | + iimg[:,0,::strh,::strw] |= 1 |
| 44 | + |
| 45 | + # ravel the image |
| 46 | + n,c,h,w = np.array(core.shape) |
| 47 | + shp = ((0,0),(0,0),(h//2,h//2),(w//2,w//2)) |
| 48 | + pdimg = np.pad(iimg, shp, 'constant', constant_values=0) |
| 49 | + nbs = neighbors(pdimg.shape[1:], core.shape[1:]) |
| 50 | + fill_col(pdimg.ravel(), nbs, col_img) |
| 51 | + col_img = col_img.view(np.float32) |
| 52 | + col_img = col_img.reshape((cimg_h, cimg_w)) |
| 53 | + |
| 54 | + # dot |
| 55 | + col_core = core.reshape((core.shape[0],-1)) |
| 56 | + rst = col_core.dot(col_img.T) |
| 57 | + ni, ci, hi, wi = img.shape |
| 58 | + return rst.reshape((ni, n, hi//strh, wi//strw)) |
| 59 | + |
| 60 | +if __name__ == '__main__': |
| 61 | + from skimage.data import camera |
| 62 | + import matplotlib.pyplot as plt |
| 63 | + from scipy.ndimage import convolve |
| 64 | + img = np.zeros((1, 3, 512, 512), dtype=np.float32) |
| 65 | + img.ravel()[:] = np.arange(3*512*512) |
| 66 | + |
| 67 | + iimg = img.view(dtype=np.int32) |
| 68 | + iimg.view(dtype=np.uint8).ravel()[::4]&=0xfe |
| 69 | + iimg[:,0,:,:] |= 1 |
| 70 | + |
| 71 | + core = np.zeros((32, 3, 3, 3), dtype=np.float32) |
| 72 | + core.ravel()[:] = np.arange(3*3*3*32) |
| 73 | + |
| 74 | + start = time() |
| 75 | + rst1 = img2col(img, core, (1,1)) |
| 76 | + print(time()-start) |
| 77 | + |
| 78 | + start = time() |
| 79 | + rst2 = img2col(img, core, (2,2)) |
| 80 | + print('jit cost: x10', time()-start) |
0 commit comments