Determine scale factor difference between edge triangulation and world
163553
static double determineScale(View base , Motion edge )
		throws Exception
	{
		View viewA = edge.viewSrc;
		View viewB = edge.viewDst;

		boolean baseIsA = base == viewA;

		// determine the scale factor difference
		Point3D_F64 worldInBase3D = new Point3D_F64();
		Point3D_F64 localInBase3D = new Point3D_F64();
		GrowQueue_F64 scales = new GrowQueue_F64();
		for (int i = 0; i < edge.stereoTriangulations.size(); i++) {
			// get the feature as triangulated in this edge.
			Feature3D edge3D = edge.stereoTriangulations.get(i);

			int indexSrc = edge3D.obsIdx.get(0);
			int indexDst = edge3D.obsIdx.get(1);

			Feature3D world3D = baseIsA ? viewA.features3D[indexSrc] : viewB.features3D[indexDst];
			if( world3D == null )
				continue;
			// Find the world point in the local coordinate system
			SePointOps_F64.transformReverse(base.viewToWorld,world3D.worldPt,worldInBase3D);

			// put this point into the base frame
			if( !baseIsA ) {
				SePointOps_F64.transform(edge.a_to_b,edge3D.worldPt,localInBase3D);
			} else {
				localInBase3D.set(edge3D.worldPt);
			}

			scales.add(worldInBase3D.z / localInBase3D.z);
		}

		if( scales.size < 20 ) {
			throw new Exception("Not enough matches with known points");
		}
		// Get the scale offset as the median value to make it robust to noise
		scales.sort();
		return scales.getFraction(0.5);
	}
421281
private float calculateScale(int x, int y) {
		final float centerX = getWidth() / 2f;
		final float centerY = getHeight() / 2f;
		final float maxDistance = (float) Math.sqrt(centerX * centerX + centerY * centerY);

		final float deltaX = centerX - x;
		final float deltaY = centerY - y;
		final float distance = (float) Math.sqrt(deltaX * deltaX + deltaY * deltaY);
		final float scale = 0.5f + (distance / maxDistance) * 0.5f;
		return scale;
	}
743708
def get_scaled_cutout_basic_view(shp, p1, p2, scales):
    

    x1, y1 = p1[:2]
    x2, y2 = p2[:2]
    scale_x, scale_y = scales[:2]
    # calculate dimensions of NON-scaled cutout
    old_wd = x2 - x1 + 1
    old_ht = y2 - y1 + 1
    # calculate dimensions of scaled cutout
    new_wd = int(round(scale_x * old_wd))
    new_ht = int(round(scale_y * old_ht))

    if len(scales) == 2:
        return get_scaled_cutout_wdht_view(shp, x1, y1, x2, y2, new_wd, new_ht)

    z1, z2, scale_z = p1[2], p2[2], scales[2]
    old_dp = z2 - z1 + 1
    new_dp = int(round(scale_z * old_dp))
    return get_scaled_cutout_wdhtdp_view(shp, p1, p2, (new_wd, new_ht, new_dp))
261028
def _adjacent_tri(self, edge, i):
        
        if not np.isscalar(i):
            i = [x for x in i if x not in edge][0]

        try:
            pt1 = self._edges_lookup[edge]
            pt2 = self._edges_lookup[(edge[1], edge[0])]
        except KeyError:
            return None
            
        if pt1 == i:
            return (edge[1], edge[0], pt2)
        elif pt2 == i:
            return (edge[1], edge[0], pt1)
        else:
            raise RuntimeError("Edge %s and point %d do not form a triangle "
                               "in this mesh." % (edge, i))
603292
def edge_refine_triangulation_by_vertices(self, vertices):
        

        triangles = self.identify_vertex_triangles(vertices)

        return self.edge_refine_triangulation_by_triangles(triangles)
507618
def distance_to_edge(geom, skydir):
    """"""

    # FIXME: We should add a pixel_size property in gammapy.maps
    # FIXME: We should make this into a MapGeom method 
    
    xpix, ypix = skydir.to_pixel(geom.wcs, origin=0)
    deltax = np.array((xpix - geom.center_pix[0]) * geom._cdelt[0],
                      ndmin=1)
    deltay = np.array((ypix - geom.center_pix[1]) * geom._cdelt[1],
                      ndmin=1)

    deltax = np.abs(deltax) - 0.5 * geom.width[0]
    deltay = np.abs(deltay) - 0.5 * geom.width[1]

    m0 = (deltax < 0) & (deltay < 0)
    m1 = (deltax > 0) & (deltay < 0)
    m2 = (deltax < 0) & (deltay > 0)
    m3 = (deltax > 0) & (deltay > 0)
    mx = np.abs(deltax) <= np.abs(deltay)
    my = np.abs(deltay) < np.abs(deltax)

    delta = np.zeros(len(deltax))
    delta[(m0 & mx) | (m3 & my) | m1] = deltax[(m0 & mx) | (m3 & my) | m1]
    delta[(m0 & my) | (m3 & mx) | m2] = deltay[(m0 & my) | (m3 & mx) | m2]
    return delta
706626
def wt(u, v, dfs_data):
    """"""
    # Determine the edge_id
    edge_id = dfs_data['graph'].get_first_edge_id_by_node_ids(u, v)
    # Pull the weight of that edge
    return dfs_data['edge_weights'][edge_id]
1122478
def __msgc_step3_discontinuity_localization(self):
        
        import scipy

        start = self._start_time
        seg = 1 - self.segmentation.astype(np.int8)
        self.stats["low level object voxels"] = np.sum(seg)
        self.stats["low level image voxels"] = np.prod(seg.shape)
        # in seg is now stored low resolution segmentation
        # back to normal parameters
        # step 2: discontinuity localization
        # self.segparams = sparams_hi
        seg_border = scipy.ndimage.filters.laplace(seg, mode="constant")
        logger.debug("seg_border: %s", scipy.stats.describe(seg_border, axis=None))
        # logger.debug(str(np.max(seg_border)))
        # logger.debug(str(np.min(seg_border)))
        seg_border[seg_border != 0] = 1
        logger.debug("seg_border: %s", scipy.stats.describe(seg_border, axis=None))
        # scipy.ndimage.morphology.distance_transform_edt
        boundary_dilatation_distance = self.segparams["boundary_dilatation_distance"]
        seg = scipy.ndimage.morphology.binary_dilation(
            seg_border,
            # seg,
            np.ones(
                [
                    (boundary_dilatation_distance * 2) + 1,
                    (boundary_dilatation_distance * 2) + 1,
                    (boundary_dilatation_distance * 2) + 1,
                ]
            ),
        )
        if self.keep_temp_properties:
            self.temp_msgc_lowres_discontinuity = seg
        else:
            self.temp_msgc_lowres_discontinuity = None

        if self.debug_images:
            import sed3

            pd = sed3.sed3(seg_border)  # ), contour=seg)
            pd.show()
            pd = sed3.sed3(seg)  # ), contour=seg)
            pd.show()
        # segzoom = scipy.ndimage.interpolation.zoom(seg.astype('float'), zoom,
        #                                                order=0).astype('int8')
        self.stats["t3"] = time.time() - start
        return seg
89758
def triangulate(self):
        
        self._initialize()
        
        pts = self.pts
        front = self._front
        
        ## Begin sweep (sec. 3.4)
        for i in range(3, pts.shape[0]):
            pi = pts[i]
            #debug("========== New point %d: %s ==========" % (i, pi))
            
            # First, triangulate from front to new point
            # This applies to both "point events" (3.4.1) 
            # and "edge events" (3.4.2).

            # get index along front that intersects pts[i]
            l = 0
            while pts[front[l+1], 0] <= pi[0]:
                l += 1
            pl = pts[front[l]]
            
            # "(i) middle case"
            if pi[0] > pl[0]:  
                #debug("  mid case")
                # Add a single triangle connecting pi,pl,pr
                self._add_tri(front[l], front[l+1], i)
                front.insert(l+1, i)
            # "(ii) left case"
            else:
                #debug("  left case")
                # Add triangles connecting pi,pl,ps and pi,pl,pr
                self._add_tri(front[l], front[l+1], i)
                self._add_tri(front[l-1], front[l], i)
                front[l] = i
            
            #debug(front)
                
            # Continue adding triangles to smooth out front
            # (heuristics shown in figs. 9, 10)
            #debug("Smoothing front...")
            for direction in -1, 1:
                while True:
                    # Find point connected to pi
                    ind0 = front.index(i)
                    ind1 = ind0 + direction
                    ind2 = ind1 + direction
                    if ind2 < 0 or ind2 >= len(front):
                        break
                    
                    # measure angle made with front
                    p1 = pts[front[ind1]]
                    p2 = pts[front[ind2]]
                    err = np.geterr()
                    np.seterr(invalid='ignore')
                    try:
                        angle = np.arccos(self._cosine(pi, p1, p2))
                    finally:
                        np.seterr(**err)
                    
                    # if angle is < pi/2, make new triangle
                    #debug("Smooth angle:", pi, p1, p2, angle)
                    if angle > np.pi/2. or np.isnan(angle):
                        break
                    
                    assert (i != front[ind1] and 
                            front[ind1] != front[ind2] and 
                            front[ind2] != i)
                    self._add_tri(i, front[ind1], front[ind2],
                                  source='smooth1')
                    front.pop(ind1)
            #debug("Finished smoothing front.")
            
            # "edge event" (sec. 3.4.2)
            # remove any triangles cut by completed edges and re-fill 
            # the holes.
            if i in self._tops:
                for j in self._bottoms[self._tops == i]:
                    # Make sure edge (j, i) is present in mesh
                    # because edge event may have created a new front list
                    self._edge_event(i, j)  
                    front = self._front 
                
        self._finalize()
        
        self.tris = np.array(list(self.tris.keys()), dtype=int)
692288
function getGeoBounds({face, ij, level}) {
  const result = [];
  const offsets = [[0, 0], [0, 1], [1, 1], [1, 0], [0, 0]];

  // The S2 cell edge is curved: http://s2geometry.io/
  // This is more prominent at lower levels
  // resolution is the number of segments to generate per edge.
  // We exponentially reduce resolution as level increases so it doesn't affect perf
  // when there are a large number of cells
  const resolution = Math.max(1, MAX_RESOLUTION * Math.pow(2, -level));

  for (let i = 0; i < 4; i++) {
    const offset = offsets[i].slice(0);
    const nextOffset = offsets[i + 1];
    const stepI = (nextOffset[0] - offset[0]) / resolution;
    const stepJ = (nextOffset[1] - offset[1]) / resolution;

    for (let j = 0; j < resolution; j++) {
      offset[0] += stepI;
      offset[1] += stepJ;
      // Cell can be represented by coordinates IJ, ST, UV, XYZ
      // http://s2geometry.io/devguide/s2cell_hierarchy#coordinate-systems
      const st = S2.IJToST(ij, level, offset);
      const uv = S2.STToUV(st);
      const xyz = S2.FaceUVToXYZ(face, uv);

      result.push(XYZToLngLat(xyz));
    }
  }
  return result;
}
1233540
def distance_to_edge(self, skydir):
        """"""

        xpix, ypix = skydir.to_pixel(self.wcs, origin=0)
        deltax = np.array((xpix - self._pix_center[0]) * self._pix_size[0],
                          ndmin=1)
        deltay = np.array((ypix - self._pix_center[1]) * self._pix_size[1],
                          ndmin=1)

        deltax = np.abs(deltax) - 0.5 * self._width[0]
        deltay = np.abs(deltay) - 0.5 * self._width[1]

        m0 = (deltax < 0) & (deltay < 0)
        m1 = (deltax > 0) & (deltay < 0)
        m2 = (deltax < 0) & (deltay > 0)
        m3 = (deltax > 0) & (deltay > 0)
        mx = np.abs(deltax) <= np.abs(deltay)
        my = np.abs(deltay) < np.abs(deltax)

        delta = np.zeros(len(deltax))
        delta[(m0 & mx) | (m3 & my) | m1] = deltax[(m0 & mx) | (m3 & my) | m1]
        delta[(m0 & my) | (m3 & mx) | m2] = deltay[(m0 & my) | (m3 & mx) | m2]
        return delta
1233316
def get_edge_mark(ttree):
    
    
    ## tree style
    if ttree._kwargs["tree_style"] in ["c", "cladogram"]:
        a=ttree.edges
        vcoordinates=ttree.verts
    else:
        a=ttree._lines               
        vcoordinates=ttree._coords    
   
    ## fixed args
    along='x'
    vmarker='o'
    vcolor=None
    vlshow=False            
    vsize=0.         
    estyle=ttree._kwargs["edge_style"]

    ## get axes
    layout = toyplot.layout.graph(a, vcoordinates=vcoordinates)
    along = toyplot.require.value_in(along, ["x", "y"])
    if along == "x":
        coordinate_axes = ["x", "y"]
    elif along == "y":
        coordinate_axes = ["y", "x"]
        
    ## broadcast args along axes
    vlabel = layout.vids
    vmarker = toyplot.broadcast.pyobject(vmarker, layout.vcount)
    vsize = toyplot.broadcast.scalar(vsize, layout.vcount)
    estyle = toyplot.style.require(estyle, allowed=toyplot.style.allowed.line)

    ## fixed args
    vcolor = toyplot.color.broadcast(colors=None, shape=layout.vcount, default=toyplot.color.black)
    vopacity = toyplot.broadcast.scalar(1.0, layout.vcount)
    vtitle = toyplot.broadcast.pyobject(None, layout.vcount)
    vstyle = None
    vlstyle = None
    
    ## this could be modified in the future to allow diff color edges
    ecolor = toyplot.color.broadcast(colors=None, shape=layout.ecount, default=toyplot.color.black)
    ewidth = toyplot.broadcast.scalar(1.0, layout.ecount)
    eopacity = toyplot.broadcast.scalar(1.0, layout.ecount)
    hmarker = toyplot.broadcast.pyobject(None, layout.ecount)
    mmarker = toyplot.broadcast.pyobject(None, layout.ecount)
    mposition = toyplot.broadcast.scalar(0.5, layout.ecount)
    tmarker = toyplot.broadcast.pyobject(None, layout.ecount)
    
    ## tables are required if I don't want to edit the class
    vtable = toyplot.data.Table()
    vtable["id"] = layout.vids
    for axis, coordinates in zip(coordinate_axes, layout.vcoordinates.T):
        vtable[axis] = coordinates
        #_mark_exportable(vtable, axis)
    vtable["label"] = vlabel
    vtable["marker"] = vmarker
    vtable["size"] = vsize
    vtable["color"] = vcolor
    vtable["opacity"] = vopacity
    vtable["title"] = vtitle

    etable = toyplot.data.Table()
    etable["source"] = layout.edges.T[0]
    #_mark_exportable(etable, "source")
    etable["target"] = layout.edges.T[1]
    #_mark_exportable(etable, "target")
    etable["shape"] = layout.eshapes
    etable["color"] = ecolor
    etable["width"] = ewidth
    etable["opacity"] = eopacity
    etable["hmarker"] = hmarker
    etable["mmarker"] = mmarker
    etable["mposition"] = mposition
    etable["tmarker"] = tmarker
    
    edge_mark = toyplot.mark.Graph(
        coordinate_axes=['x', 'y'],
        ecolor=["color"],
        ecoordinates=layout.ecoordinates,
        efilename=None,
        eopacity=["opacity"],
        eshape=["shape"],
        esource=["source"],
        estyle=estyle,
        etable=etable,
        etarget=["target"],
        ewidth=["width"],
        hmarker=["hmarker"],
        mmarker=["mmarker"],
        mposition=["mposition"],
        tmarker=["tmarker"],
        vcolor=["color"],
        vcoordinates=['x', 'y'],
        vfilename=None,
        vid=["id"],
        vlabel=["label"],
        vlshow=False,
        vlstyle=None,
        vmarker=["marker"],
        vopacity=["opacity"],
        vsize=["size"],
        vstyle=None,
        vtable=vtable,
        vtitle=["title"],
        )
    return edge_mark
1176157
def get_sampleEdge(self, res, DS=None, resMode='abs', offsetIn=0.):
        
        pts, dlr, ind = _comp._Ves_get_sampleEdge(self.Poly, res, DS=DS,
                                                  dLMode=resMode, DIn=offsetIn,
                                                  VIn=self.dgeom['VIn'],
                                                  margin=1.e-9)
        return pts, dlr, ind
726682
private void boundaryAreaExteriorArea_(int half_edge, int id_a, int id_b,
			int predicate) {
		if (m_matrix[predicate] == 1)
			return;

		int faceParentage = m_topo_graph.getHalfEdgeFaceParentage(half_edge);
		int twinFaceParentage = m_topo_graph
				.getHalfEdgeFaceParentage(m_topo_graph
						.getHalfEdgeTwin(half_edge));

		if ((faceParentage & id_b) == 0 && (twinFaceParentage & id_b) == 0)
			m_matrix[predicate] = 1;
	}
1040092
def _edge_opposite_point(self, tri, i):
        
        ind = tri.index(i)
        return (tri[(ind+1) % 3], tri[(ind+2) % 3])
550623
def midpoint_refine_triangulation_by_vertices(self, vertices):
        

        mlons, mlats = self.segment_midpoints_by_vertices(vertices=vertices)

        lonv1 = np.concatenate((self.lons, mlons), axis=0)
        latv1 = np.concatenate((self.lats, mlats), axis=0)

        return lonv1, latv1
478071
def _scale_mesh(self, scale):
        
        pos_ks = ['vertices', 'centers']

        # TODO: scale velocities???

        # handle scale
        self.update_columns_dict({k: self[k]*scale for k in pos_ks})

        self.update_columns(areas=self.areas*(scale**2))
        self._volume *= scale**3
        if self._area is not None:
            # self._area is None for wd meshes
            self._area += scale**2