Python numpy.arctan2() Examples
The following are 30
code examples of numpy.arctan2().
You can vote up the ones you like or vote down the ones you don't like,
and go to the original project or source file by following the links above each example.
You may also want to check out all available functions/classes of the module
numpy
, or try the search function
.
Example #1
Source File: point_cloud.py From FRIDA with MIT License | 7 votes |
def doa(self, receiver, source): ''' Computes the direction of arrival wrt a source and receiver ''' s_ind = self.key2ind(source) r_ind = self.key2ind(receiver) # vector from receiver to source v = self.X[:,s_ind] - self.X[:,r_ind] azimuth = np.arctan2(v[1], v[0]) elevation = np.arctan2(v[2], la.norm(v[:2])) azimuth = azimuth + 2*np.pi if azimuth < 0. else azimuth elevation = elevation + 2*np.pi if elevation < 0. else elevation return np.array([azimuth, elevation])
Example #2
Source File: util.py From neuropythy with GNU Affero General Public License v3.0 | 7 votes |
def vector_angle(u, v, direction=None): ''' vector_angle(u, v) yields the angle between the two vectors u and v. The optional argument direction is by default None, which specifies that the smallest possible angle between the vectors be reported; if the vectors u and v are 2D vectors and direction parameters True and False specify the clockwise or counter-clockwise directions, respectively; if the vectors are 3D vectors, then direction may be a 3D point that is not in the plane containing u, v, and the origin, and it specifies around which direction (u x v or v x u) the the counter-clockwise angle from u to v should be reported (the cross product vector that has a positive dot product with the direction argument is used as the rotation axis). ''' if direction is None: return np.arccos(vector_angle_cos(u, v)) elif direction is True: return np.arctan2(v[1], v[0]) - np.arctan2(u[1], u[0]) elif direction is False: return np.arctan2(u[1], u[0]) - np.arctan2(v[1], v[0]) else: axis1 = normalize(u) axis2 = normalize(np.cross(u, v)) if np.dot(axis2, direction) < 0: axis2 = -axis2 return np.arctan2(np.dot(axis2, v), np.dot(axis1, v))
Example #3
Source File: test_pose.py From SfmLearner-Pytorch with MIT License | 6 votes |
def compute_pose_error(gt, pred): RE = 0 snippet_length = gt.shape[0] scale_factor = np.sum(gt[:,:,-1] * pred[:,:,-1])/np.sum(pred[:,:,-1] ** 2) ATE = np.linalg.norm((gt[:,:,-1] - scale_factor * pred[:,:,-1]).reshape(-1)) for gt_pose, pred_pose in zip(gt, pred): # Residual matrix to which we compute angle's sin and cos R = gt_pose[:,:3] @ np.linalg.inv(pred_pose[:,:3]) s = np.linalg.norm([R[0,1]-R[1,0], R[1,2]-R[2,1], R[0,2]-R[2,0]]) c = np.trace(R) - 1 # Note: we actually compute double of cos and sin, but arctan2 is invariant to scale RE += np.arctan2(s,c) return ATE/snippet_length, RE/snippet_length
Example #4
Source File: geometry.py From connecting_the_dots with MIT License | 6 votes |
def xyz_from_rotm(R): R = R.reshape(-1,3,3) xyz = np.empty((R.shape[0],3), dtype=R.dtype) for bidx in range(R.shape[0]): if R[bidx,0,2] < 1: if R[bidx,0,2] > -1: xyz[bidx,1] = np.arcsin(R[bidx,0,2]) xyz[bidx,0] = np.arctan2(-R[bidx,1,2], R[bidx,2,2]) xyz[bidx,2] = np.arctan2(-R[bidx,0,1], R[bidx,0,0]) else: xyz[bidx,1] = -np.pi/2 xyz[bidx,0] = -np.arctan2(R[bidx,1,0],R[bidx,1,1]) xyz[bidx,2] = 0 else: xyz[bidx,1] = np.pi/2 xyz[bidx,0] = np.arctan2(R[bidx,1,0], R[bidx,1,1]) xyz[bidx,2] = 0 return xyz.squeeze()
Example #5
Source File: LSDMappingTools.py From LSDMappingTools with MIT License | 6 votes |
def Hillshade(raster_file, azimuth, angle_altitude): array = ReadRasterArrayBlocks(raster_file,raster_band=1) x, y = np.gradient(array) slope = np.pi/2. - np.arctan(np.sqrt(x*x + y*y)) aspect = np.arctan2(-x, y) azimuthrad = np.azimuth*np.pi / 180. altituderad = np.angle_altitude*np.pi / 180. shaded = np.sin(altituderad) * np.sin(slope)\ + np.cos(altituderad) * np.cos(slope)\ * np.cos(azimuthrad - aspect) return 255*(shaded + 1)/2 #==============================================================================
Example #6
Source File: geometry.py From connecting_the_dots with MIT License | 6 votes |
def zyx_from_rotm(R): R = R.reshape(-1,3,3) zyx = np.empty((R.shape[0],3), dtype=R.dtype) for bidx in range(R.shape[0]): if R[bidx,2,0] < 1: if R[bidx,2,0] > -1: zyx[bidx,1] = np.arcsin(-R[bidx,2,0]) zyx[bidx,0] = np.arctan2(R[bidx,1,0], R[bidx,0,0]) zyx[bidx,2] = np.arctan2(R[bidx,2,1], R[bidx,2,2]) else: zyx[bidx,1] = np.pi / 2 zyx[bidx,0] = -np.arctan2(-R[bidx,1,2], R[bidx,1,1]) zyx[bidx,2] = 0 else: zyx[bidx,1] = -np.pi / 2 zyx[bidx,0] = np.arctan2(-R[bidx,1,2], R[bidx,1,1]) zyx[bidx,2] = 0 return zyx.squeeze()
Example #7
Source File: pyoptic.py From tf-pose with Apache License 2.0 | 6 votes |
def propagateRay(self, ray): """Refract, reflect, absorb, and/or scatter ray. This function may create and return new rays""" surface = self.surfaces[0] p1, ai = surface.intersectRay(ray) if p1 is not None: p1 = surface.mapToItem(ray, p1) rd = ray['dir'] a1 = np.arctan2(rd[1], rd[0]) ar = a1 + np.pi - 2*ai ray.setEnd(p1) dp = Point(np.cos(ar), np.sin(ar)) ray = Ray(parent=ray, dir=dp) else: ray.setEnd(None) return [ray]
Example #8
Source File: geometry.py From connecting_the_dots with MIT License | 6 votes |
def axisangle_from_rotm(R): # logarithm of rotation matrix # R = R.reshape(-1,3,3) # tr = np.trace(R, axis1=1, axis2=2) # phi = np.arccos(np.clip((tr - 1) / 2, -1, 1)) # scale = np.zeros_like(phi) # div = 2 * np.sin(phi) # np.divide(phi, div, out=scale, where=np.abs(div) > 1e-6) # A = (R - R.transpose(0,2,1)) * scale.reshape(-1,1,1) # aa = np.stack((A[:,2,1], A[:,0,2], A[:,1,0]), axis=1) # return aa.squeeze() R = R.reshape(-1,3,3) omega = np.empty((R.shape[0], 3), dtype=R.dtype) omega[:,0] = R[:,2,1] - R[:,1,2] omega[:,1] = R[:,0,2] - R[:,2,0] omega[:,2] = R[:,1,0] - R[:,0,1] r = np.linalg.norm(omega, axis=1).reshape(-1,1) t = np.trace(R, axis1=1, axis2=2).reshape(-1,1) omega = np.arctan2(r, t-1) * omega aa = np.zeros_like(omega) np.divide(omega, r, out=aa, where=r != 0) return aa.squeeze()
Example #9
Source File: SRTTransform.py From tf-pose with Apache License 2.0 | 6 votes |
def setFromQTransform(self, tr): p1 = Point(tr.map(0., 0.)) p2 = Point(tr.map(1., 0.)) p3 = Point(tr.map(0., 1.)) dp2 = Point(p2-p1) dp3 = Point(p3-p1) ## detect flipped axes if dp2.angle(dp3) > 0: #da = 180 da = 0 sy = -1.0 else: da = 0 sy = 1.0 self._state = { 'pos': Point(p1), 'scale': Point(dp2.length(), dp3.length() * sy), 'angle': (np.arctan2(dp2[1], dp2[0]) * 180. / np.pi) + da } self.update()
Example #10
Source File: ROI.py From tf-pose with Apache License 2.0 | 6 votes |
def generateShape(self): dt = self.deviceTransform() if dt is None: self._shape = self.path return None v = dt.map(QtCore.QPointF(1, 0)) - dt.map(QtCore.QPointF(0, 0)) va = np.arctan2(v.y(), v.x()) dti = fn.invertQTransform(dt) devPos = dt.map(QtCore.QPointF(0,0)) tr = QtGui.QTransform() tr.translate(devPos.x(), devPos.y()) tr.rotate(va * 180. / 3.1415926) return dti.map(tr.map(self.path))
Example #11
Source File: MovementMap.py From poeai with MIT License | 6 votes |
def GetNeighborCells(self, p, nr, dp = None): ''' Returns all cells no more than a given distance in any direction from a specified cell p: The cell of which to get the neighbors nr: Neighbor radius dp: Direction preference ''' pi, pj, pk = p tqm = self.qm * self.qp nc = [(pi - i * tqm, pj - j * tqm, pk) for i in range(-nr, nr + 1) for j in range(-nr, nr + 1)] if dp is not None: #Sort points based on direction preference dpa = np.arctan2(dp[1], dp[0]) #Get angle of direction prefered #Prefer directions in the direction of dp; sort based on magnitude of angle from last direction nc = sorted(nc, key = lambda t : np.abs(np.arctan2(t[1], t[0]) - dpa)) return nc #Gets the current 3d position of the player
Example #12
Source File: sunrgbd_utils.py From H3DNet with MIT License | 6 votes |
def __init__(self, line): data = line.split(' ') data[1:] = [float(x) for x in data[1:]] self.classname = data[0] self.xmin = data[1] self.ymin = data[2] self.xmax = data[1]+data[3] self.ymax = data[2]+data[4] self.box2d = np.array([self.xmin,self.ymin,self.xmax,self.ymax]) self.centroid = np.array([data[5],data[6],data[7]]) self.unused_dimension = np.array([data[8],data[9],data[10]]) self.w = data[8] self.l = data[9] self.h = data[10] self.orientation = np.zeros((3,)) self.orientation[0] = data[11] self.orientation[1] = data[12] self.heading_angle = -1 * np.arctan2(self.orientation[1], self.orientation[0])
Example #13
Source File: analyses.py From PySpice with GNU General Public License v3.0 | 6 votes |
def do_ac_analysis(): circuit, n = simple_bjt_amp() simulator = circuit.simulator(temperature=25, nominal_temperature=25) analysis = simulator.ac(start_frequency=10, stop_frequency=1e6, number_of_points=100, variation='dec') gain = np.array(analysis[n.n5]) figure = plt.figure(1, (20, 10)) axe = plt.subplot(211) axe.grid(True) axe.set_xlabel("Frequency [Hz]") axe.set_ylabel("dB gain.") axe.semilogx(analysis.frequency, 20*np.log10(np.abs(gain))) axe = plt.subplot(212) axe.grid(True) axe.set_xlabel("Frequency [Hz]") axe.set_ylabel("Phase.") axe.semilogx(analysis.frequency, np.arctan2(gain.imag, gain.real)) plt.show()
Example #14
Source File: stanley_controller.py From RacingRobot with MIT License | 6 votes |
def stanleyControl(state, cx, cy, cyaw, last_target_idx): """ :param state: (State object) :param cx: ([float]) :param cy: ([float]) :param cyaw: ([float]) :param last_target_idx: (int) :return: (float, int, float) """ # Cross track error current_target_idx, error_front_axle = calcTargetIndex(state, cx, cy) if last_target_idx >= current_target_idx: current_target_idx = last_target_idx # theta_e corrects the heading error theta_e = normalizeAngle(cyaw[current_target_idx] - state.yaw) # theta_d corrects the cross track error theta_d = np.arctan2(K_STANLEY_CONTROL * error_front_axle, state.v) # Steering control delta = theta_e + theta_d return delta, current_target_idx, error_front_axle
Example #15
Source File: stanley_controller.py From RacingRobot with MIT License | 6 votes |
def calcTargetIndex(state, cx, cy): """ :param state: (State object) :param cx: [float] :param cy: [float] :return: (int, float) """ # Calc front axle position fx = state.x + CAR_LENGTH * np.cos(state.yaw) fy = state.y + CAR_LENGTH * np.sin(state.yaw) # Search nearest point index dx = [fx - icx for icx in cx] dy = [fy - icy for icy in cy] d = [np.sqrt(idx ** 2 + idy ** 2) for (idx, idy) in zip(dx, dy)] error_front_axle = min(d) target_idx = d.index(error_front_axle) target_yaw = normalizeAngle(np.arctan2(fy - cy[target_idx], fx - cx[target_idx]) - state.yaw) if target_yaw > 0.0: error_front_axle = - error_front_axle return target_idx, error_front_axle
Example #16
Source File: rpe.py From pyGSTi with Apache License 2.0 | 6 votes |
def raw_angles(self, measured): """ Determine the raw angles from the count data. This corresponds to the angle of U^N, i.e., it is N times the phase of U. """ angles = OrderedDict() # The ordering here is chosen to maintain compatibility. for N, (Cp_Ns, Cm_Ns, Cp_Nc, Cm_Nc) in measured.items(): # See the description of RobustPhaseEstimationDesign. # We estimate P^+_{Ns} and P^-_{Nc} from the similarly named counts. # The MLE for these probabilities is: Pp_Ns = Cp_Ns / (Cp_Ns + Cm_Ns) Pp_Nc = Cp_Nc / (Cp_Nc + Cm_Nc) angles[N] = numpy.arctan2(2 * Pp_Ns - 1, 2 * Pp_Nc - 1) % (2 * numpy.pi) return angles
Example #17
Source File: kincar-flatsys.py From python-control with BSD 3-Clause "New" or "Revised" License | 6 votes |
def vehicle_flat_reverse(zflag, params={}): # Get the parameter values b = params.get('wheelbase', 3.) # Create a vector to store the state and inputs x = np.zeros(3) u = np.zeros(2) # Given the flat variables, solve for the state x[0] = zflag[0][0] # x position x[1] = zflag[1][0] # y position x[2] = np.arctan2(zflag[1][1], zflag[0][1]) # tan(theta) = ydot/xdot # And next solve for the inputs u[0] = zflag[0][1] * np.cos(x[2]) + zflag[1][1] * np.sin(x[2]) thdot_v = zflag[1][2] * np.cos(x[2]) - zflag[0][2] * np.sin(x[2]) u[1] = np.arctan2(thdot_v, u[0]**2 / b) return x, u # Function to compute the RHS of the system dynamics
Example #18
Source File: kincar-flatsys.py From python-control with BSD 3-Clause "New" or "Revised" License | 6 votes |
def vehicle_flat_reverse(zflag, params={}): # Get the parameter values b = params.get('wheelbase', 3.) # Create a vector to store the state and inputs x = np.zeros(3) u = np.zeros(2) # Given the flat variables, solve for the state x[0] = zflag[0][0] # x position x[1] = zflag[1][0] # y position x[2] = np.arctan2(zflag[1][1], zflag[0][1]) # tan(theta) = ydot/xdot # And next solve for the inputs u[0] = zflag[0][1] * np.cos(x[2]) + zflag[1][1] * np.sin(x[2]) thdot_v = zflag[1][2] * np.cos(x[2]) - zflag[0][2] * np.sin(x[2]) u[1] = np.arctan2(thdot_v, u[0]**2 / b) return x, u # Function to compute the RHS of the system dynamics
Example #19
Source File: tools_fri_doa_plane.py From FRIDA with MIT License | 6 votes |
def mtx_freq2visi(M, p_mic_x, p_mic_y): """ build the matrix that maps the Fourier series to the visibility :param M: the Fourier series expansion is limited from -M to M :param p_mic_x: a vector that constains microphones x coordinates :param p_mic_y: a vector that constains microphones y coordinates :return: """ num_mic = p_mic_x.size ms = np.reshape(np.arange(-M, M + 1, step=1), (1, -1), order='F') G = np.zeros((num_mic * (num_mic - 1), 2 * M + 1), dtype=complex, order='C') count_G = 0 for q in range(num_mic): p_x_outer = p_mic_x[q] p_y_outer = p_mic_y[q] for qp in range(num_mic): if not q == qp: p_x_qqp = p_x_outer - p_mic_x[qp] p_y_qqp = p_y_outer - p_mic_y[qp] norm_p_qqp = np.sqrt(p_x_qqp ** 2 + p_y_qqp ** 2) phi_qqp = np.arctan2(p_y_qqp, p_x_qqp) G[count_G, :] = (-1j) ** ms * sp.special.jv(ms, norm_p_qqp) * \ np.exp(1j * ms * phi_qqp) count_G += 1 return G
Example #20
Source File: transform_util.py From lingvo with Apache License 2.0 | 5 votes |
def TransformHeading(transform, heading): """Compute 'heading' given transform. The heading provided as input is assumed to be in the original coordinate space. When the coordinate space undergoes a transformation (e.g., with CarToImageTransform), the heading in the new coordinate space must be recomputed. We compute this by deriving the formula for the angle of transformed unit vector defined by 'heading'. Args: transform: 4x4 numpy matrix used to convert from car to image coordinates. heading: Floating point scalar heading. Returns: Heading in the transformed coordinate system. """ x1, y1 = np.cos(heading), np.sin(heading) # Transform the unit ray. unit_ray = np.array([x1, y1, 0.0, 1.0]) transform_no_shift = CopyTransform(transform) transform_no_shift[0, 3] = 0 transform_no_shift[1, 3] = 0 transformed_ray = np.matmul(transform_no_shift, unit_ray) x2, y2 = transformed_ray[0:2] # Use arctan2 to compute the new rotation angle; note that arctan2 takes 'y' # and then 'x'. new_heading = np.arctan2(y2, x2) return new_heading
Example #21
Source File: test_core.py From recruit with Apache License 2.0 | 5 votes |
def test_testUfuncRegression(self): # Tests new ufuncs on MaskedArrays. for f in ['sqrt', 'log', 'log10', 'exp', 'conjugate', 'sin', 'cos', 'tan', 'arcsin', 'arccos', 'arctan', 'sinh', 'cosh', 'tanh', 'arcsinh', 'arccosh', 'arctanh', 'absolute', 'fabs', 'negative', 'floor', 'ceil', 'logical_not', 'add', 'subtract', 'multiply', 'divide', 'true_divide', 'floor_divide', 'remainder', 'fmod', 'hypot', 'arctan2', 'equal', 'not_equal', 'less_equal', 'greater_equal', 'less', 'greater', 'logical_and', 'logical_or', 'logical_xor', ]: try: uf = getattr(umath, f) except AttributeError: uf = getattr(fromnumeric, f) mf = getattr(numpy.ma.core, f) args = self.d[:uf.nin] ur = uf(*args) mr = mf(*args) assert_equal(ur.filled(0), mr.filled(0), f) assert_mask_equal(ur.mask, mr.mask, err_msg=f)
Example #22
Source File: _footprint.py From buzzard with Apache License 2.0 | 5 votes |
def angle(self): """Angle in degree: rotation used in the affine transformation, (0 is north-up)""" lrvec = self.lrvec return float(np.arctan2(lrvec[1], lrvec[0]) * 180. / np.pi)
Example #23
Source File: cmap.py From ehtplot with GNU General Public License v3.0 | 5 votes |
def gethue(color): """Get the hue of a color""" if isinstance(color, float): return color if is_color_like(color): RGB = to_rgba(color)[:3] else: raise ValueError("color is not color like") Jp, ap, bp = transform(np.array([RGB]))[0] hp = np.arctan2(bp, ap) * 180 / np.pi print("Decode color \"{}\"; h' = {}".format(color, hp)) return hp
Example #24
Source File: cmap.py From ehtplot with GNU General Public License v3.0 | 5 votes |
def _JChp(ax1, cmap): """Plot J', C', and h' of a colormap as function of the mapped value Args: ax1 (matplotlib.axes.Axes): The matplotlib Axes to be plot on. cmap (string or matplotlib.colors.Colormap): The colormap to be plotted. """ ctab = get_ctab(cmap) # get the colormap as a color table in sRGB Jpapbp = transform(ctab) # transform color table into CAM02-UCS colorspace Jp = Jpapbp[:,0] ap = Jpapbp[:,1] bp = Jpapbp[:,2] Cp = np.sqrt(ap * ap + bp * bp) hp = np.arctan2(bp, ap) * 180 / np.pi v = np.linspace(0.0, 1.0, len(Jp)) ax1.set_title(cmap if isinstance(cmap, basestring) else cmap.name) ax1.set_xlabel("Value") ax2 = ax1.twinx() ax1.set_ylim(0, 100) ax1.set_ylabel("J' & C' (0-100)") ax2.set_ylim(-180,180) ax2.set_ylabel("h' (degrees)") ax1.scatter(v, Jp, color=ctab) ax1.plot (v, Cp, c='k', linestyle='--') ax2.scatter(v[::15], hp[::15], s=12, c='k')
Example #25
Source File: create_table_pose.py From ObjectPoseEstimationSummary with MIT License | 5 votes |
def compute_angles_from_pose(R, t): #c = - np.dot(np.transpose(R), t) direction = -np.dot(np.transpose(R), np.array([0., 0., 1.]).reshape(3, 1)) azimuth = np.arctan2(direction[0], direction[1])[0] * 180./ np.pi elevation = np.arcsin(direction[-1])[0] * 180./ np.pi inplane = compute_inplane_from_rotation(R, direction.reshape(-1)) * 180./ np.pi return azimuth, elevation, inplane
Example #26
Source File: atlas.py From ibllib with MIT License | 5 votes |
def cart2sph(x, y, z): """ Converts cartesian to spherical Coordinates theta: polar angle, phi: azimuth """ r = np.sqrt(x ** 2 + y ** 2 + z ** 2) phi = np.arctan2(y, x) * 180 / np.pi theta = np.zeros_like(r) iok = r != 0 theta[iok] = np.arccos(z[iok] / r[iok]) * 180 / np.pi return r, theta, phi
Example #27
Source File: basic_observables.py From pytim with GNU General Public License v3.0 | 5 votes |
def _cartesian_to_spherical(self, cartesian): spherical = np.empty(cartesian.shape) spherical[:, 0] = np.linalg.norm(cartesian, axis=1) if len(self.dirmask) > 1: # 2d or 3d spherical[:, 1] = np.arctan2(cartesian[:, 1], cartesian[:, 0]) if len(self.dirmask) == 3: # 3d xy = np.sum(cartesian[:, [0, 1]]**2, axis=1) spherical[:, 2] = np.arctan2(np.sqrt(xy), cartesian[:, 2]) return spherical
Example #28
Source File: sasa.py From pytim with GNU General Public License v3.0 | 5 votes |
def _overlap(self, Ri, Rj, pij, dzi): dzj = pij[:, 2] - dzi ri2 = Ri**2 - dzi**2 rj2 = Rj**2 - dzj**2 cond = np.where(rj2 >= 0.0)[0] if len(cond) == 0: return [], [] rj2 = rj2[cond] pij = pij[cond] ri, rj = np.sqrt(ri2), np.sqrt(rj2) pij2 = pij**2 dij2 = pij2[:, 0] + pij2[:, 1] dij = np.sqrt(dij2) # slice within neighboring one if np.any(dij <= rj - ri): return np.array([2. * np.pi]), np.array([0.0]) # c1 => no overlap; c2 => neighboring slice enclosed c1, c2 = dij < ri + rj, dij > ri - rj cond = np.where(np.logical_and(c1, c2))[0] if len(cond) == 0: return [], [] arg = (ri2 + dij2 - rj2)[cond] / (ri * dij * 2.)[cond] alpha = 2. * np.arccos(arg) argx, argy = pij[:, 0], pij[:, 1] beta = np.arctan2(argx[cond], argy[cond]) return alpha, beta
Example #29
Source File: test_eval.py From recruit with Apache License 2.0 | 5 votes |
def test_df_use_case(self): df = DataFrame({'a': np.random.randn(10), 'b': np.random.randn(10)}) df.eval("e = arctan2(sin(a), b)", engine=self.engine, parser=self.parser, inplace=True) got = df.e expect = np.arctan2(np.sin(df.a), df.b) tm.assert_series_equal(got, expect, check_names=False)
Example #30
Source File: beyeler2019.py From pulse2percept with BSD 3-Clause "New" or "Revised" License | 5 votes |
def calc_bundle_tangent(self, xc, yc): """Calculates orientation of fiber bundle tangent at (xc, yc) Parameters ---------- xc, yc: float (x, y) location of point at which to calculate bundle orientation in microns. """ # Check for scalar: if isinstance(xc, (list, np.ndarray)): raise TypeError("xc must be a scalar") if isinstance(yc, (list, np.ndarray)): raise TypeError("yc must be a scalar") # Find the fiber bundle closest to (xc, yc): bundles = self.grow_axon_bundles() bundle = self.find_closest_axon(bundles, xret=xc, yret=yc)[0] # For that bundle, find the bundle segment closest to (xc, yc): idx = np.argmin((bundle[:, 0] - xc) ** 2 + (bundle[:, 1] - yc) ** 2) # Calculate orientation from atan2(dy, dx): if idx == 0: # Bundle index 0: there's no index -1 dx = bundle[1, :] - bundle[0, :] elif idx == bundle.shape[0] - 1: # Bundle index -1: there's no index len(bundle) dx = bundle[-1, :] - bundle[-2, :] else: # Else: Look at previous and subsequent segments: dx = (bundle[idx + 1, :] - bundle[idx - 1, :]) / 2 dx[1] *= -1 tangent = np.arctan2(*dx[::-1]) # Confine to (-pi/2, pi/2): if tangent < np.deg2rad(-90): tangent += np.deg2rad(180) if tangent > np.deg2rad(90): tangent -= np.deg2rad(180) return tangent