C++ CvRTrees类代码示例

void Model::Train_rtrees( const SampleSet& samples )
{
	CvRTrees* model = (CvRTrees*)m_pModel;
	CvRTParams* para = (CvRTParams*)m_trainPara;
	model->train(samples.Samples(), CV_ROW_SAMPLE, samples.Labels(), 
		cv::Mat(), cv::Mat(), cv::Mat(), cv::Mat(), *para);
}

int main()
{
    const int train_sample_count = 300;

//#define LEPIOTA
#ifdef LEPIOTA
    const char* filename = "../../../OpenCV_SVN/samples/c/agaricus-lepiota.data";
#else
    const char* filename = "../../../OpenCV_SVN/samples/c/waveform.data";
#endif

    CvDTree dtree;
    CvBoost boost;
    CvRTrees rtrees;
    CvERTrees ertrees;

    CvMLData data;

    CvTrainTestSplit spl( train_sample_count );
    
    data.read_csv( filename );

#ifdef LEPIOTA
    data.set_response_idx( 0 );     
#else
    data.set_response_idx( 21 );     
    data.change_var_type( 21, CV_VAR_CATEGORICAL );
#endif

    data.set_train_test_split( &spl );
    
    printf("======DTREE=====\n");
    dtree.train( &data, CvDTreeParams( 10, 2, 0, false, 16, 0, false, false, 0 ));
    print_result( dtree.calc_error( &data, CV_TRAIN_ERROR), dtree.calc_error( &data ), dtree.get_var_importance() );

#ifdef LEPIOTA
    printf("======BOOST=====\n");
    boost.train( &data, CvBoostParams(CvBoost::DISCRETE, 100, 0.95, 2, false, 0));
    print_result( boost.calc_error( &data, CV_TRAIN_ERROR ), boost.calc_error( &data ), 0 );
#endif

    printf("======RTREES=====\n");
    rtrees.train( &data, CvRTParams( 10, 2, 0, false, 16, 0, true, 0, 100, 0, CV_TERMCRIT_ITER ));
    print_result( rtrees.calc_error( &data, CV_TRAIN_ERROR), rtrees.calc_error( &data ), rtrees.get_var_importance() );

    printf("======ERTREES=====\n");
    ertrees.train( &data, CvRTParams( 10, 2, 0, false, 16, 0, true, 0, 100, 0, CV_TERMCRIT_ITER ));
    print_result( ertrees.calc_error( &data, CV_TRAIN_ERROR), ertrees.calc_error( &data ), ertrees.get_var_importance() );

    return 0;
}

void Model::Predict_rtrees( const SampleSet& samples, SampleSet& outError )
{	
	int true_resp = 0;
	CvRTrees *model = (CvRTrees*)m_pModel;
	
	for (int i = 0; i < samples.N(); i++)
	{
		float ret = model->predict(samples.GetSampleAt(i), cv::Mat());
		if (ret != samples.GetLabelAt(i))
		{
			outError.Add(samples.GetSampleAt(i), samples.GetLabelAt(i));
		}
		else
		{
			true_resp++;
		}
	}
	printf("%d %d",samples.N(), true_resp);
}

CvRTrees* train_rf(CvMat* predictors, CvMat* labels)
{
	int stat[2];
	get_stat(labels, stat);
	printf("%d negative samples, %d positive samples\n", stat[0], stat[1]);
	
	const int tree_count = 500;
	const float priors[] = {0.25f,0.75f};
	CvRTrees* rtrees = new CvRTrees();
	CvRTParams rtparams = CvRTParams(5, 10, 0, false, 2, priors, true, 
									 (int)sqrt((float)predictors->cols), tree_count, 1e-6, 
									 CV_TERMCRIT_ITER + CV_TERMCRIT_EPS);
	CvMat* var_type = cvCreateMat(predictors->cols + 1, 1, CV_8UC1);
	for(int i = 0; i < predictors->cols; i++)
	{
		*(int*)(var_type->data.ptr + i*var_type->step) = CV_VAR_NUMERICAL;
	}
	*(int*)(var_type->data.ptr + predictors->cols*var_type->step) = CV_VAR_CATEGORICAL;
	rtrees->train(predictors, CV_ROW_SAMPLE, labels, 0, 0, var_type, 0, rtparams);
	return rtrees;
}

/* Examines the values at each leaf node in order to see what the distribution of data
  we put in is doing */
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
	ASSERT_NUM_RHS_ARGS_EQUALS(1);
	
	const mxArray* forest_ptr = prhs[0];
	ASSERT_IS_POINTER(forest_ptr);
	CvRTrees *forest = (CvRTrees *) unpack_pointer(forest_ptr);
	
	// We are going to return a cell array with one cell per tree, so need this number
	int num_trees = forest->get_tree_count();
	mexPrintf("Loaded forest of %d trees, retrieving leave node values.\n", num_trees);

	mxArray *output_cell_array = mxCreateCellMatrix(1, num_trees);
	ASSERT_NON_NULL(output_cell_array);
	
	for (unsigned int t = 0; t < num_trees; t++) {
		mxArray* tree_struct = mxCreateStructArray(num_dims, dims, tree_num_fields, tree_field_names);
		ASSERT_NON_NULL(tree_struct);
		mxSetCell(output_cell_array, t, make_matlab_tree_struct(forest->get_tree(t)));
	}
	plhs[0] = output_cell_array;
}

void find_decision_boundary_RF()
{
    img.copyTo( imgDst );

    Mat trainSamples, trainClasses;
    prepare_train_data( trainSamples, trainClasses );

    // learn classifier
    CvRTrees  rtrees;
    CvRTParams  params( 4, // max_depth,
                        2, // min_sample_count,
                        0.f, // regression_accuracy,
                        false, // use_surrogates,
                        16, // max_categories,
                        0, // priors,
                        false, // calc_var_importance,
                        1, // nactive_vars,
                        5, // max_num_of_trees_in_the_forest,
                        0, // forest_accuracy,
                        CV_TERMCRIT_ITER // termcrit_type
                       );

    rtrees.train( trainSamples, CV_ROW_SAMPLE, trainClasses, Mat(), Mat(), Mat(), Mat(), params );

    Mat testSample(1, 2, CV_32FC1 );
    for( int y = 0; y < img.rows; y += testStep )
    {
        for( int x = 0; x < img.cols; x += testStep )
        {
            testSample.at<float>(0) = (float)x;
            testSample.at<float>(1) = (float)y;

            int response = (int)rtrees.predict( testSample );
            circle( imgDst, Point(x,y), 2, classColors[response], 1 );
        }
    }
}

static
int build_rtrees_classifier( char* data_filename,
    char* filename_to_save, char* filename_to_load )
{
    CvMat* data = 0;
    CvMat* responses = 0;
    CvMat* var_type = 0;
    CvMat* sample_idx = 0;

    int ok = read_num_class_data( data_filename, 16, &data, &responses );
    int nsamples_all = 0, ntrain_samples = 0;
    int i = 0;
    double train_hr = 0, test_hr = 0;
    CvRTrees forest;
    CvMat* var_importance = 0;

    if( !ok )
    {
        printf( "Could not read the database %s\n", data_filename );
        return -1;
    }

    printf( "The database %s is loaded.\n", data_filename );
    nsamples_all = data->rows;
    ntrain_samples = (int)(nsamples_all*0.8);

    // Create or load Random Trees classifier
    if( filename_to_load )
    {
        // load classifier from the specified file
        forest.load( filename_to_load );
        ntrain_samples = 0;
        if( forest.get_tree_count() == 0 )
        {
            printf( "Could not read the classifier %s\n", filename_to_load );
            return -1;
        }
        printf( "The classifier %s is loaded.\n", data_filename );
    }
    else
    {
        // create classifier by using <data> and <responses>
        printf( "Training the classifier ...\n");

        // 1. create type mask
        var_type = cvCreateMat( data->cols + 1, 1, CV_8U );
        cvSet( var_type, cvScalarAll(CV_VAR_ORDERED) );
        cvSetReal1D( var_type, data->cols, CV_VAR_CATEGORICAL );

        // 2. create sample_idx
        sample_idx = cvCreateMat( 1, nsamples_all, CV_8UC1 );
        {
            CvMat mat;
            cvGetCols( sample_idx, &mat, 0, ntrain_samples );
            cvSet( &mat, cvRealScalar(1) );

            cvGetCols( sample_idx, &mat, ntrain_samples, nsamples_all );
            cvSetZero( &mat );
        }

        // 3. train classifier
        forest.train( data, CV_ROW_SAMPLE, responses, 0, sample_idx, var_type, 0,
            CvRTParams(10,10,0,false,15,0,true,4,100,0.01f,CV_TERMCRIT_ITER));
        printf( "\n");
    }

    // compute prediction error on train and test data
    for( i = 0; i < nsamples_all; i++ )
    {
        double r;
        CvMat sample;
        cvGetRow( data, &sample, i );

        r = forest.predict( &sample );
        r = fabs((double)r - responses->data.fl[i]) <= FLT_EPSILON ? 1 : 0;

        if( i < ntrain_samples )
            train_hr += r;
        else
            test_hr += r;
    }

    test_hr /= (double)(nsamples_all-ntrain_samples);
    train_hr /= (double)ntrain_samples;
    printf( "Recognition rate: train = %.1f%%, test = %.1f%%\n",
            train_hr*100., test_hr*100. );

    printf( "Number of trees: %d\n", forest.get_tree_count() );

    // Print variable importance
    var_importance = (CvMat*)forest.get_var_importance();
    if( var_importance )
    {
        double rt_imp_sum = cvSum( var_importance ).val[0];
        printf("var#\timportance (in %%):\n");
        for( i = 0; i < var_importance->cols; i++ )
            printf( "%-2d\t%-4.1f\n", i,
            100.f*var_importance->data.fl[i]/rt_imp_sum);
    }

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

int main()
{
    const int train_sample_count = 300;
    bool is_regression = false;

    const char* filename = "data/waveform.data";
    int response_idx = 21;

    CvMLData data;

    CvTrainTestSplit spl( train_sample_count );
    
    if(data.read_csv(filename) != 0)
    {
        printf("couldn't read %s\n", filename);
        exit(0);
    }

    data.set_response_idx(response_idx);
    data.change_var_type(response_idx, CV_VAR_CATEGORICAL);
    data.set_train_test_split( &spl );

    const CvMat* values = data.get_values();
    const CvMat* response = data.get_responses();
    const CvMat* missing = data.get_missing();
    const CvMat* var_types = data.get_var_types();
    const CvMat* train_sidx = data.get_train_sample_idx();
    const CvMat* var_idx = data.get_var_idx();
    CvMat*response_map;
    CvMat*ordered_response = cv_preprocess_categories(response, var_idx, response->rows, &response_map, NULL);
    int num_classes = response_map->cols;
    
    CvDTree dtree;
    printf("======DTREE=====\n");
    CvDTreeParams cvd_params( 10, 1, 0, false, 16, 0, false, false, 0);
    dtree.train( &data, cvd_params);
    print_result( dtree.calc_error( &data, CV_TRAIN_ERROR), dtree.calc_error( &data, CV_TEST_ERROR ), dtree.get_var_importance() );

#if 0
    /* boosted trees are only implemented for two classes */
    printf("======BOOST=====\n");
    CvBoost boost;
    boost.train( &data, CvBoostParams(CvBoost::DISCRETE, 100, 0.95, 2, false, 0));
    print_result( boost.calc_error( &data, CV_TRAIN_ERROR ), boost.calc_error( &data, CV_TEST_ERROR), 0 );
#endif

    printf("======RTREES=====\n");
    CvRTrees rtrees;
    rtrees.train( &data, CvRTParams( 10, 2, 0, false, 16, 0, true, 0, 100, 0, CV_TERMCRIT_ITER ));
    print_result( rtrees.calc_error( &data, CV_TRAIN_ERROR), rtrees.calc_error( &data, CV_TEST_ERROR ), rtrees.get_var_importance() );

    printf("======ERTREES=====\n");
    CvERTrees ertrees;
    ertrees.train( &data, CvRTParams( 10, 2, 0, false, 16, 0, true, 0, 100, 0, CV_TERMCRIT_ITER ));
    print_result( ertrees.calc_error( &data, CV_TRAIN_ERROR), ertrees.calc_error( &data, CV_TEST_ERROR ), ertrees.get_var_importance() );

    printf("======GBTREES=====\n");
    CvGBTrees gbtrees;
    CvGBTreesParams gbparams;
    gbparams.loss_function_type = CvGBTrees::DEVIANCE_LOSS; // classification, not regression
    gbtrees.train( &data, gbparams);
    
    //gbt_print_error(&gbtrees, values, response, response_idx, train_sidx);
    print_result( gbtrees.calc_error( &data, CV_TRAIN_ERROR), gbtrees.calc_error( &data, CV_TEST_ERROR ), 0);

    printf("======KNEAREST=====\n");
    CvKNearest knearest;
    //bool CvKNearest::train( const Mat& _train_data, const Mat& _responses,
    //                const Mat& _sample_idx, bool _is_regression,
    //                int _max_k, bool _update_base )
    bool is_classifier = var_types->data.ptr[var_types->cols-1] == CV_VAR_CATEGORICAL;
    assert(is_classifier);
    int max_k = 10;
    knearest.train(values, response, train_sidx, is_regression, max_k, false);

    CvMat* new_response = cvCreateMat(response->rows, 1, values->type);
    //print_types();

    //const CvMat* train_sidx = data.get_train_sample_idx();
    knearest.find_nearest(values, max_k, new_response, 0, 0, 0);

    print_result(knearest_calc_error(values, response, new_response, train_sidx, is_regression, CV_TRAIN_ERROR),
                 knearest_calc_error(values, response, new_response, train_sidx, is_regression, CV_TEST_ERROR), 0);

    printf("======== RBF SVM =======\n");
    //printf("indexes: %d / %d, responses: %d\n", train_sidx->cols, var_idx->cols, values->rows);
    CvMySVM svm1;
    CvSVMParams params1 = CvSVMParams(CvSVM::C_SVC, CvSVM::RBF,
                                     /*degree*/0, /*gamma*/1, /*coef0*/0, /*C*/1,
                                     /*nu*/0, /*p*/0, /*class_weights*/0,
                                     cvTermCriteria(CV_TERMCRIT_ITER+CV_TERMCRIT_EPS, 1000, FLT_EPSILON));
    //svm1.train(values, response, train_sidx, var_idx, params1);
    svm1.train_auto(values, response, var_idx, train_sidx, params1);
    svm_print_error(&svm1, values, response, response_idx, train_sidx);

    printf("======== Linear SVM =======\n");
    CvMySVM svm2;
    CvSVMParams params2 = CvSVMParams(CvSVM::C_SVC, CvSVM::LINEAR,
                                     /*degree*/0, /*gamma*/1, /*coef0*/0, /*C*/1,
                                     /*nu*/0, /*p*/0, /*class_weights*/0,
//.........这里部分代码省略.........

int RandomTrees::train(const char* samples_filename, const char* model_filename, const double ratio, double &train_error, double &test_error)
{
	CvMat* data = 0;
	CvMat* responses = 0;
	CvMat* var_type = 0;
	CvMat* sample_idx = 0;

	this->tree_parameters_.nactive_vars = (int)sqrt(this->number_of_features_);

	int ok = read_num_class_data( samples_filename, this->number_of_features_, &data, &responses );
	int nsamples_all = 0, ntrain_samples = 0;
	int i = 0;
	double train_hr = 0, test_hr = 0;
	CvRTrees forest;
	CvMat* var_importance = 0;

	if( !ok )
	{
		cout << "Could not read the sample in" << samples_filename << endl;;
		return -1;
	}

	cout << "The sample file " << samples_filename << " is loaded." << endl;
	nsamples_all = data->rows;
	ntrain_samples = (int)(nsamples_all * ratio);


	// create classifier by using <data> and <responses>
	cout << "Training the classifier ..." << endl;

	// 1. create type mask
	var_type = cvCreateMat( data->cols + 1, 1, CV_8U );
	cvSet( var_type, cvScalarAll(CV_VAR_ORDERED) );
	cvSetReal1D( var_type, data->cols, CV_VAR_CATEGORICAL );

	// 2. create sample_idx
	sample_idx = cvCreateMat( 1, nsamples_all, CV_8UC1 );
	{
		CvMat mat;
		cvGetCols( sample_idx, &mat, 0, ntrain_samples );
		cvSet( &mat, cvRealScalar(1) );

		cvGetCols( sample_idx, &mat, ntrain_samples, nsamples_all );
		cvSetZero( &mat );
	}

	// 3. train classifier
	forest.train( data, CV_ROW_SAMPLE, responses, 0, sample_idx, var_type, 0, this->tree_parameters_);
	cout << endl;


	// compute prediction error on train and test data
	for( i = 0; i < nsamples_all; i++ )
	{
		double r;
		CvMat sample;
		cvGetRow( data, &sample, i );

		r = forest.predict( &sample );
		r = fabs((double)r - responses->data.fl[i]) <= FLT_EPSILON ? 1 : 0;

		if( i < ntrain_samples )
			train_hr += r;
		else
			test_hr += r;
	}

	test_hr /= (double)(nsamples_all-ntrain_samples);
	train_hr /= (double)ntrain_samples;

	train_error = 1 - train_hr;
	test_error = 1 - test_hr;

	cout << "Recognition rate: train = " << train_hr*100 << ", test = " << test_hr*100 << endl;
	cout << "Number of trees: " << forest.get_tree_count() << endl;

	// Print variable importance
	var_importance = (CvMat*)forest.get_var_importance();
	if( var_importance )
	{
		double rt_imp_sum = cvSum( var_importance ).val[0];
		printf("var#\timportance (in %%):\n");
		for( i = 0; i < var_importance->cols; i++ )
			printf( "%-2d\t%-4.1f\n", i,100.f*var_importance->data.fl[i]/rt_imp_sum);
	}

	// Save Random Trees classifier to file if needed
	if( model_filename )
		forest.save( model_filename );

	//cvReleaseMat( &var_importance );		//causes a segmentation fault
	cvReleaseMat( &sample_idx );
	cvReleaseMat( &var_type );
	cvReleaseMat( &data );
	cvReleaseMat( &responses );

	return 0;
}

int main(int argc, char** argv)
{
  // std::cout<<FLT_EPSILON<<std::endl; 
  cv::Mat training_data, training_labels,testing_data, testing_labels;
  
  training_data = read_rgbd_data_cv(argv[1],NUMBER_OF_TRAINING_SAMPLES);
  training_labels = read_rgbd_data_cv(argv[2], NUMBER_OF_TRAINING_SAMPLES);
  testing_data = read_rgbd_data_cv(argv[3],NUMBER_OF_TESTING_SAMPLES);
  testing_labels = read_rgbd_data_cv(argv[4], NUMBER_OF_TESTING_SAMPLES);
  
 
  printf("dataset specs: %d samples with %d features\n", training_data.rows, training_data.cols);

  // define all the attributes as numerical
  // alternatives are CV_VAR_CATEGORICAL or CV_VAR_ORDERED(=CV_VAR_NUMERICAL)
  // that can be assigned on a per attribute basis

  cv::Mat var_type = cv::Mat(training_data.cols + 1, 1, CV_8U );
  var_type.setTo(cv::Scalar(CV_VAR_NUMERICAL) ); // all inputs are numerical
  var_type.at<uchar>(training_data.cols, 0) = CV_VAR_CATEGORICAL; // the labels are categorical

  /********************************步骤1：定义初始化Random Trees的参数******************************/
  float priors[] = {1,1,1,1,1};  // weights of each classification for classes
  CvRTParams params = CvRTParams(25, // max depth
				 50, // min sample count
				 0, // regression accuracy: N/A here
				 false, // compute surrogate split, no missing data
				 15, // max number of categories (use sub-optimal algorithm for larger numbers)
				 priors, // the array of priors
				 false,  // calculate variable importance
				 20,       // number of variables randomly selected at node and used to find the best split(s).
				 NUMBER_OF_TREES,	 // max number of trees in the forest
				 0.01f,				// forrest accuracy
				 CV_TERMCRIT_ITER |	CV_TERMCRIT_EPS // termination cirteria
				 );
  
  /****************************步骤2：训练 Random Decision Forest(RDF)分类器*********************/
  // printf( "\nUsing training database: %s\n\n", argv[1]);
  CvRTrees* rtree = new CvRTrees;
  rtree->train(training_data, CV_ROW_SAMPLE, training_labels,
	       cv::Mat(), cv::Mat(), var_type, cv::Mat(), params);
  
  // perform classifier testing and report results
  cv::Mat test_sample, train_sample;
  int correct_class = 0;
  int wrong_class = 0;
  int result;
  int label;
  int false_positives [NUMBER_OF_CLASSES] = {0,0,0,0,0};
  int false_negatives [NUMBER_OF_CLASSES] = {0,0,0,0,0};

  CvDTreeNode* leaf_nodes [training_data.rows];

  for (int tsample = 0; tsample < training_data.rows; tsample++)
    {
      train_sample = training_data.row(tsample);
      CvForestTree* tree = rtree->get_tree(1);
      CvDTreeNode* leaf_node = tree->predict(train_sample, cv::Mat());
      leaf_nodes[tsample] = leaf_node; 
    }

  // printf( "\nUsing testing database: %s\n\n", argv[2]);

  for (int tsample = 0; tsample < testing_data.rows; tsample++)
    {	       
      // extract a row from the testing matrix
      test_sample = testing_data.row(tsample);
      // train on the testing data:
      // test_sample = training_data.row(tsample);
      /********************************步骤3：预测*********************************************/

      result = (int) rtree->predict(test_sample, cv::Mat());
      label = (int) testing_labels.at<float>(tsample, 0);

      printf("Testing Sample %i -> class result (digit %d) - label (digit %d)\n", tsample, result, label);

      // get the leaf nodes of the first tree in the forest
      /*CvForestTree* tree = rtree->get_tree(0);
      std::list<const CvDTreeNode*> leaf_list;
      leaf_list = get_leaf_node( tree );
      printf("Number of Leaf nodes: %ld\n", leaf_list.size());*/

      // if the prediction and the (true) testing classification are the same
      // (N.B. openCV uses a floating point decision tree implementation!)
      if (fabs(result - label)
	  >= FLT_EPSILON)
	{
	  // if they differ more than floating point error => wrong class
	  wrong_class++;
	  false_positives[(int) result]++;
	  false_negatives[(int) testing_labels.at<float>(tsample, 0)]++;
	}
      else
	{
	  // otherwise correct
	  correct_class++;
	}
    }

  printf( // "\nResults on the testing database: %s\n"
//.........这里部分代码省略.........

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
    ASSERT_NUM_RHS_ARGS_GTE(2);
    ASSERT_NUM_LHS_ARGS_LT(3);

    const mxArray* dataMtx = prhs[0];
    const mxArray* targetValueVec = prhs[1];
    
    //see if we have been provided a struct containing options for the training. 
    //if not, then use defaults provided by opencv
    CvRTParams* rtParams;
    if (nrhs > 2) {
        mexPrintf("Parsing struct argument for parameters\n");
        rtParams = parse_struct_to_forest_config(prhs[2]);
    }
    else {
        mexPrintf("Using default parameters\n");
        rtParams = parse_struct_to_forest_config(NULL);
    }

    mexPrintf("Parameters:\n");
    print_forest_params(rtParams);
    
    unsigned int numSamples, numVariables;

    CvMat* dataCvMtx = matlab_matrix_to_opencv_matrix(dataMtx);
    numSamples = dataCvMtx->rows;
    numVariables = dataCvMtx->cols;
    mexPrintf("training data converted to opencv format. %d samples, each with %d variables\n",
              numSamples, numVariables);
#ifdef PRINT_INPUTS
    print_opencv_matrix(dataCvMtx);
#endif

    CvMat* targetCvMtx = matlab_array_to_opencv_array(targetValueVec);
    if (targetCvMtx->rows != numSamples) {
		MEX_ERR_PRINTF("training data had %d samples, labels contain %d values.", 
		               numSamples, targetCvMtx->rows);
    }
    mexPrintf("training labels converted to opencv format.\n");
#ifdef PRINT_INPUTS
    print_opencv_matrix(targetCvMtx);
#endif

    //specify the type of our variables. In this case, all our variables are 
    CvMat* var_type = cvCreateMat(dataCvMtx->cols + 1, 1, CV_8U);
    cvSet(var_type, cvScalarAll(CV_VAR_ORDERED));

    //actually make the forest and do the training
    clock_t start_time, end_time;
    mexPrintf("training now...");
    start_time = clock();
    CvRTrees *forest = new CvRTrees;
    forest->train(dataCvMtx, CV_ROW_SAMPLE, targetCvMtx, NULL, NULL, var_type, NULL, *rtParams);
    end_time = clock();
	clock_t diff_time = end_time - start_time;
	double seconds_passed = ((float)diff_time) / CLOCKS_PER_SEC;
    mexPrintf("training done in %fs\n", seconds_passed);

    //pack the pointer and return it to matlab
    plhs[0] = pack_pointer((void *)forest);

	// If the user supplied a second lhs argument, return them the time taken to train
	if (nlhs > 1) {
		plhs[1] = mxCreateDoubleScalar(seconds_passed);
	}
    
    cvReleaseMat(&var_type);
    cvReleaseMat(&dataCvMtx);
    cvReleaseMat(&targetCvMtx);
} 

void evaluation(CvRTrees& forest, DataSet& data, Mat& sampleIdx, TrainResult& result)
{

	int numTrainSamples = (cv::sum( sampleIdx ))[0];

	// retrieve variable_importance
	result.var_importance = forest.get_var_importance();
//	result.var_importance = forest.get_subtree_weights();
//	cout << result.var_importance << endl;

	double min,max;
	Point minLoc,maxLoc;

	minMaxLoc(result.var_importance,&min,&max,&minLoc,&maxLoc);
//	printf("variable importance (max:%.2f%%):\n\n",max*100.f);

	// compute prediction error on train and test data
	result.train_hr = 0; result.test_hr = 0; result.fpRate = 0; result.fnRate = 0;


	Mat responses_new = Mat(data.numSamples,1,CV_32F,9.0);

	for(int i = 0; i < data.numSamples; i++ )
	{
		double r;
		Mat sample = data.data.row(i);

		// do prediction with trained forest
		r = forest.predict(sample);
		responses_new.at<float>(i,0) = r;
		float respo = data.responses.at<float>(i,0);

		// prediction correct ?
		r = fabs(r - respo) <= FLT_EPSILON ? 1 : 0;


		if( sampleIdx.at<char>(0,i) )
			result.train_hr += r;
		else
			result.test_hr += r;

		// false prediction, increase appropriate counter
		if(!r)
		{
			if(respo)
				result.fnRate += 1;
			else
				result.fpRate += 1;
		}
	}
//	cout << sampleIdx << endl;
//	cout << data.responses << endl;
//	cout << responses_new << endl;

	result.test_hr /= (double)(data.numSamples-numTrainSamples);
	result.train_hr /= (double)numTrainSamples;

	result.fpRate /= (double) data.numNeg;
	result.fnRate /= (double) data.numPos;

}

//.........这里部分代码省略.........
				// this is a classification problem (i.e. predict a discrete number of class
				// outputs) so reset the last (+1) output var_type element to CV_VAR_CATEGORICAL

				var_type.at<uchar>(ATTRIBUTES_PER_SAMPLE, 0) = CV_VAR_CATEGORICAL;

				double result; // value returned from a prediction

				// load training and testing data sets

				if (read_data_from_csv(argv[1], training_data, training_classifications, NUMBER_OF_TRAINING_SAMPLES) &&
						read_data_from_csv(argv[2], testing_data, testing_classifications, NUMBER_OF_TESTING_SAMPLES))
				{
					// define the parameters for training the random forest (trees)

					float priors[] = {1,1,1,1,1,1,1,1,1,1,1,1,1,1,1};  // weights of each classification for classes
					// (all equal as equal samples of each digit)

					CvRTParams params = CvRTParams(20, // max depth
												5, // min sample count
												0, // regression accuracy: N/A here
												false, // compute surrogate split, no missing data
												15, // max number of categories (use sub-optimal algorithm for larger numbers)
												priors, // the array of priors
												false,  // calculate variable importance
												40,       // number of variables randomly selected at node and used to find the best split(s).
												100,	 // max number of trees in the forest
												0.01f,				// forrest accuracy
												CV_TERMCRIT_ITER |	CV_TERMCRIT_EPS // termination cirteria
												);

					// train random forest classifier (using training data)

					printf( "\nUsing training database: %s\n\n", argv[1]);
					CvRTrees* rtree = new CvRTrees;

					rtree->train(training_data, CV_ROW_SAMPLE, training_classifications,
								 Mat(), Mat(), var_type, Mat(), params);

					// perform classifier testing and report results

					Mat test_sample;
					int correct_class = 0;
					int wrong_class = 0;
					int false_positives [NUMBER_OF_CLASSES] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};

					printf( "\nUsing testing database: %s\n\n", argv[2]);

					for (int tsample = 0; tsample < NUMBER_OF_TESTING_SAMPLES; tsample++)
					{

						// extract a row from the testing matrix

						test_sample = testing_data.row(tsample);

						// run random forest prediction

						result = rtree->predict(test_sample, Mat());

						printf("Testing Sample %i -> class result (digit %d)\n", tsample, (int) result);

						// if the prediction and the (true) testing classification are the same
						// (N.B. openCV uses a floating point decision tree implementation!)

						if (fabs(result - testing_classifications.at<float>(tsample, 0))
								>= FLT_EPSILON)
						{

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

		tex_imgs[i - START][0] = flair_tex;
		tex_imgs[i - START][1] =t1_tex;
		tex_imgs[i - START][2] = t1c_tex;
		tex_imgs[i - START][3] = t2_tex;
	}
	//----------------------------------------------------------
	cout << "read training data............" << endl;
	Mat train_datas(HEIGHT*WIDTH*(END - START + 1), ATTRIBUTES_PER_SAMPLE, CV_32FC1);    
	Mat responses(HEIGHT*WIDTH*(END - START + 1), 1, CV_32SC1);
	//---读取训练数据----
	int dataline=read_training_data(imgs,tex_imgs, train_datas, responses);
	Mat _train_datas(dataline, ATTRIBUTES_PER_SAMPLE, CV_32FC1);
	Mat _responses(dataline, 1, CV_32SC1);
	//减少训练数据为dataline个
	for (int i = 0; i < dataline; i++)
	{
		float* float_data = train_datas.ptr<float>(i);
		int* int_data = responses.ptr<int>(i);

		_train_datas.at<float>(i, 0) = float_data[0];
		_train_datas.at<float>(i, 1) = float_data[1];
		_train_datas.at<float>(i, 2) = float_data[2];
		_train_datas.at<float>(i, 3) = float_data[3];

		_train_datas.at<float>(i, 4) = float_data[4];
		_train_datas.at<float>(i, 5) = float_data[5];
		_train_datas.at<float>(i, 6) = float_data[6];
		_train_datas.at<float>(i, 7) = float_data[7];

		_train_datas.at<float>(i, 8) = float_data[8];

		_responses.at<int>(i, 0) = int_data[0];
	}
	//----设置输入类型---
	Mat var_type = Mat(ATTRIBUTES_PER_SAMPLE+1, 1, CV_8U);
	var_type.setTo(Scalar(CV_VAR_NUMERICAL)); // all inputs are numerical  
	var_type.at<uchar>(ATTRIBUTES_PER_SAMPLE, 0) = CV_VAR_CATEGORICAL;
	//---训练数据---
	cout << "training......." << endl;
	float priors[NUMBER_OF_CLASSES] = { 1, 1 };
	CvRTParams params = CvRTParams(25, // max depth  
	                       4, // min sample count  
	                       0, // regression accuracy: N/A here  
                           false, // compute surrogate split, no missing data  
		                   5, // max number of categories (use sub-optimal algorithm for larger numbers)  
	                        priors, // the array of priors  
	                        false,  // calculate variable importance  
		                    3,       // number of variables randomly selected at node and used to find the best split(s).  
							3,  // max number of trees in the forest  
                            0.01f,               // forrest accuracy  
	                        CV_TERMCRIT_ITER | CV_TERMCRIT_EPS // termination cirteria  
		);
	CvRTrees* rtree = new CvRTrees;
	bool train_result = rtree->train(_train_datas, CV_ROW_SAMPLE, _responses,
		Mat(), Mat(), var_type, Mat(), params);
	if (train_result == false)
		cout << "random trees train failed!" << endl;
	
	cout << "predicting.........." << endl;
	//-------预测数据生成图片并存储---------
	for (int k = 0; k < END - START + 1; k++)
	{
		IplImage* img_dst = cvCreateImage(cvGetSize(imgs[k][0]), IPL_DEPTH_8U, 1);
		uchar* ptr;
		for (int i = 0; i <HEIGHT ; i++)
		{
			ptr = (uchar*)img_dst->imageData + i*img_dst->widthStep;
			for (int j = 0; j < WIDTH; j++)
			{//读一行数据
				Mat test_data(1, ATTRIBUTES_PER_SAMPLE, CV_32FC1);
				test_data.at<float>(0, 0) = cvGet2D(imgs[k][0], i, j).val[0];
				test_data.at<float>(0, 1) = cvGet2D(imgs[k][1], i, j).val[0];
				test_data.at<float>(0, 2) = cvGet2D(imgs[k][2], i, j).val[0];
				test_data.at<float>(0, 3) = cvGet2D(imgs[k][3], i, j).val[0];
			
				test_data.at<float>(0, 4) = cvGet2D(tex_imgs[k][0], i, j).val[0];
				test_data.at<float>(0, 5) = cvGet2D(tex_imgs[k][1], i, j).val[0];
				test_data.at<float>(0, 6) = cvGet2D(tex_imgs[k][2], i, j).val[0];
				test_data.at<float>(0, 7) = cvGet2D(tex_imgs[k][3], i, j).val[0];

				test_data.at<float>(0, 8) = k;
				//产生结果
				int result = rtree->predict(test_data, Mat());
				*(ptr + j) = result * 255;
			}
		}
		IplConvKernel* strel = cvCreateStructuringElementEx(5, 5, 2, 2, CV_SHAPE_ELLIPSE);
		cvErode(img_dst, img_dst, strel, 1);
		//cvDilate(img_dst, img_dst, strel, 1);
		
		cout << "save image  " << k + START << endl;
		char result_name[100];
		memset(result_name, 0, 100);
		sprintf(result_name, "BRATS_HG0005_RESULT/BRATS_HG0005_RESULT_%d.png", k+START);
		cvSaveImage(result_name, img_dst);
	}
	cout << "complete!!!" << endl;
	cvWaitKey(0);
}

//.........这里部分代码省略.........
    for (Int_t i = 0; i < TRAIN; i++)
    {
	for (Int_t j = 0; j < VARS; j++)
	    cout << train[i][j] << "\t";
	cout << endl;
    }
*/
    // Create type mask
    Int_t var[VARS + 1];
    for (Int_t i = 0; i < VARS; i++)
    	var[i] = CV_VAR_ORDERED;
    var[VARS] = CV_VAR_CATEGORICAL;
    Mat var_type(VARS + 1, 1, CV_32SC1, var);
    var_type.convertTo(var_type, CV_8SC1);	// Convert to 8-bit ints
    
    // Create missing data mask
    Int_t miss_t[TRAIN][VARS];
    for (Int_t i = 0; i < TRAIN; i++)
    {
	for (Int_t j = 0; j < VARS; j++)
	    miss_t[i][j] = 0;
    }
    Mat missing_data_mask(TRAIN, VARS, CV_32SC1, miss_t);
    missing_data_mask.convertTo(missing_data_mask, CV_8UC1);

    // Create indices
    Mat var_idx = Mat::ones(VARS, 1, CV_8UC1);
    Mat sample_idx = Mat::ones(TRAIN, 1, CV_8UC1);

    // Train forest, print variable importance (if used)
    cout << "Trees: " << max_trees << endl;
    cout << "Depth: " << max_depth << endl;
    cout << "m: " << nactive_vars << endl;
    CvRTrees forest;
    forest.train(train_data, tflag, responses, var_idx, sample_idx, var_type, missing_data_mask, CvRTParams(max_depth, min_sample_count, regression_accuracy, use_surrogates, max_categories, priors, calc_var_importance, nactive_vars, max_trees, forest_accuracy, termcrit_type));
    if (calc_var_importance)
    {
	Mat imp = forest.getVarImportance();
	cout << endl << imp << endl << endl;
    }

    // Create solving array and data mask
    Int_t solve_good = event[0] - train_good;
    Int_t solve_bad = event[1] - train_bad;
    const Int_t SOLVE = solve_good + solve_bad;
    Int_t flag[SOLVE];
    Float_t solve[SOLVE][VARS];
    for (Int_t i = 0; i < solve_good; i++)
    {
	tree[0]->GetEvent(i + train_good);
	flag[i] = 0;
//	solve[i][0] = tibHits[0];
//	solve[i][1] = tidHits[1];
//	solve[i][2] = tobHits[1];
//	solve[i][3] = tecHits[0];
	solve[i][0] = ptR[0];
	solve[i][1] = etaR[0];
	solve[i][2] = phiR[0];
	solve[i][3] = foundR[0];
    }
    for (Int_t i = 0; i < solve_bad; i++)
    {
	tree[1]->GetEvent(i + train_bad);
	flag[i + solve_good] = 1;
//	solve[i + solve_good][0] = tibHits[1];
//	solve[i + solve_good][1] = tidHits[1];