1+ import numpy as np
2+ from time import time
3+ try : from numba import njit
4+ except : njit = None
5+
6+ if njit is None : print ('install numba can double speed!' )
7+
8+ def neighbors (shape , core , offset = 0 ):
9+ shp = [slice (0 ,i ) for i in core ]
10+ idx = np .mgrid [tuple (shp )]
11+ idx = idx .reshape ((len (core ),- 1 ))
12+ offset = np .array (core )// 2 * offset
13+ idx -= offset .reshape ((- 1 ,1 ))
14+ acc = np .cumprod ((1 ,)+ shape [::- 1 ][:- 1 ])
15+ return np .dot (idx .T , acc [::- 1 ])
16+
17+ def jit_fill_col (pdimg , idx , colimg ):
18+ s = 0
19+ for i in range (len (pdimg )):
20+ if pdimg [i ]& 1 == 0 : continue
21+ for j in idx :
22+ colimg [s ] = pdimg [i + j ]
23+ s += 1
24+ return colimg
25+
26+ def fill_col (pdimg , idx , colimg ):
27+ rc = np .where (pdimg & 1 )[0 ]
28+ rc = rc .reshape ((- 1 ,1 ))+ idx
29+ colimg [:] = pdimg [rc .ravel ()]
30+ return colimg
31+
32+ if not njit is None : fill_col = njit (jit_fill_col )
33+
34+ def conv (img , core , stride = (1 ,1 ), buf = ['' ]):
35+ # new the col_img, if needed
36+ strh , strw = stride
37+ cimg_w = np .cumprod (core .shape [1 :])[- 1 ]
38+ n ,c ,h ,w = img .shape
39+ cimg_h = n * (h // strh )* (w // strw )
40+ if len (buf [0 ])< cimg_h * cimg_w :
41+ col_img = np .zeros (cimg_h * cimg_w , dtype = np .int32 )
42+ buf [0 ] = col_img
43+ else :
44+ col_img = buf [0 ][:cimg_h * cimg_w ]
45+ col_img [:] = 0
46+ # mark where need
47+ iimg = img .view (dtype = np .int32 )
48+ iimg &= 0xfffffffe
49+ iimg [:,0 ,::strh ,::strw ] |= 1
50+ # ravel the image
51+ n ,c ,h ,w = np .array (core .shape )
52+ shp = ((0 ,0 ),(0 ,0 ),(h // 2 ,h // 2 ),(w // 2 ,w // 2 ))
53+ pdimg = np .pad (iimg , shp , 'constant' , constant_values = 0 )
54+ nbs = neighbors (pdimg .shape [1 :], core .shape [1 :], (0 ,1 ,1 ))
55+ fill_col (pdimg .ravel (), nbs , col_img )
56+ col_img = col_img .view (np .float32 )
57+ col_img = col_img .reshape ((cimg_h , cimg_w ))
58+ # dot
59+ col_core = core .reshape ((core .shape [0 ],- 1 ))
60+ rst = col_core .dot (col_img .T )
61+ ni , ci , hi , wi = img .shape
62+ return rst .reshape ((ni , n , hi // strh , wi // strw ))
63+
64+
65+ def jit_fill_max (pdimg , idx , colimg ):
66+ s = 0
67+ for i in range (len (pdimg )):
68+ if pdimg [i ]& 1 == 0 : continue
69+ for j in idx :
70+ colimg [s ] = max (colimg [s ], pdimg [i + j ])
71+ s += 1
72+ return colimg
73+
74+ def fill_max (pdimg , idx , colimg ):
75+ rc = np .where (pdimg & 1 )[0 ]
76+ rc = rc .reshape ((- 1 ,1 ))+ idx
77+ vs = pdimg [rc .ravel ()].reshape ((- 1 , len (idx )))
78+ np .max (vs , axis = - 1 , out = colimg )
79+
80+ if not njit is None : fill_max = njit (jit_fill_max )
81+
82+ def maxpool (img , stride = (2 ,2 )):
83+ strh , strw = stride
84+ n ,c ,h ,w = img .shape
85+ cimg_h = n * c * (h // strh )* (w // strw )
86+
87+ iimg = img .view (dtype = np .int32 )
88+ iimg &= 0xfffffffe
89+ iimg [:,:,::strh ,::strw ] |= 1
90+
91+ nbs = neighbors (img .shape [1 :], (1 ,)+ stride )
92+ shp = (n , c , h // strh , w // strw )
93+ colimg = np .zeros (shp , dtype = np .int32 )
94+ fill_max (iimg .ravel (), nbs , colimg .ravel ())
95+ return colimg .view (np .float32 )
96+
97+ def jit_resize (img , k , ra , rb , rs , _rs , ca , cb , cs , _cs , out ):
98+ h , w = img .shape
99+ for r in range (h * k ):
100+ rar = ra [r ]
101+ rbr = rar + 1
102+ rsr = rs [r ]
103+ _rsr = _rs [r ]
104+ for c in range (w * k ):
105+ cac = ca [c ]
106+ cbc = cac + 1
107+ rra = img [rar ,cac ]* _rsr
108+ rra += img [rbr ,cac ]* rsr
109+ rrb = img [rar ,cbc ]* _rsr
110+ rrb += img [rbr ,cbc ]* rsr
111+ rcab = rra * _cs [c ] + rrb * cs [c ]
112+ out [r ,c ] = rcab
113+
114+ def resize (img , k , ra , rb , rs , _rs , ca , cb , cs , _cs , out ):
115+ out [:img .shape [0 ]] = img [:,ca ]* _cs + img [:,cb ]* cs
116+ out [:] = (out [ra ].T * _rs + out [rb ].T * rs ).T
117+
118+ if not njit is None : resize = njit (jit_resize )
119+
120+ def upsample (img , k , out = None ):
121+ nc , (h , w ) = img .shape [:- 2 ], img .shape [- 2 :]
122+ if out is None :
123+ out = np .zeros (nc + (h * k , w * k ), dtype = img .dtype )
124+ rs = np .linspace (- 0.5 + 0.5 / k , h - 0.5 - 0.5 / k , h * k , dtype = np .float32 )
125+ cs = np .linspace (- 0.5 + 0.5 / k , w - 0.5 - 0.5 / k , w * k , dtype = np .float32 )
126+ np .clip (rs , 0 , h - 1 , out = rs )
127+ np .clip (cs , 0 , w - 1 , out = cs )
128+ ra = np .floor (rs ).astype (np .uint32 )
129+ ca = np .floor (cs ).astype (np .uint32 )
130+ np .clip (ra , 0 , h - 1.5 , out = ra )
131+ np .clip (ca , 0 , w - 1.5 , out = ca )
132+ rs -= ra ; cs -= ca ;
133+ outcol = out .reshape ((- 1 , h * k , w * k ))
134+ imgcol = img .reshape ((- 1 , h , w ))
135+ for i , o in zip (imgcol , outcol ):
136+ resize (i , k , ra , ra + 1 , rs , 1 - rs , ca , ca + 1 , cs , 1 - cs , o )
137+ return out
138+
139+ if __name__ == '__main__' :
140+ from skimage .data import camera
141+ import matplotlib .pyplot as plt
142+ from scipy .ndimage import convolve
143+ img = np .zeros ((1 , 3 , 512 , 512 ), dtype = np .float32 )
144+ #img.ravel()[:] = np.arange(3*512*512)
145+ core = np .zeros ((32 , 3 , 3 , 3 ), dtype = np .float32 )
146+ #core.ravel()[:] = np.arange(3*3*3*32)
147+
148+ rst1 = conv (img , core , (1 ,1 ))
149+ start = time ()
150+ rst1 = conv (img , core , (1 ,1 ))
151+ print ('jit cost:' , time ()- start )
152+
153+ start = time ()
154+ rst2 = conv (img , core , (1 ,1 ))
155+ print ('numpy cost:' , time ()- start )
0 commit comments