我们从Python开源项目中，提取了以下50个代码示例，用于说明如何使用numpy.gradient()。
def get_mag_ang(img): """ Gets image gradient (magnitude) and orientation (angle) Args: img Returns: Gradient, orientation """ img = np.sqrt(img) gx = cv2.Sobel(np.float32(img), cv2.CV_32F, 1, 0) gy = cv2.Sobel(np.float32(img), cv2.CV_32F, 0, 1) mag, ang = cv2.cartToPolar(gx, gy) return mag, ang, gx, gy
def calculate_angle_quality_thread(self, voi, gantry, couch, calculate_from=0, stepsize=1.0, q=None, avoid=[], voi_cube=None, gradient=True): """ TODO: Documentation """ os.nice(1) for couch_angle in couch: qual = self.calculate_angle_quality(voi, gantry, couch_angle, calculate_from, stepsize, avoid, voi_cube, gradient) q.put({"couch": couch_angle, "gantry": gantry, "data": qual})
def test_model_monovariate(func, var, par, k): model = Model(func, var, par) x, dx = np.linspace(0, 10, 100, retstep=True, endpoint=False) U = np.cos(x * 2 * np.pi / 10) fields = model.fields_template(x=x, U=U) parameters = dict(periodic=True, k=k) F = model.F(fields, parameters) J_sparse = model.J(fields, parameters) J_dense = model.J(fields, parameters, sparse=False) J_approx = model.F.diff_approx(fields, parameters) dxU = np.gradient(np.pad(U, 2, mode='wrap')) / dx dxxU = np.gradient(dxU) / dx dxU = dxU[2:-2] dxxU = dxxU[2:-2] assert np.isclose(F, k * dxxU, atol=1E-3).all() assert np.isclose(J_approx, J_sparse.todense(), atol=1E-3).all() assert np.isclose(J_approx, J_dense, atol=1E-3).all()
def test_model_bivariate(func, var, par): model = Model(func, var, par) x, dx = np.linspace(0, 10, 100, retstep=True, endpoint=False) u = np.cos(x * 2 * np.pi / 10) v = np.sin(x * 2 * np.pi / 10) fields = model.fields_template(x=x, u=u, v=v) parameters = dict(periodic=True, k1=1, k2=1) F = model.F(fields, parameters) J_sparse = model.J(fields, parameters) J_dense = model.J(fields, parameters, sparse=False) J_approx = model.F.diff_approx(fields, parameters) dxu = np.gradient(np.pad(u, 2, mode='wrap')) / dx dxxu = np.gradient(dxu) / dx dxu = dxu[2:-2] dxxu = dxxu[2:-2] dxv = np.gradient(np.pad(v, 2, mode='wrap')) / dx dxxv = np.gradient(dxv) / dx dxv = dxv[2:-2] dxxv = dxxv[2:-2] assert np.isclose(F, np.vstack([dxxv, dxxu]).flatten('F'), atol=1E-3).all() assert np.isclose(J_approx, J_sparse.todense(), atol=1E-3).all() assert np.isclose(J_approx, J_dense, atol=1E-3).all()
def inflection_points(points, axis, span): ''' Find the list of vertices that preceed inflection points in a curve. The curve is differentiated with respect to the coordinate system defined by axis and span. axis: A vector representing the vertical axis of the coordinate system. span: A vector representing the the horiztonal axis of the coordinate system. returns: a list of points in space corresponding to the vertices that immediately preceed inflection points in the curve ''' coords_on_span = points.dot(span) dx = np.gradient(coords_on_span) coords_on_axis = points.dot(axis) # Take the second order finite difference of the curve with respect to the # defined coordinate system finite_difference_2 = np.gradient(np.gradient(coords_on_axis, dx), dx) # Compare the product of all neighboring pairs of points in the second derivative # If a pair of points has a negative product, then the second derivative changes sign # at one of those points, signalling an inflection point is_inflection_point = [finite_difference_2[i] * finite_difference_2[i + 1] <= 0 for i in range(len(finite_difference_2) - 1)] inflection_point_indices = [i for i, b in enumerate(is_inflection_point) if b] if len(inflection_point_indices) == 0: # pylint: disable=len-as-condition return [] return points[inflection_point_indices]
def calculate_sift_grid(self, image, gridH, gridW): H, W = image.shape Npatches = gridH.size feaArr = np.zeros((Npatches, Nsamples * Nangles)) # calculate gradient GH, GW = gen_dgauss(self.sigma) IH = signal.convolve2d(image, GH, mode='same') IW = signal.convolve2d(image, GW, mode='same') Imag = np.sqrt(IH ** 2 + IW ** 2) Itheta = np.arctan2(IH, IW) Iorient = np.zeros((Nangles, H, W)) for i in range(Nangles): Iorient[i] = Imag * np.maximum(np.cos(Itheta - angles[i]) ** alpha, 0) for i in range(Npatches): currFeature = np.zeros((Nangles, Nsamples)) for j in range(Nangles): currFeature[j] = np.dot(self.weights, \ Iorient[j, gridH[i]:gridH[i] + self.pS, gridW[i]:gridW[i] + self.pS].flatten()) feaArr[i] = currFeature.flatten() return feaArr
def extract_sift_patches(self, image, gridH, gridW): # extracts the sift descriptor of patches # in positions (gridH,gridW) in the image H, W = image.shape Npatches = gridH.size feaArr = np.zeros((Npatches, Nsamples * Nangles)) # calculate gradient GH, GW = gen_dgauss(self.sigma) IH = signal.convolve2d(image, GH, mode='same') IW = signal.convolve2d(image, GW, mode='same') Imag = np.sqrt(IH ** 2 + IW ** 2) Itheta = np.arctan2(IH, IW) Iorient = np.zeros((Nangles, H, W)) for i in range(Nangles): Iorient[i] = Imag * np.maximum(np.cos(Itheta - angles[i]) ** alpha, 0) for i in range(Npatches): currFeature = np.zeros((Nangles, Nsamples)) for j in range(Nangles): currFeature[j] = np.dot(self.weights, \ Iorient[j, gridH[i]:gridH[i] + self.pS, gridW[i]:gridW[i] + self.pS].flatten()) feaArr[i] = currFeature.flatten() # feaArr contains each descriptor in a row feaArr = self.normalize_sift(feaArr) return feaArr
def test_specific_axes(self): # Testing that gradient can work on a given axis only v = [[1, 1], [3, 4]] x = np.array(v) dx = [np.array([[2., 3.], [2., 3.]]), np.array([[0., 0.], [1., 1.]])] assert_array_equal(gradient(x, axis=0), dx[0]) assert_array_equal(gradient(x, axis=1), dx[1]) assert_array_equal(gradient(x, axis=-1), dx[1]) assert_array_equal(gradient(x, axis=(1, 0)), [dx[1], dx[0]]) # test axis=None which means all axes assert_almost_equal(gradient(x, axis=None), [dx[0], dx[1]]) # and is the same as no axis keyword given assert_almost_equal(gradient(x, axis=None), gradient(x)) # test vararg order assert_array_equal(gradient(x, 2, 3, axis=(1, 0)), [dx[1]/2.0, dx[0]/3.0]) # test maximal number of varargs assert_raises(SyntaxError, gradient, x, 1, 2, axis=1) assert_raises(ValueError, gradient, x, axis=3) assert_raises(ValueError, gradient, x, axis=-3) assert_raises(TypeError, gradient, x, axis=[1,])
def contour_from_array_at_label_l(im_arr, l, thr=0.3, omit_axis=None, verbose=0): if verbose > 0: print(l) array_label_l = im_arr == l assert isinstance(array_label_l, np.ndarray) gra = np.gradient(array_label_l.astype(np.bool).astype(np.float64)) if omit_axis is None: norm_gra = np.sqrt(gra[0] ** 2 + gra[1] ** 2 + gra[2] ** 2) > thr elif omit_axis == 'x': norm_gra = np.sqrt(gra[1] ** 2 + gra[2] ** 2) > thr elif omit_axis == 'y': norm_gra = np.sqrt(gra[0] ** 2 + gra[2] ** 2) > thr elif omit_axis == 'z': norm_gra = np.sqrt(gra[0] ** 2 + gra[1] ** 2) > thr else: raise IOError return norm_gra
def computeWeightsLocallyNormalized(I, centered_gradient=True, norm_radius=45): h,w = I.shape[:2] if centered_gradient: gy,gx = np.gradient(I)[:2] gysq = (gy**2).mean(axis=2) if gy.ndim > 2 else gy**2 gxsq = (gx**2).mean(axis=2) if gx.ndim > 2 else gx**2 gxsq_local_mean = cv2.blur(gxsq, ksize=(norm_radius, norm_radius)) gysq_local_mean = cv2.blur(gysq, ksize=(norm_radius, norm_radius)) w_horizontal = np.exp( - gxsq * 1.0/(2*np.maximum(1e-6, gxsq_local_mean))) w_vertical = np.exp( - gysq * 1.0/(2*np.maximum(1e-6, gysq_local_mean))) else: raise Exception("NotImplementedYet") return w_horizontal, w_vertical
def rel_vort(u,v,dx,dy): #get the gradient of the wind in the u- and v-directions du = np.gradient(u) dv = np.gradient(v) #compute the relative vorticity (units : 10^-5 s^-1) vort = ((dv[-1]/dx) - (du[-2]/dy))*pow(10,5) #return the vorticity (Units: 10^-5 s-1) return vort ############################################################################### ########################## End function rel_vort() ########################## ############################################################################### #################### 25.Begin function of wind_chill() ###################### ## Required libraries: numpy # ## # ## Inputs: t2m = ndarray of 2-meter temperature values (Units: F) # ## # ## wind = ndarray of 10-meter bulk wind values (Units: MPH) # ## # ###############################################################################
def richardson(self): u = self.dataSet.readNCVariable('U') v = self.dataSet.readNCVariable('V') u_corr = wrf.unstaggerX(u) v_corr = wrf.unstaggerY(v) wind = (u_corr*u_corr + v_corr*v_corr)**(0.5) height = wrf.unstaggerZ(self.height) #compute the vertical gradient of theta dtheta = np.gradient(self.theta)[0] du = np.gradient(wind)[0] dz = np.gradient(height)[0] #compute the richardson number self.u10 = u_corr self.v10 = v_corr self.var = ((self.g/self.theta)*(dtheta/dz))/(du/dz)**(2) self.var2 = self.var self.varTitle = "Richardson Number\n" + self.dataSet.getTime() #Set short variable title for time series self.sTitle = "Richardson Number"
def smooth_abs(x, dx=0.01): """smoothed version of absolute vaue function, with quadratic instead of sharp bottom. Derivative w.r.t. variable of interest. Width of quadratic can be controlled""" x, n = _checkIfFloat(x) y = np.abs(x) idx = np.logical_and(x > -dx, x < dx) y[idx] = x[idx]**2/(2.0*dx) + dx/2.0 # gradient dydx = np.ones_like(x) dydx[x <= -dx] = -1.0 dydx[idx] = x[idx]/dx if n == 1: y = y[0] dydx = dydx[0] return y, dydx
def main(): """Load image, compute derivatives, plot.""" img = data.camera() imgy_np, imgx_np = np.gradient(img) imgx_ski, imgy_ski = _compute_derivatives(img) # dx util.plot_images_grayscale( [img, dx_symmetric(img), dx_forward(img), dx_backward(img), imgx_np, imgx_ski], ["Image", "dx (symmetric)", "dx (forward)", "dx (backward)", "Ground Truth (numpy)", "Ground Truth (scikit-image)"] ) # dy util.plot_images_grayscale( [img, dy_symmetric(img), dy_forward(img), dy_backward(img), imgy_np, imgy_ski], ["Image", "dy (symmetric)", "dy (forward)", "dy (backward)", "Ground Truth (numpy)", "Ground Truth (scikit-image)"] )
def direction_map(dmap: DistanceMap): r"""Computes normalized gradient of distance map. Not defined when length of the gradient is zero. .. math:: \hat{\mathbf{e}}_{S} = -\frac{\nabla S(\mathbf{x})}{\| \nabla S(\mathbf{x}) \|} Args: dmap (numpy.ndarray): Distance map. Returns: DirectionMap: Direction map. """ u, v = np.gradient(dmap) l = np.hypot(u, v) # Avoids zero division l[l == 0] = np.nan # Flip order from (row, col) to (x, y) return v / l, u / l # Potentials
def __entropy(self, spectrum, m): """Calculates get m-th derivative of the real part of spectrum and returns entropy of its absolute value. """ assert type(m) is int, 'm should be a (possitive) integer' assert m > 0, 'need m > 0' #get the real part of the spectrum spect = np.array(spectrum) spect = spect.real # calculate the m-th derivative of the real part of the spectrum spectrumDerivative = spect for i in range(m): spectrumDerivative = np.gradient(spectrumDerivative) # now get the entropy of the abslolute value of the m-th derivative: entropy = sp.stats.entropy(np.abs(spectrumDerivative)) return entropy
def threshold_gradients(self, grad_thresh): """Creates a new DepthImage by zeroing out all depths where the magnitude of the gradient at that point is greater than grad_thresh. Parameters ---------- grad_thresh : float A threshold for the gradient magnitude. Returns ------- :obj:`DepthImage` A new DepthImage created from the thresholding operation. """ data = np.copy(self._data) gx, gy = self.gradients() gradients = np.zeros([gx.shape[0], gx.shape[1], 2]) gradients[:, :, 0] = gx gradients[:, :, 1] = gy gradient_mags = np.linalg.norm(gradients, axis=2) ind = np.where(gradient_mags > grad_thresh) data[ind[0], ind[1]] = 0.0 return DepthImage(data, self._frame)
def normal_cloud_im(self): """Generate a NormalCloudImage from the PointCloudImage. Returns ------- :obj:`NormalCloudImage` The corresponding NormalCloudImage. """ gx, gy, _ = np.gradient(self.data) gx_data = gx.reshape(self.height * self.width, 3) gy_data = gy.reshape(self.height * self.width, 3) pc_grads = np.cross(gx_data, gy_data) # default to point toward camera pc_grad_norms = np.linalg.norm(pc_grads, axis=1) normal_data = pc_grads / np.tile(pc_grad_norms[:, np.newaxis], [1, 3]) normal_im_data = normal_data.reshape(self.height, self.width, 3) return NormalCloudImage(normal_im_data, frame=self.frame)
def visualize_derivatives(image): ''' Plot gradient on left and Laplacian on right. Only tested on 2D 1-channel float imags ''' dx1,dy1 = np.gradient(image) gradient = dx1 + 1j*dy1 a1 = np.abs(gradient) plt.figure(None,(12,6)) plt.subplot(121) a1 = mean_center(blur(exposure.equalize_hist(unitize(a1)),1)) plt.imshow(a1, origin='lower',interpolation='nearest',cmap='gray',extent=(0,64,)*2) plt.title('Gradient Magnitude') plt.subplot(122) laplacian = scipy.ndimage.filters.laplace(image) lhist = mean_center( blur(exposure.equalize_hist(unitize(laplacian)),1)) plt.imshow(lhist, origin='lower',interpolation='nearest',cmap='gray',extent=(0,64,)*2) plt.title('Laplacian') return gradient, laplacian
def _get_displacement_function(self, f): """ Getter function for calculating the displacement. :param f: field that is used for the displacement, mainly velocity components :returns: function of the Taylor expanded field to first order """ dx = self._distance f_x, f_y = np.gradient(f , dx) f_xx, f_xy = np.gradient(f_x, dx) f_yx, f_yy = np.gradient(f_y, dx) return lambda i, j, x, y : (f[i, j] + x*f_x[i, j] + y*f_y[i, j] + 0.5*(f_xx[i, j]*x**2 + 2*f_xy[i, j]*x*y + f_yy[i, j]*y**2)) #For the bilinear method the build in scipy method `map_coordinates <https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.ndimage.interpolation.map_coordinates.html>`_ is used with *order* set to 1.
def psf_shift(lenslet_efl, dx, dy, mag=1): ''' Computes the shift of a PSF, in microns. Args: lenslet_efl (`microns`): EFL of lenslets. dx (`np.ndarray`): dx gradient of wavefront. dy (`np.ndarray`): dy gradient of wavefront. mag (`float`): magnification of the collimation system. Returns: `numpy.ndarray`: array of PSF shifts. Notes: see eq. 12 of Chanan, "Principles of Wavefront Sensing and Reconstruction" delta = m * fl * grad(z) m is magnification of SH system fl is lenslet focal length grad(z) is the x, y gradient of the opd, z, which is expressed in radians. ''' coef = -mag * lenslet_efl return coef * dx, coef * dy
def test_basic(self): v = [[1, 1], [3, 4]] x = np.array(v) dx = [np.array([[2., 3.], [2., 3.]]), np.array([[0., 0.], [1., 1.]])] assert_array_equal(gradient(x), dx) assert_array_equal(gradient(v), dx)
def test_badargs(self): # for 2D array, gradient can take 0, 1, or 2 extra args x = np.array([[1, 1], [3, 4]]) assert_raises(SyntaxError, gradient, x, np.array([1., 1.]), np.array([1., 1.]), np.array([1., 1.]))
def test_masked(self): # Make sure that gradient supports subclasses like masked arrays x = np.ma.array([[1, 1], [3, 4]], mask=[[False, False], [False, False]]) out = gradient(x)[0] assert_equal(type(out), type(x)) # And make sure that the output and input don't have aliased mask # arrays assert_(x.mask is not out.mask) # Also check that edge_order=2 doesn't alter the original mask x2 = np.ma.arange(5) x2[2] = np.ma.masked np.gradient(x2, edge_order=2) assert_array_equal(x2.mask, [False, False, True, False, False])
def test_datetime64(self): # Make sure gradient() can handle special types like datetime64 x = np.array( ['1910-08-16', '1910-08-11', '1910-08-10', '1910-08-12', '1910-10-12', '1910-12-12', '1912-12-12'], dtype='datetime64[D]') dx = np.array( [-5, -3, 0, 31, 61, 396, 731], dtype='timedelta64[D]') assert_array_equal(gradient(x), dx) assert_(dx.dtype == np.dtype('timedelta64[D]'))
def test_timedelta64(self): # Make sure gradient() can handle special types like timedelta64 x = np.array( [-5, -3, 10, 12, 61, 321, 300], dtype='timedelta64[D]') dx = np.array( [2, 7, 7, 25, 154, 119, -21], dtype='timedelta64[D]') assert_array_equal(gradient(x), dx) assert_(dx.dtype == np.dtype('timedelta64[D]'))
def calc_v_vec(self, tp): v_vec_all_bands = [] for ib in range(self.num_bands[tp]): v_vec_k_ordered = self.velocity_signed[tp][ib][self.pos_idx[tp]] v_vec_all_bands.append(self.grid_from_ordered_list(v_vec_k_ordered, tp, none_missing=True)) return np.array(v_vec_all_bands) # def calc_v_vec(self, tp): # # TODO: Take into account the fact that this gradient is found in three directions specified by the lattice, not # # the x, y, and z directions. It must be corrected to account for this. # energy_grid = self.array_from_kgrid('energy', tp) # # print('energy:') # # np.set_printoptions(precision=3) # # print(energy_grid[0,:,:,:,0]) # N = self.kgrid_array[tp].shape # k_grid = self.kgrid_array[tp] # v_vec_result = [] # for ib in range(self.num_bands[tp]): # v_vec = np.gradient(energy_grid[ib][:,:,:,0], k_grid[:,0,0,0] * self._rec_lattice.a, k_grid[0,:,0,1] * self._rec_lattice.b, k_grid[0,0,:,2] * self._rec_lattice.c) # v_vec_rearranged = np.zeros((N[0], N[1], N[2], 3)) # for i in range(N[0]): # for j in range(N[1]): # for k in range(N[2]): # v_vec_rearranged[i,j,k,:] = np.array([v_vec[0][i,j,k], v_vec[1][i,j,k], v_vec[2][i,j,k]]) # v_vec_rearranged *= A_to_m * m_to_cm / hbar # v_vec_result.append(v_vec_rearranged) # return np.array(v_vec_result) # turns a kgrid property into a list of grid arrays of that property for k integration
def compute_beta(TT, minT): """ This function computes the volumetric thermal expansion as a numerical derivative of the volume as a function of temperature V(T) given in the input array *minT*. This array can obtained from the free energy minimization which should be done before. """ grad=np.gradient(np.array(minT)) # numerical derivatives with numpy betaT = np.array(grad) # grad contains the derivatives with respect to T # also convert to np.array format return betaT/minT
def compute_alpha(minT,ibrav): """ This function calculates the thermal expansion alphaT at different temperatures from the input minT matrix by computing the numerical derivatives with numpy. The input matrix minT has shape nT*6, where the first index is the temperature and the second the lattice parameter. For example, minT[i,0] and minT[i,2] are the lattice parameters a and c at the temperature i. More ibrav types must be implemented """ grad=np.gradient(np.array(minT)) # numerical derivatives with numpy alphaT = np.array(grad[0]) # grad[0] contains the derivatives with respect to T, which is the first axis in minT # also convert to np.array format # now normalize the alpha properly. It must be different for different ibrav # to avoid a divide by 0 error (minT is zero for lattice parameters not defined # in the system) if ibrav==4: alphaT[:,0] = alphaT[:,0]/minT[:,0] alphaT[:,2] = alphaT[:,2]/minT[:,2] return alphaT ################################################################################
def calc_eud(dvh, a): v = -np.gradient(dvh) dose_bins = np.linspace(1, np.size(dvh), np.size(dvh)) dose_bins = np.round(dose_bins, 3) bin_centers = dose_bins - 0.5 eud = np.power(np.sum(np.multiply(v, np.power(bin_centers, a))), 1 / float(a)) eud = np.round(eud, 2) * 0.01 return eud
def _compute_gradients(self): """Computes the gradients of the SDF. Returns ------- :obj:`list` of :obj:`numpy.ndarray` of float A list of ndarrays of the same dimension as the SDF. The arrays are in axis order and specify the gradients for that axis at each point. """ self.gradients_ = np.gradient(self.data_)
def curvature(self, coords, delta=0.001): """ Returns an approximation to the local SDF curvature (Hessian) at the given coordinate in grid basis. Parameters --------- coords : numpy 3-vector the grid coordinates at which to get the curvature Returns ------- curvature : 3x3 ndarray of the curvature at the surface points """ # perturb local coords coords_x_up = coords + np.array([delta, 0, 0]) coords_x_down = coords + np.array([-delta, 0, 0]) coords_y_up = coords + np.array([0, delta, 0]) coords_y_down = coords + np.array([0, -delta, 0]) coords_z_up = coords + np.array([0, 0, delta]) coords_z_down = coords + np.array([0, 0, -delta]) # get gradient grad_x_up = self.gradient(coords_x_up) grad_x_down = self.gradient(coords_x_down) grad_y_up = self.gradient(coords_y_up) grad_y_down = self.gradient(coords_y_down) grad_z_up = self.gradient(coords_z_up) grad_z_down = self.gradient(coords_z_down) # finite differences curvature_x = (grad_x_up - grad_x_down) / (4 * delta) curvature_y = (grad_y_up - grad_y_down) / (4 * delta) curvature_z = (grad_z_up - grad_z_down) / (4 * delta) curvature = np.c_[curvature_x, np.c_[curvature_y, curvature_z]] curvature = curvature + curvature.T return curvature
def call_tad_borders(ii_results, cutoff=0): """ Calls TAD borders using the first derivative of the insulation index. :param ii_results: (raw) insulation index results, numpy vector :param cutoff: raw insulation index threshold for "true" TAD boundaries """ tad_borders = [] g = np.gradient(ii_results) for i in range(len(ii_results)): border_type = _border_type(i, g) if border_type == 1 and ii_results[i] <= cutoff: tad_borders.append(i) return tad_borders
def orientation_field(bmp, sigma=3): # Author: Shaun Harker, 2016 # Based on algorithm by Bazen and Gerez from "Systematic methods for the # computation of the directional fields and singular points of fingerprints," 2002. """ Computes orientation field (result everywhere between -pi/2 and pi/2) from the given vector field. """ u = bmp.astype(float) du = np.gradient(u) [ux, uy] = du Y = scipy.ndimage.filters.gaussian_filter(2.0*ux*uy, sigma=sigma) X = scipy.ndimage.filters.gaussian_filter(ux**2.0 - uy**2.0, sigma=sigma) return .5 * np.arctan2(Y, X)
def calculate_quality_grid(self, voi, gantry, couch, calculate_from=0, stepsize=1.0, avoid=[], gradient=True): """ TODO: Documentation """ result = self.calculate_quality_list(voi, gantry, couch, calculate_from, stepsize, avoid=avoid, gradient=gradient) result = sorted(result, key=cmp_to_key(cmp_sort)) grid_data = [] for x in result: grid_data.append(x["data"][0]) result = np.reshape(grid_data, (len(gantry), len(couch))) return result
def calculate_quality_list(self, voi, gantry, couch, calculate_from=0, stepsize=1.0, avoid=[], gradient=True): """ TODO: Documentation """ q = Queue(32767) process = [] d = voi.get_voi_cube() d.cube = np.array(d.cube, dtype=np.float32) voi_cube = DensityProjections(d) result = [] for gantry_angle in gantry: p = Process( target=self.calculate_angle_quality_thread, args=(voi, gantry_angle, couch, calculate_from, stepsize, q, avoid, voi_cube, gradient)) p.start() p.deamon = True process.append(p) if len(process) > 2: tmp = q.get() result.append(tmp) for p in process: if not p.is_alive(): process.remove(p) while not len(result) == len(gantry) * len(couch): tmp = q.get() result.append(tmp) return result
def calculate_angle_quality(self, voi, gantry, couch, calculate_from=0, stepsize=1.0, avoid=[], voi_cube=None, gradient=True): """ Calculates a quality index for a given gantry/couch combination. """ voi_min, voi_max = voi.get_min_max() for v in avoid: v_min, v_max = v.get_min_max() if voi_cube is None: d = voi.get_voi_cube() d.cube = np.array(d.cube, dtype=np.float32) voi_cube = DensityProjections(d) data, start, basis = self.calculate_projection(voi, gantry, couch, calculate_from, stepsize) voi_proj, t1, t2 = voi_cube.calculate_projection(voi, gantry, couch, 1, stepsize) if gradient: gradient = np.gradient(data) data = (gradient[0]**2 + gradient[1]**2)**0.5 a = data * (voi_proj > 0.0) quality = sum(a) area = sum(voi_proj > 0.0) # ~ area = sum(data>0.0)/10 if gradient: mean_quality = 10 - abs(quality / area) else: mean_quality = abs(quality / area) return mean_quality, quality, area
def get_direction(asa, *, sigma=None): """Return epochs during which an animal was running left to right, or right to left. Parameters ---------- asa : AnalogSignalArray 1D AnalogSignalArray containing the 1D position data. sigma : float, optional Smoothing to apply to position (x) before computing gradient estimate. Default is 0. Returns ------- l2r, r2l : EpochArrays EpochArrays corresponding to left-to-right and right-to-left movement. """ if sigma is None: sigma = 0 if not isinstance(asa, core.AnalogSignalArray): raise TypeError('AnalogSignalArray expected!') assert asa.n_signals == 1, "1D AnalogSignalArray expected!" direction = dxdt_AnalogSignalArray(asa.smooth(sigma=sigma), rectify=False).ydata direction[direction>=0] = 1 direction[direction<0] = -1 direction = direction.squeeze() l2r = get_contiguous_segments(np.argwhere(direction>0).squeeze(), step=1) l2r[:,1] -= 1 # change bounds from [inclusive, exclusive] to [inclusive, inclusive] l2r = core.EpochArray(asa.time[l2r]) r2l = get_contiguous_segments(np.argwhere(direction<0).squeeze(), step=1) r2l[:,1] -= 1 # change bounds from [inclusive, exclusive] to [inclusive, inclusive] r2l = core.EpochArray(asa.time[r2l]) return l2r, r2l
def hessian(array): (dy, dx) = np.gradient(array) (dydy, dxdy) = np.gradient(dy) (dydx, dxdx) = np.gradient(dx) return np.dstack((dxdx, dydx, dxdy, dydy))
def gen_dgauss(sigma): ''' gradient of the gaussian on both directions. ''' fwid = np.int(2 * np.ceil(sigma)) G = np.array(range(-fwid, fwid + 1)) ** 2 G = G.reshape((G.size, 1)) + G G = np.exp(- G / 2.0 / sigma / sigma) G /= np.sum(G) GH, GW = np.gradient(G) GH *= 2.0 / np.sum(np.abs(GH)) GW *= 2.0 / np.sum(np.abs(GW)) return GH, GW
def gradient_pdf(self, x): """Return the gradient of density of the joint prior at x.""" raise NotImplementedError
def gradient_logpdf(self, x, stepsize=None): """Return the gradient of log density of the joint prior at x. Parameters ---------- x : float or np.ndarray stepsize : float or list Stepsize or stepsizes for the dimensions """ x = np.asanyarray(x) ndim = x.ndim x = x.reshape((-1, self.dim)) grads = np.zeros_like(x) for i in range(len(grads)): xi = x[i] grads[i] = numgrad(self.logpdf, xi, h=stepsize) grads[np.isinf(grads)] = 0 grads[np.isnan(grads)] = 0 if ndim == 0 or (ndim == 1 and self.dim > 1): grads = grads[0] return grads
