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 vote down vote up
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 vote down vote up
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 vote down vote up
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 vote down vote up
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 vote down vote up
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 vote down vote up
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 vote down vote up
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 vote down vote up
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 vote down vote up
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 vote down vote up
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 vote down vote up
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 vote down vote up
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 vote down vote up
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 vote down vote up
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 vote down vote up
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 vote down vote up
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 vote down vote up
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 vote down vote up
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 vote down vote up
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 vote down vote up
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 vote down vote up
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 vote down vote up
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 vote down vote up
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 vote down vote up
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 vote down vote up
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 vote down vote up
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 vote down vote up
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 vote down vote up
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 vote down vote up
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 vote down vote up
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