Python numpy.cross() Examples

The following are 30 code examples of numpy.cross(). 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: 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 #2
Source File: dicomwrappers.py    From me-ica with GNU Lesser General Public License v2.1 6 votes vote down vote up
def slice_normal(self):
        #The std_slice_normal comes from the cross product of the directions 
        #in the ImageOrientationPatient
        std_slice_normal = super(SiemensWrapper, self).slice_normal
        csa_slice_normal = csar.get_slice_normal(self.csa_header)
        if std_slice_normal is None and csa_slice_normal is None:
            return None
        elif std_slice_normal is None:
            return np.array(csa_slice_normal)
        elif csa_slice_normal is None:
            return std_slice_normal
        else:
            #Make sure the two normals are very close to parallel unit vectors
            dot_prod = np.dot(csa_slice_normal, std_slice_normal)
            assert np.allclose(np.fabs(dot_prod), 1.0, atol=1e-5)
            #Use the slice normal computed with the cross product as it will 
            #always be the most orthogonal, but take the sign from the CSA
            #slice normal
            if dot_prod < 0:
                return -std_slice_normal
            else:
                return std_slice_normal 
Example #3
Source File: retinotopy.py    From neuropythy with GNU Affero General Public License v3.0 6 votes vote down vote up
def _retinotopic_field_sign_triangles(m, retinotopy):
    t = m.tess if isinstance(m, geo.Mesh) or isinstance(m, geo.Topology) else m
    # get the polar angle and eccen data as a complex number in degrees
    if pimms.is_str(retinotopy):
        (x,y) = as_retinotopy(retinotopy_data(m, retinotopy), 'geographical')
    elif retinotopy is Ellipsis:
        (x,y) = as_retinotopy(retinotopy_data(m, 'any'),      'geographical')
    else:
        (x,y) = as_retinotopy(retinotopy,                     'geographical')
    # Okay, now we want to make some coordinates...
    coords = np.asarray([x, y])
    us = coords[:, t.indexed_faces[1]] - coords[:, t.indexed_faces[0]]
    vs = coords[:, t.indexed_faces[2]] - coords[:, t.indexed_faces[0]]
    (us,vs) = [np.concatenate((xs, np.full((1, t.face_count), 0.0))) for xs in [us,vs]]
    xs = np.cross(us, vs, axis=0)[2]
    xs[np.isclose(xs, 0)] = 0
    return np.sign(xs) 
Example #4
Source File: PolyMesh.py    From laplacian-meshes with GNU General Public License v3.0 6 votes vote down vote up
def updateNormalBuffer(self):
        #First compute cross products of all face triangles
        V1 = self.VPos[self.ITris[:, 1], :] - self.VPos[self.ITris[:, 0], :]
        V2 = self.VPos[self.ITris[:, 2], :] - self.VPos[self.ITris[:, 0], :]
        FNormals = np.cross(V1, V2)
        FAreas = np.reshape(np.sqrt(np.sum(FNormals**2, 1)), (FNormals.shape[0], 1))
        FAreas[FAreas == 0] = 1
        FNormals = FNormals/FAreas
        self.FNormals = FNormals
        self.FCentroid = 0*FNormals
        self.VNormals = 0*self.VPos
        VAreas = np.zeros((self.VPos.shape[0], 1))
        for k in range(3):
            self.VNormals[self.ITris[:, k], :] += FAreas*FNormals
            VAreas[self.ITris[:, k]] += FAreas
            self.FCentroid += self.VPos[self.ITris[:, k], :]
        self.FCentroid /= 3
        #Normalize normals
        VAreas[VAreas == 0] = 1
        self.VNormals = self.VNormals / VAreas
    
    #buffersOnly: True if there is no mesh structure other than VPos
    #VColors and ITris 
Example #5
Source File: csys.py    From StructEngPy with MIT License 6 votes vote down vote up
def __init__(self,origin, pt1, pt2, name=None):
        """
        origin: 3x1 vector
        pt1: 3x1 vector
        pt2: 3x1 vector
        """
        self.__origin=origin    
        vec1 = np.array([pt1[0] - origin[0] , pt1[1] - origin[1] , pt1[2] - origin[2]])
        vec2 = np.array([pt2[0] - origin[0] , pt2[1] - origin[1] , pt2[2] - origin[2]])
        cos = np.dot(vec1, vec2)/np.linalg.norm(vec1)/np.linalg.norm(vec2)
        if  cos == 1 or cos == -1:
            raise Exception("Three points should not in a line!!")        
        self.__x = vec1/np.linalg.norm(vec1)
        z = np.cross(vec1, vec2)
        self.__z = z/np.linalg.norm(z)
        self.__y = np.cross(self.z, self.x)
        self.__name=uuid.uuid1() if name==None else name 
Example #6
Source File: rotation_utils.py    From DOTA_models with Apache License 2.0 6 votes vote down vote up
def rotate_camera_to_point_at(up_from, lookat_from, up_to, lookat_to):
  inputs = [up_from, lookat_from, up_to, lookat_to]
  for i in range(4):
    inputs[i] = normalize(np.array(inputs[i]).reshape((-1,)))
  up_from, lookat_from, up_to, lookat_to = inputs
  r1 = r_between(lookat_from, lookat_to)

  new_x = np.dot(r1, np.array([1, 0, 0]).reshape((-1, 1))).reshape((-1))
  to_x = normalize(np.cross(lookat_to, up_to))
  angle = np.arccos(np.dot(new_x, to_x))
  if angle > ANGLE_EPS:
    if angle < np.pi - ANGLE_EPS:
      ax = normalize(np.cross(new_x, to_x))
      flip = np.dot(lookat_to, ax)
      if flip > 0:
        r2 = get_r_matrix(lookat_to, angle)
      elif flip < 0:
        r2 = get_r_matrix(lookat_to, -1. * angle)
    else:
      # Angle of rotation is too close to 180 degrees, direction of rotation
      # does not matter.
      r2 = get_r_matrix(lookat_to, angle)
  else:
    r2 = np.eye(3)
  return np.dot(r2, r1) 
Example #7
Source File: densepose_methods.py    From Parsing-R-CNN with MIT License 6 votes vote down vote up
def barycentric_coordinates_exists(self, P0, P1, P2, P):
        u = P1 - P0
        v = P2 - P0
        w = P - P0

        vCrossW = np.cross(v, w)
        vCrossU = np.cross(v, u)
        if np.dot(vCrossW, vCrossU) < 0:
            return False

        uCrossW = np.cross(u, w)
        uCrossV = np.cross(u, v)

        if np.dot(uCrossW, uCrossV) < 0:
            return False

        denom = np.sqrt((uCrossV ** 2).sum())
        r = np.sqrt((vCrossW ** 2).sum()) / denom
        t = np.sqrt((uCrossW ** 2).sum()) / denom

        return (r <= 1) & (t <= 1) & (r + t <= 1) 
Example #8
Source File: common.py    From PyOptiX with MIT License 6 votes vote down vote up
def calculate_camera_variables(eye, lookat, up, fov, aspect_ratio, fov_is_vertical=False):
    import numpy as np
    import math

    W = np.array(lookat) - np.array(eye)
    wlen = np.linalg.norm(W)
    U = np.cross(W, np.array(up))
    U /= np.linalg.norm(U)
    V = np.cross(U, W)
    V /= np.linalg.norm(V)

    if fov_is_vertical:
        vlen = wlen * math.tan(0.5 * fov * math.pi / 180.0)
        V *= vlen
        ulen = vlen * aspect_ratio
        U *= ulen
    else:
        ulen = wlen * math.tan(0.5 * fov * math.pi / 180.0)
        U *= ulen
        vlen = ulen * aspect_ratio
        V *= vlen

    return U, V, W 
Example #9
Source File: motion.py    From 2D-Motion-Retargeting with MIT License 6 votes vote down vote up
def get_local3d(motion3d, angles=None):
    """
    Get the unit vectors for local rectangular coordinates for given 3D motion
    :param motion3d: numpy array. 3D motion from 3D joints positions, shape (nr_joints, 3, nr_frames).
    :param angles: tuple of length 3. Rotation angles around each axis.
    :return: numpy array. unit vectors for local rectangular coordinates's , shape (3, 3).
    """
    # 2 RightArm 5 LeftArm 9 RightUpLeg 12 LeftUpLeg
    horizontal = (motion3d[2] - motion3d[5] + motion3d[9] - motion3d[12]) / 2
    horizontal = np.mean(horizontal, axis=1)
    horizontal = horizontal / np.linalg.norm(horizontal)
    local_z = np.array([0, 0, 1])
    local_y = np.cross(horizontal, local_z)  # bugs!!!, horizontal and local_Z may not be perpendicular
    local_y = local_y / np.linalg.norm(local_y)
    local_x = np.cross(local_y, local_z)
    local = np.stack([local_x, local_y, local_z], axis=0)

    if angles is not None:
        local = rotate_coordinates(local, angles)

    return local 
Example #10
Source File: base.py    From pymesh with MIT License 6 votes vote down vote up
def update_normals(self):
        v0 = self.vectors[:, 0, :3]
        v1 = self.vectors[:, 1, :3]
        v2 = self.vectors[:, 2, :3]
        _normals = numpy.cross(v1 - v0, v2 - v0)

        for i in range(len(_normals)):
            norm = numpy.linalg.norm(_normals[i])
            if norm != 0:
                _normals[i] /= numpy.linalg.norm(_normals[i])

        self.normals[:] = _normals
        return self

    #####################################################################
    # Analyze functions
    # 
Example #11
Source File: transformations.py    From rai-python with MIT License 6 votes vote down vote up
def shear_matrix(angle, direction, point, normal):
    """Return matrix to shear by angle along direction vector on shear plane.
    The shear plane is defined by a point and normal vector. The direction
    vector must be orthogonal to the plane's normal vector.
    A point P is transformed by the shear matrix into P" such that
    the vector P-P" is parallel to the direction vector and its extent is
    given by the angle of P-P'-P", where P' is the orthogonal projection
    of P onto the shear plane.
    >>> angle = (random.random() - 0.5) * 4*math.pi
    >>> direct = numpy.random.random(3) - 0.5
    >>> point = numpy.random.random(3) - 0.5
    >>> normal = numpy.cross(direct, numpy.random.random(3))
    >>> S = shear_matrix(angle, direct, point, normal)
    >>> numpy.allclose(1.0, numpy.linalg.det(S))
    True
    """
    normal = unit_vector(normal[:3])
    direction = unit_vector(direction[:3])
    if abs(numpy.dot(normal, direction)) > 1e-6:
        raise ValueError("direction and normal vectors are not orthogonal")
    angle = math.tan(angle)
    M = numpy.identity(4)
    M[:3, :3] += angle * numpy.outer(direction, normal)
    M[:3, 3] = -angle * numpy.dot(point[:3], normal) * direction
    return M 
Example #12
Source File: utils_perspective.py    From DIB-R with MIT License 6 votes vote down vote up
def camera_info(param):
    theta = np.deg2rad(param[0])
    phi = np.deg2rad(param[1])

    camY = param[3] * np.sin(phi)
    temp = param[3] * np.cos(phi)
    camX = temp * np.cos(theta)    
    camZ = temp * np.sin(theta)        
    cam_pos = np.array([camX, camY, camZ])        

    axisZ = cam_pos.copy()
    axisY = np.array([0, 1, 0], dtype=np.float32)
    axisX = np.cross(axisY, axisZ)
    axisY = np.cross(axisZ, axisX)
    
    # cam_mat = np.array([axisX, axisY, axisZ])
    cam_mat = np.array([unit(axisX), unit(axisY), unit(axisZ)])
    
    return cam_mat, cam_pos


##################################################### 
Example #13
Source File: geometry.py    From connecting_the_dots with MIT License 6 votes vote down vote up
def mesh_adjust_winding_order(verts, faces, normals):
  n0 = normals[faces[:,0]]
  n1 = normals[faces[:,1]]
  n2 = normals[faces[:,2]]
  fnormals = (n0 + n1 + n2) / 3

  v0 = verts[faces[:,0]]
  v1 = verts[faces[:,1]]
  v2 = verts[faces[:,2]]

  e0 = v1 - v0
  e1 = v2 - v0
  fn = np.cross(e0, e1)

  dot = np.sum(fnormals * fn, axis=1)
  ma = dot < 0

  nfaces = faces.copy()
  nfaces[ma,1], nfaces[ma,2] = nfaces[ma,2], nfaces[ma,1]

  return nfaces 
Example #14
Source File: test_fem.py    From simnibs with GNU General Public License v3.0 6 votes vote down vote up
def tms_sphere(sphere3_msh):
    m = sphere3_msh.crop_mesh(elm_type=4)
    dipole_pos = np.array([0., 0., 300])
    dipole_moment = np.array([1., 0., 0.])
    didt = 1e6
    r = (m.nodes.node_coord - dipole_pos) * 1e-3
    dAdt = 1e-7 * didt * np.cross(dipole_moment, r) / (np.linalg.norm(r, axis=1)[:, None] ** 3)
    dAdt = mesh_io.NodeData(dAdt, mesh=m)
    dAdt.field_name = 'dAdt'
    dAdt.mesh = m
    pos = m.elements_baricenters().value
    E_analytical = analytical_solutions.tms_E_field(dipole_pos * 1e-3,
                                                    dipole_moment, didt,
                                                    pos * 1e-3)
    cond = mesh_io.ElementData(np.ones(m.elm.nr))
    cond.mesh = m
    return m, cond, dAdt, E_analytical 
Example #15
Source File: pywannier90.py    From pyscf with Apache License 2.0 6 votes vote down vote up
def transform(x_vec, z_vec):
    '''
    Construct a transformation matrix to transform r_vec to the new coordinate system defined by x_vec and z_vec
    '''

    x_vec = x_vec/np.linalg.norm(np.asarray(x_vec))
    z_vec = z_vec/np.linalg.norm(np.asarray(z_vec))
    assert x_vec.dot(z_vec) == 0
    y_vec = np.cross(x_vec,z_vec)
    new = np.asarray([x_vec, y_vec, z_vec])
    original = np.asarray([[1,0,0],[0,1,0],[0,0,1]])

    tran_matrix = np.empty([3,3])
    for row in range(3):
        for col in range(3):
            tran_matrix[row,col] = np.cos(angle(original[row],new[col]))

    return tran_matrix.T 
Example #16
Source File: coil_numpy.py    From simnibs with GNU General Public License v3.0 6 votes vote down vote up
def _calculate_dadt_ccd(msh, ccd_file, coil_matrix, didt, geo_fn):
    """ auxiliary function to calculate the dA/dt field from a ccd file """
    # read ccd file
    d_position, d_moment = read_ccd(ccd_file)
    # transfrom positions to mm
    d_position *= 1e3
    # add a column to the position in order to apply the transformation matrix
    d_position = np.hstack([d_position, np.ones((d_position.shape[0], 1))])
    d_position = coil_matrix.dot(d_position.T).T[:, :3]
    # rotate the moment
    d_moment = coil_matrix[:3, :3].dot(d_moment.T).T
    A = np.zeros((msh.nodes.nr, 3), dtype=float)
    for p, m in zip(d_position, d_moment):
        # get distance of point to dipole, transform back to meters
        r = (msh.nodes.node_coord - p) * 1e-3
        A += 1e-7 * didt * np.cross(m, r) / (np.linalg.norm(r, axis=1)[:, None] ** 3)
    node_data = mesh_io.NodeData(A)

    if geo_fn is not None:
        mesh_io.write_geo_spheres(d_position, geo_fn,
                               np.linalg.norm(d_moment, axis=1),
                               'coil_dipoles')

    return node_data 
Example #17
Source File: csys.py    From StructEngPy with MIT License 6 votes vote down vote up
def set_by_3pts(self,origin, pt1, pt2):
        """
        origin: tuple 3
        pt1: tuple 3
        pt2: tuple 3
        """
        self.origin=origin    
        vec1 = np.array([pt1[0] - origin[0] , pt1[1] - origin[1] , pt1[2] - origin[2]])
        vec2 = np.array([pt2[0] - origin[0] , pt2[1] - origin[1] , pt2[2] - origin[2]])
        cos = np.dot(vec1, vec2)/np.linalg.norm(vec1)/np.linalg.norm(vec2)
        if  cos == 1 or cos == -1:
            raise Exception("Three points should not in a line!!")        
        self.x = vec1/np.linalg.norm(vec1)
        z = np.cross(vec1, vec2)
        self.z = z/np.linalg.norm(z)
        self.y = np.cross(self.z, self.x) 
Example #18
Source File: mesh_affine_equ.py    From pyscf with Apache License 2.0 6 votes vote down vote up
def __init__(self, **kw):
    """  
      Constructor of affine, equidistant 3d mesh class
      ucell : unit cell vectors (in coordinate space)
      Ecut  : Energy cutoff to parametrize the discretization 
    """
    from scipy.fftpack import next_fast_len
    
    self.ucell = kw['ucell'] if 'ucell' in kw else 30.0*np.eye(3) # Not even unit cells vectors are required by default
    self.Ecut = Ecut = kw['Ecut'] if 'Ecut' in kw else 50.0 # 50.0 Hartree by default
    luc = np.sqrt(np.einsum('ix,ix->i', self.ucell, self.ucell))
    self.shape = nn = np.array([next_fast_len( int(np.rint(l * np.sqrt(Ecut)/2))) for l in luc], dtype=int)
    self.size  = np.prod(self.shape)
    gc = self.ucell/(nn) # This is probable the best for finite systems, for PBC use nn, not (nn-1)
    self.dv = np.abs(np.dot(gc[0], np.cross(gc[1], gc[2] )))
    rr = [np.array([gc[i]*j for j in range(nn[i])]) for i in range(3)]
    self.rr = rr
    self.origin = kw['origin'] if 'origin' in kw else np.zeros(3) 
Example #19
Source File: coordinates.py    From ratcave with MIT License 5 votes vote down vote up
def rotation_matrix_between_vectors(from_vec, to_vec):
    """
    Returns a rotation matrix to rotate from 3d vector "from_vec" to 3d vector "to_vec".
    Equation from https://math.stackexchange.com/questions/180418/calculate-rotation-matrix-to-align-vector-a-to-vector-b-in-3d
    """
    a, b = (vec/np.linalg.norm(vec) for vec in (from_vec, to_vec))

    v = np.cross(a, b)
    cos = np.dot(a, b)
    if cos == -1.:
        raise ValueError("Orientation in complete opposite direction")
    v_cpm = cross_product_matrix(v)
    rot_mat = np.identity(3) + v_cpm + np.dot(v_cpm, v_cpm) * (1. / 1. + cos)
    return rot_mat 
Example #20
Source File: coordinates.py    From ratcave with MIT License 5 votes vote down vote up
def cross_product_matrix(vec):
    """Returns a 3x3 cross-product matrix from a 3-element vector."""
    return np.array([[0, -vec[2], vec[1]],
                     [vec[2], 0, -vec[0]],
                     [-vec[1], vec[0], 0]]) 
Example #21
Source File: commons.py    From connecting_the_dots with MIT License 5 votes vote down vote up
def get_rotation_matrix(v0, v1):
  v0 = v0/np.linalg.norm(v0)
  v1 = v1/np.linalg.norm(v1)
  v = np.cross(v0,v1)
  c = np.dot(v0,v1)
  s = np.linalg.norm(v)
  I = np.eye(3)
  vXStr = '{} {} {}; {} {} {}; {} {} {}'.format(0, -v[2], v[1], v[2], 0, -v[0], -v[1], v[0], 0)
  k = np.matrix(vXStr)
  r = I + k + k @ k * ((1 -c)/(s**2))
  return np.asarray(r.astype(np.float32)) 
Example #22
Source File: test_sphere.py    From simnibs with GNU General Public License v3.0 5 votes vote down vote up
def test_radial_source(self, quadrant_points):
        """ The field of a radial source is given by a 3.10 (Heller and Hulsteyn) """
        dipole_pos = np.array([3, 0, 0])
        dipole_moment = np.array([1, 0, 0])
        E = sphere.tms_E_field(dipole_pos, dipole_moment, 1, quadrant_points)
        E_simple = -1e-7 * np.cross(dipole_pos - quadrant_points, dipole_moment) \
                / np.linalg.norm(quadrant_points - dipole_pos, axis=1)[:, None]**3  
        # Explanation for the minus: see appendix 1 in Heller and Hulsteyn
        assert np.allclose(E, E_simple) 
Example #23
Source File: mesh_io.py    From simnibs with GNU General Public License v3.0 5 votes vote down vote up
def nodes_normals(self, smooth=0):
        ''' Normals for all nodes in a surface

        Parameters
        ------------
        smooth: int (optional)
            Number of smoothing cycles to perform. Default: 0

        Returns
        ---------
        nd: NodeData
            NodeData structure with normals for each surface node

        '''
        nodes = np.unique(self.elm[self.elm.triangles, :3])
        elements = self.elm.triangles

        nd = np.zeros((self.nodes.nr, 3))

        node_tr = self.nodes[self.elm.node_number_list[elements - 1, :3]]
        sideA = node_tr[:, 1] - node_tr[:, 0]

        sideB = node_tr[:, 2] - node_tr[:, 0]
        normals = np.cross(sideA, sideB)

        triangle_nodes = self.elm.node_number_list[elements - 1, :3] - 1
        for s in range(smooth + 1):
            for i in range(3):
                nd[:, i] = \
                    np.bincount(triangle_nodes.reshape(-1),
                                np.repeat(normals[:, i], 3),
                                self.nodes.nr)

            normals = np.sum(nd[self.elm.node_number_list[elements - 1, :3] - 1], axis=1)
            normals /= np.linalg.norm(normals, axis=1)[:, None]

        nd[nodes - 1] = nd[nodes-1] / \
            np.linalg.norm(nd[nodes-1], axis=1)[:, None]
 
        return NodeData(nd, 'normals') 
Example #24
Source File: mesh_io.py    From simnibs with GNU General Public License v3.0 5 votes vote down vote up
def triangle_normals(self, smooth=False):
        """ Calculates the normals of triangles

        Parameters
        ------------
        smooth (optional): int
            Number of smoothing steps to perform. Default: 0
        Returns
        --------
        normals: ElementData
            normals of triangles, zero at the tetrahedra

        """
        normals = ElementData(
            np.zeros((self.elm.nr, 3), dtype=float),
            'normals'
        )
        tr_indexes = self.elm.triangles
        node_tr = self.nodes[self.elm[tr_indexes, :3]]
        sideA = node_tr[:, 1] - node_tr[:, 0]
        sideB = node_tr[:, 2] - node_tr[:, 0]
        n = np.cross(sideA, sideB)
        normals[tr_indexes] = n / np.linalg.norm(n, axis=1)[:, None]
        if smooth == 0:
            node_tr = self.nodes[self.elm[tr_indexes, :3]]
            sideA = node_tr[:, 1] - node_tr[:, 0]

            sideB = node_tr[:, 2] - node_tr[:, 0]
            n = np.cross(sideA, sideB)
            normals[tr_indexes] = n / np.linalg.norm(n, axis=1)[:, None]
        elif smooth > 0:
            normals_nodes = self.nodes_normals(smooth)
            tr = self.elm[tr_indexes, :3]
            n = np.mean(normals_nodes[tr], axis=1)
            normals[tr_indexes] = n / np.linalg.norm(n, axis=1)[:, None]
        else:
            raise ValueError('smooth parameter must be >= 0')
        return normals 
Example #25
Source File: mesh_io.py    From simnibs with GNU General Public License v3.0 5 votes vote down vote up
def elements_volumes_and_areas(self):
        """ Calculates the volumes of tetrahedra and areas of triangles

        Returns
        ----------
        v: simnibs.Msh.ElementData
            Volume/areas of tetrahedra/triangles

        Note
        ------
            In the mesh's unit (normally mm)
        """
        vol = ElementData(np.zeros(self.elm.nr, dtype=float),
                          'volumes_and_areas')
        th_indexes = self.elm.tetrahedra
        tr_indexes = self.elm.triangles

        node_tr = self.nodes[self.elm[tr_indexes, :3]]
        sideA = node_tr[:, 1] - node_tr[:, 0]

        sideB = node_tr[:, 2] - node_tr[:, 0]

        n = np.cross(sideA, sideB)

        vol[tr_indexes] = np.linalg.norm(n, axis=1) * 0.5

        node_th = self.nodes[self.elm[th_indexes]]
        M = node_th[:, 1:] - node_th[:, 0, None]
        vol[th_indexes] = np.abs(np.linalg.det(M)) / 6.

        return vol 
Example #26
Source File: surface.py    From simnibs with GNU General Public License v3.0 5 votes vote down vote up
def calculateMatSimnibs(self, p1, p2, skin_distance=4):
        p1_np = np.array(p1, 'float')
        p2_np = np.array(p2, 'float')
        center, z_axis = self.projectPoint(p1_np, smooth=True)
        z_axis = -1 * z_axis
        center = center - skin_distance * z_axis

        # y direction: user input, orthogonalized
        y_axis = p2_np - p1_np
        y_axis = y_axis - y_axis.dot(z_axis) * z_axis
        if np.linalg.norm(y_axis) > 1e-6:
            y_axis /= np.linalg.norm(y_axis)
        else:
            y_axis = np.random.rand(3)
            y_axis = y_axis - y_axis.dot(z_axis) * z_axis
            y_axis /= np.linalg.norm(y_axis)


        x_axis = np.cross(y_axis, z_axis)

        matsimnibs = [[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]]

        (matsimnibs[0][0], matsimnibs[1][0],
         matsimnibs[2][0]) = x_axis.tolist()
        (matsimnibs[0][1], matsimnibs[1][1],
         matsimnibs[2][1]) = y_axis.tolist()
        (matsimnibs[0][2], matsimnibs[1][2],
         matsimnibs[2][2]) = z_axis.tolist()
        (matsimnibs[0][3], matsimnibs[1][3],
         matsimnibs[2][3]) = center.tolist()
        matsimnibs[3][3] = 1

        return matsimnibs 
Example #27
Source File: sim_struct.py    From simnibs with GNU General Public License v3.0 5 votes vote down vote up
def add_positions_from_csv(self, fn_csv):
        ''' Reads a csv file and adds the positions defined to the tmslist

        Parameters
        ----------
        fn_csv: str
            CSV file with the fields
                Type, pos_x, pos_y, pos_z, ez_x, ez_y, ez_z, ey_x, ey_y, ey_z, dist, name, ...
                "Type" needs to be CoilPos. The positions are in subject space. The
                transformations module can transfrom from MNI to subject space
        '''
        type_, coordinates, extra, name, _, _ = transformations._read_csv(fn_csv)
        for t, c, e, n in zip(type_, coordinates, extra, name):
            if t == 'CoilPos':
                p = POSITION()
                vz = e[:3]
                vy = e[3:6]
                vx = np.cross(vy, vz)
                mat = np.eye(4)
                mat[:3, 0] = vx
                mat[:3, 1] = vy
                mat[:3, 2] = vz
                mat[:3, 3] = c
                p.matsimnibs = mat.tolist()
                p.name = n
                self.pos.append(p) 
Example #28
Source File: Cameras3D.py    From laplacian-meshes with GNU General Public License v3.0 5 votes vote down vote up
def gotoCameraFrame(self):
        gotoCameraFrame(self.towards, self.up, np.cross(self.towards, self.up), self.eye) 
Example #29
Source File: Cameras3D.py    From laplacian-meshes with GNU General Public License v3.0 5 votes vote down vote up
def translate(self, dx, dy):
        length = self.center - self.eye
        length = np.sqrt(np.sum(length**2))*math.tan(self.yfov)
        dx = length*dx / float(self.pixWidth)
        dy = length*dy / float(self.pixHeight)
        r = np.cross(self.towards, self.up)
        self.center = self.center - dx*r - dy*self.up
        self.eye = self.eye - dx*r - dy*self.up
        self.updateVecsFromPolar() 
Example #30
Source File: dicomwrappers.py    From me-ica with GNU Lesser General Public License v2.1 5 votes vote down vote up
def slice_normal(self):
        iop = self.image_orient_patient
        if iop is None:
            return None
        return np.cross(*iop.T[:])