Interpolates along a jagged path in 3D
    
    Args:
        xyz = section path specified in cartesian coordinates
        nseg = number of segment paths in section path
        
    Returns:
        interp_xyz = interpolated path
174016
def interpolate_jagged(xyz,nseg):
    
    
    # Spherical coordinates specifying the angles of all line
    # segments that make up the section path
    (r,theta,phi) = sequential_spherical(xyz)
    
    # cumulative length of section path at each coordinate
    rcum = np.append(0,np.cumsum(r))

    # breakpoints for segment paths along section path
    breakpoints = np.linspace(0,rcum[-1],nseg+1)
    np.delete(breakpoints,0)
    
    # Find segment paths
    seg_paths = []
    for a in range(nseg):
        path = []
        
        # find (x,y,z) starting coordinate of path
        if a == 0:
            start_coord = xyz[0,:]
        else:
            start_coord = end_coord # start at end of last path
        path.append(start_coord)

        # find all coordinates between the start and end points
        start_length = breakpoints[a]
        end_length = breakpoints[a+1]
        mid_boolean = (rcum > start_length) & (rcum < end_length)
        mid_indices = np.nonzero(mid_boolean)[0]
        for mi in mid_indices:
            path.append(xyz[mi,:])

        # find (x,y,z) ending coordinate of path
        end_coord = find_coord(end_length,xyz,rcum,theta,phi)
        path.append(end_coord)

        # Append path to list of segment paths
        seg_paths.append(np.array(path))
    
    # Return all segment paths
    return seg_paths
505101
def interpolate_cubic(self, lons, lats, data):
        
        return self.interpolate(lons, lats, data, order=3)
514772
private int segdir(double xend, double yend, int[] ij) {
        if (YMATCH(yend, y[ij[1]])) {
            if (ij[1] == 0) {
                return 0;
            }
            ij[1] -= 1;
            return 3;
        }

        if (XMATCH(xend, x[ij[0]])) {
            if (ij[0] == 0) {
                return 0;
            }
            ij[0] -= 1;
            return 4;
        }

        if (YMATCH(yend, y[ij[1] + 1])) {
            if (ij[1] >= y.length - 1) {
                return 0;
            }
            ij[1] += 1;
            return 1;
        }

        if (XMATCH(xend, x[ij[0] + 1])) {
            if (ij[0] >= x.length - 1) {
                return 0;
            }
            ij[0] += 1;
            return 2;
        }

        return 0;
    }
255201
def interpolate(self, lon, lat, egy=None, interp_log=True):
        

        if self.data.ndim == 1:
            theta = np.pi / 2. - np.radians(lat)
            phi = np.radians(lon)
            return hp.pixelfunc.get_interp_val(self.counts, theta,
                                               phi, nest=self.hpx.nest)
        else:
            return self._interpolate_cube(lon, lat, egy, interp_log)
1281742
def get_xyz_1d(self, xcoord, x, ycoord, y, u, v):
        """"""
        xclose = xcoord.indexes[xcoord.name].get_loc(x, method='nearest')
        yclose = ycoord.indexes[ycoord.name].get_loc(y, method='nearest')
        uval = u[yclose, xclose].values
        vval = v[yclose, xclose].values
        return xcoord[xclose].values, ycoord[yclose].values, uval, vval
958265
public SVGPath relativeSmoothCubicTo(double[] c2xy, double[] xy) {
    return append(PATH_SMOOTH_CUBIC_TO_RELATIVE).append(c2xy[0]).append(c2xy[1]).append(xy[0]).append(xy[1]);
  }
486555
def window_xy x,y
      s = ""
      s << "+" if x > 0
      s << x.to_s
      s << "+" if y > 0
      s << y.to_s
      @root.geometry s
    end
538758
def clamped_interpolation(skimage_range, sinogram):
    
    min_x = skimage_range.domain.min()[1]
    max_x = skimage_range.domain.max()[1]

    def interpolation_wrapper(x):
        x = (x[0], np.maximum(min_x, np.minimum(max_x, x[1])))

        return sinogram.interpolation(x)
    return interpolation_wrapper
756783
public SVGPath quadTo(double c1x, double c1y, double x, double y) {
    return append(PATH_QUAD_TO).append(c1x).append(c1y).append(x).append(y);
  }
683043
def xy2geom(x, y, t_srs=None):
    
    geom_wkt = 'POINT({0} {1})'.format(x, y)
    geom = ogr.CreateGeometryFromWkt(geom_wkt)
    if t_srs is not None and not wgs_srs.IsSame(t_srs):
        ct = osr.CoordinateTransformation(t_srs, wgs_srs)
        geom.Transform(ct)
        geom.AssignSpatialReference(t_srs)
    return geom
554458
public static void srgbToXyz( double r , double g , double b , double xyz[] ) {
		xyz[0] = 0.412453*r + 0.35758*g + 0.180423*b;
		xyz[1] = 0.212671*r + 0.71516*g + 0.072169*b;
		xyz[2] = 0.019334*r + 0.119193*g + 0.950227*b;
	}
872116
def smoothpolygon(document, coords):
	""
	
	element = document.createElement('path')
	path    = []
	points  = [(coords[i], coords[i+1]) for i  in range(0, len(coords), 2)]
	def pt(points):
		p = points
		n = len(points)
		for i in range(0, len(points)):
			a = p[(i-1) % n]
			b = p[i]
			c = p[(i+1) % n]

			yield lerp(a, b, 0.5), b, lerp(b, c, 0.5)
		
	for i, (A, B, C) in enumerate(pt(points)):
		if i == 0:
			path.append("M%s,%s Q%s,%s %s,%s" % (A[0], A[1], B[0], B[1], C[0], C[1]))
		else:
			path.append("T%s,%s" % (C[0], C[1]))
	
	path.append("z")

	element.setAttribute('d', ' '.join(path))
	return element
68285
def xyzinterp(self, x, y, z, titletext):
    
    # Help from http://wiki.scipy.org/Cookbook/Matplotlib/Gridding_irregularly_spaced_data
    
    if self.Verbose:
      print("Starting to interpolate grid for plotting -- can be a slow process!")
    
    from scipy.interpolate import griddata
    import numpy.ma as ma
    
    # define grid.
    xmin = np.min(self.xw)
    xmean = np.mean(self.xw) # not used right now
    xmax = np.max(self.xw)
    ymin = np.min(self.yw)
    ymean = np.mean(self.yw) # not used right now
    ymax = np.max(self.yw)
    x_range = xmax - xmin
    y_range = ymax - ymin
    
    # x and y grids
    # 100 cells on each side -- just for plotting, not so important
    # to optimize with how many points are plotted
    #xi = np.linspace(xmin-.05*x_range, xmax+.05*x_range, 200)
    #yi = np.linspace(ymin-.05*y_range, ymax+.05*y_range, 200)
    xi = np.linspace(xmin, xmax, 200)
    yi = np.linspace(ymin, ymax, 200)
    # grid the z-axis
    zi = griddata((x, y), z, (xi[None,:], yi[:,None]), method='cubic')
    # turn nan into 0 -- this will just be outside computation area for q
    zi[np.isnan(zi)] = 0
    # contour the gridded outputs, plotting dots at the randomly spaced data points.
    #CS = plt.contour(xi,yi,zi,15,linewidths=0.5,colors='k') -- don't need lines
    if self.latlon:
      CS = plt.contourf(xi, yi, zi, 100, cmap=plt.cm.jet)
    else:
      CS = plt.contourf(xi/1000., yi/1000., zi, 100, cmap=plt.cm.jet)
    plt.colorbar() # draw colorbar
    # plot model points.
    # Computed at
    if self.latlon:
      plt.plot(x, y, 'o', markerfacecolor='.6', markeredgecolor='.6', markersize=1)
      plt.plot(self.x, self.y, 'o', markerfacecolor='.2', markeredgecolor='.2', markersize=1)
    else:
      plt.plot(x/1000., y/1000., 'o', markerfacecolor='.6', markeredgecolor='.6', markersize=1)
      # Load sources (overlay computed at)
      plt.plot(self.x/1000., self.y/1000., 'o', markerfacecolor='.2', markeredgecolor='.2', markersize=1)
      #plt.hexbin(self.x, self.y, C=self.w) -- show colors on points -- harder to see
    if self.latlon:
      plt.xlabel('longitude [deg E]', fontsize=12, fontweight='bold')
      plt.ylabel('latitude [deg N]', fontsize=12, fontweight='bold')
    else:
      plt.xlabel('x [km]', fontsize=12, fontweight='bold')
      plt.ylabel('y [km]', fontsize=12, fontweight='bold')
    # Limits -- to not get messed up by points (view wants to be wider so whole point visible)
    if self.latlon:
      plt.xlim( (xi[0], xi[-1]) )
      plt.ylim( (yi[0], yi[-1]) )
    else:
      plt.xlim( (xi[0]/1000., xi[-1]/1000.) )
      plt.ylim( (yi[0]/1000., yi[-1]/1000.) )
    # Title
    plt.title(titletext, fontsize=16)
625217
void interpolateAttributes(int from_path_index, int from_point_index,
			int to_path_index, int to_point_index) {
		for (int ipath = from_path_index; ipath < to_path_index - 1; ipath++) {
			if (isClosedPath(ipath))
				throw new IllegalArgumentException(
						"cannot interpolate across closed paths");
		}

		int nattr = m_description.getAttributeCount();

		if (nattr == 1)
			return; // only has position

		double sub_length = calculateSubLength2D(from_path_index,
				from_point_index, to_path_index, to_point_index);

		if (sub_length == 0.0)
			return;

		for (int iattr = 1; iattr < nattr; iattr++) {
			int semantics = m_description.getSemantics(iattr);

			int interpolation = VertexDescription.getInterpolation(semantics);
			if (interpolation == VertexDescription.Interpolation.ANGULAR)
				continue;

			int components = VertexDescription.getComponentCount(semantics);

			for (int ordinate = 0; ordinate < components; ordinate++)
				interpolateAttributes_(semantics, from_path_index,
						from_point_index, to_path_index, to_point_index,
						sub_length, ordinate);
		}
	}

	// TODO: Add code for interpolation type (none and angular)
	void interpolateAttributesForSemantics(int semantics, int from_path_index,
			int from_point_index, int to_path_index, int to_point_index) {
		if (semantics == VertexDescription.Semantics.POSITION)
			return;

		if (!hasAttribute(semantics))
			throw new IllegalArgumentException(
					"does not have the given attribute");

		int interpolation = VertexDescription.getInterpolation(semantics);
		if (interpolation == VertexDescription.Interpolation.ANGULAR)
			throw new IllegalArgumentException(
					"not implemented for the given semantics");

		for (int ipath = from_path_index; ipath < to_path_index - 1; ipath++) {
			if (isClosedPath(ipath))
				throw new IllegalArgumentException(
						"cannot interpolate across closed paths");
		}

		double sub_length = calculateSubLength2D(from_path_index,
				from_point_index, to_path_index, to_point_index);

		if (sub_length == 0.0)
			return;

		int components = VertexDescription.getComponentCount(semantics);

		for (int ordinate = 0; ordinate < components; ordinate++)
			interpolateAttributes_(semantics, from_path_index,
					from_point_index, to_path_index, to_point_index,
					sub_length, ordinate);
	}

	void interpolateAttributes(int path_index, int from_point_index,
			int to_point_index) {
		int nattr = m_description.getAttributeCount();

		if (nattr == 1)
			return; // only has position

		double sub_length = calculateSubLength2D(path_index, from_point_index,
				to_point_index);

		if (sub_length == 0.0)
			return;

		for (int iattr = 1; iattr < nattr; iattr++) {
			int semantics = m_description.getSemantics(iattr);

			int interpolation = VertexDescription.getInterpolation(semantics);
			if (interpolation == VertexDescription.Interpolation.ANGULAR)
				continue;

			int components = VertexDescription.getComponentCount(semantics);

			for (int ordinate = 0; ordinate < components; ordinate++)
				interpolateAttributes_(semantics, path_index, from_point_index,
						to_point_index, sub_length, ordinate);
		}
	}

	void interpolateAttributesForSemantics(int semantics, int path_index,
			int from_point_index, int to_point_index) {
		if (semantics == VertexDescription.Semantics.POSITION)
			return;

		if (!hasAttribute(semantics))
			throw new IllegalArgumentException(
					"does not have the given attribute");

		int interpolation = VertexDescription.getInterpolation(semantics);
		if (interpolation == VertexDescription.Interpolation.ANGULAR)
			throw new IllegalArgumentException(
					"not implemented for the given semantics");

		double sub_length = calculateSubLength2D(path_index, from_point_index,
				to_point_index);

		if (sub_length == 0.0)
			return;

		int components = VertexDescription.getComponentCount(semantics);

		for (int ordinate = 0; ordinate < components; ordinate++)
			interpolateAttributes_(semantics, path_index, from_point_index,
					to_point_index, sub_length, ordinate);
	}

	// 
	void interpolateAttributes_(int semantics, int from_path_index,
			int from_point_index, int to_path_index, int to_point_index,
			double sub_length, int ordinate) {
		SegmentIteratorImpl seg_iter = querySegmentIterator();

		int absolute_from_index = getPathStart(from_path_index)
				+ from_point_index;
		int absolute_to_index = getPathStart(to_path_index) + to_point_index;

		double from_attribute = getAttributeAsDbl(semantics,
				absolute_from_index, ordinate);
		double to_attribute = getAttributeAsDbl(semantics, absolute_to_index,
				ordinate);
		double interpolated_attribute = from_attribute;
		double cumulative_length = 0.0;

		seg_iter.resetToVertex(absolute_from_index);

		do {
			if (seg_iter.hasNextSegment()) {
				seg_iter.nextSegment();

				if (seg_iter.getStartPointIndex() == absolute_to_index)
					return;

				setAttribute(semantics, seg_iter.getStartPointIndex(),
						ordinate, interpolated_attribute);

				seg_iter.previousSegment();

				do {
					Segment segment = seg_iter.nextSegment();

					if (seg_iter.getEndPointIndex() == absolute_to_index)
						return;

					double segment_length = segment.calculateLength2D();
					cumulative_length += segment_length;
					double t = cumulative_length / sub_length;
					interpolated_attribute = MathUtils.lerp(from_attribute,  to_attribute, t);

					if (!seg_iter.isClosingSegment())
						setAttribute(semantics, seg_iter.getEndPointIndex(),
								ordinate, interpolated_attribute);

				} while (seg_iter.hasNextSegment());
			}

		} while (seg_iter.nextPath());
	}

	void interpolateAttributes_(int semantics, int path_index,
			int from_point_index, int to_point_index, double sub_length,
			int ordinate) {
		assert (m_bPolygon);
		SegmentIteratorImpl seg_iter = querySegmentIterator();

		int absolute_from_index = getPathStart(path_index) + from_point_index;
		int absolute_to_index = getPathStart(path_index) + to_point_index;

		if (absolute_to_index == absolute_from_index)
			return;

		double from_attribute = getAttributeAsDbl(semantics,
				absolute_from_index, ordinate);
		double to_attribute = getAttributeAsDbl(semantics, absolute_to_index,
				ordinate);
		double cumulative_length = 0.0;

		seg_iter.resetToVertex(absolute_from_index);
		seg_iter.setCirculator(true);

		double prev_interpolated_attribute = from_attribute;

		do {
			Segment segment = seg_iter.nextSegment();
			setAttribute(semantics, seg_iter.getStartPointIndex(), ordinate,
					prev_interpolated_attribute);

			double segment_length = segment.calculateLength2D();
			cumulative_length += segment_length;
			double t = cumulative_length / sub_length;
			prev_interpolated_attribute = MathUtils.lerp(from_attribute, to_attribute, t);

		} while (seg_iter.getEndPointIndex() != absolute_to_index);
	}

	@Override
	public void setEmpty() {
		m_curveParamwritePoint = 0;
		m_bPathStarted = false;
		m_paths = null;
		m_pathFlags = null;
		m_segmentParamIndex = null;
		m_segmentFlags = null;
		m_segmentParams = null;
		_setEmptyImpl();
	}

	@Override
	public void applyTransformation(Transformation2D transform) {
		applyTransformation(transform, -1);
	}

	public void applyTransformation(Transformation2D transform, int pathIndex) {
		if (isEmpty())
			return;

		if (transform.isIdentity())
			return;

		_verifyAllStreams();
		AttributeStreamOfDbl points = (AttributeStreamOfDbl) m_vertexAttributes[0];
		Point2D ptStart = new Point2D();
		Point2D ptControl = new Point2D();

		boolean bHasNonLinear;
		int fistIdx;
		int lastIdx;
		if (pathIndex < 0) {
			bHasNonLinear = hasNonLinearSegments();
			fistIdx = 0;
			lastIdx = m_pointCount;
		} else {
			bHasNonLinear = hasNonLinearSegments(pathIndex);
			fistIdx = getPathStart(pathIndex);
			lastIdx = getPathEnd(pathIndex);
		}

		for (int ipoint = fistIdx; ipoint < lastIdx; ipoint++) {
			ptStart.x = points.read(ipoint * 2);
			ptStart.y = points.read(ipoint * 2 + 1);

			if (bHasNonLinear) {
				int segIndex = m_segmentParamIndex.read(ipoint);
				if (segIndex >= 0) {
					int segmentType = (int) m_segmentFlags.read(ipoint);
					int type = segmentType & SegmentFlags.enumSegmentMask;
					switch (type) {
					case SegmentFlags.enumBezierSeg: {
						ptControl.x = m_segmentParams.read(segIndex);
						ptControl.y = m_segmentParams.read(segIndex + 1);
						transform.transform(ptControl, ptControl);
						m_segmentParams.write(segIndex, ptControl.x);
						m_segmentParams.write(segIndex + 1, ptControl.y);

						ptControl.x = m_segmentParams.read(segIndex + 3);
						ptControl.y = m_segmentParams.read(segIndex + 4);
						transform.transform(ptControl, ptControl);
						m_segmentParams.write(segIndex + 3, ptControl.x);
						m_segmentParams.write(segIndex + 4, ptControl.y);
					}
						break;
					case SegmentFlags.enumArcSeg:
						throw GeometryException.GeometryInternalError();

					}
				}
			}

			transform.transform(ptStart, ptStart);
			points.write(ipoint * 2, ptStart.x);
			points.write(ipoint * 2 + 1, ptStart.y);
		}

		notifyModified(DirtyFlags.DirtyCoordinates);
		// REFACTOR: reset the exact envelope only and transform the loose
		// envelope
	}

	@Override
	public void applyTransformation(Transformation3D transform) {
		if (isEmpty())
			return;

		addAttribute(VertexDescription.Semantics.Z);
		_verifyAllStreams();
		AttributeStreamOfDbl points = (AttributeStreamOfDbl) m_vertexAttributes[0];
		AttributeStreamOfDbl zs = (AttributeStreamOfDbl) m_vertexAttributes[1];
		Point3D ptStart = new Point3D();
		Point3D ptControl = new Point3D();
		boolean bHasNonLinear = hasNonLinearSegments();
		for (int ipoint = 0; ipoint < m_pointCount; ipoint++) {
			ptStart.x = points.read(ipoint * 2);
			ptStart.y = points.read(ipoint * 2 + 1);
			ptStart.z = zs.read(ipoint);

			if (bHasNonLinear) {
				int segIndex = m_segmentParamIndex.read(ipoint);
				if (segIndex >= 0) {
					int segmentType = (int) m_segmentFlags.read(ipoint);
					int type = segmentType & (int) SegmentFlags.enumSegmentMask;
					switch (type) {
					case SegmentFlags.enumBezierSeg: {
						ptControl.x = m_segmentParams.read(segIndex);
						ptControl.y = m_segmentParams.read(segIndex + 1);
						ptControl.z = m_segmentParams.read(segIndex + 2);
						ptControl = transform.transform(ptControl);
						m_segmentParams.write(segIndex, ptControl.x);
						m_segmentParams.write(segIndex + 1, ptControl.y);
						m_segmentParams.write(segIndex + 1, ptControl.z);

						ptControl.x = m_segmentParams.read(segIndex + 3);
						ptControl.y = m_segmentParams.read(segIndex + 4);
						ptControl.z = m_segmentParams.read(segIndex + 5);
						ptControl = transform.transform(ptControl);
						m_segmentParams.write(segIndex + 3, ptControl.x);
						m_segmentParams.write(segIndex + 4, ptControl.y);
						m_segmentParams.write(segIndex + 5, ptControl.z);
					}
						break;
					case SegmentFlags.enumArcSeg:
						throw GeometryException.GeometryInternalError();

					}
				}
			}

			ptStart = transform.transform(ptStart);
			points.write(ipoint * 2, ptStart.x);
			points.write(ipoint * 2 + 1, ptStart.y);
			zs.write(ipoint, ptStart.z);
		}

		// REFACTOR: reset the exact envelope only and transform the loose
		// envelope

		notifyModified(DirtyFlags.DirtyCoordinates);
	}

	@Override
	protected void _verifyStreamsImpl() {
		if (m_paths == null) {
			m_paths = (AttributeStreamOfInt32) AttributeStreamBase
					.createIndexStream(1, 0);
			m_pathFlags = (AttributeStreamOfInt8) AttributeStreamBase
					.createByteStream(1, (byte) 0);
		}

		if (m_segmentFlags != null) {
			m_segmentFlags.resize(m_reservedPointCount,
					(byte) SegmentFlags.enumLineSeg);
			m_segmentParamIndex.resize(m_reservedPointCount, -1);
		}
	}

	@Override
	void _copyToImpl(MultiVertexGeometryImpl dst) {
		MultiPathImpl dstPoly = (MultiPathImpl) dst;
		dstPoly.m_bPathStarted = false;
		dstPoly.m_curveParamwritePoint = m_curveParamwritePoint;
		dstPoly.m_fill_rule = m_fill_rule;
		
		if (m_paths != null)
			dstPoly.m_paths = new AttributeStreamOfInt32(m_paths);
		else
			dstPoly.m_paths = null;

		if (m_pathFlags != null)
			dstPoly.m_pathFlags = new AttributeStreamOfInt8(m_pathFlags);
		else
			dstPoly.m_pathFlags = null;

		if (m_segmentParamIndex != null)
			dstPoly.m_segmentParamIndex = new AttributeStreamOfInt32(
					m_segmentParamIndex);
		else
			dstPoly.m_segmentParamIndex = null;

		if (m_segmentFlags != null)
			dstPoly.m_segmentFlags = new AttributeStreamOfInt8(m_segmentFlags);
		else
			dstPoly.m_segmentFlags = null;

		if (m_segmentParams != null)
			dstPoly.m_segmentParams = new AttributeStreamOfDbl(m_segmentParams);
		else
			dstPoly.m_segmentParams = null;

		dstPoly.m_cachedLength2D = m_cachedLength2D;
		dstPoly.m_cachedArea2D = m_cachedArea2D;

		if (!_hasDirtyFlag(DirtyFlags.DirtyRingAreas2D)) {
			dstPoly.m_cachedRingAreas2D = (AttributeStreamOfDbl) m_cachedRingAreas2D;
		} else
			dstPoly.m_cachedRingAreas2D = null;

	}

	@Override
	public double calculateLength2D() {
		if (!_hasDirtyFlag(DirtyFlags.DirtyLength2D)) {
			return m_cachedLength2D;
		}

		SegmentIteratorImpl segIter = querySegmentIterator();
		MathUtils.KahanSummator len = new MathUtils.KahanSummator(0);
		while (segIter.nextPath()) {
			while (segIter.hasNextSegment()) {
				len.add(segIter.nextSegment().calculateLength2D());
			}
		}

		m_cachedLength2D = len.getResult();
		_setDirtyFlag(DirtyFlags.DirtyLength2D, false);

		return len.getResult();
	}

	@Override
	public boolean equals(Object other) {
		if (other == this)
			return true;

		if (!(other instanceof MultiPathImpl))
			return false;

		if (!super.equals(other))
			return false;

		MultiPathImpl otherMultiPath = (MultiPathImpl) other;

		int pathCount = getPathCount();
		int pathCountOther = otherMultiPath.getPathCount();

		if (pathCount != pathCountOther)
			return false;

		if (pathCount > 0 && m_paths != null
				&& !m_paths.equals(otherMultiPath.m_paths, 0, pathCount + 1))
			return false;

		if (m_fill_rule != otherMultiPath.m_fill_rule)
			return false;

		{
			// Note: OGC flags do not participate in the equals operation by
			// design.
			// Because for the polygon pathFlags will have all enum_closed set,
			// we do not need to compare this stream. Only for polyline.
			// Polyline does not have OGC flags set.
			if (!m_bPolygon) {
				if (m_pathFlags != null
						&& !m_pathFlags.equals(otherMultiPath.m_pathFlags, 0,
								pathCount))
					return false;
			}
		}
	      
		return super.equals(other);
	}

	/**
	 * Returns a SegmentIterator that set to a specific vertex of the
	 * MultiPathImpl. The call to NextSegment will return the segment that
	 * starts at the vertex. Call to PreviousSegment will return the segment
	 * that starts at the previous vertex.
	 */
	public SegmentIteratorImpl querySegmentIteratorAtVertex(int startVertexIndex) {
		if (startVertexIndex < 0 || startVertexIndex >= getPointCount())
			throw new IndexOutOfBoundsException();

		SegmentIteratorImpl iter = new SegmentIteratorImpl(this,
				startVertexIndex);
		return iter;
	}

	// void QuerySegmentIterator(int fromVertex, SegmentIterator iterator);
	public SegmentIteratorImpl querySegmentIterator() {
		return new SegmentIteratorImpl(this);
	}

	@Override
	public void _updateXYImpl(boolean bExact) {
		super._updateXYImpl(bExact);
		boolean bHasCurves = hasNonLinearSegments();
		if (bHasCurves) {
			SegmentIteratorImpl segIter = querySegmentIterator();
			while (segIter.nextPath()) {
				while (segIter.hasNextSegment()) {
					Segment curve = segIter.nextCurve();
					if (curve != null) {
						Envelope2D env2D = new Envelope2D();
						curve.queryEnvelope2D(env2D);
						m_envelope.merge(env2D);
					} else
						break;
				}
			}
		}
	}

	@Override
	void calculateEnvelope2D(Envelope2D env, boolean bExact) {
		super.calculateEnvelope2D(env, bExact);
		boolean bHasCurves = hasNonLinearSegments();
		if (bHasCurves) {
			SegmentIteratorImpl segIter = querySegmentIterator();
			while (segIter.nextPath()) {
				while (segIter.hasNextSegment()) {
					Segment curve = segIter.nextCurve();
					if (curve != null) {
						Envelope2D env2D = new Envelope2D();
						curve.queryEnvelope2D(env2D);
						env.merge(env2D);
					} else
						break;
				}
			}
		}
	}

	@Override
	public void _notifyModifiedAllImpl() {
		if (m_paths == null || m_paths.size() == 0)// if (m_paths == null ||
													// !m_paths.size())
			m_pointCount = 0;
		else
			m_pointCount = m_paths.read(m_paths.size() - 1);
	}

	@Override
	public double calculateArea2D() {
		if (!m_bPolygon)
			return 0.0;

		_updateRingAreas2D();

		return m_cachedArea2D;
	}

	/**
	 * Returns True if the ring is an exterior ring. Valid only for simple
	 * polygons.
	 */
	public boolean isExteriorRing(int ringIndex) {
		if (!m_bPolygon)
			return false;

		if (!_hasDirtyFlag(DirtyFlags.DirtyOGCFlags))
			return (m_pathFlags.read(ringIndex) & (byte) PathFlags.enumOGCStartPolygon) != 0;

		_updateRingAreas2D();
		return m_cachedRingAreas2D.read(ringIndex) > 0;
		// Should we make a function called _UpdateHasNonLinearSegmentsFlags and
		// call it here?
	}

	public double calculateRingArea2D(int pathIndex) {
		if (!m_bPolygon)
			return 0.0;

		_updateRingAreas2D();

		return m_cachedRingAreas2D.read(pathIndex);
	}

	public void _updateRingAreas2D() {
		if (_hasDirtyFlag(DirtyFlags.DirtyRingAreas2D)) {
			int pathCount = getPathCount();

			if (m_cachedRingAreas2D == null)
				m_cachedRingAreas2D = new AttributeStreamOfDbl(pathCount);
			else if (m_cachedRingAreas2D.size() != pathCount)
				m_cachedRingAreas2D.resize(pathCount);

			MathUtils.KahanSummator totalArea = new MathUtils.KahanSummator(0);
			MathUtils.KahanSummator pathArea = new MathUtils.KahanSummator(0);
			Point2D pt = new Point2D();
			int ipath = 0;
			SegmentIteratorImpl segIter = querySegmentIterator();
			while (segIter.nextPath()) {
				pathArea.reset();
				getXY(getPathStart(segIter.getPathIndex()), pt);// get the area
																// calculation
																// origin to be
																// the origin of
																// the ring.
				while (segIter.hasNextSegment()) {
					pathArea.add(segIter.nextSegment()._calculateArea2DHelper(
							pt.x, pt.y));
				}

				totalArea.add(pathArea.getResult());

				int i = ipath++;
				m_cachedRingAreas2D.write(i, pathArea.getResult());
			}

			m_cachedArea2D = totalArea.getResult();
			_setDirtyFlag(DirtyFlags.DirtyRingAreas2D, false);
		}
	}

	int getOGCPolygonCount() {
		if (!m_bPolygon)
			return 0;

		_updateOGCFlags();

		int polygonCount = 0;
		int partCount = getPathCount();
		for (int ipart = 0; ipart < partCount; ipart++) {
			if (((int) m_pathFlags.read(ipart) & (int) PathFlags.enumOGCStartPolygon) != 0)
				polygonCount++;
		}

		return polygonCount;
	}

	protected void _updateOGCFlags() {
		if (_hasDirtyFlag(DirtyFlags.DirtyOGCFlags)) {
			_updateRingAreas2D();

			int pathCount = getPathCount();
			if (pathCount > 0 && (m_pathFlags == null || m_pathFlags.size() < pathCount))
				m_pathFlags = (AttributeStreamOfInt8) AttributeStreamBase
						.createByteStream(pathCount + 1);

			int firstSign = 1;
			for (int ipath = 0; ipath < pathCount; ipath++) {
				double area = m_cachedRingAreas2D.read(ipath);
				if (ipath == 0)
					firstSign = area > 0 ? 1 : -1;
				if (area * firstSign > 0.0)
					m_pathFlags.setBits(ipath,
							(byte) PathFlags.enumOGCStartPolygon);
				else
					m_pathFlags.clearBits(ipath,
							(byte) PathFlags.enumOGCStartPolygon);
			}
			_setDirtyFlag(DirtyFlags.DirtyOGCFlags, false);
		}
	}

	public int getPathIndexFromPointIndex(int pointIndex) {
		int positionHint = m_currentPathIndex;// in case of multithreading
												// thiswould simply produce an
												// invalid value
		int pathCount = getPathCount();

		// Try using the hint position first to get the path index.
		if (positionHint >= 0 && positionHint < pathCount) {
			if (pointIndex < getPathEnd(positionHint)) {
				if (pointIndex >= getPathStart(positionHint))
					return positionHint;
				positionHint--;
			} else {
				positionHint++;
			}

			if (positionHint >= 0 && positionHint < pathCount) {
				if (pointIndex >= getPathStart(positionHint)
						&& pointIndex < getPathEnd(positionHint)) {
					m_currentPathIndex = positionHint;
					return positionHint;
				}
			}
		}

		if (pathCount < 5) {// TODO: time the performance to choose when to use
							// linear search.
			for (int i = 0; i < pathCount; i++) {
				if (pointIndex < getPathEnd(i)) {
					m_currentPathIndex = i;
					return i;
				}
			}
			throw new GeometryException("corrupted geometry");
		}

		// Do binary search:
		int minPathIndex = 0;
		int maxPathIndex = pathCount - 1;
		while (maxPathIndex > minPathIndex) {
			int mid = minPathIndex + ((maxPathIndex - minPathIndex) >> 1);
			int pathStart = getPathStart(mid);
			if (pointIndex < pathStart)
				maxPathIndex = mid - 1;
			else {
				int pathEnd = getPathEnd(mid);
				if (pointIndex >= pathEnd)
					minPathIndex = mid + 1;
				else {
					m_currentPathIndex = mid;
					return mid;
				}
			}
		}

		m_currentPathIndex = minPathIndex;
		return minPathIndex;
	}

	int getHighestPointIndex(int path_index) {
		assert (path_index >= 0 && path_index < getPathCount());

		AttributeStreamOfDbl position = (AttributeStreamOfDbl) (getAttributeStreamRef(VertexDescription.Semantics.POSITION));
		AttributeStreamOfInt32 paths = (AttributeStreamOfInt32) (getPathStreamRef());

		int path_end = getPathEnd(path_index);
		int path_start = getPathStart(path_index);
		int max_index = -1;
		Point2D max_point = new Point2D(), pt = new Point2D();
		max_point.y = NumberUtils.negativeInf();
		max_point.x = NumberUtils.negativeInf();

		for (int i = path_start + 0; i < path_end; i++) {
			position.read(2 * i, pt);
			if (max_point.compare(pt) == -1) {
				max_index = i;
				max_point.setCoords(pt);
			}
		}

		return max_index;
	}

	/**
	 * Returns total segment count in the MultiPathImpl.
	 */
	public int getSegmentCount() {
		int segCount = getPointCount();
		if (!m_bPolygon) {
			segCount -= getPathCount();
			for (int i = 0, n = getPathCount(); i < n; i++)
				if (isClosedPath(i))
					segCount++;
		}

		return segCount;
	}

	public int getSegmentCount(int path_index) {
		int segCount = getPathSize(path_index);
		if (!isClosedPath(path_index))
			segCount--;
		return segCount;
	}

	// HEADER defintions
	@Override
	public Geometry createInstance() {
		return new MultiPathImpl(m_bPolygon, getDescription());
	}

	@Override
	public int getDimension() {
		return m_bPolygon ? 2 : 1;
	}

	@Override
	public Geometry.Type getType() {
		return m_bPolygon ? Type.Polygon : Type.Polyline;
	}

	/**
	 * Returns True if the class is envelope. THis is not an exact method. Only
	 * addEnvelope makes this true.
	 */
	public boolean isEnvelope() {
		return !_hasDirtyFlag(DirtyFlags.DirtyIsEnvelope);
	}

	/**
	 * Returns a reference to the AttributeStream of MultiPathImpl parts
	 * (Paths).
	 * 
	 * For the non empty MultiPathImpl, that stream contains start points of the
	 * MultiPathImpl curves. In addition, the last element is the total point
	 * count. The number of vertices in a given part is parts[i + 1] - parts[i].
	 */
	public AttributeStreamOfInt32 getPathStreamRef() {
		throwIfEmpty();
		return m_paths;
	}

	/**
	 * sets a reference to an AttributeStream of MultiPathImpl paths (Paths).
	 */
	public void setPathStreamRef(AttributeStreamOfInt32 paths) {
		m_paths = paths;
		notifyModified(DirtyFlags.DirtyAll);
	}

	/**
	 * Returns a reference to the AttributeStream of Segment flags (SegmentFlags
	 * flags). Can be NULL when no non-linear segments are present.
	 * 
	 * Segment flags indicate what kind of segment originates (starts) on the
	 * given point. The last vertices of open Path parts has enumNone flag.
	 */
	public AttributeStreamOfInt8 getSegmentFlagsStreamRef() {
		throwIfEmpty();
		return m_segmentFlags;
	}

	/**
	 * Returns a reference to the AttributeStream of Path flags (PathFlags
	 * flags).
	 * 
	 * Each start point of a path has a flag set to indicate if the Path is open
	 * or closed.
	 */
	public AttributeStreamOfInt8 getPathFlagsStreamRef() {
		throwIfEmpty();
		return m_pathFlags;
	}

	/**
	 * sets a reference to an AttributeStream of Path flags (PathFlags flags).
	 */
	public void setPathFlagsStreamRef(AttributeStreamOfInt8 pathFlags) {
		m_pathFlags = pathFlags;
		notifyModified(DirtyFlags.DirtyAll);
	}

	public AttributeStreamOfInt32 getSegmentIndexStreamRef() {
		throwIfEmpty();
		return m_segmentParamIndex;
	}

	public AttributeStreamOfDbl getSegmentDataStreamRef() {
		throwIfEmpty();
		return m_segmentParams;
	}

	public int getPathCount() {
		return (m_paths != null) ? m_paths.size() - 1 : 0;
	}

	public int getPathEnd(int partIndex) {
		return m_paths.read(partIndex + 1);
	}

	public int getPathSize(int partIndex) {
		return m_paths.read(partIndex + 1) - m_paths.read(partIndex);
	}

	public int getPathStart(int partIndex) {
		return m_paths.read(partIndex);
	}

	@Override
	public Object _getImpl() {
		return this;
	}

	public void setDirtyOGCFlags(boolean bYesNo) {
		_setDirtyFlag(DirtyFlags.DirtyOGCFlags, bYesNo);
	}

	public boolean hasDirtyOGCStartFlags() {
		return _hasDirtyFlag(DirtyFlags.DirtyOGCFlags);
	}

	public void setDirtyRingAreas2D(boolean bYesNo) {
		_setDirtyFlag(DirtyFlags.DirtyRingAreas2D, bYesNo);
	}

	public boolean hasDirtyRingAreas2D() {
		return _hasDirtyFlag(DirtyFlags.DirtyRingAreas2D);
	}

	public void setRingAreasStreamRef(AttributeStreamOfDbl ringAreas) {
		m_cachedRingAreas2D = ringAreas;
		_setDirtyFlag(DirtyFlags.DirtyRingAreas2D, false);
	}

	// HEADER defintions

	// // TODO check this against current implementation in native
	// public void notifyModified(int flags)
	// {
	// if(flags == DirtyFlags.DirtyAll)
	// {
	// m_reservedPointCount = -1;
	// _notifyModifiedAllImpl();
	// }
	// m_flagsMask |= flags;
	// _clearAccelerators();
	//
	//
	// // ROHIT's implementation
	// // if (m_paths == null || 0 == m_paths.size())
	// // m_pointCount = 0;
	// // else
	// // m_pointCount = m_paths.read(m_paths.size() - 1);
	// //
	// // super.notifyModified(flags);
	// }

	@Override
	public boolean _buildRasterizedGeometryAccelerator(double toleranceXY,
			GeometryAccelerationDegree accelDegree) {
		if (m_accelerators == null)// (!m_accelerators)
		{
			m_accelerators = new GeometryAccelerators();
		}

		int rasterSize = RasterizedGeometry2D
				.rasterSizeFromAccelerationDegree(accelDegree);
		RasterizedGeometry2D rgeom = m_accelerators.getRasterizedGeometry();
		if (rgeom != null) {
			if (rgeom.getToleranceXY() < toleranceXY
					|| rasterSize > rgeom.getRasterSize()) {
				m_accelerators._setRasterizedGeometry(null);
			} else
				return true;
		}

		rgeom = RasterizedGeometry2D.create(this, toleranceXY, rasterSize);
		m_accelerators._setRasterizedGeometry(rgeom);
		//rgeom.dbgSaveToBitmap("c:/temp/ddd.bmp");
		return true;
	}

	@Override
	public int hashCode() {
		int hashCode = super.hashCode();

		if (!isEmptyImpl()) {
			int pathCount = getPathCount();

			if (m_paths != null)
				m_paths.calculateHashImpl(hashCode, 0, pathCount + 1);

			if (m_pathFlags != null)
				m_pathFlags.calculateHashImpl(hashCode, 0, pathCount);
		}

		return hashCode;
	}

	public byte getSegmentFlags(int ivertex) {
		if (m_segmentFlags != null)
			return m_segmentFlags.read(ivertex);
		else
			return (byte) SegmentFlags.enumLineSeg;
	}

	public void getSegment(int startVertexIndex, SegmentBuffer segBuffer,
			boolean bStripAttributes) {
		int ipath = getPathIndexFromPointIndex(startVertexIndex);
		if (startVertexIndex == getPathEnd(ipath) - 1 && !isClosedPath(ipath))
			throw new GeometryException("index out of bounds");

		_verifyAllStreams();
		AttributeStreamOfInt8 segFlagStream = getSegmentFlagsStreamRef();
		int segFlag = SegmentFlags.enumLineSeg;
		if (segFlagStream != null)
			segFlag = segFlagStream.read(startVertexIndex)
					& SegmentFlags.enumSegmentMask;

		switch (segFlag) {
		case SegmentFlags.enumLineSeg:
			segBuffer.createLine();
			break;
		case SegmentFlags.enumBezierSeg:
			throw GeometryException.GeometryInternalError();
		case SegmentFlags.enumArcSeg:
			throw GeometryException.GeometryInternalError();
		default:
			throw GeometryException.GeometryInternalError();
		}

		Segment currentSegment = segBuffer.get();
		if (!bStripAttributes)
			currentSegment.assignVertexDescription(m_description);
		else
			currentSegment
					.assignVertexDescription(VertexDescriptionDesignerImpl
							.getDefaultDescriptor2D());

		int endVertexIndex;
		if (startVertexIndex == getPathEnd(ipath) - 1 && isClosedPath(ipath)) {
			endVertexIndex = getPathStart(ipath);
		} else
			endVertexIndex = startVertexIndex + 1;

		Point2D pt = new Point2D();
		getXY(startVertexIndex, pt);
		currentSegment.setStartXY(pt);
		getXY(endVertexIndex, pt);
		currentSegment.setEndXY(pt);

		if (!bStripAttributes) {
			for (int i = 1, nattr = m_description.getAttributeCount(); i < nattr; i++) {
				int semantics = m_description._getSemanticsImpl(i);
				int ncomp = VertexDescription.getComponentCount(semantics);
				for (int ord = 0; ord < ncomp; ord++) {
					double vs = getAttributeAsDbl(semantics, startVertexIndex,
							ord);
					currentSegment.setStartAttribute(semantics, ord, vs);
					double ve = getAttributeAsDbl(semantics, endVertexIndex,
							ord);
					currentSegment.setEndAttribute(semantics, ord, ve);
				}
			}
		}
	}

	void queryPathEnvelope2D(int path_index, Envelope2D envelope) {
		if (path_index >= getPathCount())
			throw new IllegalArgumentException();

		if (isEmpty()) {
			envelope.setEmpty();
			return;
		}

		if (hasNonLinearSegments(path_index)) {
			throw new GeometryException("not implemented");
		} else {
			AttributeStreamOfDbl stream = (AttributeStreamOfDbl) getAttributeStreamRef(VertexDescription.Semantics.POSITION);
			Point2D pt = new Point2D();
			Envelope2D env = new Envelope2D();
			env.setEmpty();
			for (int i = getPathStart(path_index), iend = getPathEnd(path_index); i < iend; i++) {
				stream.read(2 * i, pt);
				env.merge(pt);
			}
			envelope.setCoords(env);
		}
	}

	public void queryLoosePathEnvelope2D(int path_index, Envelope2D envelope) {
		if (path_index >= getPathCount())
			throw new IllegalArgumentException();

		if (isEmpty()) {
			envelope.setEmpty();
			return;
		}

		if (hasNonLinearSegments(path_index)) {
			throw new GeometryException("not implemented");
		} else {
			AttributeStreamOfDbl stream = (AttributeStreamOfDbl) getAttributeStreamRef(VertexDescription.Semantics.POSITION);
			Point2D pt = new Point2D();
			Envelope2D env = new Envelope2D();
			env.setEmpty();
			for (int i = getPathStart(path_index), iend = getPathEnd(path_index); i < iend; i++) {
				stream.read(2 * i, pt);
				env.merge(pt);
			}
			envelope.setCoords(env);
		}
	}
	
	@Override
	public boolean _buildQuadTreeAccelerator(GeometryAccelerationDegree d) {
		if (m_accelerators == null)// (!m_accelerators)
		{
			m_accelerators = new GeometryAccelerators();
		}

		if (d == GeometryAccelerationDegree.enumMild || getPointCount() < 16)
			return false;

		QuadTreeImpl quad_tree_impl = InternalUtils.buildQuadTree(this);
		m_accelerators._setQuadTree(quad_tree_impl);

		return true;
	}

	boolean _buildQuadTreeForPathsAccelerator(GeometryAccelerationDegree degree) {
		if (m_accelerators == null) {
			m_accelerators = new GeometryAccelerators();
		}

		// TODO: when less than two envelopes - no need to this.

		if (m_accelerators.getQuadTreeForPaths() != null)
			return true;

		m_accelerators._setQuadTreeForPaths(null);
		QuadTreeImpl quad_tree_impl = InternalUtils.buildQuadTreeForPaths(this);
		m_accelerators._setQuadTreeForPaths(quad_tree_impl);

		return true;
	}

	void setFillRule(int rule) {
		assert (m_bPolygon);
		m_fill_rule = rule;
	}

	int getFillRule() {
		return m_fill_rule;
	}

	void clearDirtyOGCFlags() { 
		_setDirtyFlag(DirtyFlags.DirtyOGCFlags, false);
	}
}
871280
def from_path(kls, vertices):
    
    if vertices.shape[0] == 0:
      return PrecomputedSkeleton()

    skel = PrecomputedSkeleton(vertices)
    edges = np.zeros(shape=(skel.vertices.shape[0] - 1, 2), dtype=np.uint32)
    edges[:,0] = np.arange(skel.vertices.shape[0] - 1)
    edges[:,1] = np.arange(1, skel.vertices.shape[0])
    skel.edges = edges
    return skel
549536
protected function interpolate($x, $y, $nw, $ne, $sw, $se)
    {
        list($r0, $g0, $b0) = $this->getRGB($nw);
        list($r1, $g1, $b1) = $this->getRGB($ne);
        list($r2, $g2, $b2) = $this->getRGB($sw);
        list($r3, $g3, $b3) = $this->getRGB($se);

        $cx = 1.0 - $x;
        $cy = 1.0 - $y;

        $m0 = $cx * $r0 + $x * $r1;
        $m1 = $cx * $r2 + $x * $r3;
        $r  = (int) ($cy * $m0 + $y * $m1);

        $m0 = $cx * $g0 + $x * $g1;
        $m1 = $cx * $g2 + $x * $g3;
        $g  = (int) ($cy * $m0 + $y * $m1);

        $m0 = $cx * $b0 + $x * $b1;
        $m1 = $cx * $b2 + $x * $b3;
        $b  = (int) ($cy * $m0 + $y * $m1);

        return ($r << 16) | ($g << 8) | $b;
    }
861319
def xyz2lonlat(x,y,z):
    

    xs = np.array(x)
    ys = np.array(y)
    zs = np.array(z)

    lons = np.arctan2(ys, xs)
    lats = np.arcsin(zs)

    return lons, lats