matlab 图像插interpn,在两个图像之间插值

fromscipy.ndimage.morphologyimportdistance_transform_edtfromscipy.interpolateimportinterpndefndgrid(*args,**kwargs):"""Same as calling ``meshgrid`` with *indexing* = ``'ij'`` (see``meshgrid`` for documentation)."""kwargs['indexing']='ij'returnnp.meshgrid(*args,**kwargs)defbwperim(bw,n=4):"""perim = bwperim(bw, n=4)Find the perimeter of objects in binary images.A pixel is part of an object perimeter if its value is one and thereis at least one zero-valued pixel in its neighborhood.By default the neighborhood of a pixel is 4 nearest pixels, butif `n` is set to 8 the 8 nearest pixels will be considered.Parameters----------bw : A black-and-white imagen : Connectivity. Must be 4 or 8 (default: 8)Returns-------perim : A boolean imageFrom Mahotas: http://nullege.com/codes/search/mahotas.bwperim"""ifnnotin(4,8):raiseValueError('mahotas.bwperim: n must be 4 or 8')rows,cols=bw.shape# Translate image by one pixel in all directionsnorth=np.zeros((rows,cols))south=np.zeros((rows,cols))west=np.zeros((rows,cols))east=np.zeros((rows,cols))north[:-1,:]=bw[1:,:]south[1:,:]=bw[:-1,:]west[:,:-1]=bw[:,1:]east[:,1:]=bw[:,:-1]idx=(north==bw)&\(south==bw)&\(west==bw)&\(east==bw)ifn==8:north_east=np.zeros((rows,cols))north_west=np.zeros((rows,cols))south_east=np.zeros((rows,cols))south_west=np.zeros((rows,cols))north_east[:-1,1:]=bw[1:,:-1]north_west[:-1,:-1]=bw[1:,1:]south_east[1:,1:]=bw[:-1,:-1]south_west[1:,:-1]=bw[:-1,1:]idx&=(north_east==bw)&\(south_east==bw)&\(south_west==bw)&\(north_west==bw)return~idx*bwdefsigned_bwdist(im):'''Find perim and return masked image (signed/reversed)'''im=-bwdist(bwperim(im))*np.logical_not(im)+bwdist(bwperim(im))*imreturnimdefbwdist(im):'''Find distance map of image'''dist_im=distance_transform_edt(1-im)returndist_imdefinterp_shape(top,bottom,num):ifnum<0andround(num)==num:print("Error: number of slices to be interpolated must be integer>0")top=signed_bwdist(top)bottom=signed_bwdist(bottom)r,c=top.shapet=num+2print("Rows - Cols - Slices")print(r,c,t)print("")# rejoin top, bottom into a single array of shape (2, r, c)# MATLAB: cat(3,bottom,top)top_and_bottom=np.r_['0,3',top,bottom]#top_and_bottom = np.rollaxis(top_and_bottom, 0, 3)# create ndgridsx,y,z=np.mgrid[0:r,0:c,0:t-1]# existing datax1,y1,z1=np.mgrid[0:r,0:c,0:t]# including new sliceprint("Shape x y z:",x.shape,y.shape,z.shape)print("Shape x1 y1 z1:",x1.shape,y1.shape,z1.shape)print(top_and_bottom.shape,len(x),len(y),len(z))# Do interpolationout=interpn((x,y,z),top_and_bottom,(x1,y1,z1))# MATLAB: out = out(:,:,2:end-1)>=0;array_lim=out[-1]-1out[out[:,:,2:out]>=0]=1returnout