parallel.cc

主从式并行遗传算法

	{
		
		point_k.x = X[i];
		point_k.y = X[i+1];
		Ak = CalAk(point_k);                         //Ak为考虑点源之间影响的情况下单点的覆盖面积
		//Ak = CalSingleArea(point_k);
		//Fa+=Ak;
		//Fill(point_k);
		Ak_Ac += 2 / ( 1 + ( Ak / Ac ) ) - 1;
		i++;
	}
	F = a_temp * Fa / Global_area + ( 1 - a_temp ) * ( 1 - ( Ak_Ac / Dem ) );
	//F =  1 - ( Ak_Ac / Dem ) ;
	return F;
}                

int CalAk(CPoint pt)             //考虑前面已经有点布置的情况计算单点覆盖
{
	int kk = 0;
	float dlta = pi/36;
  	float dltas=0,alfa,r=0;
   	int len,wide,x,y,maxx=0,minx=260,maxy=0,miny=130;
  	int **polygon;
   	bool flag;
   	//bool flag1;
   
   	polygon = new int*[73];
   	for(int m=0;m<73;m++)
    		 polygon[m] = new int[2];

   	polygon[0][0] = 73;
   	polygon[0][1] = 0;
   	
   	//下面三项好像可以预处理得到每个点源的    m_maxdis[72]
                                                         // 		m_shademax[72];
       //主要得到pt所对应的 m_maxdis[72]
                       // 		m_shademax[72];
   	Shade(pt);
   	MaxHDis(pt);//已全部换算到256*128分辨率下
   	AdjustDis(pt);
  /* /////////////////////////
	FILE *fp;

	fp=fopen("f:/lastdata/mresult.txt","w");
    fprintf(fp,"%d%c%d\n",pt.x,' ',pt.y);
	  
	   for(int h=0;h<72;h++)
	   {
	      fprintf(fp,"%f\n",m_maxdis[h]);
	   }

    fclose(fp);
	//////////////////////////////////*/
	//以下得到点源pt的polygon矩阵，同时得到多边形最左和最右的点的坐标(minx,miny)  (maxx,maxy)
   	for(int i=1;i<73;i++)
   	{
      		alfa = (float)(i-1)*dlta;
      		polygon[i][0] = (int)(pt.x + m_maxdis[i-1]*cos(alfa));
	  	if(maxx<polygon[i][0])
			  maxx = polygon[i][0];
	 	if(minx>polygon[i][0])
		  	minx = polygon[i][0];
	  	polygon[i][1] = (int)(pt.y + m_maxdis[i-1]*sin(alfa)); 
	  	if(maxy<polygon[i][1])
		  	maxy = polygon[i][1];
	  	if(miny>polygon[i][1])
		  	miny = polygon[i][1];
   	}

   	len = maxx - minx;
   	wide = maxy - miny;
   	for(int l=0;l<len;l++)
	  	 for(int k=0;k<wide;k++)
	  	 {
	   	 	x = minx + l;
		  	y = miny + k;
		  	CPoint mpt;
			mpt.x = x;
			mpt.y = y;
		  	//if(Region_Flag[x][y])           这里也不判断点是否在陆地上
		  	//flag1 = IsPointInBoundary(mpt);
		  	//if(flag1==true)
		  	//{
		     	flag = IsPointInPolygon(x,y,polygon,73);
		     	if(flag)
			 	{
		       		if(Flag[x][y]!=1)              //在没有确定在此放置点时,先不把Flag[x][y]设置为1
					{
						kk++;
						Flag[x][y] = 1;
					}
				 // Flag[x][y] = 1;
			  
			 //}
		  }
	   }
   
   	for(int n=0;n<73;n++)
	   	delete[] polygon[n];
   	delete[] polygon;
  	return kk;
}


void Shade(CPoint pt)               //是否可以先计算出所有点的m_shademax[72]
{
	float delta = pi/36;
	float dis,ran,mean,ratio;
	float alfa=0,maxratio=0,deltaz,tmpy,tmpx;
	int x,y;
	CPoint point;
	int R = 25;

  	for(int i=0;i<72;i++)
  	{
      		maxratio=0;
	  	alfa = (float)i*delta;
	   
      		for(int dd=0;dd<R;dd++)
	  	{
		   	tmpx = pt.x + dd*cos(alfa);
		   	tmpy = pt.y + dd*sin(alfa);//将极坐标化为直角坐标
		   	x = (int)tmpx;
		   	y = (int)tmpy;
		   
		   	if((x>0)&&(x<255)&&(y>0)&&(y<127))
		   	{
               		double mean1 = (tmpy-y)*(high[x][y+1]-high[x][y])+high[x][y];
		       		double mean2 = (tmpy-y)*(high[x+1][y+1]-high[x+1][y])+high[x+1][y];
		       		mean = (tmpx-x)*(mean2-mean1)+mean1;
               		deltaz = mean - high[pt.x][pt.y];
               		point.x = x;
			   	point.y = y;
			   	dis = DistanceOfTwoPoints(pt.x,pt.y,point.x,point.y);
               		ratio = deltaz/dis;
			   	if(maxratio<ratio)
				   	maxratio = ratio; 
		   }
	 	 }
	  	ran = atan(maxratio);
	  	m_shademax[i] = (ran*180)/pi;
  	}
}

void MaxHDis(CPoint pt)     //是否也可以将 m_maxdis[i]计算保存到文件，每次读次文件即可                  
{
	float dlta = pi/36;
	float alfa;
	float height,aimheight=1020; //4km； 256*128分辨率下1象素对应171.7952m 
  //水平面上1个象素对应于竖直方向4.38单位,现取25象素为10km
  /***************************************************************************************
    255~1km,aimheight是目标高度(高度单位)，这里取绝对高度，因为在高度方向上和平面上的分辨率
	不同,高度方向上1单位对应3.92m，平面上1像素对应于42.9488m,因此平面上1单位对应于高度方向
	上11个单位。对于10km而言，在平面上对应于233象素，高度方向上对应于255单位。现在令100象素
    对应10km，则高度方向上109.5单位对应1km
  *****************************************************************************************/
	double maxR,temp;
//	CPoint crtpoint;
	int hR=2550;  //最大半径10km对应的高度单位数
 
  //int hR=1530;  //最大半径10km对应的高度单位数
	height = aimheight - high[pt.x][pt.y];  //换成现在比例尺下的高度单位数
	maxR = sqrt(hR*hR-height*height);
  
	for(int i=0;i<72;i++)
  	{
	  	alfa = (float)i*dlta;
	  	if(m_shademax[i]>0)
	 	{
	    		float a = m_shademax[i]*pi/180;
			temp = height/tan(a);
			if(temp>maxR)
				m_maxdis[i] = maxR/102;
			else
				m_maxdis[i] = temp/102;
      
	  	}
	  	else
		 	 m_maxdis[i] = maxR/102;  //将高度单位换算成象素单位
  	}   
}

void AdjustDis(CPoint pt)
{
	float dlta = pi/36;
	float alfa,r=0;
	CPoint crtpoint,tpt;
	int flag,flag1;
	//int x;
	//int y;
  
	for(int i=0;i<72;i++)
   	{
		alfa = (float)i*dlta;
      	crtpoint.x = (int)(pt.x + m_maxdis[i]*cos(alfa));
	  	crtpoint.y = (int)(pt.y + m_maxdis[i]*sin(alfa));
      		//flag = IsPointInBoundary(crtpoint);
      	flag = Region_Flag[crtpoint.x][crtpoint.y];
	  	if(!flag)
	 	{
	     	for(float r=m_maxdis[i];r>0;r-=0.1)
		 	{
				tpt.x = (int)(pt.x+r*cos(alfa));
				tpt.y = (int)(pt.y+r*sin(alfa));
				//flag1 = IsPointInBoundary(tpt);
				flag1 = Region_Flag[tpt.x][tpt.y];
				if(flag1)
				{
			  		m_maxdis[i] = r;
			   		break;
				}
	    			else
			   		m_maxdis[i] = 0; //有问题，死循环??
			}	
	  	}
   	}
}

int MaxFit_index(float *fit, int Size)
{
	int size = Size;
	int max = 0;
       for( int i = 1; i < size; i++ )
	{
		if(fit[i] > fit[max])
			max = i;
	} 
	return max;
}

int MinFit_index(float *fit, int Size)
{
	int size = Size;
	int min = 0;
       for( int i = 1; i < size; i++ )
	{
		if(fit[i] < fit[min])
			min = i;
	} 
	return min;
}

bool IsPointInPolygon(float x, float y, int **polygon, int PointNumber)
{
	//判断点(x,y)是否在多边形polygon的内部.
	//对多边形的要求，多边形是一些离散的点序列组成，可以不闭合
	//这里对程序进行了改动是因为我的多边形序列存储时是从下标为1时开始的，polygon[0][0]中存放的是多边形序列的长度
	//即就是PointNumber-1的值，polygon[0][1]中没有放有用的值
	int i;
	float *Alpha;
	float a1,b1;
	float d1,d2,d3;

	Alpha = (float*)malloc(sizeof(float)*PointNumber);

	for(i = 1; i < PointNumber -1;i++)
	{
		d1 = DistanceOfTwoPoints(x,y,polygon[i][0],polygon[i][1]);
		d2 = DistanceOfTwoPoints(x,y,polygon[i+1][0],polygon[i+1][1]);
		d3 = DistanceOfTwoPoints(polygon[i][0],polygon[i][1],polygon[i+1][0],polygon[i+1][1]);
		Alpha[i] = (double)acos((double)(d1*d1+d2*d2-d3*d3)/(2.0*d1*d2));

		a1 = CrossProductOfTwoLine(x,y,polygon[i][0],polygon[i][1],polygon[i+1][0],polygon[i+1][1]);
		if(a1 > 0.0)	
			Alpha[i] = Alpha[i]*(-1.0f);
	}
	d1 = DistanceOfTwoPoints(x,y,polygon[PointNumber-1][0],polygon[PointNumber-1][1]);
	d2 = DistanceOfTwoPoints(x,y,polygon[1][0],polygon[1][1]);
	d3 = DistanceOfTwoPoints(polygon[PointNumber-1][0],polygon[PointNumber-1][1],polygon[1][0],polygon[1][1]);
	Alpha[PointNumber-1] = (double)acos((double)(d1*d1+d2*d2-d3*d3)/(2.0*d1*d2));

	a1 = CrossProductOfTwoLine(x,y,polygon[PointNumber-1][0],polygon[PointNumber-1][1],polygon[1][0],polygon[1][1]);
	if(a1 > 0.0)	
		Alpha[PointNumber-1] = Alpha[PointNumber-1]*(-1.0f);

	b1 = 0.0f;
	for(i = 1; i < PointNumber; i++)
	{
		//	TRACE("%f\n",Alpha[i]);
		b1 += Alpha[i];
	}
	free(Alpha);
	if(fabs((float)b1) <= 0.001) 
		return false;
	else
		return true;

}

float CrossProductOfTwoLine(float x1, float y1, float x2, float y2, float x3, float y3)
{
	float CrossProduct;
	CrossProduct = (x2-x1)*(y3-y1)-(y2-y1)*(x3-x1);
	return CrossProduct;
}

float DistanceOfTwoPoints(float x1, float y1, float x2, float y2)
{
	float dis;
	dis = float(sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1)));
	return dis;
}

 void InitialArray()
 {
    individual_vector = new int* [size];      //int*类型的数组
	//newindividual_vector = new int* [size*2]; 
    fit = new float[size];          //存放适应度值的数组，一个适应度值有相应的individual_vector  
    binary = new int*[size];
    for(int i=0;i<size;i++)
    {
      individual_vector[i] = new int[2*num];
      binary[i] = new int[15*num];                  //binary是否可以用char
    }
	new_fit = new float[size*2];                 //存放40个个体的适应度  
/*	for( i = 0; i < 2*size; i++)
	{
		newindividual_vector[i] = new int[2*num];
	} 
	*/
 }

int FitFunction(int *X, int Dem)
{
   int Fa=0,Ak;
   CPoint p;
   int m=0;
   for(int h=0;h<256;h++)
	  for(int l=0;l<128;l++)
	  {
	    Flag[h][l] = 0;
	  }
   /*while(m<2*Dem)
   {
	  p.x = X[m];
      m++;
	  p.y = X[m];
	  CalArea(p);
	  m++;
   }

   for(int k=0;k<256;k++)
	   for(int j=0;j<128;j++)
	   {
	      S+=Flag[k][j];
	   }
	   */
	for( int i = 0; i < Dem*2; i++ )
	{
		
		p.x = X[i];
		p.y = X[i+1];
		Ak = CalAk(p);               //CalAk(p)中已经完成了Fill(p)的功能
		Fa+=Ak;
		//Fill(p);
		i++;
	}

   return Fa;
}

bool WithDeadArea(int cs)
{

   if(Global_area-cs<=E_CUT)
	   return true;// no deadarea
   else
	   return false; // have deadarea

}

void Initial()
{
	//int size = Size;
	//int num = Num;
	CPoint p;
	
	//每次初始化之前把标志为置为0 ！！！！！！！！！！！！！！！！！！！
	for( int i = 0; i < size; i++ )
	{
		for(int h=0;h<256;h++)
			for(int l=0;l<128;l++)