Reconstruct point set using S(quared)R(ange)L(east)S(quares) method.
968432
def reconstruct_srls(edm, all_points, W=None, print_out=False, n=1, rescale=False,
                     z=None):
    
    from .lateration import SRLS, get_lateration_parameters
    Y = all_points.copy()
    indices = range(n)
    for index in indices:
        anchors, w, r2 = get_lateration_parameters(all_points, indices, index,
                                                   edm, W)
        if print_out:
            print('SRLS parameters:')
            print('anchors', anchors)
            print('w', w)
            print('r2', r2)

        srls = SRLS(anchors, w, r2, rescale, z, print_out)
        if rescale:
            srls = srls[0]  # second element of output is the scale
        Y[index, :] = srls
    return Y
353548
def _reshape_G(G):
    """"""
    G = np.moveaxis(G, (1, 3), (3, 4))
    (nsrc, nchan, filters_len) = G.shape[0:3]
    G = np.reshape(
        G, (nsrc * nchan * filters_len, nsrc * nchan * filters_len)
    )
    return G
967455
def warp_image_by_corner_points_projection(corner_points, image):
    
    # Clarify by storing in named variables.
    top_left, top_right, bottom_left, bottom_right = np.array(corner_points)

    top_edge = np.linalg.norm(top_right - top_left)
    bottom_edge = np.linalg.norm(bottom_right - bottom_left)
    left_edge = np.linalg.norm(top_left - bottom_left)
    right_edge = np.linalg.norm(top_right - bottom_right)

    L = int(np.ceil(max([top_edge, bottom_edge, left_edge, right_edge])))
    src = np.array([top_left, top_right, bottom_left, bottom_right])
    dst = np.array([[0, 0], [L - 1, 0], [0, L - 1], [L - 1, L - 1]])

    tr = ProjectiveTransform()
    tr.estimate(dst, src)
    warped_image = warp(image, tr, output_shape=(L, L))
    out = resize(warped_image, (500, 500))

    return out
916564
def reconstructImage(self):
        '''
        
        '''
        f_ishift = np.fft.ifftshift(self.fshift)
        return np.real(np.fft.ifft2(f_ishift))
456862
def _polygon_from_coords(coords, fix_geom=False, swap=True, dims=2):
    
    assert len(coords) % dims == 0
    number_of_points = len(coords)/dims
    coords_as_array = np.array(coords)
    reshaped = coords_as_array.reshape(number_of_points, dims)
    points = [
        (float(i[1]), float(i[0])) if swap else ((float(i[0]), float(i[1])))
        for i in reshaped.tolist()
        ]
    polygon = Polygon(points).buffer(0)
    try:
        assert polygon.is_valid
        return polygon
    except AssertionError:
        if fix_geom:
            return polygon.buffer(0)
        else:
            raise RuntimeError("Geometry is not valid.")
1181603
void constructLinearSystem(List<DMatrixRMaj> homographies,
							   List<List<PointIndex2D_F64>> observations) {
		// parameters are encoded points first then the
		int startView = totalFeatures*3;

		A.reshape(numEquations,numUnknown);

		DMatrixRMaj H = new DMatrixRMaj(3,3);
		Point2D_F64 p = new Point2D_F64();

		int row = 0;
		for (int viewIdx = 0; viewIdx < homographies.size(); viewIdx++) {
			N.apply(homographies.get(viewIdx),H);

			int colView = startView + viewIdx*3;

			List<PointIndex2D_F64> obs = observations.get(viewIdx);
			for (int i = 0; i < obs.size(); i++) {
				PointIndex2D_F64 p_pixel = obs.get(i);

				N.apply(p_pixel,p);

				// column this feature is at
				int col = p_pixel.index*3;

				// x component of pixel
				// A(row,colView) =  ...
				A.data[row*A.numCols + col    ] = p.x*H.get(2,0) - H.get(0,0);
				A.data[row*A.numCols + col + 1] = p.x*H.get(2,1) - H.get(0,1);
				A.data[row*A.numCols + col + 2] = p.x*H.get(2,2) - H.get(0,2);

				A.data[row*A.numCols + colView     ] = -1;
				A.data[row*A.numCols + colView + 1 ] = 0;
				A.data[row*A.numCols + colView + 2 ] = p.x;

				// y component of pixel
				row += 1;
				A.data[row*A.numCols + col    ] = p.y*H.get(2,0) - H.get(1,0);
				A.data[row*A.numCols + col + 1] = p.y*H.get(2,1) - H.get(1,1);
				A.data[row*A.numCols + col + 2] = p.y*H.get(2,2) - H.get(1,2);

				A.data[row*A.numCols + colView     ] = 0;
				A.data[row*A.numCols + colView + 1 ] = -1;
				A.data[row*A.numCols + colView + 2 ] = p.y;

				row += 1;
			}
		}
	}
1310374
def ds_cT(ds, x, y, xy_srs=wgs_srs):
    
    #Convert lat/lon to projected srs
    ds_srs = get_ds_srs(ds)
    #If xy_srs is undefined, assume it is the same as ds_srs
    mX = x
    mY = y
    if xy_srs is not None:
        if not ds_srs.IsSame(xy_srs):
            mX, mY, mZ = cT_helper(x, y, 0, xy_srs, ds_srs)
    return mX, mY
519200
def undistortPoints(self, points, keepSize=False):
        '''
        
        '''
        s = self.img.shape
        cam = self.coeffs['cameraMatrix']
        d = self.coeffs['distortionCoeffs']

        pts = np.asarray(points, dtype=np.float32)
        if pts.ndim == 2:
            pts = np.expand_dims(pts, axis=0)

        (newCameraMatrix, roi) = cv2.getOptimalNewCameraMatrix(cam,
                                                               d, s[::-1], 1,
                                                               s[::-1])
        if not keepSize:
            xx, yy = roi[:2]
            pts[0, 0] -= xx
            pts[0, 1] -= yy

        return cv2.undistortPoints(pts,
                                   cam, d, P=newCameraMatrix)
450117
def linearToPolar(img, center=None,
                  final_radius=None,
                  initial_radius=None,
                  phase_width=None,
                  interpolation=cv2.INTER_AREA, maps=None,
                  borderValue=0, borderMode=cv2.BORDER_REFLECT, **opts):
    '''
    
    '''
    if maps is None:
        mapY, mapX = linearToPolarMaps(img.shape[:2], center, final_radius,
                                       initial_radius, phase_width)
    else:
        mapY, mapX = maps

    o = {'interpolation': interpolation,
         'borderValue': borderValue,
         'borderMode': borderMode}
    o.update(opts)

    return cv2.remap(img, mapY, mapX, **o)
753378
def to_contour(grid_size=@npoints)
		m = Matrix.alloc(grid_size, grid_size)
		raise TypeError("Must be 3d data") unless @gridpoints.shape[1] == 2
# 			p @gridpoints.shape
# 			p @gridpoints
		xmax, xmin = @gridpoints.col(0).max, @gridpoints.col(0).min
		ymax, ymin = @gridpoints.col(1).max, @gridpoints.col(1).min
		p 'x', xmax, 'y', ymax
		xvec = GSL::Vector.alloc((0...grid_size).to_a) * (xmax - xmin) / (grid_size - 1.0).to_f + xmin
		yvec = GSL::Vector.alloc((0...grid_size).to_a) * (ymax - ymin) / (grid_size - 1.0).to_f + ymin	
		p 'x', xvec.max, 'y', yvec.max

		for i in 0...grid_size
			for j in 0...grid_size
# 					p xvec[i], yvec[j]
				m[i,j] = eval(xvec[i], yvec[j])
			end
		end
		p 'm', m.max
		Contour.alloc(xvec, yvec, m)
	end
608292
def _get_cartesian_edge_set(self):
        
        # Get projection space for cartesian projection
        edge_sets = numpy.vstack(self.edge_set)
        west, east, north, south = utils.get_spherical_bounding_box(
            edge_sets[:, 0],
            edge_sets[:, 1])
        self.proj = utils.OrthographicProjection(west, east, north, south)

        for edges in self.edge_set:
            # Project edges into cartesian space
            px, py = self.proj(edges[:, 0], edges[:, 1])
            # Store the two end-points of the trace
            self.cartesian_endpoints.append(
                numpy.array([[px[0], py[0], edges[0, 2]],
                             [px[-1], py[-1], edges[-1, 2]]]))
            self.cartesian_edges.append(numpy.column_stack([px, py,
                                                            edges[:, 2]]))
            # Get surface length vector for the trace - easier in cartesian
            lengths = numpy.sqrt((px[:-1] - px[1:]) ** 2. +
                                 (py[:-1] - py[1:]) ** 2.)
            self.length_set.append(lengths)
            # Get cumulative surface length vector
            self.cum_length_set.append(
                numpy.hstack([0., numpy.cumsum(lengths)]))
        return edge_sets
975442
def convex_hull(features):
    """"""
    points = sorted([s.point() for s in features])
    l = reduce(_keep_left, points, [])
    u = reduce(_keep_left, reversed(points), [])
    return l.extend(u[i] for i in xrange(1, len(u) - 1)) or l
100378
def recoSurface(points, bins=256):
    

    if isinstance(points, vtk.vtkActor):
        points = points.coordinates()
    N = len(points)
    if N < 50:
        print("recoSurface: Use at least 50 points.")
        return None
    points = np.array(points)

    ptsSource = vtk.vtkPointSource()
    ptsSource.SetNumberOfPoints(N)
    ptsSource.Update()
    vpts = ptsSource.GetOutput().GetPoints()
    for i, p in enumerate(points):
        vpts.SetPoint(i, p)
    polyData = ptsSource.GetOutput()

    distance = vtk.vtkSignedDistance()
    f = 0.1
    x0, x1, y0, y1, z0, z1 = polyData.GetBounds()
    distance.SetBounds(x0-(x1-x0)*f, x1+(x1-x0)*f,
                       y0-(y1-y0)*f, y1+(y1-y0)*f,
                       z0-(z1-z0)*f, z1+(z1-z0)*f)
    if polyData.GetPointData().GetNormals():
        distance.SetInputData(polyData)
    else:
        normals = vtk.vtkPCANormalEstimation()
        normals.SetInputData(polyData)
        normals.SetSampleSize(int(N / 50))
        normals.SetNormalOrientationToGraphTraversal()
        distance.SetInputConnection(normals.GetOutputPort())
        print("Recalculating normals for", N, "Points, sample size=", int(N / 50))

    b = polyData.GetBounds()
    diagsize = np.sqrt((b[1] - b[0]) ** 2 + (b[3] - b[2]) ** 2 + (b[5] - b[4]) ** 2)
    radius = diagsize / bins * 5
    distance.SetRadius(radius)
    distance.SetDimensions(bins, bins, bins)
    distance.Update()

    print("Calculating mesh from points with R =", radius)
    surface = vtk.vtkExtractSurface()
    surface.SetRadius(radius * 0.99)
    surface.HoleFillingOn()
    surface.ComputeNormalsOff()
    surface.ComputeGradientsOff()
    surface.SetInputConnection(distance.GetOutputPort())
    surface.Update()
    return Actor(surface.GetOutput(), "gold", 1, 0, "tomato")
438989
def qlaplace_dd(d):
    """"""
    res = _matrix.matrix()
    d0 = d[::-1]
    D = len(d0)
    I = _np.eye(2)
    J = _np.array([[0, 1], [0, 0]])
    cr = []
    if D is 1:
        for k in xrange(1, d0[0] + 1):
            if k is 1:
                cur_core = _np.zeros((1, 2, 2, 3))
                cur_core[:, :, :, 0] = 2 * I - J - J.T
                cur_core[:, :, :, 1] = -J
                cur_core[:, :, :, 2] = -J.T
            elif k is d0[0]:
                cur_core = _np.zeros((3, 2, 2, 1))
                cur_core[0, :, :, 0] = I
                cur_core[1, :, :, 0] = J.T
                cur_core[2, :, :, 0] = J
            else:
                cur_core = _np.zeros((3, 2, 2, 3))
                cur_core[0, :, :, 0] = I
                cur_core[1, :, :, 1] = J
                cur_core[2, :, :, 2] = J.T
                cur_core[1, :, :, 0] = J.T
                cur_core[2, :, :, 0] = J
            cr.append(cur_core)
    else:
        for k in xrange(D):
            for kappa in xrange(1, d0[k] + 1):
                if kappa is 1:
                    if k is 0:
                        cur_core = _np.zeros((1, 2, 2, 4))
                        cur_core[:, :, :, 0] = 2 * I - J - J.T
                        cur_core[:, :, :, 1] = -J
                        cur_core[:, :, :, 2] = -J.T
                        cur_core[:, :, :, 3] = I
                    elif k is D - 1:
                        cur_core = _np.zeros((2, 2, 2, 3))
                        cur_core[0, :, :, 0] = 2 * I - J - J.T
                        cur_core[0, :, :, 1] = -J
                        cur_core[0, :, :, 2] = -J.T
                        cur_core[1, :, :, 0] = I
                    else:
                        cur_core = _np.zeros((2, 2, 2, 4))
                        cur_core[0, :, :, 0] = 2 * I - J - J.T
                        cur_core[0, :, :, 1] = -J
                        cur_core[0, :, :, 2] = -J.T
                        cur_core[0, :, :, 3] = I
                        cur_core[1, :, :, 0] = I
                elif kappa is d0[k]:
                    if k is D - 1:
                        cur_core = _np.zeros((3, 2, 2, 1))
                        cur_core[0, :, :, 0] = I
                        cur_core[1, :, :, 0] = J.T
                        cur_core[2, :, :, 0] = J
                    else:
                        cur_core = _np.zeros((4, 2, 2, 2))
                        cur_core[3, :, :, 0] = I
                        cur_core[0, :, :, 1] = I
                        cur_core[1, :, :, 1] = J.T
                        cur_core[2, :, :, 1] = J
                else:
                    if k is D - 1:
                        cur_core = _np.zeros((3, 2, 2, 3))
                        cur_core[0, :, :, 0] = I
                        cur_core[1, :, :, 1] = J
                        cur_core[2, :, :, 2] = J.T
                        cur_core[1, :, :, 0] = J.T
                        cur_core[2, :, :, 0] = J
                    else:
                        cur_core = _np.zeros((4, 2, 2, 4))
                        cur_core[0, :, :, 0] = I
                        cur_core[1, :, :, 1] = J
                        cur_core[2, :, :, 2] = J.T
                        cur_core[1, :, :, 0] = J.T
                        cur_core[2, :, :, 0] = J
                        cur_core[3, :, :, 3] = I
                cr.append(cur_core)
    return _matrix.matrix.from_list(cr)
318022
def _cdist_scipy(x, y, exponent=1):
    """"""
    metric = 'euclidean'

    if exponent != 1:
        metric = 'sqeuclidean'

    distances = _spatial.distance.cdist(x, y, metric=metric)

    if exponent != 1:
        distances **= exponent / 2

    return distances
1346130
def pixel_kernel(point_source_kernel, subgrid_res=7):
    
    kernel_subgrid = subgrid_kernel(point_source_kernel, subgrid_res, num_iter=10)
    kernel_size = len(point_source_kernel)
    kernel_pixel = np.zeros((kernel_size*subgrid_res, kernel_size*subgrid_res))
    for i in range(subgrid_res):
        k_x = int((kernel_size-1) / 2 * subgrid_res + i)
        for j in range(subgrid_res):
            k_y = int((kernel_size-1) / 2 * subgrid_res + j)
            kernel_pixel = image_util.add_layer2image(kernel_pixel, k_x, k_y, kernel_subgrid)
    kernel_pixel = util.averaging(kernel_pixel, numGrid=kernel_size*subgrid_res, numPix=kernel_size)
    return kernel_norm(kernel_pixel)
902691
def to_scanner(self, x, y, z):
        
        if self.transform is None:
            raise ValueError("No transform set for MRSData object {}".format(self))

        transformed_point = self.transform * numpy.matrix([x, y, z, 1]).T

        return numpy.squeeze(numpy.asarray(transformed_point))[0:3]