Foam::point Foam::CatmullRomSpline::position
(
    const label segment,
    const scalar mu
) const
{
    // out-of-bounds
    if (segment < 0)
    {
        return points().first();
    }
    else if (segment > nSegments())
    {
        return points().last();
    }

    const point& p0 = points()[segment];
    const point& p1 = points()[segment+1];

    // special cases - no calculation needed
    if (mu <= 0.0)
    {
        return p0;
    }
    else if (mu >= 1.0)
    {
        return p1;
    }


    // determine the end points
    point e0;
    point e1;

    if (segment == 0)
    {
        // end: simple reflection
        e0 = 2*p0 - p1;
    }
    else
    {
        e0 = points()[segment-1];
    }

    if (segment+1 == nSegments())
    {
        // end: simple reflection
        e1 = 2*p1 - p0;
    }
    else
    {
        e1 = points()[segment+2];
    }


    return 0.5 *
    (
        ( 2*p0 )
      + mu *
        (
            ( -e0 + p1 )
          + mu *
            (
                ( 2*e0 - 5*p0 + 4*p1 - e1 )
              + mu *
                ( -e0 + 3*p0 - 3*p1 + e1 )
            )
        )
    );
}

// virtual
bool MayaNurbsCurveWriter::writeNurbsCurveAttrs(const UsdTimeCode &usdTime, UsdGeomNurbsCurves &primSchema)
{
    MStatus status = MS::kSuccess;

    // Write parent class attrs
    writeTransformAttrs(usdTime, primSchema);

    // Return if usdTime does not match if shape is animated
    if (usdTime.IsDefault() == isShapeAnimated() ) {
        // skip shape as the usdTime does not match if shape isAnimated value
        return true; 
    }

    MFnDependencyNode fnDepNode(getDagPath().node(), &status);
    MString name = fnDepNode.name();

    MFnNurbsCurve curveFn( getDagPath(), &status );
    if (!status) {
        MGlobal::displayError("MFnNurbsCurve() failed for MayaNurbsCurveWriter");
        return false;
    }

    // Get curve attrs ======
    unsigned int numCurves = 1; // Assuming only 1 curve for now
    VtArray<int> curveOrder(numCurves);
    VtArray<int> curveVertexCounts(numCurves);
    VtArray<float> curveWidths(numCurves);
    VtArray<GfVec2d> ranges(numCurves);

    curveOrder[0] = curveFn.degree()+1;
    curveVertexCounts[0] = curveFn.numCVs();
    TF_AXIOM(curveOrder[0] <= curveVertexCounts[0] );
    curveWidths[0] = 1.0; // TODO: Retrieve from custom attr

    double mayaKnotDomainMin;
    double mayaKnotDomainMax;
    status = curveFn.getKnotDomain(mayaKnotDomainMin, mayaKnotDomainMax);
    TF_AXIOM(status == MS::kSuccess);
    ranges[0][0] = mayaKnotDomainMin;
    ranges[0][1] = mayaKnotDomainMax;

    MPointArray mayaCurveCVs;
    status = curveFn.getCVs(mayaCurveCVs, MSpace::kObject);
    TF_AXIOM(status == MS::kSuccess);
    VtArray<GfVec3f> points(mayaCurveCVs.length()); // all CVs batched together
    for (unsigned int i=0; i < mayaCurveCVs.length(); i++) {
        points[i].Set(mayaCurveCVs[i].x, mayaCurveCVs[i].y, mayaCurveCVs[i].z);
    }

    MDoubleArray mayaCurveKnots;
    status = curveFn.getKnots(mayaCurveKnots);
    TF_AXIOM(status == MS::kSuccess);
    VtArray<double> curveKnots(mayaCurveKnots.length()); // all knots batched together
    for (unsigned int i=0; i < mayaCurveKnots.length(); i++) {
        curveKnots[i] = mayaCurveKnots[i];
    }

    // Gprim
    VtArray<GfVec3f> extent(2);
    UsdGeomCurves::ComputeExtent(points, curveWidths, &extent);
    primSchema.CreateExtentAttr().Set(extent, usdTime);

    // Curve
    primSchema.GetOrderAttr().Set(curveOrder);   // not animatable
    primSchema.GetCurveVertexCountsAttr().Set(curveVertexCounts); // not animatable
    primSchema.GetWidthsAttr().Set(curveWidths); // not animatable
    primSchema.GetKnotsAttr().Set(curveKnots);   // not animatable
    primSchema.GetRangesAttr().Set(ranges); // not animatable
    primSchema.GetPointsAttr().Set(points, usdTime); // CVs

    // TODO: Handle periodic and non-periodic cases

    return true;
}

void mexFunction(int nlhs, mxArray  *plhs[], int nrhs, const mxArray  *prhs[])
{	
	int ref_or_direct = 1;			// ref_or_direct = 1 means interpret second input argument as reference indices
									// = 0 means interpret second input argument as reference vectors
	
	long past = 0;
	double epsilon = 0;
	
	/* check input args */
	
	if (nrhs < 4)
	{
		mexErrMsgTxt("Fast nearest neighbour searcher : Data set of points (row vectors), preprocessing data, reference indices or reference points \nand number of nearest neighbours must be given, epsilon (default = 0) is optional");
		return;
	}
	
	/* handle matrix I/O */
	
	const long N 		= mxGetM(prhs[0]);	
	const long dim  	= mxGetN(prhs[0]);
	const double* p 	= (double *)mxGetPr(prhs[0]);
	
	double* ref 	= (double *)mxGetPr(prhs[2]);
	long R; 										// number of query (reference) points
	
	const long NNR 	= (long) *((double *)mxGetPr(prhs[3]));
	
	if (N < 1) {
		mexErrMsgTxt("Data set must consist of at least two points (row vectors)");
		return;
	}		
	if (dim < 1) {
		mexErrMsgTxt("Data points must be at least of dimension one");
		return;
	}	
	if (NNR<1) {
		mexErrMsgTxt("At least one nearest neighbour must be requested");
		return;
	}	
	if ((mxGetN(prhs[2]) == 0) || (mxGetN(prhs[2]) == 0)) {
		mexErrMsgTxt("Wrong reference indices or reference points given");
		return;
	}	

	
	if (mxGetN(prhs[2]) == 1) {
		R = mxGetM(prhs[2]);
		ref_or_direct = 1;	
	} 
	else if ((mxGetM(prhs[2]) == 1) && (mxGetN(prhs[2]) != dim)) {
		R = mxGetN(prhs[2]);
		ref_or_direct = 1;	
	}
	else if (mxGetN(prhs[2]) == dim) {
		R = mxGetM(prhs[2]);
		ref_or_direct = 0;	
	} else  {
		mexErrMsgTxt("Cannot determine if second argument are reference indices or reference points");
		return;
	}	
	
/*******************************/
/*** Quick&Dirty patch, to manage onedimensional "direct" query-points. I. Wedekind, 16.9.2002 */
	if((mxGetN(prhs[2])==dim) && (ref_or_direct==1) && (nrhs<5))
	    {
		ref_or_direct=0;
	    }
/*******************************/

	if (R < 1) {
		mexErrMsgTxt("At least one reference index or point must be given");
		return;
	}	

	if (ref_or_direct) {		// interpret second argument as list of indices into the data set given as first argument 
		if (nrhs < 5)
		{
			mexErrMsgTxt("Fast nearest neighbour searcher : Data set of points (row vectors), preprocessing data, reference indices,\nnumber of nearest neighbours and past must be given");
			return;
		}
		past	= (long) *((double *)mxGetPr(prhs[4]));
		for (long i=0; i < R; i++) {
			if ((ref[i] < 1) || (ref[i]>N)) {
				mexErrMsgTxt("Reference indices out of range");
				return;
			}	
		}
		if ((N - (2*past)-1) < NNR)
		{
			mexErrMsgTxt("Fast nearest neighbour searcher : To many neighbors for each query point are requested");
			return;
		}
		if (nrhs > 5) 
			epsilon = (double) *((double *)mxGetPr(prhs[5]));		// support approximative queries		
	}  else {
		if (nrhs > 4) 
			epsilon = (double) *((double *)mxGetPr(prhs[4]));		// support approximative queries	
	}
	
	plhs[0] = mxCreateDoubleMatrix(R, NNR, mxREAL);
//.........这里部分代码省略.........

void TestExportHighMesh(const char *filename, int frame_count)
{
    auto *ctx = usdiCreateContext();

    usdi::ExportSettings settings;
    settings.instanceable_by_default = true;
    usdiSetExportSettings(ctx, &settings);

    usdiCreateStage(ctx, filename);
    auto *root = usdiGetRoot(ctx);

    auto *xf = usdiCreateXform(ctx, root, "WaveMeshRoot");
    usdiPrimSetInstanceable(xf, true);
    {
        usdi::XformData data;
        usdiXformWriteSample(xf, &data);
    }

    auto *mesh = usdiCreateMesh(ctx, xf, "WaveMesh");
    {
        std::vector<std::vector<int>> counts(frame_count);
        std::vector<std::vector<int>> indices(frame_count);
        std::vector<std::vector<float3>> points(frame_count);
        std::vector<std::vector<float2>> uv(frame_count);
        std::vector<std::vector<float4>> colors(frame_count);

        usdi::Time t = 0.0;

        tbb::parallel_for(0, frame_count, [frame_count, &counts, &indices, &points, &uv, &colors](int i) {
            usdi::Time t = i;
            int resolution = 8;
            if (i < 30)      { resolution = 8; }
            else if (i < 60) { resolution = 16; }
            else if (i < 90) { resolution = 32; }
            else if (i < 120) { resolution = 64; }
            else if (i < 150) { resolution = 128; }
            else { resolution = 256; }
            GenerateWaveMesh(counts[i], indices[i], points[i], uv[i], 1.0f, 0.5f, resolution, (360.0 * 5 * DegToRad / frame_count) * i);

            auto& color = colors[i];
            color.resize(points[i].size());
            for (size_t ci = 0; ci < color.size(); ++ci) {
                color[ci] = {
                    std::sin((float)i * DegToRad * 2.1f) * 0.5f + 0.5f,
                    std::cos((float)i * DegToRad * 6.8f) * 0.5f + 0.5f,
                    std::sin((float)i * DegToRad * 11.7f) * 0.5f + 0.5f,
                    1.0f
                };
            }
        });
        for (int i = 0; i < frame_count; ++i) {
            auto& vertices = points[i];

            usdi::Time t = i;
            usdi::MeshData data;
            data.counts = counts[i].data();
            data.num_counts = counts[i].size();
            data.indices = indices[i].data();
            data.num_indices = indices[i].size();
            data.points = points[i].data();
            data.num_points = points[i].size();
            data.uvs = uv[i].data();
            data.colors = colors[i].data();
            usdiMeshWriteSample(mesh, &data, t);
        }

        usdiMeshPreComputeNormals(mesh, false);
        //usdiMeshPreComputeNormals(mesh, true);
    }

    {
        auto *ref1 = usdiCreateXform(ctx, root, "WaveMeshRef1");
        usdi::XformData data;
        data.position.x = 1.5f;
        usdiXformWriteSample(ref1, &data);

        auto *ref = usdiCreateOverride(ctx, "/WaveMeshRef1/Ref");
        usdiPrimAddReference(ref, nullptr, "/WaveMeshRoot");
    }
    {
        auto *ref2 = usdiCreateXform(ctx, root, "WaveMeshRef2");
        usdi::XformData data;
        data.position.x = -1.5f;
        usdiXformWriteSample(ref2, &data);

        auto *ref = usdiCreateOverride(ctx, "/WaveMeshRef2/Ref");
        usdiPrimAddReference(ref, nullptr, "/WaveMeshRoot");
    }

    usdiSave(ctx);
    usdiDestroyContext(ctx);
}

// virtual
bool MayaMeshWriter::writeMeshAttrs(const UsdTimeCode &usdTime, UsdGeomMesh &primSchema)
{

    MStatus status = MS::kSuccess;

    // Write parent class attrs
    writeTransformAttrs(usdTime, primSchema);

    // Return if usdTime does not match if shape is animated
    if (usdTime.IsDefault() == isShapeAnimated() ) {
        // skip shape as the usdTime does not match if shape isAnimated value
        return true; 
    }

    MFnMesh lMesh( getDagPath(), &status );
    if ( !status )
    {
        MGlobal::displayError( "MFnMesh() failed for MayaMeshWriter" );
        return false;
    }
    unsigned int numVertices = lMesh.numVertices();
    unsigned int numPolygons = lMesh.numPolygons();

    // Set mesh attrs ==========
    // Get points
    // TODO: Use memcpy()
    const float* mayaRawPoints = lMesh.getRawPoints(&status);
    VtArray<GfVec3f> points(numVertices);
    for (unsigned int i = 0; i < numVertices; i++) {
        unsigned int floatIndex = i*3;
        points[i].Set(mayaRawPoints[floatIndex],
                      mayaRawPoints[floatIndex+1],
                      mayaRawPoints[floatIndex+2]);
    }
    primSchema.GetPointsAttr().Set(points, usdTime); // ANIMATED

    // Compute the extent using the raw points
    VtArray<GfVec3f> extent(2);
    UsdGeomPointBased::ComputeExtent(points, &extent);
    primSchema.CreateExtentAttr().Set(extent, usdTime);

    // Get faceVertexIndices
    unsigned int numFaceVertices = lMesh.numFaceVertices(&status);
    VtArray<int>     faceVertexCounts(numPolygons);
    VtArray<int>     faceVertexIndices(numFaceVertices);
    MIntArray mayaFaceVertexIndices; // used in loop below
    unsigned int curFaceVertexIndex = 0;
    for (unsigned int i = 0; i < numPolygons; i++) {
        lMesh.getPolygonVertices(i, mayaFaceVertexIndices);
        faceVertexCounts[i] = mayaFaceVertexIndices.length();
        for (unsigned int j=0; j < mayaFaceVertexIndices.length(); j++) {
            faceVertexIndices[ curFaceVertexIndex ] = mayaFaceVertexIndices[j]; // push_back
            curFaceVertexIndex++;
        }
    }
    primSchema.GetFaceVertexCountsAttr().Set(faceVertexCounts);   // not animatable
    primSchema.GetFaceVertexIndicesAttr().Set(faceVertexIndices); // not animatable

    // Read usdSdScheme attribute. If not set, we default to defaultMeshScheme
    // flag that can be user defined and initialized to catmullClark
    TfToken sdScheme = PxrUsdMayaMeshUtil::getSubdivScheme(lMesh, getArgs().defaultMeshScheme);    
    primSchema.CreateSubdivisionSchemeAttr(VtValue(sdScheme), true);

    // Polygonal Mesh Case
    if (sdScheme==UsdGeomTokens->none) {
        // Support for standard USD bool and with Mojito bool tag
        TfToken normalInterp=PxrUsdMayaMeshUtil::getEmitNormals(lMesh, UsdGeomTokens->none);
        
        if (normalInterp==UsdGeomTokens->faceVarying) {
            // Get References to members of meshData object
            MFloatVectorArray normalArray;
            MFloatVectorArray vertexNormalArray;
 
            lMesh.getNormals(normalArray, MSpace::kObject);

            // Iterate through each face in the mesh.
            vertexNormalArray.setLength(lMesh.numFaceVertices());
            VtArray<GfVec3f> meshNormals(lMesh.numFaceVertices());
            size_t faceVertIdx = 0;
            for (MItMeshPolygon faceIter(getDagPath()); !faceIter.isDone(); faceIter.next()) {
                // Iterate through each face-vertex.
                for (size_t locVertIdx = 0; locVertIdx < faceIter.polygonVertexCount();
                        ++locVertIdx, ++faceVertIdx) {
                    int index=faceIter.normalIndex(locVertIdx);
                    for (int j=0;j<3;j++) {
                        meshNormals[faceVertIdx][j]=normalArray[index][j];
                    }
                }
            }
            primSchema.GetNormalsAttr().Set(meshNormals, usdTime);
            primSchema.SetNormalsInterpolation(normalInterp);
        }
    } else {
        TfToken sdInterpBound = PxrUsdMayaMeshUtil::getSubdivInterpBoundary(
            lMesh, UsdGeomTokens->edgeAndCorner);

        primSchema.CreateInterpolateBoundaryAttr(VtValue(sdInterpBound), true);
        
        TfToken sdFVInterpBound = PxrUsdMayaMeshUtil::getSubdivFVInterpBoundary(
//.........这里部分代码省略.........

void QgsMapToolSplitParts::cadCanvasReleaseEvent( QgsMapMouseEvent *e )
{
  //check if we operate on a vector layer
  QgsVectorLayer *vlayer = qobject_cast<QgsVectorLayer *>( mCanvas->currentLayer() );

  if ( !vlayer )
  {
    notifyNotVectorLayer();
    return;
  }

  if ( !vlayer->isEditable() )
  {
    notifyNotEditableLayer();
    return;
  }

  bool split = false;

  //add point to list and to rubber band
  if ( e->button() == Qt::LeftButton )
  {
    //If we snap the first point on a vertex of a line layer, we directly split the feature at this point
    if ( vlayer->geometryType() == QgsWkbTypes::LineGeometry && points().isEmpty() )
    {
      QgsPointLocator::Match m = mCanvas->snappingUtils()->snapToCurrentLayer( e->pos(), QgsPointLocator::Vertex );
      if ( m.isValid() )
      {
        split = true;
      }
    }

    int error = addVertex( e->mapPoint() );
    if ( error == 1 )
    {
      //current layer is not a vector layer
      return;
    }
    else if ( error == 2 )
    {
      //problem with coordinate transformation
      QgisApp::instance()->messageBar()->pushMessage(
        tr( "Coordinate transform error" ),
        tr( "Cannot transform the point to the layers coordinate system" ),
        QgsMessageBar::INFO,
        QgisApp::instance()->messageTimeout() );
      return;
    }

    startCapturing();
  }
  else if ( e->button() == Qt::RightButton )
  {
    split = true;
  }
  if ( split )
  {
    deleteTempRubberBand();

    //bring up dialog if a split was not possible (polygon) or only done once (line)
    bool topologicalEditing = QgsProject::instance()->topologicalEditing();
    vlayer->beginEditCommand( tr( "Parts split" ) );
    QgsGeometry::OperationResult returnCode = vlayer->splitParts( points(), topologicalEditing );
    vlayer->endEditCommand();
    if ( returnCode == QgsGeometry::OperationResult::NothingHappened )
    {
      QgisApp::instance()->messageBar()->pushMessage(
        tr( "No parts were split" ),
        tr( "If there are selected parts, the split tool only applies to those. If you would like to split all parts under the split line, clear the selection." ),
        QgsMessageBar::WARNING,
        QgisApp::instance()->messageTimeout() );
    }
    else if ( returnCode == QgsGeometry::OperationResult::GeometryEngineError )
    {
      QgisApp::instance()->messageBar()->pushMessage(
        tr( "No part split done" ),
        tr( "Cut edges detected. Make sure the line splits parts into multiple parts." ),
        QgsMessageBar::WARNING,
        QgisApp::instance()->messageTimeout() );
    }
    else if ( returnCode == QgsGeometry::OperationResult::InvalidBaseGeometry )
    {
      QgisApp::instance()->messageBar()->pushMessage(
        tr( "No part split done" ),
        tr( "The geometry is invalid. Please repair before trying to split it." ),
        QgsMessageBar::WARNING,
        QgisApp::instance()->messageTimeout() );
    }
    else if ( returnCode != QgsGeometry::OperationResult::Success )
    {
      //several intersections but only one split (most likely line)
      QgisApp::instance()->messageBar()->pushMessage(
        tr( "Split error" ),
        tr( "An error occurred during splitting." ),
        QgsMessageBar::WARNING,
        QgisApp::instance()->messageTimeout() );
    }

    stopCapturing();
  }
//.........这里部分代码省略.........

int main() {

    const unsigned oversample = 2;
    const unsigned w = 2048*oversample, h = 2048*oversample;
    cout << "allocating memory for a " << w << " x " << h << " image buffer..." << flush;
    simg img( w, h, {0, 0, 0, 255} );
    cout << " done." << endl;

    // list of points
    const int numpoints = 5;
    uniform_int_distribution<> dist_int( 0, numpoints-1 );
    vector< double[2] > points( numpoints );

    // radius: 0.5 will just touch the edges of the image
    const double rad = 0.99 * 0.5;
    // initial rotation: not important but makes each image unique
    const double rot = rand01();

    // distribute points evenly around a circle (not necessary, but makes it easy to compare)
    for( int p = 0; p < numpoints; p++ ) {
        double angle = (rot + ((double)p / numpoints)) * 2. * 3.14159265358979323846;
        points[p][0] = 0.5 + rad*cos(angle);
        points[p][1] = 0.5 + rad*sin(angle);
    }

    // additive RGBA (may be negative, usually want to leave alpha = 0)
    const srgba_light light = {10, 10, 10, 0};

    // a touchy constant: too big and the render appears fuzzy
    const double beta = 1. / 2.;


    // note: may need to try many different starting points to get a good image
    // * note: so far has worked with a single starting point
    double point[2] = { rand01(), rand01() };

    // how many times to plot the point: too big = slow
    const unsigned long numplots = 100000000;

    // do a render
    cout << "rendering " << numplots << " points..." << flush;
    for( unsigned long i = 0; i < numplots; i++ ) {
        
        // pick a point to move toward
        int p = dist_int( rand_gen );

        // compute new coordinates
        point[0] = points[p][0] + beta * ( point[0] - points[p][0] );
        point[1] = points[p][1] + beta * ( point[1] - points[p][1] );

        // plot point
        if( point[0] >= 0. && point[0] < 1. && point[1] >= 0. && point[1] < 1. ) {
            unsigned x = floor(point[0] * w);
            unsigned y = floor(point[1] * h);
            img.add( x, y, light );
        }
    }
    cout << " done." << endl;

    char filename[1024];
    sprintf( filename, "%d-points_%ux%u_%lu-samples.png", numpoints, w, h, numplots );

    img.save( string(filename) );

    return 0;
}

    int
    process(const tendrils& in, const tendrils& out)
    {
      cv::Mat depth = in.get<cv::Mat>("depth");
      if (depth.empty())
        return ecto::OK;

      cv::Mat R, T, K;
      in.get<cv::Mat>("R").convertTo(R, CV_64F);
      in.get<cv::Mat>("T").convertTo(T, CV_64F);
      in.get<cv::Mat>("K").convertTo(K, CV_64F);

      if (R.empty() || T.empty() || K.empty())
        return ecto::OK;

      cv::Mat mask = cv::Mat::zeros(depth.size(), CV_8UC1);

      box_mask.create(depth.size());
      box_mask.setTo(cv::Scalar(0));

      std::vector<cv::Point3f> box(8);
      box[0] = cv::Point3f(*x_crop, *y_crop, *z_min);
      box[1] = cv::Point3f(-*x_crop, *y_crop, *z_min);
      box[2] = cv::Point3f(-*x_crop, -*y_crop, *z_min);
      box[3] = cv::Point3f(*x_crop, -*y_crop, *z_min);
      box[4] = cv::Point3f(*x_crop, *y_crop, *z_crop);
      box[5] = cv::Point3f(-*x_crop, *y_crop, *z_crop);
      box[6] = cv::Point3f(-*x_crop, -*y_crop, *z_crop);
      box[7] = cv::Point3f(*x_crop, -*y_crop, *z_crop);

      std::vector<cv::Point2f> projected, hull;
      cv::projectPoints(box, R, T, K, cv::Mat(4, 1, CV_64FC1, cv::Scalar(0)), projected);

      cv::convexHull(projected, hull, true);
      std::vector<cv::Point> points(hull.size());
      std::copy(hull.begin(), hull.end(), points.begin());
      cv::fillConvexPoly(box_mask, points.data(), points.size(), cv::Scalar::all(255));

      int width = mask.size().width;
      int height = mask.size().height;
      cv::Mat_<uint8_t>::iterator it = mask.begin<uint8_t>();
      cv::Mat_<uint8_t>::iterator mit = box_mask.begin();

      cv::Mat_<cv::Vec3f> points3d;
      depthTo3dMask(K, depth, box_mask, points3d);
      if (points3d.empty())
        return ecto::OK;

      cv::Matx<double, 3, 1> p, p_r, Tx(T); //Translation
      cv::Matx<double, 3, 3> Rx; //inverse Rotation
      Rx = cv::Mat(R.t());
      cv::Mat_<cv::Vec3f>::iterator point = points3d.begin();

//      std::cout << cv::Mat(Rx) << "\n" << cv::Mat(Tx) << std::endl;
//      std::cout << fx << " " << fy << " " << cx << " " << cy << " " << std::endl;
      double z_min_ = *z_min, z_max_ = *z_crop, x_min_ = -*x_crop, x_max_ = *x_crop, y_min_ = -*y_crop,
          y_max_ = *y_crop;
      for (int v = 0; v < height; v++)
      {
        for (int u = 0; u < width; u++, ++it, ++mit)
        {
          if (*mit == 0)
            continue;
          //calculate the point based on the depth
          p(0) = (*point).val[0];
          p(1) = (*point).val[1];
          p(2) = (*point).val[2];
          ++point;
          p_r = Rx * (p - Tx);
//          std::cout <<"p=" << cv::Mat(p) << ",p_r="<< cv::Mat(p_r) << std::endl;
          if (p_r(2) > z_min_ && p_r(2) < z_max_ && p_r(0) > x_min_ && p_r(0) < x_max_ && p_r(1) > y_min_
              && p_r(1) < y_max_)
            *it = 255;
        }
      }
      out["mask"] << mask;
      return ecto::OK;
    }

ON_BOOL32 ON_Light::GetBBox( // returns true if successful
       double* boxmin,    // boxmin[dim]
       double* boxmax,    // boxmax[dim]
       ON_BOOL32 bGrowBox
       ) const
{
  bool rc = true;
  ON_3dPointArray points(16);

  switch(m_style)
  {
  case ON::camera_directional_light:
  case ON::world_directional_light:
    points.Append(m_location);
    points.Append(m_location+m_direction);
    break;

  case ON::camera_point_light:
  case ON::world_point_light:
    points.Append(m_location);
    break;

  case ON::camera_spot_light:
  case ON::world_spot_light:
    if ( m_spot_angle > 0.0 && m_spot_angle < 90.0 )
    {
      double r = m_direction.Length()*tan(ON_PI*m_spot_angle/180.0);
      ON_Circle c(ON_Plane(m_location+m_direction,m_direction),r);
      ON_BoundingBox cbox = c.BoundingBox();
      cbox.GetCorners( points );
    }
    else
    {
      points.Append(m_location+m_direction);
    }
    points.Append(m_location);
    break;

  case ON::ambient_light:
    points.Append(m_location);
    rc = false;
    break;
  
  case ON::world_linear_light:
    points.Append(m_location);
    points.Append(m_location+m_length);
    break;

  case ON::world_rectangular_light:
    points.Append(m_location);
    points.Append(m_location+m_length);
    points.Append(m_location+m_width);
    points.Append(m_location+m_width+m_length);
    {
      // include target and direction marker to avoid display clipping
      ON_3dPoint center(m_location+(m_width+m_length)*0.5);
      points.Append(center+m_direction);
      ON_3dVector marker(m_direction); 
      marker.Unitize();
      marker *= (m_width+m_length).Length()/12.0; // from GetRectangularLightSegments
      points.Append(center+marker);
    }
    break;

  default:
    rc = false;
    break;
  }

  if ( rc && points.Count() > 0 )
  {
     rc = ON_GetPointListBoundingBox( 3, 0, points.Count(), 3, 
                                      (double*)points.Array(), 
                                      boxmin, boxmax, 
                                      bGrowBox?true:false )
        ? true 
        : false;
  }

  return rc;
}

int Delta::execute()
{

    Options sourceOptions;
    {
        sourceOptions.add<std::string>("filename", m_sourceFile);
        sourceOptions.add<bool>("debug", isDebug());
        sourceOptions.add<boost::uint32_t>("verbose", getVerboseLevel());
    }
    Stage* source = AppSupport::makeReader(sourceOptions);
    source->initialize();
    
    boost::uint32_t totalPointCount(source->getNumPoints());
    
    PointBuffer source_data(source->getSchema(), totalPointCount);
    StageSequentialIterator* source_iter = source->createSequentialIterator(source_data);

    boost::uint32_t  numRead = source_iter->read(source_data);
    assert(numRead == source_data.getNumPoints());

    delete source_iter;
    delete source;



    Options candidateOptions;
    {
        candidateOptions.add<std::string>("filename", m_candidateFile);
        candidateOptions.add<bool>("debug", isDebug());
        candidateOptions.add<boost::uint32_t>("verbose", getVerboseLevel());
    }

    Stage* candidate = AppSupport::makeReader(candidateOptions);

    candidate->initialize();    


    IndexedPointBuffer candidate_data(candidate->getSchema(), totalPointCount);
    StageSequentialIterator* candidate_iter = candidate->createSequentialIterator(candidate_data);

    numRead = candidate_iter->read(candidate_data);
    assert(numRead == candidate_data.getNumPoints());
        
    delete candidate_iter;    


    if (source_data.getNumPoints() != candidate_data.getNumPoints())
    {
        std::cerr << "Source and candidate files do not have the same point count, testing each source point only!" << std::endl;
    }
    

    // m_summary_x(xd);
    // m_summary_y(yd);
    // m_summary_z(zd);

    if (m_outputFileName.size())
    {
        m_outputStream = FileUtils::createFile(m_outputFileName);
    }

    candidate_data.build(m_3d);
    boost::uint32_t count(std::min(source_data.getNumPoints(), candidate_data.getNumPoints()));
    


    boost::scoped_ptr<std::map<Point, Point> > points(cumulatePoints(source_data, candidate_data));
    if (m_OutputDetail)
    {
        outputDetail(source_data, candidate_data, points.get());
        return 0;
    }
    
    std::map<Point, Point>::const_iterator i;
    for(i = points->begin(); i != points->end(); ++i)
    {
        Point const& s = i->first;
        Point const& c = i->second;

        double xd = s.x - c.x;
        double yd = s.y - c.y;
        double zd = s.z - c.z;        
        m_summary_x(xd);
        m_summary_y(yd);
        m_summary_z(zd);        
    }
    
    std::string headline("------------------------------------------------------------------------------------------");
    std::cout << headline << std::endl;
    std::cout << " Delta summary for source '" << m_sourceFile << "' and candidate '" << m_candidateFile <<"'" << std::endl;
    std::cout << headline << std::endl;
    std::cout << std::endl;
    
    std::string thead("----------- --------------- --------------- --------------");
    std::cout << thead << std::endl;
    std::cout << " Dimension       X             Y                  Z    " << std::endl;
    std::cout << thead << std::endl;
    
    boost::format fmt("%.4f");
    double sminx  = (boost::accumulators::min)(m_summary_x);
//.........这里部分代码省略.........

std::vector<Triple> Seek::getVel(unsigned int ticks, unsigned int delta_ticks) {
    Triple cp, tp;

    std::tie(cp, tp) = points(this->character, this->target);
    return std::vector<Triple>(1, (tp - cp).normalized() * maxSpeed);
}

void ElementTests::setup()
{
  /**********************************

   _________________________________
   |               |               |
   |               |               |
   |               |               |
   |       1       |       3       |
   |               |               |
   |               |               |
   |               |               |
   ---------------------------------
   |       |       |               |
   |   7   |   6   |               |
   |       |       |               |
   |-------0-------|       2       |
   |       | 11| 10|               |
   |   4   |---5---|               |
   |       | 8 | 9 |               |
   ---------------------------------

   in the code:
   _sw: 0
   _nw: 1
   _se: 2
   _ne: 3
   _sw_se: 5
   _sw_ne: 6
   _sw_se_se:  9
   _sw_se_ne: 10

   *********************************/

  // first, build a simple mesh
  FieldContainer<double> quadPoints(4,2);

  quadPoints(0,0) = 0.0; // x1
  quadPoints(0,1) = 0.0; // y1
  quadPoints(1,0) = 1.0;
  quadPoints(1,1) = 0.0;
  quadPoints(2,0) = 1.0;
  quadPoints(2,1) = 1.0;
  quadPoints(3,0) = 0.0;
  quadPoints(3,1) = 1.0;

  int H1Order = 2;
  int testOrder = H1Order;
  int horizontalCells = 2;
  int verticalCells = 2;

  double eps = 1.0; // not really testing for sharp gradients right now--just want to see if things basically work
  double beta_x = 1.0;
  double beta_y = 1.0;

  BFPtr confusionBF = ConfusionBilinearForm::confusionBF(eps,beta_x,beta_y);

  _mesh = MeshFactory::buildQuadMesh(quadPoints, horizontalCells, verticalCells, confusionBF, H1Order, testOrder);

  // the right way to determine the southwest element, etc. is as follows:
  FieldContainer<double> points(4,2);
  // southwest center:
  points(0,0) = 0.25;
  points(0,1) = 0.25;
  // southeast center:
  points(1,0) = 0.75;
  points(1,1) = 0.25;
  // northwest center:
  points(2,0) = 0.25;
  points(2,1) = 0.75;
  // northeast center:
  points(3,0) = 0.75;
  points(3,1) = 0.75;
  vector<ElementPtr> elements = _mesh->elementsForPoints(points, false);

  _sw = elements[0]; // as presently implemented, cellID = 0
  _se = elements[1]; // as presently implemented, cellID = 2
  _nw = elements[2]; // as presently implemented, cellID = 1
  _ne = elements[3]; // as presently implemented, cellID = 3

  vector<GlobalIndexType> cellIDsToRefine;
  cellIDsToRefine.push_back(_sw->cellID());

  _mesh->hRefine(cellIDsToRefine,RefinementPattern::regularRefinementPatternQuad());

  // get the new elements; these are all within the original SW quadrant
  // southwest center:
  points(0,0) = 0.125;
  points(0,1) = 0.125;
  // southeast center:
  points(1,0) = 0.375;
  points(1,1) = 0.125;
  // northwest center:
  points(2,0) = 0.125;
  points(2,1) = 0.375;
  // northeast center:
  points(3,0) = 0.375;
  points(3,1) = 0.375;

  elements = _mesh->elementsForPoints(points, false);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    }

    if ( e->button() == Qt::LeftButton )
    {
      startCapturing();
    }
    else if ( e->button() == Qt::RightButton )
    {
      // End of string

      //lines: bail out if there are not at least two vertices
      if ( mode() == CaptureLine && size() < 2 )
      {
        stopCapturing();
        return;
      }

      //polygons: bail out if there are not at least two vertices
      if ( mode() == CapturePolygon && size() < 3 )
      {
        stopCapturing();
        return;
      }

      //create QgsFeature with wkb representation
      QgsFeature* f = new QgsFeature( 0, "WKBLineString" );

      QgsGeometry *g;

      if ( mode() == CaptureLine )
      {
        if ( layerWKBType == QGis::WKBLineString || layerWKBType == QGis::WKBLineString25D )
        {
          g = QgsGeometry::fromPolyline( points().toVector() );
        }
        else if ( layerWKBType == QGis::WKBMultiLineString || layerWKBType == QGis::WKBMultiLineString25D )
        {
          g = QgsGeometry::fromMultiPolyline( QgsMultiPolyline() << points().toVector() );
        }
        else
        {
          QMessageBox::critical( 0, tr( "Error" ), tr( "Cannot add feature. Unknown WKB type" ) );
          stopCapturing();
          return; //unknown wkbtype
        }

        f->setGeometry( g );
      }
      else // polygon
      {
        QgsGeometry *g;
        if ( layerWKBType == QGis::WKBPolygon ||  layerWKBType == QGis::WKBPolygon25D )
        {
          g = QgsGeometry::fromPolygon( QgsPolygon() << points().toVector() );
        }
        else if ( layerWKBType == QGis::WKBMultiPolygon ||  layerWKBType == QGis::WKBMultiPolygon25D )
        {
          g = QgsGeometry::fromMultiPolygon( QgsMultiPolygon() << ( QgsPolygon() << points().toVector() ) );
        }
        else
        {
          QMessageBox::critical( 0, tr( "Error" ), tr( "Cannot add feature. Unknown WKB type" ) );
          stopCapturing();
          return; //unknown wkbtype
        }

int points(int start, int end, long long x, long long y) {
    int c = (start + end)/2;
    if (start == end) return start;
    if (r[c] < (x*x)+(y*y)) points(c+1, end, x, y);
    else points(start, c, x, y);
}

    // transfer to normal lists
    storedPoints().transfer(dynPoints);

    pointId.shrink();
    dynEdges.shrink();

    // Build inverse mapping (NASTRAN pointId -> index)
    Map<label> mapPointId(2*pointId.size());
    forAll(pointId, i)
    {
        mapPointId.insert(pointId[i], i);
    }

    // note which points were really used and which can be culled
    PackedBoolList usedPoints(points().size());


    // Pass1: relabel edges
    // ~~~~~~~~~~~~~~~~~~~~
    forAll(dynEdges, i)
    {
        edge& e = dynEdges[i];
        e[0] = mapPointId[e[0]];
        e[1] = mapPointId[e[1]];

        usedPoints.set(e[0]);
        usedPoints.set(e[1]);
    }
    pointId.clearStorage();
    mapPointId.clear();

void Foam::processorPolyPatch::calcGeometry(PstreamBuffers& pBufs)
{
    if (Pstream::parRun())
    {
        {
            UIPstream fromNeighbProc(neighbProcNo(), pBufs);

            fromNeighbProc
                >> neighbFaceCentres_
                >> neighbFaceAreas_
                >> neighbFaceCellCentres_;
        }

        // My normals
        vectorField faceNormals(size());

        // Neighbour normals
        vectorField nbrFaceNormals(neighbFaceAreas_.size());

        // Face match tolerances
        scalarField tols = calcFaceTol(*this, points(), faceCentres());

        // Calculate normals from areas and check
        forAll(faceNormals, facei)
        {
            scalar magSf = mag(faceAreas()[facei]);
            scalar nbrMagSf = mag(neighbFaceAreas_[facei]);
            scalar avSf = (magSf + nbrMagSf)/2.0;

            // For small face area calculation the results of the area
            // calculation have been found to only be accurate to ~1e-20
            if (magSf < SMALL || nbrMagSf < SMALL)
            {
                // Undetermined normal. Use dummy normal to force separation
                // check.
                faceNormals[facei] = point(1, 0, 0);
                nbrFaceNormals[facei] = -faceNormals[facei];
                tols[facei] = GREAT;
            }
            else if (mag(magSf - nbrMagSf) > matchTolerance()*sqr(tols[facei]))
            {
                fileName nm
                (
                    boundaryMesh().mesh().time().path()
                   /name()+"_faces.obj"
                );

                Pout<< "processorPolyPatch::calcGeometry : Writing my "
                    << size()
                    << " faces to OBJ file " << nm << endl;

                writeOBJ(nm, *this, points());

                OFstream ccStr
                (
                    boundaryMesh().mesh().time().path()
                    /name() + "_faceCentresConnections.obj"
                );

                Pout<< "processorPolyPatch::calcGeometry :"
                    << " Dumping cell centre lines between"
                    << " corresponding face centres to OBJ file" << ccStr.name()
                    << endl;

                label vertI = 0;

                forAll(faceCentres(), faceI)
                {
                    const point& c0 = neighbFaceCentres_[faceI];
                    const point& c1 = faceCentres()[faceI];

                    writeOBJ(ccStr, c0, c1, vertI);
                }

                FatalErrorInFunction
                    << "face " << facei << " area does not match neighbour by "
                    << 100*mag(magSf - nbrMagSf)/avSf
                    << "% -- possible face ordering problem." << endl
                    << "patch:" << name()
                    << " my area:" << magSf
                    << " neighbour area:" << nbrMagSf
                    << " matching tolerance:"
                    << matchTolerance()*sqr(tols[facei])
                    << endl
                    << "Mesh face:" << start()+facei
                    << " vertices:"
                    << UIndirectList<point>(points(), operator[](facei))()
                    << endl
                    << "If you are certain your matching is correct"
                    << " you can increase the 'matchTolerance' setting"
                    << " in the patch dictionary in the boundary file."
                    << endl
                    << "Rerun with processor debug flag set for"
                    << " more information." << exit(FatalError);
            }
            else
            {

v