C++ Real函数代码示例

	Real Math<Real>::Log(Real fValue)
	{
		return Real(log((long double) fValue));
	}

    void FdSimpleExtOUJumpSwingEngine::calculate() const {

        // 1. Exercise
        boost::shared_ptr<SwingExercise> swingExercise(
            boost::dynamic_pointer_cast<SwingExercise>(arguments_.exercise));

        QL_REQUIRE(swingExercise, "Swing exercise supported only");

        // 2. Mesher
        const std::vector<Time> exerciseTimes
            = swingExercise->exerciseTimes(rTS_->dayCounter(),
                                           rTS_->referenceDate());

        const Time maturity = exerciseTimes.back();
        const boost::shared_ptr<StochasticProcess1D> ouProcess(
                              process_->getExtendedOrnsteinUhlenbeckProcess());
        const boost::shared_ptr<Fdm1dMesher> xMesher(
                     new FdmSimpleProcess1dMesher(xGrid_, ouProcess,maturity));

        const boost::shared_ptr<Fdm1dMesher> yMesher(
                        new ExponentialJump1dMesher(yGrid_,
                                                    process_->beta(),
                                                    process_->jumpIntensity(),
                                                    process_->eta()));
        const boost::shared_ptr<Fdm1dMesher> exerciseMesher(
                       new Uniform1dMesher(
                           0, static_cast<Real>(arguments_.maxExerciseRights),
                           arguments_.maxExerciseRights+1));

        const boost::shared_ptr<FdmMesher> mesher(
            new FdmMesherComposite(xMesher, yMesher, exerciseMesher));

        // 3. Calculator
        boost::shared_ptr<FdmInnerValueCalculator> calculator(
                                                    new FdmZeroInnerValue());
        // 4. Step conditions
        std::list<boost::shared_ptr<StepCondition<Array> > > stepConditions;
        std::list<std::vector<Time> > stoppingTimes;

        // 4.1 Bermudan step conditions
        stoppingTimes.push_back(exerciseTimes);

        const boost::shared_ptr<StrikedTypePayoff> payoff =
            boost::dynamic_pointer_cast<StrikedTypePayoff>(arguments_.payoff);
        boost::shared_ptr<FdmInnerValueCalculator> exerciseCalculator(
                      new FdmExtOUJumpModelInnerValue(payoff, mesher, shape_));

        stepConditions.push_back(boost::shared_ptr<StepCondition<Array> >(
                new FdmSimpleSwingCondition(exerciseTimes, mesher,
                                            exerciseCalculator, 2)));

        boost::shared_ptr<FdmStepConditionComposite> conditions(
                new FdmStepConditionComposite(stoppingTimes, stepConditions));


        // 5. Boundary conditions
        const FdmBoundaryConditionSet boundaries;

        // 6. set-up solver
        FdmSolverDesc solverDesc = { mesher, boundaries, conditions,
                                     calculator, maturity, tGrid_, 0 };

        const boost::shared_ptr<FdmSimple3dExtOUJumpSolver> solver(
            new FdmSimple3dExtOUJumpSolver(
                                    Handle<ExtOUWithJumpsProcess>(process_),
                                    rTS_, solverDesc, schemeDesc_));

        const Real x = process_->initialValues()[0];
        const Real y = process_->initialValues()[1];

        std::vector< std::pair<Real, Real> > exerciseValues;
        for (Size i=arguments_.minExerciseRights;
             i <= arguments_.maxExerciseRights; ++i) {
            const Real z = Real(i);
            exerciseValues.push_back(std::pair<Real, Real>(
                                       solver->valueAt(x, y, z), z));
        }
        const Real z = std::max_element(exerciseValues.begin(),
                                        exerciseValues.end())->second;

        results_.value = solver->valueAt(x, y, z);
    }

int Execute(int argc,char* argv[])
{
	int i;
	cmdLineString In,Out;
	cmdLineReadable Binary,Verbose,NoResetSamples,NoClipTree,Confidence,Manifold,PolygonMesh;
	cmdLineInt Depth(8),SolverDivide(8),IsoDivide(8),Refine(3);
	cmdLineInt KernelDepth;
	cmdLineFloat SamplesPerNode(1.0f),Scale(1.1f);
	char* paramNames[]=
	{
		"in","depth","out","refine","noResetSamples","noClipTree",
		"binary","solverDivide","isoDivide","scale","verbose",
		"kernelDepth","samplesPerNode","confidence","manifold","polygonMesh"
	};
	cmdLineReadable* params[]=
	{
		&In,&Depth,&Out,&Refine,&NoResetSamples,&NoClipTree,
		&Binary,&SolverDivide,&IsoDivide,&Scale,&Verbose,
		&KernelDepth,&SamplesPerNode,&Confidence,&Manifold,&PolygonMesh
	};
	int paramNum=sizeof(paramNames)/sizeof(char*);
	int commentNum=0;
	char **comments;

	comments=new char*[paramNum+7];
	for(i=0;i<paramNum+7;i++){comments[i]=new char[1024];}

	const char* Rev = "Rev: V2 ";
	const char* Date = "Date: 2006-11-09 (Thur, 09 Nov 2006) ";

	cmdLineParse(argc-1,&argv[1],paramNames,paramNum,params,0);

	if(Verbose.set){echoStdout=1;}

	DumpOutput2(comments[commentNum++],"Running Multi-Grid Octree Surface Reconstructor (degree %d). Version 3\n", Degree);
	if(In.set)				{DumpOutput2(comments[commentNum++],"\t--in %s\n",In.value);}
	if(Out.set)				{DumpOutput2(comments[commentNum++],"\t--out %s\n",Out.value);}
	if(Binary.set)			{DumpOutput2(comments[commentNum++],"\t--binary\n");}
	if(Depth.set)			{DumpOutput2(comments[commentNum++],"\t--depth %d\n",Depth.value);}
	if(SolverDivide.set)	{DumpOutput2(comments[commentNum++],"\t--solverDivide %d\n",SolverDivide.value);}
	if(IsoDivide.set)		{DumpOutput2(comments[commentNum++],"\t--isoDivide %d\n",IsoDivide.value);}
	if(Refine.set)			{DumpOutput2(comments[commentNum++],"\t--refine %d\n",Refine.value);}
	if(Scale.set)			{DumpOutput2(comments[commentNum++],"\t--scale %f\n",Scale.value);}
	if(KernelDepth.set)		{DumpOutput2(comments[commentNum++],"\t--kernelDepth %d\n",KernelDepth.value);}
	if(SamplesPerNode.set)	{DumpOutput2(comments[commentNum++],"\t--samplesPerNode %f\n",SamplesPerNode.value);}
	if(NoResetSamples.set)	{DumpOutput2(comments[commentNum++],"\t--noResetSamples\n");}
	if(NoClipTree.set)		{DumpOutput2(comments[commentNum++],"\t--noClipTree\n");}
	if(Confidence.set)		{DumpOutput2(comments[commentNum++],"\t--confidence\n");}
	if(Manifold.set)		{DumpOutput2(comments[commentNum++],"\t--manifold\n");}
	if(PolygonMesh.set)		{DumpOutput2(comments[commentNum++],"\t--polygonMesh\n");}

	double t;
	double tt=Time();
	Point3D<float> center;
	Real scale=1.0;
	Real isoValue=0;
	//////////////////////////////////
	// Fix courtesy of David Gallup //
	TreeNodeData::UseIndex = 1;     //
	//////////////////////////////////
	Octree<Degree> tree;
	PPolynomial<Degree> ReconstructionFunction=PPolynomial<Degree>::GaussianApproximation();

	center.coords[0]=center.coords[1]=center.coords[2]=0;
	if(!In.set || !Out.set)
	{
		ShowUsage(argv[0]);
		return 0;
	}
	
	TreeOctNode::SetAllocator(MEMORY_ALLOCATOR_BLOCK_SIZE);

	t=Time();
	int kernelDepth=Depth.value-2;
	if(KernelDepth.set){kernelDepth=KernelDepth.value;}

	tree.setFunctionData(ReconstructionFunction,Depth.value,0,Real(1.0)/(1<<Depth.value));
	DumpOutput("Function Data Set In: %lg\n",Time()-t);
	DumpOutput("Memory Usage: %.3f MB\n",float(MemoryInfo::Usage())/(1<<20));
	if(kernelDepth>Depth.value){
		fprintf(stderr,"KernelDepth can't be greater than Depth: %d <= %d\n",kernelDepth,Depth.value);
		return EXIT_FAILURE;
	}


	t=Time();
#if 1
	tree.setTree(In.value,Depth.value,Binary.set,kernelDepth,Real(SamplesPerNode.value),Scale.value,center,scale,!NoResetSamples.set,Confidence.set);
#else
if(Confidence.set){
	tree.setTree(In.value,Depth.value,Binary.set,kernelDepth,Real(SamplesPerNode.value),Scale.value,center,scale,!NoResetSamples.set,0,1);
}
else{
	tree.setTree(In.value,Depth.value,Binary.set,kernelDepth,Real(SamplesPerNode.value),Scale.value,center,scale,!NoResetSamples.set,0,0);
}
#endif
	DumpOutput2(comments[commentNum++],"#             Tree set in: %9.1f (s), %9.1f (MB)\n",Time()-t,tree.maxMemoryUsage);
	DumpOutput("Leaves/Nodes: %d/%d\n",tree.tree.leaves(),tree.tree.nodes());
	DumpOutput("   Tree Size: %.3f MB\n",float(sizeof(TreeOctNode)*tree.tree.nodes())/(1<<20));
	DumpOutput("Memory Usage: %.3f MB\n",float(MemoryInfo::Usage())/(1<<20));
//.........这里部分代码省略.........

 double modelSpaceScale() { return Real(_params[12]); }

 double minUserResolution() { return Real(_params[18]); }

 Real PiecewiseConstantVariance::totalVariance(Size i) const {
     QL_REQUIRE(i<variances().size(),
                "invalid step index");
     return std::accumulate(variances().begin(),
                            variances().begin()+i+1, Real(0.0));
 }

template<class Real> Real
Random (void)
{
  return Real (rand ()) / RAND_MAX;
}

void* cGraphicsMFC::_getScaledImage(cSpriteIcon *picon, CRect &targetRect)
{
	int i;
	BOOL resourceidwasfound = FALSE;
	cMapEntry_ResourceIDToScaledDC mapentry; /* There are four possible cases here: (1), (2.1), (2.2) and (2.3).
		(1) If this is the first time this view is drawing an icon with this resourceID, then
		we need to intialize all the field of mapentry and add it to cGraphicsMFC::_map_ResourceIDToScaledDC. 
		(2) We have already done step (a) once in some instance of cGraphicsMFC during this program
		run, and mapentry is a copy of this entry with a valid _pResourceImage.  Now three subcases
		arise.  (2.1) This is the first time picon has called _getScaledImage, and !picon->imageloaded()
		holds.  In this case, the targetRect passed in was not in fact of the correct aspect, that is
		it doesn't match the aspect of the _pResourceImage.  In this case, we recalculate targerRect,
		and once again search for a mapentry with matching resourceID and targetRect.  We will
		at least find one with matching resourceID.  At this point whether or not we had to do
		case (2.1) we now move on to (2.2) or (2.3).  (2.2) None of	the mapentry with the correct
		_resourceID has the correct _targetRect. In this case we make a new _pScaledImage, put it
		in the mapentry, and add the new mapentry to  cGraphicsMFC::_map_ResourceIDToScaledDC. (2.3) We find
		a mapentry with a _targetRect so that it matchesTargetRect(targetRect). In	this case we can return
		mapentry._pScaledImage. 
		As a further complication, we will make separate mapentries with separate _pResourceImage
		for the same _resoruceID in the case where we have two presetaspect icons with different
		aspect.*/
	for (i=0; i< cGraphicsMFC::_map_ResourceIDToScaledDC.GetSize(); i++)
		if (cGraphicsMFC::_map_ResourceIDToScaledDC[i].matchesResourceAndAspect(picon))
		{
			resourceidwasfound = TRUE;
			mapentry = cGraphicsMFC::_map_ResourceIDToScaledDC.GetAt(i); //Make a copy.
			if (mapentry.matchesTargetRect(targetRect))
			{//This is Case 2.3, the easy case where you find just what you need.
				cGraphicsMFC::_map_ResourceIDToScaledDC.ElementAt(i)._lifespan = 
					cMapEntry_ResourceIDToScaledDC::FRESHLIFESPAN;
				return mapentry._pScaledImage; //You're done!
			}
		}
	/* If we're here, we're in case (1) or case (2.2) or we may be in case (2.1) where the
		picon fed in an incorrect targetrect.  We assume the worst for now and start to
		we are going to fill in some or all of the fields of mapentry to later add to
		 cGraphicsMFC::_map_ResourceIDToScaledDC. */
	if (!resourceidwasfound) /* In this case (1), mapentry is the default constructor entry; no
		corresponding resourceID was found. Initialize the mapentry _resourceID and _pResourceImage
		to start with. */
	{		// First set the dummy mapentry's _resourceID
		mapentry._resourceID = picon->resourceID();
			//Second construct the dummy mapentry's _pResourceImage
		if (picon->usesTransparentMask())
			mapentry._pResourceImage = new cTransparentMemoryDC();
		else
			mapentry._pResourceImage = new cMemoryDC();
			/* The default constructors make the cMemoryDC 1 x 1 and set 
			the _resize_on_load flag to TRUE so we will make the bitmap match
			the image being loaded in loadResourceBitmap. */
		ASSERT(mapentry._pResourceImage);
			//Third, load the resource bitmap into the mapentry._pResourceImage
		mapentry._pResourceImage->loadResourceBitmap(picon->resourceID());
			//If necessary, adjust the size of the mapentry._pResourceImage if picon->presetaspect()
		if (picon->presetaspect())
		{
			int sourcex = mapentry._pResourceImage->cx();
			int sourcey = mapentry._pResourceImage->cy();
			ASSERT(sourcey);
			ASSERT(picon->aspect());
			Real tempaspect = Real(sourcex)/sourcey;
			// We want to make the source bitmap match the _aspect of the picon.
			if (tempaspect > picon->aspect()) //sourcex too big, make it smaller than tempaspect*sourcey.
				sourcex = picon->aspect() * sourcey;
			if (tempaspect < picon->aspect()) //sourcey too big, make it smaller than sourcex/tempaspect.
				sourcey = sourcex / picon->aspect();
			cMemoryDC *pMemDC_bitmap = new cMemoryDC();
			ASSERT(pMemDC_bitmap);
			pMemDC_bitmap->setSize(sourcex, sourcey);
			mapentry._pResourceImage->copyTo(pMemDC_bitmap, CRect(0,0,sourcex, sourcey));
			delete mapentry._pResourceImage;
			mapentry._pResourceImage = NULL;
			mapentry._pResourceImage = pMemDC_bitmap;
		}
		//Initialize mapentry._aspect to match the dimensions of _pResourceImage.
		mapentry._aspect = (Real)(mapentry._pResourceImage->cx())/ mapentry._pResourceImage->cy();
		/* What about putting the newly initialized pmapentry into cGraphicsMFC::_map_ResourceIDToScaledDC?
			We'll do that in the next block of code, after we have initialized ALL of 
			its fields. */
		picon->setAspect(mapentry._aspect); //Do this here to be sure in case roundoff error made not match
	}

		/* At this point, we know that mapentry has a valid _resourceID and _pResourceImage, 
		and that in case of picon->presetaspect, the _apsect of the _pResourceImage as been 
		adjusted to match the picon->aspect. But since
		we didn't manage to do a case 2.3 bail in the initial loop, we know the _targetRect 
		doesn't match the size of the _pScaledImage.  There are two possibilities.  Either (2.2)
		targetRect is a legitimate size and we actually need a new _pScaledImage.
			Or maybe (2.1) targetRect is wrong because picon didn't know what it's shaped like yet.
		Tell the picon the correct aspect of the basic bitmap in case this is the very first time
		the picon is calling draw.  Also correct the targetRect to match the correct aspect. Note
		that the default starting aspect of picon was 1.0. 
			Of course in the case of picon->presetaspect, we shouldn't need to adjust pcion->aspect,
		but there could be some discrepancy due to round off, so back-correcting it here can't
		hurt. */
	if (!picon->imageloaded()) /* This means the picon has an _aspect
		that may not match the aspect of the image you loaded, and this  means that the aspect
		of the	targetRect doens't necessarily match the image.  Fix targetRect in here. */
	{
//.........这里部分代码省略.........

Ogre::Vector3* Maze::translatePos( Ogre::Vector3* pos )
{
	return new Ogre::Vector3(Real(pos->x - (this->mWidth / 2.0f * 100)),Real(0),Real(pos->y - (this->mHeight / 2.0f * 100)));
}

    //-----------------------------------------------------------------------
    void InstanceBatch::_notifyCurrentCamera( Camera* cam )
    {
        mCurrentCamera = cam;

        //See DistanceLodStrategy::getValueImpl()
        //We use our own because our SceneNode is just filled with zeroes, and updating it
        //with real values is expensive, plus we would need to make sure it doesn't get to
        //the shader
		Real depth = Math::Sqrt(getSquaredViewDepth(cam)) - getBoundingRadius();          
        depth = std::max( depth, Real(0) );

        Real lodValue = depth * cam->_getLodBiasInverse();

        //Now calculate Material LOD
        /*const LodStrategy *materialStrategy = m_material->getLodStrategy();
        
        //Calculate LOD value for given strategy
        Real lodValue = materialStrategy->getValue( this, cam );*/

        //Get the index at this depth
        unsigned short idx = mMaterial->getLodIndex( lodValue );

        //TODO: Replace subEntity for MovableObject
        // Construct event object
        /*EntityMaterialLodChangedEvent subEntEvt;
        subEntEvt.subEntity = this;
        subEntEvt.camera = cam;
        subEntEvt.lodValue = lodValue;
        subEntEvt.previousLodIndex = m_materialLodIndex;
        subEntEvt.newLodIndex = idx;

        //Notify LOD event listeners
        cam->getSceneManager()->_notifyEntityMaterialLodChanged(subEntEvt);*/

        // Change LOD index
        mMaterialLodIndex = idx;

        mBeyondFarDistance = false;

        if (cam->getUseRenderingDistance() && mUpperDistance > 0)
        {
            if (depth > mUpperDistance)
                mBeyondFarDistance = true;
        }

        if (!mBeyondFarDistance && cam->getUseMinPixelSize() && mMinPixelSize > 0)
        {
            Real pixelRatio = cam->getPixelDisplayRatio();

            Ogre::Vector3 objBound =
                getBoundingBox().getSize() * getParentNode()->_getDerivedScale();
            objBound.x = Math::Sqr(objBound.x);
            objBound.y = Math::Sqr(objBound.y);
            objBound.z = Math::Sqr(objBound.z);
            float sqrObjMedianSize = std::max(
                std::max(std::min(objBound.x, objBound.y), std::min(objBound.x, objBound.z)),
                std::min(objBound.y, objBound.z));

            // If we have a perspective camera calculations are done relative to distance
            Real sqrDistance = 1;

            if (cam->getProjectionType() == PT_PERSPECTIVE)
                sqrDistance = getSquaredViewDepth(cam->getLodCamera()); // it's ok

            mBeyondFarDistance =
                sqrObjMedianSize < sqrDistance * Math::Sqr(pixelRatio * mMinPixelSize);
        }

        if (mParentNode)
        {
            MovableObjectLodChangedEvent evt;
            evt.movableObject = this;
            evt.camera = cam;

            cam->getSceneManager()->_notifyMovableObjectLodChanged(evt);
        }

        mRenderingDisabled = mListener && !mListener->objectRendering(this, cam);

        // MovableObject::_notifyCurrentCamera( cam ); // it does not suit
    }

/* Implementation *************************************************************/
void CSyncUsingPil::ProcessDataInternal(CParameter& ReceiverParam)
{
	int i;

	/**************************************************************************\
	* Frame synchronization detection										   *
	\**************************************************************************/
	_BOOLEAN bSymbolIDHasChanged = FALSE;

	if ((bSyncInput == FALSE) && (bAquisition == TRUE))
	{
#ifdef USE_DRM_FRAME_SYNC_IR_BASED
		/* DRM frame synchronization using impulse response ----------------- */
		/* We assume that the current received OFDM symbol is the first symbol
		   in a DRM frame and estimate the channel transfer function at the
		   pilot positions (the positions of pilots in the first OFDM symbol in
		   a DRM frame). Then we calculate an FFT to get the impulse response of
		   the channel. If the assumption was correct and this really was the
		   correct OFDM symbol, we will get something which looks like an
		   impulse response (contains peaks -> peak-to-average ratio is high).
		   If not, we will certainly get only noise -> no peaks -> peak to
		   average ratio is small. This is because the transmitted values at
		   the pilot positions are different from the values at the pilot cells
		   when transmitting the correct OFDM symbol (which we assumed) */

		/* Pick pilot positions and calculate "test" channel estimation */
		int iCurIndex = 0;
		for (i = 0; i < iNumCarrier; i++)
		{
			if (_IsScatPil(ReceiverParam.matiMapTab[0][i]))
			{
				/* Get channel estimate */
				veccChan[iCurIndex] =
					(*pvecInputData)[i] / ReceiverParam.matcPilotCells[0][i];

				/* We have to introduce a new index because not on all carriers
				   is a pilot */
				iCurIndex++;
			}
		}

		/* Calculate abs(IFFT) for getting estimate of impulse response */
		vecrTestImpResp = Abs(Ifft(veccChan, FftPlan));

		/* Calculate peak to average */
		const CReal rResultIREst = Max(vecrTestImpResp) / Sum(vecrTestImpResp);

		/* Store correlation results in a shift register for finding the peak */
		vecrCorrHistory.AddEnd(rResultIREst);

#else

		/* DRM frame synchronization based on time pilots ------------------- */
		/* Calculate correlation of received cells with pilot pairs */
		CReal rResultPilPairCorr = (CReal) 0.0;
		for (i = 0; i < iNumPilPairs; i++)
		{
			/* Actual correlation */
			const CComplex cCorrRes = (*pvecInputData)[vecPilCorr[i].iIdx1] *
				Conj(vecPilCorr[i].cPil1) *
				Conj((*pvecInputData)[vecPilCorr[i].iIdx2]) *
				vecPilCorr[i].cPil2 * cR_HH;

			rResultPilPairCorr += Real(cCorrRes);
		}

		/* Store correlation results in a shift register for finding the peak */
		vecrCorrHistory.AddEnd(rResultPilPairCorr);
#endif


		/* Finding beginning of DRM frame in results ------------------------ */
		/* Wait until history is filled completly */
		if (iInitCntFraSy > 0)
			iInitCntFraSy--;
		else
		{
			/* Search for maximum */
			int iMaxIndex = 0;
			CReal rMaxValue = -_MAXREAL;
			for (i = 0; i < iNumSymPerFrame; i++)
			{
				if (vecrCorrHistory[i] > rMaxValue)
				{
					rMaxValue = vecrCorrHistory[i];
					iMaxIndex = i;
				}
			}

			/* For initial frame synchronization, use maximum directly */
			if (bInitFrameSync == TRUE)
			{
				/* Reset init flag */
				bInitFrameSync = FALSE;

				/* Set symbol ID index according to received data */
				iSymbCntFraSy = iNumSymPerFrame - iMaxIndex - 1;
			}
			else
//.........这里部分代码省略.........

int Execute( int argc , char* argv[] )
{
	int i;
	int paramNum = sizeof(params)/sizeof(cmdLineReadable*);
	int commentNum=0;
	char **comments;

	comments = new char*[paramNum+7];
	for( i=0 ; i<paramNum+7 ; i++ ) comments[i] = new char[1024];

	if( Verbose.set ) echoStdout=1;

	XForm4x4< Real > xForm , iXForm;
	if( XForm.set )
	{
		FILE* fp = fopen( XForm.value , "r" );
		if( !fp )
		{
			fprintf( stderr , "[WARNING] Could not read x-form from: %s\n" , XForm.value );
			xForm = XForm4x4< Real >::Identity();
		}
		else
		{
			for( int i=0 ; i<4 ; i++ ) for( int j=0 ; j<4 ; j++ ) fscanf( fp , " %f " , &xForm( i , j ) );
			fclose( fp );
		}
	}
	else xForm = XForm4x4< Real >::Identity();
	iXForm = xForm.inverse();

	DumpOutput2( comments[commentNum++] , "Running Screened Poisson Reconstruction (Version 5.5)\n" , Degree );
	char str[1024];
	for( int i=0 ; i<paramNum ; i++ )
		if( params[i]->set )
		{
			params[i]->writeValue( str );
			if( strlen( str ) ) DumpOutput2( comments[commentNum++] , "\t--%s %s\n" , params[i]->name , str );
			else                DumpOutput2( comments[commentNum++] , "\t--%s\n" , params[i]->name );
		}

	double t;
	double tt=Time();
	Real isoValue = 0;

	POctree< Degree , OutputDensity > tree;
	tree.threads = Threads.value;
	if( !In.set )
	{
		ShowUsage(argv[0]);
		return 0;
	}
	if( !MaxSolveDepth.set ) MaxSolveDepth.value = Depth.value;
	if( SolverDivide.value<MinDepth.value )
	{
		fprintf( stderr , "[WARNING] %s must be at least as large as %s: %d>=%d\n" , SolverDivide.name , MinDepth.name , SolverDivide.value , MinDepth.value );
		SolverDivide.value = MinDepth.value;
	}
	if( IsoDivide.value<MinDepth.value )
	{
		fprintf( stderr , "[WARNING] %s must be at least as large as %s: %d>=%d\n" , IsoDivide.name , MinDepth.name , IsoDivide.value , IsoDivide.value );
		IsoDivide.value = MinDepth.value;
	}
	
	OctNode< TreeNodeData< OutputDensity > , Real >::SetAllocator( MEMORY_ALLOCATOR_BLOCK_SIZE );

	t=Time();
	int kernelDepth = KernelDepth.set ?  KernelDepth.value : Depth.value-2;

	tree.setBSplineData( Depth.value , BoundaryType.value );
	if( kernelDepth>Depth.value )
	{
		fprintf( stderr,"[ERROR] %s can't be greater than %s: %d <= %d\n" , KernelDepth.name , Depth.name , KernelDepth.value , Depth.value );
		return EXIT_FAILURE;
	}

	double maxMemoryUsage;
	t=Time() , tree.maxMemoryUsage=0;
	int pointCount = tree.setTree( In.value , Depth.value , MinDepth.value , kernelDepth , Real(SamplesPerNode.value) , Scale.value , Confidence.set , PointWeight.value , AdaptiveExponent.value , xForm );
	tree.ClipTree();
	tree.finalize( IsoDivide.value );

	DumpOutput2( comments[commentNum++] , "#             Tree set in: %9.1f (s), %9.1f (MB)\n" , Time()-t , tree.maxMemoryUsage );
	DumpOutput( "Input Points: %d\n" , pointCount );
	DumpOutput( "Leaves/Nodes: %d/%d\n" , tree.tree.leaves() , tree.tree.nodes() );
	DumpOutput( "Memory Usage: %.3f MB\n" , float( MemoryInfo::Usage() )/(1<<20) );

	maxMemoryUsage = tree.maxMemoryUsage;
	t=Time() , tree.maxMemoryUsage=0;
	tree.SetLaplacianConstraints();
	DumpOutput2( comments[commentNum++] , "#      Constraints set in: %9.1f (s), %9.1f (MB)\n" , Time()-t , tree.maxMemoryUsage );
	DumpOutput( "Memory Usage: %.3f MB\n" , float( MemoryInfo::Usage())/(1<<20) );
	maxMemoryUsage = std::max< double >( maxMemoryUsage , tree.maxMemoryUsage );

	t=Time() , tree.maxMemoryUsage=0;
	tree.LaplacianMatrixIteration( SolverDivide.value, ShowResidual.set , MinIters.value , SolverAccuracy.value , MaxSolveDepth.value , FixedIters.value );
	DumpOutput2( comments[commentNum++] , "# Linear system solved in: %9.1f (s), %9.1f (MB)\n" , Time()-t , tree.maxMemoryUsage );
	DumpOutput( "Memory Usage: %.3f MB\n" , float( MemoryInfo::Usage() )/(1<<20) );
	maxMemoryUsage = std::max< double >( maxMemoryUsage , tree.maxMemoryUsage );

	CoredFileMeshData< Vertex > mesh;
//.........这里部分代码省略.........

inline void
Snapshot
( const Matrix<Int>& preimage,
  const Matrix<Real>& estimates, 
  const Matrix<Int>& itCounts,
        Int numIts,
        bool deflate,
        SnapshotCtrl& snapCtrl )
{
    DEBUG_ONLY(CSE cse("pspec::Snapshot"));
    auto logMap = []( Real alpha ) { return Log(alpha); };
    if( snapCtrl.realSize != 0 && snapCtrl.imagSize != 0 )
    {
        const bool numSave = 
            ( snapCtrl.numSaveFreq > 0 && 
              snapCtrl.numSaveCount >= snapCtrl.numSaveFreq );
        const bool imgSave = 
            ( snapCtrl.imgSaveFreq > 0 && 
              snapCtrl.imgSaveCount >= snapCtrl.imgSaveFreq );
        const bool imgDisp = 
            ( snapCtrl.imgDispFreq > 0 &&
              snapCtrl.imgDispCount >= snapCtrl.imgDispFreq );
        Matrix<Real> invNorms, estMap;
        Matrix<Int> itCountsReord, itCountMap;
        if( numSave || imgSave || imgDisp )
        {
            invNorms = estimates;
            if( deflate )
                RestoreOrdering( preimage, invNorms );
            ReshapeIntoGrid
            ( snapCtrl.realSize, snapCtrl.imagSize, invNorms, estMap );
            if( snapCtrl.itCounts )
            {
                itCountsReord = itCounts;
                if( deflate )
                    RestoreOrdering( preimage, itCountsReord );
                ReshapeIntoGrid
                ( snapCtrl.realSize, snapCtrl.imagSize, itCountsReord, 
                  itCountMap );
            }
        }
        if( numSave )
        {
            ostringstream os;
            os << snapCtrl.numBase << "_" << numIts;
            Write( estMap, os.str(), snapCtrl.numFormat );
            if( snapCtrl.itCounts )
                Write( itCountMap, os.str()+"_counts", snapCtrl.numFormat );
            snapCtrl.numSaveCount = 0;
        }
        if( imgSave || imgDisp )
            EntrywiseMap( estMap, function<Real(Real)>(logMap) );
        if( imgSave )
        {
            ostringstream os;
            os << snapCtrl.imgBase << "_" << numIts;
            Write( estMap, os.str(), snapCtrl.imgFormat );
            if( snapCtrl.itCounts )
                Write( itCountMap, os.str()+"_counts", snapCtrl.imgFormat );
            auto colorMap = GetColorMap();
            SetColorMap( GRAYSCALE_DISCRETE );
            Write( estMap, os.str()+"_discrete", snapCtrl.imgFormat );
            SetColorMap( colorMap );
            snapCtrl.imgSaveCount = 0;
        }
        if( imgDisp )
        {
            ostringstream os;
            os << snapCtrl.imgBase << "_" << numIts;
            Display( estMap, os.str() );       
            if( snapCtrl.itCounts )
                Display( itCountMap, os.str()+"_counts" );
            auto colorMap = GetColorMap();
            SetColorMap( GRAYSCALE_DISCRETE );
            Display( estMap, os.str()+"_discrete" );
            SetColorMap( colorMap );
            snapCtrl.imgDispCount = 0;
        }
    }
}

inline void
PanelLQ( DistMatrix<Real>& A )
{
#ifndef RELEASE
    PushCallStack("internal::PanelLQ");
#endif
    const Grid& g = A.Grid();

    // Matrix views
    DistMatrix<Real>
        ATL(g), ATR(g),  A00(g), a01(g),     A02(g),  aTopRow(g), ABottomPan(g),
        ABL(g), ABR(g),  a10(g), alpha11(g), a12(g),
                         A20(g), a21(g),     A22(g);

    // Temporary distributions
    DistMatrix<Real,STAR,MR> aTopRow_STAR_MR(g);
    DistMatrix<Real,MC,STAR> z_MC_STAR(g);

    PushBlocksizeStack( 1 );
    PartitionDownLeftDiagonal
    ( A, ATL, ATR,
         ABL, ABR, 0 );
    while( ATL.Height() < A.Height() && ATL.Width() < A.Width() )
    {
        RepartitionDownDiagonal
        ( ATL, /**/ ATR,  A00, /**/ a01,     A02,
         /*************/ /**********************/
               /**/       a10, /**/ alpha11, a12,
          ABL, /**/ ABR,  A20, /**/ a21,     A22 );

        View1x2( aTopRow, alpha11, a12 );
        View1x2( ABottomPan, a21, A22 );

        aTopRow_STAR_MR.AlignWith( ABottomPan );
        z_MC_STAR.AlignWith( ABottomPan );
        Zeros( ABottomPan.Height(), 1, z_MC_STAR );
        //--------------------------------------------------------------------//
        const Real tau = Reflector( alpha11, a12 );

        const bool myDiagonalEntry = ( g.Row() == alpha11.ColAlignment() &&
                                       g.Col() == alpha11.RowAlignment() );
        Real alpha = 0;
        if( myDiagonalEntry )
        {
            alpha = alpha11.GetLocal(0,0);
            alpha11.SetLocal(0,0,1);
        }

        aTopRow_STAR_MR = aTopRow;

        Gemv
        ( NORMAL,
          Real(1), ABottomPan.LockedLocalMatrix(),
                   aTopRow_STAR_MR.LockedLocalMatrix(),
          Real(0), z_MC_STAR.LocalMatrix() );
        z_MC_STAR.SumOverRow();

        Ger
        ( -tau,
          z_MC_STAR.LockedLocalMatrix(),
          aTopRow_STAR_MR.LockedLocalMatrix(),
          ABottomPan.LocalMatrix() );

        if( myDiagonalEntry )
            alpha11.SetLocal(0,0,alpha);
        //--------------------------------------------------------------------//
        aTopRow_STAR_MR.FreeAlignments();
        z_MC_STAR.FreeAlignments();

        SlidePartitionDownDiagonal
        ( ATL, /**/ ATR,  A00, a01,     /**/ A02,
               /**/       a10, alpha11, /**/ a12,
         /*************/ /**********************/
          ABL, /**/ ABR,  A20, a21,     /**/ A22 );
    }
    PopBlocksizeStack();
#ifndef RELEASE
    PopCallStack();
#endif
}

template<class Real> void
TriangleCollapse (const Real& edgeRatio, std::vector<TriangleIndex>& triangles,
                  std::vector<Point3D<Real> >& positions, std::vector<Point3D<Real> >* normals)
{
  int i, j, *remapTable, *pointCount, idx[3];
  Point3D<Real> p[3], q[2], c;
  double d[3], a;
  double Ratio = 12.0 / sqrt (3.0); // (Sum of Squares Length / Area) for and equilateral triangle

  remapTable = new int[positions.size ()];
  pointCount = new int[positions.size ()];
  for (i = 0; i < int (positions.size ()); i++)
  {
    remapTable[i] = i;
    pointCount[i] = 1;
  }
  for (i = int (triangles.size () - 1); i >= 0; i--)
  {
    for (j = 0; j < 3; j++)
    {
      idx[j] = triangles[i].idx[j];
      while (remapTable[idx[j]] < idx[j])
      {
        idx[j] = remapTable[idx[j]];
      }
    }
    if (idx[0] == idx[1] || idx[0] == idx[2] || idx[1] == idx[2])
    {
      triangles[i] = triangles[triangles.size () - 1];
      triangles.pop_back ();
      continue;
    }
    for (j = 0; j < 3; j++)
    {
      p[j].coords[0] = positions[idx[j]].coords[0] / pointCount[idx[j]];
      p[j].coords[1] = positions[idx[j]].coords[1] / pointCount[idx[j]];
      p[j].coords[2] = positions[idx[j]].coords[2] / pointCount[idx[j]];
    }
    for (j = 0; j < 3; j++)
    {
      q[0].coords[j] = p[1].coords[j] - p[0].coords[j];
      q[1].coords[j] = p[2].coords[j] - p[0].coords[j];
      d[j] = SquareDistance (p[j], p[(j + 1) % 3]);
    }
    CrossProduct (q[0], q[1], c);
    a = Length (c) / 2;

    if ((d[0] + d[1] + d[2]) * edgeRatio > a * Ratio)
    {
      // Find the smallest edge
      j = 0;
      if (d[1] < d[j])
      {
        j = 1;
      }
      if (d[2] < d[j])
      {
        j = 2;
      }

      int idx1, idx2, idx3;
      if (idx[0] < idx[1])
      {
        if (idx[0] < idx[2])
        {
          idx1 = idx[0];
          idx2 = idx[2];
          idx3 = idx[1];
        }
        else
        {
          idx1 = idx[2];
          idx2 = idx[0];
          idx3 = idx[1];
        }
      }
      else
      {
        if (idx[1] < idx[2])
        {
          idx1 = idx[1];
          idx2 = idx[2];
          idx3 = idx[0];
        }
        else
        {
          idx1 = idx[2];
          idx2 = idx[1];
          idx3 = idx[0];
        }
      }
      positions[idx1].coords[0] += positions[idx2].coords[0] + positions[idx3].coords[0];
      positions[idx1].coords[1] += positions[idx2].coords[1] + positions[idx3].coords[1];
      positions[idx1].coords[2] += positions[idx2].coords[2] + positions[idx3].coords[2];
      if (normals)
      {
        (*normals)[idx1].coords[0] += (*normals)[idx2].coords[0] + (*normals)[idx3].coords[0];
        (*normals)[idx1].coords[1] += (*normals)[idx2].coords[1] + (*normals)[idx3].coords[1];
        (*normals)[idx1].coords[2] += (*normals)[idx2].coords[2] + (*normals)[idx3].coords[2];
      }
//.........这里部分代码省略.........