pdfbthmt.cc

matlab官方网站中的用于图像融合技术的contourlet变换源代码

	scaling[J][i] = sum;

	for (m = 0; m < M; m++) 
	    beta[J][i][m] /= sum;
    }
}

//-----------------------------------------------------------------------------
void THMT::compute_alpha()
{
    register int J, i, m, mm, n, nNode;
    register double sum;

    // Initialize the coarsest level
    for (m = 0; m < M; m++)
	alpha[0][0][m] = model_trans[0][0][m][0];

    for (J = 1; J < nLev; J++) 
    {
	for (i = 0, nNode = alpha[J].size(); i < nNode; i++) 
	{
	  mm = subbandtree[J][i];
	  for (m = 0; m < M; m++) 
	    {
	      sum = 0.0;

	      for (n = 0; n < M; n++) 
		sum += model_trans[J][mm][m][n] *
		  alpha[J-1][i/nCh][n] *
		  beta[J-1][i/nCh][n] /
		  beta_par[J][i][n];

		alpha[J][i][m] = sum;
	    } // end for m
	} // end for i
    } // end for J
}


//-----------------------------------------------------------------------------
double THMT::compute_likelihood()
{
    register int J, i, nNode, m;
    register double f;
    register double log_scale;

    for (m = 0, f = 0.0; m < M; m++)
	f += beta[0][0][m] * model_trans[0][0][m][0];

    // Re-scale back using saved scaling factors
    for (J = 0, log_scale = 0.0; J < nLev; J++)
	for (i = 0, nNode = scaling[J].size(); i < nNode; i++) 
	    log_scale += log(scaling[J][i]);

    // Final result
    return (log_scale + log(f));
}


//-----------------------------------------------------------------------------
void THMT::update_model(double& delta, double& avf)
{
    register int J, MM, i, m, mm, n, nNode;
    register int ob_ind;           // observation index
    register double o, denom, prob1, prob2, newval;

    // Reinitialize
    delta = avf = 0.0;

    for (J = 0; J < nLev; J++)
      for (mm=0; mm<sum_prob[J].size(); mm++)
	for (m = 0; m < M; m++) 
	{
	    sum_prob[J][mm][m] = 0.0;
	    sum_mean[J][mm][m] = 0.0;
	    sum_stdv[J][mm][m] = 0.0;

	    for (n = 0; n < M; n++)
	      sum_trans[J][mm][m][n] = 0.0;
	}
  
    // For each training tree
    for (ob_ind = 0; ob_ind < nObs; ob_ind++) 
    {
	// Compute $\alpha, \beta$
	compute_beta(ob_ind);
	compute_alpha();

	// Denominator for restimated probs
	denom = 0.0;
	for (m = 0; m < M; m++)
	    denom += beta[0][0][m] * alpha[0][0][m];

	// Compute state probabilities for denoising
	for (J = 0; J < nLev; J++)
	    for (i = 0, nNode = alpha[J].size(); i < nNode; i++)
	      for (m = 0; m < M; m++)
		{
		  state_prob[J][ob_ind*ipow(nCh,J)+i][m] = 
		    alpha[J][i][m]*beta[J][i][m]/denom;
		}

	// Update total log-likelihood
	avf += compute_likelihood();

	for (J = 0; J < nLev; J++)
       	    for (i = 0, nNode = alpha[J].size(); i < nNode; i++) 
      	    {
	        // o_{J,i}^{ob_ind}
       		o = (*obs)[J][ob_ind * ipow(nCh,J) + i];

       		for (m = 0; m < M; m++) 
    		{
		    // Compute $prob(S_i | O)$
       		    prob1 = alpha[J][i][m] * beta[J][i][m]
       			/ denom;

		    mm = subbandtree[J][i];
		    // Summing for all trees
       		    sum_prob[J][mm][m] += prob1;
       		    sum_mean[J][mm][m] += o * prob1;
       		    sum_stdv[J][mm][m] += (o - model_mean[J][mm][m]) *
       			(o - model_mean[J][mm][m]) *
       			prob1;
       		}
      	    }

	for (J = 1; J < nLev; J++)
	    for (i = 0, nNode = alpha[J].size(); i < nNode; i++)
		for (m = 0; m < M; m++)
		    for (n = 0; n < M; n++) 
		    {
		        mm = subbandtree[J][i];
			// Compute $prob(S_i, S_{p(i) | O)$
			prob2 = beta[J][i][m] *
			  model_trans[J][mm][m][n] *
			  alpha[J-1][i/nCh][n] *
			  beta[J-1][i/nCh][n] /
			  beta_par[J][i][n]
			  / denom;
			sum_trans[J][mm][m][n] += prob2;
		    }

    } // end ob_ind

    // Average log-likelihood
    avf = avf / double(nObs);

    // Normalize and update model parameters
    for (J = 0; J < nLev; J++)
      for (mm = 0; mm < model_stdv[J].size(); mm++)
	for (m = 0; m < M; m++) 
	{
	    if (!zeromean)	// Only update means for non-zeromean model
	    {
		newval = sum_mean[J][mm][m] / sum_prob[J][mm][m];
		delta += fabs(newval - model_mean[J][mm][m]);
		model_mean[J][mm][m] = newval;
	    }

	    newval = sqrt(sum_stdv[J][mm][m] / sum_prob[J][mm][m]);
	    delta += fabs(newval - model_stdv[J][mm][m]);
	    model_stdv[J][mm][m] = newval;
	}

    //state probs 
    for (m = 0; m < M; m++) 
    {
	newval = sum_prob[0][0][m] / nObs;
	delta += fabs(newval - model_trans[0][0][m][0]);
	model_trans[0][0][m][0] = newval;
    }

    // And transition probs 
    for (J = 1; J < nLev; J++){
      for (mm=0; mm < model_trans[J].size(); mm++)
	for (m = 0; m < M; m++)
	    for (n = 0; n < M; n++) 
	    {
	      if (sum_prob[J].size()>sum_prob[J-1].size())
		newval = sum_trans[J][mm][m][n] / sum_prob[J-1][mm/2][n]
		  /(nCh/2);
	      else
		newval =sum_trans[J][mm][m][n] / sum_prob[J-1][mm][n]/nCh;
	      delta += fabs(newval - model_trans[J][mm][m][n]);
	      model_trans[J][mm][m][n] = newval;
	    }
    }
}

//-----------------------------------------------------------------------------
void THMT::reorder_model()
{
  int level, dir, state, state2, largeststate, k;
  double largeststdv, tempdouble;

  // for each node
  for(level = 0; level < nLev; level++){
    for(dir = 0; dir < model_stdv[level].size(); dir++) {
      for(state=0; state<M-1; state++) {
	// initialize
	largeststdv = -1;
	largeststate = M;
	for(state2=state; state2<M; state2++){
	  // search for the state with the largest standard deviation
	  if (model_stdv[level][dir][state2]>largeststdv){
	    largeststdv = model_stdv[level][dir][state2];
	    largeststate = state2;
	  }
	}

	// if the current state is not the state with the largest stdv, then
	// need to swap the order
	if(largeststate != state){

	  // swap the order of the stdv
	  model_stdv[level][dir][largeststate] = model_stdv[level][dir][state];
	  model_stdv[level][dir][state] = largeststdv;

	  // swap the order of the mean
	  if (!zeromean) {
	    tempdouble = model_mean[level][dir][largeststate];
	    model_mean[level][dir][largeststate] = model_mean[level][dir]
	      [state];
	    model_mean[level][dir][state] = tempdouble;
	  }

	  // swap the order of the transition matrix with parent
	  for(k=0; k<M; k++){
	    tempdouble = model_trans[level][dir][largeststate][k];
	    model_trans[level][dir][largeststate][k] = model_trans[level][dir]
	      [state][k];
	    model_trans[level][dir][state][k] = tempdouble;
	  }

	  // swap the order of the transition matrix with children if this
	  // is not a leaf node
	  if(level != nLev-1){
	    if(model_trans[level].size()==model_trans[level+1].size())
	      for(k=0; k<M; k++){
		tempdouble = model_trans[level+1][dir][k][largeststate];
		model_trans[level+1][dir][k][largeststate] = 
		  model_trans[level+1][dir][k][state];
		model_trans[level+1][dir][k][state] = tempdouble;
	      }
	    else if (model_trans[level].size()*2==model_trans[level+1].size())
	      for (k=0; k<M; k++){
		tempdouble =model_trans[level+1][2*dir][k][largeststate];
		model_trans[level+1][2*dir][k][largeststate] = 
		  model_trans[level+1][2*dir][k][state];
		model_trans[level+1][2*dir][k][state] = tempdouble;

		tempdouble =model_trans[level+1][2*dir+1][k][largeststate];
		model_trans[level+1][2*dir+1][k][largeststate] = 
		  model_trans[level+1][2*dir+1][k][state];
		model_trans[level+1][2*dir+1][k][state] = tempdouble;
	      } 
	  }
	}
      }
    }   
  }

}

//-----------------------------------------------------------------------------
void THMT::batch_train(tree<double> *trainTree, double min_delta)
{
  
    // Assign data pointer
    obs = trainTree;
    // Train HMT model
    train_all(min_delta);
}


//-----------------------------------------------------------------------------
void THMT::train_all(double min_delta)
{
    register int count = 0, J, i, m, nNode;
    register double delta = min_delta;
    register double avf;
    register double last_avf = -10e6;
 
    if (obs->nlev() == 0)
      mexErrMsgTxt("ERROR in THMT::train_all(): empty training data");
    if ((obs->nlev() != nLev) || (obs->nch() != nCh))
      mexErrMsgTxt("ERROR in THMT::train_all(): incompatible training data");

    nObs = obs->nrt();

    // Allocate space for training
    allocate_training();

    while ((delta >= min_delta) && (count++ <= MAX_ITR)) 
    {
	update_model(delta, avf);

	//if (avf < last_avf){
	//  mexWarnMsgTxt("WARNING: Log-likelihood decreases in training!");
	//  break;
	//}

	last_avf = avf;

	//#ifdef DEBUG
	//mexPrintf("count = %d\ndelta = %f\navf = %f\n", count, delta, avf);
	//#endif
    }

    // change the model so that the state 1 always has the largest variance
    // state 2 the second, and so on.....
    //reorder_model();

    mexPrintf("done batch-train:\ncount = %d\ndelta = %f\navf = %f\n", 
	      count, delta, avf);
}


//-----------------------------------------------------------------------------
double THMT::batch_test(tree<double> *testTree)
{
    // Assign data pointer

  obs = testTree;
 
  // Compute average log-likelihood
  return test_all();
}


//-----------------------------------------------------------------------------
double THMT::batch_test(char *filename)
{
    // Read data from file
    obs = new tree<double> (filename);

    // Compute average log-likelihood
    double avf = test_all();

    // Delete data
    delete obs;

    return avf;
}


//-----------------------------------------------------------------------------
double THMT::test_all()
{
    register int ob_ind;
    register double f, sumf = 0.0;
    
    
    if (obs->nlev() == 0)
    {
      //cerr << "ERROR in THMT::test_all(): empty data"
      //     << endl;
	return 0.0;
    }
    
   
    if ((obs->nlev() != nLev) || (obs->nch() != nCh))
    {
      //cerr << "ERROR in THMT::test_all(): incompatible data"
      //	     << endl;
	return 0.0;
    }
  
    nObs = obs->nrt();
  
    // Allocate space for training
    allocate_testing();

    for (ob_ind = 0; ob_ind < nObs; ob_ind++) 
    {
	compute_beta(ob_ind);
	f = compute_likelihood();

	//#ifdef DEBUG
	//cout << "ob_ind = " << ob_ind << "\tf = " << f << endl;
	//#endif

	sumf += f;
    }

    return (sumf / double(nObs));
}


//-----------------------------------------------------------------------------
void THMT::rnd_init_model()
{
    const double MAX_MEAN = 100.0;
    const double MAX_STDV = 100.0;
    register int J, MM, m, mm, n;
    double temp;
    vector<double> vprob(M);
    int idum = -time(NULL);    // random seed

    for (J = 0; J < nLev; J++){
      for (mm=0; mm<model_mean[J].size(); mm++)
	for (m = 0; m < M; m++) 
	  {
	    if (zeromean)
	      model_mean[J][mm][m] = 0.0;
	    else
	      model_mean[J][mm][m] = (2.0*ran1(idum) - 1.0) * MAX_MEAN;
	    model_stdv[J][mm][m]  = ran1(idum) * MAX_STDV;

	    ranprobs(vprob, idum);
	    for (n = 0; n < M; n++)
	      model_trans[J][mm][n][m] = vprob[n];
	  }
    }
    
}


//-----------------------------------------------------------------------------
void THMT::generate_data(tree<double>&obs, int n)
{
  tree<double> aTree(1, nCh, nLev);
  register int ob_ind, J;
  int idum = -time(NULL);     // random seed

  // Resize output tree if necessary
  if ((obs.nlev() != nLev) || (obs.nch() != nCh) || (obs.nrt() != n))
    obs = tree<double>(n, nCh, nLev);
  
  for (ob_ind = 0; ob_ind < n; ob_ind++) {
    generate_one(aTree, idum);
    
    for (J = 0; J < nLev; J++)
      copy(aTree[J].begin(), aTree[J].end(),
	   obs[J].begin() + ob_ind * ipow(nCh,J));
  }
}


//-----------------------------------------------------------------------------
void THMT::generate_data(char *filename, int n)
{
    tree<double> genTree(n, nCh, nLev);
    tree<double> aTree(1, nCh, nLev);
    register int ob_ind, J;
    int idum = -time(NULL);     // random seed

    for (ob_ind = 0; ob_ind < n; ob_ind++) {
	generate_one(aTree, idum);

	for (J = 0; J < nLev; J++)
	    copy(aTree[J].begin(), aTree[J].end(),
		 genTree[J].begin() + ob_ind * ipow(nCh,J));
    }

    genTree.dump(filename);
}


//-----------------------------------------------------------------------------
void THMT::generate_one(tree<double> &aTree, int& idum)
{
    tree<int> states(1, nCh, nLev);
    vector<double> vprob(M);
    register int J, i, nNode, m, mm;
    register double mean, stdv;

    // Initial state
    for (m = 0; m < M; m++)
	vprob[m] = model_trans[0][0][m][0];
    
    states[0][0] = ranind(vprob, idum);

    mean = model_mean[0][0][states[0][0]];
    stdv = model_stdv[0][0][states[0][0]];
    aTree[0][0] = mean + stdv * rangas(idum);
  
    // All others
    for (J = 1; J < nLev; J++) 
    {
	for (i = 0, nNode = aTree[J].size(); i < nNode; i++) 
	{
	  mm = subbandtree[J][i];

	  // build vector prob
	  for (m = 0; m < M; m++)
	    vprob[m] = model_trans[J][mm][m][states[J-1][i/nCh]];

	  states[J][i] = ranind(vprob, idum);
	  mean = model_mean[J][mm][states[J][i]];
	  stdv = model_stdv[J][mm][states[J][i]];
	  aTree[J][i] = mean + stdv * rangas(idum);
	}
    }
}


//-----------------------------------------------------------------------------