mi_hpv_2d.cpp

一种新的方法

#define _USE_MATH_DEFINES
#include <math.h> 
#ifndef M_PI 
	#define M_PI 3.1415926535897932384626433832795 
#endif
#include <stdio.h>
#include <limits>
#include "mex.h"
#ifndef M_LOG2_E 
	#define M_LOG2_E 0.693147180559945309417 
#endif
#define log2(x) (log(x) / M_LOG2_E)
#define NMAX    20
#define ITMAX  200
typedef double MAT[NMAX][NMAX];

// quick derivitive of mi_hpv to give four corner update locations for imtransform
void corners(double xTrans, double yTrans, double theta, long nX, long nY, double *q){
	
	double hyp, theta_0, xOff, yOff;
	long x,y;

	xOff = ((double)nX-1.0)/2.0;
	yOff = ((double)nY-1.0)/2.0;

	y=0; x = 0; 
	hyp = sqrt((x-xOff)*(x-xOff)+(yOff-y)*(yOff-y));
	theta_0 = acos((yOff-y)/hyp);
	q[0] = -hyp*cos(theta+theta_0)+yTrans+yOff;
	q[1] = -hyp*sin(theta+theta_0)+xTrans+xOff;
	
	y = nY-1; x=0;
	hyp = sqrt((xOff-x)*(xOff-x)+(y-yOff)*(y-yOff));
	theta_0 = acos((xOff-x)/hyp);
	q[2] = hyp*sin(theta+theta_0)+yTrans+yOff;
	q[3] = -hyp*cos(theta+theta_0)+xTrans+xOff;

	y=nY-1; x=nX-1; 
	hyp = sqrt((x-xOff)*(x-xOff)+(y-yOff)*(y-yOff));
	theta_0 = acos((y-yOff)/hyp);
	q[4] = hyp*cos(theta+theta_0)+yTrans+yOff;
	q[5] = hyp*sin(theta+theta_0)+xTrans+xOff;

	y=0; x = nX-1;
	hyp = sqrt((x-xOff)*(x-xOff)+(yOff-y)*(yOff-y));
	theta_0 = acos((x-xOff)/hyp);
	q[6] = -hyp*sin(theta+theta_0)+yTrans+yOff;
	q[7] = hyp*cos(theta+theta_0)+xTrans+xOff;

}
// Uses HPV interpolation algoritm 
// Lu X, et al., Mutual information-based multimodal image registration using a novel joint histogram estimation, 
// Comput Med Imaging Graph (2008), doi:10.1016/j.compmedimag.2007.12.001
// Author: Matthew Sochor
// Date: 2/20/2008
// contact: [e-mail]
// Other functions (other than main control function) are modified from the below comment author. 
// Only modification is to output negative MI (so maximization problem becomes a minimization problem) and to pass 
// extra needed arrays for the mi_hpv "function"

double mi_hpv(unsigned char *R, unsigned char *F, double xTrans, double yTrans, double theta, long nX, long nY){

	// First compute where each point lies based upon first
	// rotating regular grid by theta (counter-clockwise)
	// then translating grid x and y
	double f_d_kmx, f_d_kmy, theta_0, hyp, qkx, qky, H_A=0.0, H_B=0.0, H_AB=0.0, hsum=0;
	double *P_A = new double[nX],  *P_B = new double[nY], *h = new double[256*256]; 
	long Rx, Ry, hx, hy; 
	double xOff, yOff, xTemp, yTemp, xStart, yStart;
	int j, k;

	xOff = ((double)nX-1.0)/2.0;
	yOff = ((double)nY-1.0)/2.0;

	// initialize histogram for MI measurement
	// 0-255 values correspond to unsigned char (uint8 matlab)
	// signal intensity values
	for(long xx=0; xx<256; xx++){
		for(long yy=0; yy<256; yy++){
			h[xx*256+yy] = 0;
		}
	}
	
	for(double x=0; x<nX; x++){
		for(double y=0; y<nY; y++){
			if (x==xOff){
				if (y==yOff){
					// no rotation just translation, this is the pivot point
					qkx = xTrans+xOff;
					qky = yTrans+yOff;
				}
				else{
					// on y-axis
					qkx = (y-yOff)*sin(theta)+xTrans+xOff;
					qky = (y-yOff)*cos(theta)+yTrans+yOff;
				}
			}
			else if ((x>xOff)&&(y<yOff)){
				// 1st quadrant
				hyp = sqrt((x-xOff)*(x-xOff)+(yOff-y)*(yOff-y));
				theta_0 = acos((x-xOff)/hyp);
				qkx = hyp*cos(theta+theta_0)+xTrans+xOff;
				qky = -hyp*sin(theta+theta_0)+yTrans+yOff;
			}
			else if ((x<xOff)&&(y<yOff)){
				// 2nd quadrant
				hyp = sqrt((x-xOff)*(x-xOff)+(yOff-y)*(yOff-y));
				theta_0 = acos((yOff-y)/hyp);
				qkx = -hyp*sin(theta+theta_0)+xTrans+xOff;
				qky = -hyp*cos(theta+theta_0)+yTrans+yOff;
			}
			else if ((x<xOff)&&(y>yOff)){
				// 3rd quadrant
				hyp = sqrt((xOff-x)*(xOff-x)+(y-yOff)*(y-yOff));
				theta_0 = acos((xOff-x)/hyp);
				qkx = -hyp*cos(theta+theta_0)+xTrans+xOff;
				qky = hyp*sin(theta+theta_0)+yTrans+yOff;
			}
			else if ((x>xOff)&&(y>yOff)){
				// 4th quadrant
				hyp = sqrt((x-xOff)*(x-xOff)+(y-yOff)*(y-yOff));
				theta_0 = acos((y-yOff)/hyp);
				qkx = hyp*sin(theta+theta_0)+xTrans+xOff;
				qky = hyp*cos(theta+theta_0)+yTrans+yOff;
			}
			else{
				// on x-axis
				qkx = (x-xOff)*cos(theta)+xTrans+xOff;
				qky = -(x-xOff)*sin(theta)+yTrans+yOff;
			}

			
			// Find starting point and distance of qkx,qky from convolution point
			// If distances are outside the range 0-nX-1,ny-1 then wrap around
			// Valid if you assume infinitely repeating data
			// build histogram as per HPV algorithm as referenced above
			
			xStart = floor(qkx);
			yStart = floor(qky);
			for (double l1=-1; l1<3; l1++){
				xTemp = xStart+l1;
				if (xTemp < 0){ Rx = (long)(xTemp+nX);}
				else if (xTemp > nX-1){ Rx = (long)(xTemp-nX);}
				else{ Rx=(long)xTemp;}
				f_d_kmx = (1+cos(M_PI*(qkx-xTemp)/2))/4;
				for (double l2=-1; l2<3; l2++){
					yTemp = yStart+l2;
					if (yTemp < 0){ Ry = (long)(yTemp+nY);}
					else if (yTemp > nY-1){ Ry = (long)(yTemp-nY);}
					else{ Ry = (long)yTemp;}
					f_d_kmy = (1+cos(M_PI*(qky-yTemp)/2))/4;
					// update histogram
					hx = F[(long)x*nY+(long)y];
					hy = R[Rx*nY+Ry];
					h[hx*256+hy] += f_d_kmx*f_d_kmy;
				}
			}
		}
	}

	// Determine entropies and joint entropy
	for (j=0;j<256;j++){
		for (k=0;k<256;k++){
			hsum += h[j*256+k];
		}
	}
	for (j=0;j<256;j++){
		for (k=0;k<256;k++){
			if (k==0) {P_A[j] = h[j*256]/hsum; }// init across j
			else{P_A[j]+=h[j*256+k]/hsum;} // add across j
			if (j==0) {P_B[k] = h[k]/hsum;} // init across k
			else {P_B[k]+=h[j*256+k]/hsum;} // add across k
			if (h[j*256+k] != 0){H_AB-=h[j*256+k]/hsum*log2(h[j*256+k]/hsum);} // add j-k
		}
	}
	for (j=0;j<256;j++){
		if (P_A[j] != 0){H_A-=P_A[j]*log2(P_A[j]);}
		if (P_B[j] != 0){H_B-=P_B[j]*log2(P_B[j]);}
	}
	free(P_A);
	free(P_B);
	free(h);
	return(H_A+H_B-H_AB);
}

/*************************************************************
* Minimization of a Function FUNC of N Variables By Powell's *
*    Method Discarding the Direction of Largest Decrease     *
* ---------------------------------------------------------- *
* SAMPLE RUN: Find a minimum of function F(x,y):             *
*             F=Sin(R)/R, where R = Sqrt(x*x+y*y).           *
*                                                            *
* Number of iterations: 2                                    *
*                                                            *
* Minimum Value: -0.217234                                   *
*                                                            *
* at point: 3.177320  3.177320                               *
*                                                            *
* ---------------------------------------------------------- *
* REFERENCE: "Numerical Recipes, The Art of Scientific       *
*             Computing By W.H. Press, B.P. Flannery,        *
*             S.A. Teukolsky and W.T. Vetterling,            *
*             Cambridge University Press, 1986"              *
*             [BIBLI 08].                                    *
*                                                            *
*                         C++ Release By J-P Moreau, Paris.  *
*************************************************************/
/* Code source link: http://pagesperso-orange.fr/jean-pierre.moreau/Cplus/tpowell_cpp.txt */

double P[NMAX];
MAT    XI;

double PCOM[NMAX], XICOM[NMAX];  //PCOM,XICOM,NCOM are common variables
int    ITER,N,NCOM;              //for LINMIN and F1DIM only
double FRET,FTOL;

// user defined function to minimize
double FUNC(unsigned char *R, unsigned char *F, long nX, long nY, double *P) {
	
	double MI;	
	MI = mi_hpv(R, F, P[1], P[2], P[3], nX, nY);
	return(-1.0*MI); // Return negative Mutual Information to change maximization into a minimazation
}

double MAX(double a,double b) {
  if (a>=b) return a; else return b;
}

double MIN(double a,double b) {
  if (a<=b) return a; else return b;
}

double Sign(double a, double b) {
  if (b>=0) return fabs(a);
  else return -fabs(a);
}

double Sqr(double a) {
  return a*a;
}

void   LINMIN(unsigned char *, unsigned char *, long, long, double *, double *, int, double *);
void   MNBRAK(unsigned char *, unsigned char *, long, long, double *, double *, double *, double *, double *, double *);
double BRENT (unsigned char *, unsigned char *, long, long, double *, double *, double *, double, double *);

void POWELL(unsigned char *R, unsigned char *F, long nX, long nY, double *P, MAT XI, int N, double FTOL,int *ITER, double *FRET) {
/*--------------------------------------------------------------------------
  Minimization of a function  FUNC of N variables  (FUNC is not an argument, 
  it is a fixed function name). Input consists of an initial starting point 
  P, that is a vector of length N; an initial matrix XI  whose  logical 
  dimensions are N by N, physical dimensions NMAX by NMAX, and whose columns
  contain the initial set of directions (usually the N unit vectors); and 
  FTOL, the fractional tolerance in the function value such that failure to 
  decrease by more than this amount on one iteration signals doneness. On 
  output, P is set to the best point found, XI is the then-current direction 
  set,  FRET is the returned function value at P, and ITER is the number of 
  iterations taken. The routine LINMIN is used.
---------------------------------------------------------------------------*/
// Label: e1
  double PT[NMAX], PTT[NMAX], XIT[NMAX];
  double DEL,FP,FPTT,T;
  int I,IBIG,J;
  *FRET=FUNC(R,F,nX,nY,P);
  for (J=1; J<=N; J++)
    PT[J]=P[J];            //Save initial point
  *ITER=0;
e1:*ITER=*ITER+1;
  FP=*FRET;
  IBIG=0;
  DEL=0.0;                 //Will be the biggest function decrease
  for (I=1; I<=N; I++) {   //In each iteration, loop over all directions in the set.
    for (J=1; J<=N; J++)   //Copy the direction.
      XIT[J]=XI[J][I];
    FPTT=*FRET;
    LINMIN(R, F, nX, nY,P,XIT,N,FRET);  //Minimize along it
    if (fabs(FPTT-*FRET) > DEL) {
      DEL=fabs(FPTT-*FRET);
      IBIG=I;
    }
  }
  if (2.0*fabs(FP-*FRET) <= FTOL*(fabs(FP)+fabs(*FRET))) return; //Termination criterion
  if (*ITER == ITMAX) {
    printf("\n Powell exceeding maximum iterations.\n\n");
    return;
  }
  for (J=1; J<=N; J++) {
    PTT[J]=2.0*P[J]-PT[J]; //Construct the extrapolated point and the average
    XIT[J]=P[J]-PT[J];     //direction moved. Save the old starting point
    PT[J]=P[J];
  }
  FPTT=FUNC(R,F,nX,nY,PTT);              //Function value at extrapolated point
  if  (FPTT >= FP) goto e1;    //One reason not to use new direction