算法rx
① SGM半全局匹配算法求助
#include "stdafx.h"
#include <cstdio>
#include <cstring>
#include <iostream>
#include<cv.h>
#include<highgui.h>
#include <cmath>
using namespace std;
const int Width = 1024;
const int Height = 1024;
int Lvalue[Width][Width];
uchar C(int x, int y, int d, IplImage * matchImage, IplImage * baseImage)
{
uchar * pMatchPixel = NULL;
uchar * pBasePixel = NULL;
uchar matchPixel = 0;
uchar basePixel =0;
uchar matchMax = 0;
uchar matchMin = 0;
uchar tempMatch1 = 0;
uchar tempMatch2 = 0;
uchar difPixel = 0;
if (x+d <= matchImage->width && x+d >= 0)
{
pMatchPixel = (uchar *)matchImage->imageData + y*matchImage->widthStep + (x+d);
matchPixel = *pMatchPixel;
pBasePixel= (uchar *)baseImage->imageData + y*baseImage->widthStep + x;
basePixel = *pBasePixel;
//匹配影像上的像素点与其左、右像素点的平均值,线性插值的方法
tempMatch1 = (*pMatchPixel +(*(pMatchPixel -1)))/2;
tempMatch2 = (*pMatchPixel +(*(pMatchPixel +1)))/2;
matchMax = max(max(tempMatch1,tempMatch2),matchPixel);
matchMin = min(min(tempMatch1,tempMatch2),matchPixel);
//赋予视差d时的匹配代价C
//BT法
difPixel = max(max(basePixel - matchMax, matchMin - basePixel),0);
//AD法
//difPixel = abs(basePixel - matchPixel);
return difPixel;
}
else
return 255;
}
int main()
{
IplImage * leftImage = cvLoadImage("headL.png",0);
IplImage * rightImage = cvLoadImage("headR.png",0);
int imageWidth = leftImage->width;
int imageHeight =leftImage->height;
int minLvalue = 1000;
int minL1 = 1000;
int minL2 = 1000;
int P1 = 2;
int P2 = 5;
int disparity= 0;
int minDis = 0;
int maxDis = 21;
int scale = 12;
unsigned char * pPixel = NULL;
#pragma region 实现横向DP（以右影像为准,0度方向）
IplImage * MyDPImage_0 = cvCreateImage(cvGetSize(leftImage),leftImage->depth,1);
cvZero(MyDPImage_0);
int t3 = clock();
for (int i = 0; i < imageHeight;i++)
{
for (int j = 0; j<imageWidth ;j++)
{
disparity = 0;
minL1 = 1000;
minL2 = 1000;
for (int d = minDis; d <= maxDis; d++)
{
//初始化为代价函数的值
Lvalue[j][d] = C(j, i, d, leftImage, rightImage);
if (j > 0)
{
minL1 = min(minL1, Lvalue[j-1][d]);
}
}
for (int d = minDis; d <= maxDis; d++)
{
if (j > 0)
{
minL2 = min(minL2, Lvalue[j-1][d]);
minL2 = min(minL2, (Lvalue[j-1][d+1] + P1));
minL2 = min(minL2, (Lvalue[j-1][d-1] + P1));
minL2 = min(minL2, (minL1 + P2));
Lvalue[j][d] = Lvalue[j][d] + (minL2 - minL1);
}
}
int max=Lvalue[j][0];
for(int d=minDis;d<=maxDis;++d)
{
if(Lvalue[j][d]<max)
{
disparity=d;
max=Lvalue[j][d];
}
}
disparity=disparity*scale;
//生成视差图
pPixel = (uchar *)MyDPImage_0->imageData + i*MyDPImage_0->widthStep + j;
*pPixel =disparity;
}
}
int t4 = clock();
cout<<"横向DP共用时: "<<t4-t3<<"ms"<<endl;
cvNamedWindow("MyDPImage_0", 1);
cvShowImage("MyDPImage_0", MyDPImage_0);
cvSaveImage("MyDPImage_0.jpg", MyDPImage_0);
#pragma endregion
#pragma region 实现横向DP（以左影像为准，0度方向）
IplImage * MyDPImage_0_L = cvCreateImage(cvGetSize(leftImage),leftImage->depth,1);
cvZero(MyDPImage_0_L);
for (int i = 0; i < imageHeight;i++)
{
for (int j = 0; j<imageWidth;j++)
{
disparity = 0;
minL1 = 1000;
minL2 = 1000;
for (int d = minDis; d <= maxDis; d++)
{
//初始化为代价函数的值
Lvalue[j][d] = C(j, i, -d, rightImage, leftImage);
if (j > 0)
{
minL1 = min(minL1, Lvalue[j-1][d]);
}
}
for (int d = minDis; d <= maxDis; d++)
{
if (j > 0)
{
minL2 = min(minL2, Lvalue[j-1][d]);
minL2 = min(minL2, (Lvalue[j-1][d+1] + P1));
minL2 = min(minL2, (Lvalue[j-1][d-1] + P1));
minL2 = min(minL2, (minL1 + P2));
Lvalue[j][d] = Lvalue[j][d] + minL2 - minL1;
}
}
int max=Lvalue[j][0];
for(int d=0;d<=maxDis;++d)
{
if(Lvalue[j][d]<max)
{
disparity=d;
max=Lvalue[j][d];
}
}
disparity = disparity * scale;
//生成视差图
pPixel = (uchar *)MyDPImage_0_L->imageData + i*MyDPImage_0_L->widthStep + j;
*pPixel = disparity;
}
}
cvNamedWindow("MyDPImage_0_L", 1);
cvShowImage("MyDPImage_0_L", MyDPImage_0_L);
cvSaveImage("MyDPImage_0_L.jpg", MyDPImage_0_L);
#pragma endregion
#pragma region 一致性检查
uchar * pLeftDepthPixel = NULL;
uchar * pRightDepthPixel = NULL;
uchar leftDepthPixel = 0;
uchar rightDepthPixel =0;
uchar difDepthPixel = 0;
IplImage * CheckImage_0 = cvCloneImage(MyDPImage_0);
cvZero(CheckImage_0);
for (int i = 0; i < imageHeight; i++)
{
for(int j = 0; j < imageWidth; j++)
{
pRightDepthPixel = (uchar *)MyDPImage_0->imageData + i*MyDPImage_0->widthStep + j;
rightDepthPixel = *pRightDepthPixel;
if(j + rightDepthPixel/scale < imageWidth)
{
pLeftDepthPixel = (uchar *)MyDPImage_0_L->imageData + i*MyDPImage_0_L->widthStep + j + rightDepthPixel/scale;
leftDepthPixel = *pLeftDepthPixel;
}
else
leftDepthPixel = 0;
difDepthPixel = abs((leftDepthPixel - rightDepthPixel)/scale);
pPixel = (uchar *)CheckImage_0->imageData + i * CheckImage_0->widthStep +j;
if (difDepthPixel <= 1)
{
*pPixel = rightDepthPixel;
}
else
{
//否则，当前像素的视差值为零
*pPixel = 0;
}
}
}
int lp,rp;
int lx,rx;
for(int i=0;i<imageHeight;i++)
{
for(int j=0;j<imageWidth;j++)
{
int tem = *((uchar *)CheckImage_0->imageData+i*CheckImage_0->widthStep + j);
if(tem==0)
{
lp = rp = 0;
lx = j;
rx = j;
if(lx-1<0)
lp= *((uchar *)CheckImage_0->imageData+i*CheckImage_0->widthStep + lx);
while((lp==0)&&( lx-1 >= 0 ))
lp = *((uchar *)CheckImage_0->imageData+i*CheckImage_0->widthStep + (--lx));
if(rx+1>=imageWidth)
rp = *((uchar *)CheckImage_0->imageData+i*CheckImage_0->widthStep + rx);
while((rp==0)&&(rx+1<imageWidth))
rp = *((uchar *)CheckImage_0->imageData+i*CheckImage_0->widthStep +(++rx));
if (lp > rp)
{
*((uchar *)CheckImage_0->imageData+i*CheckImage_0->widthStep + j) = rp;
}
else
{
*((uchar *)CheckImage_0->imageData+i*CheckImage_0->widthStep + j) = lp;
}
}
}
}
cvSmooth(CheckImage_0,CheckImage_0,CV_MEDIAN,3,0,0);
cvNamedWindow("CheckImage_0", 1);
cvShowImage("CheckImage_0", CheckImage_0);
cvSaveImage("CheckImage_0.jpg", CheckImage_0);
#pragma endregion
cout << "完成！"<<endl;
cvWaitKey(0);
return 0;
}
② 3个led接在10电源上需多大限流电阻
单个普通白光：电流20mA电压3.2V
3个串联是9.6V,20组并联式400mA
电源可以满足输出要求
每组限流电阻Rx算法：3.2V×3+Rx×20mA=12V
可得：Rx=120Ω
③ MATLAB中RX算法的程序~急用
function [ROUTES,PL,Tau]=ACASP(G,Tau,K,M,S,E,Alpha,Beta,Rho,Q)
%% ---------------------------------------------------------------
% ACASP.m
% 蚁群算法动态寻路算法
% ChengAihua,PLA Information Engineering University,ZhengZhou,China
% Email:[email protected]
% All rights reserved
%% ---------------------------------------------------------------
% 输入参数列表
% G 地形图为01矩阵，如果为1表示障碍物
% Tau 初始信息素矩阵（认为前面的觅食活动中有残留的信息素）
% K 迭代次数（指蚂蚁出动多少波）
% M 蚂蚁个数（每一波蚂蚁有多少个）
% S 起始点（最短路径的起始点）
% E 终止点（最短路径的目的点）
% Alpha 表征信息素重要程度的参数
% Beta 表征启发式因子重要程度的参数
% Rho 信息素蒸发系数
% Q 信息素增加强度系数
%
% 输出参数列表
% ROUTES 每一代的每一只蚂蚁的爬行路线
% PL 每一代的每一只蚂蚁的爬行路线长度
% Tau 输出动态修正过的信息素
%% --------------------变量初始化----------------------------------
%load
D=G2D(G);
N=size(D,1);%N表示问题的规模（象素个数）
MM=size(G,1);
a=1;%小方格象素的边长
Ex=a*(mod(E,MM)-0.5);%终止点横坐标
if Ex==-0.5
Ex=MM-0.5;
end
Ey=a*(MM+0.5-ceil(E/MM));%终止点纵坐标
Eta=zeros(1,N);%启发式信息，取为至目标点的直线距离的倒数
%下面构造启发式信息矩阵
for i=1:N
if ix==-0.5
ix=MM-0.5;
end
iy=a*(MM+0.5-ceil(i/MM));
if i~=E
Eta(1,i)=1/((ix-Ex)^2+(iy-Ey)^2)^0.5;
else
Eta(1,i)=100;
end
end
ROUTES=cell(K,M);%用细胞结构存储每一代的每一只蚂蚁的爬行路线
PL=zeros(K,M);%用矩阵存储每一代的每一只蚂蚁的爬行路线长度
%% -----------启动K轮蚂蚁觅食活动，每轮派出M只蚂蚁--------------------
for k=1:K
disp(k);
for m=1:M
%% 第一步：状态初始化
W=S;%当前节点初始化为起始点
Path=S;%爬行路线初始化
PLkm=0;%爬行路线长度初始化
TABUkm=ones(1,N);%禁忌表初始化
TABUkm(S)=0;%已经在初始点了，因此要排除
DD=D;%邻接矩阵初始化
%% 第二步：下一步可以前往的节点
DW=DD(W,:);
DW1=find(DW
for j=1:length(DW1)
if TABUkm(DW1(j))==0
DW(j)=inf;
end
end
LJD=find(DW
Len_LJD=length(LJD);%可选节点的个数
%% 觅食停止条件：蚂蚁未遇到食物或者陷入死胡同
while W~=E&&Len_LJD>=1
%% 第三步：转轮赌法选择下一步怎么走
PP=zeros(1,Len_LJD);
for i=1:Len_LJD
PP(i)=(Tau(W,LJD(i))^Alpha)*(Eta(LJD(i))^Beta);
end
PP=PP/(sum(PP));%建立概率分布
Pcum=cumsum(PP);
Select=find(Pcum>=rand);
%% 第四步：状态更新和记录
Path=[Path,to_visit];%路径增加
PLkm=PLkm+DD(W,to_visit);%路径长度增加
W=to_visit;%蚂蚁移到下一个节点
for kk=1:N
if TABUkm(kk)==0
DD(W,kk)=inf;
DD(kk,W)=inf;
end
end
TABUkm(W)=0;%已访问过的节点从禁忌表中删除
for j=1:length(DW1)
if TABUkm(DW1(j))==0
DW(j)=inf;
end
end
LJD=find(DW
Len_LJD=length(LJD);%可选节点的个数
end
%% 第五步：记下每一代每一只蚂蚁的觅食路线和路线长度
ROUTES{k,m}=Path;
if Path(end)==E
PL(k,m)=PLkm;
else
PL(k,m)=inf;
end
end
%% 第六步：更新信息素
Delta_Tau=zeros(N,N);%更新量初始化
for m=1:M
if PL(k,m) ROUT=ROUTES{k,m};
TS=length(ROUT)-1;%跳数
PL_km=PL(k,m);
for s=1:TS
x=ROUT(s);
Delta_Tau(x,y)=Delta_Tau(x,y)+Q/PL_km;
Delta_Tau(y,x)=Delta_Tau(y,x)+Q/PL_km;
end
end
end
Tau=(1-Rho).*Tau+Delta_Tau;%信息素挥发一部分，新增加一部分
end
%% ---------------------------绘图--------------------------------
plotif=1;%是否绘图的控制参数
if plotif==1
%绘收敛曲线
meanPL=zeros(1,K);
minPL=zeros(1,K);
for i=1:K
PLK=PL(i,:);
Nonzero=find(PLK
PLKPLK=PLK(Nonzero);
meanPL(i)=mean(PLKPLK);
minPL(i)=min(PLKPLK);
end
figure(1)
plot(minPL);
hold on
plot(meanPL);
grid on
title('收敛曲线（平均路径长度和最小路径长度）');
xlabel('迭代次数');
ylabel('路径长度');
%绘爬行图
figure(2)
axis([0,MM,0,MM])
for i=1:MM
for j=1:MM
if G(i,j)==1
x1=j-1;y1=MM-i;
x2=j;y2=MM-i;
x3=j;y3=MM-i+1;
x4=j-1;y4=MM-i+1;
fill([x1,x2,x3,x4],[y1,y2,y3,y4],[0.2,0.2,0.2]);
hold on
else
x1=j-1;y1=MM-i;
x2=j;y2=MM-i;
x3=j;y3=MM-i+1;
x4=j-1;y4=MM-i+1;
fill([x1,x2,x3,x4],[y1,y2,y3,y4],[1,1,1]);
hold on
end
end
end
hold on
ROUT=ROUTES{K,M};
LENROUT=length(ROUT);
Rx=ROUT;
Ry=ROUT;
for ii=1:LENROUT
Rx(ii)=a*(mod(ROUT(ii),MM)-0.5);
if Rx(ii)==-0.5
Rx(ii)=MM-0.5;
end
Ry(ii)=a*(MM+0.5-ceil(ROUT(ii)/MM));
end
plot(Rx,Ry)
end
plotif2=1;%绘各代蚂蚁爬行图
if plotif2==1
figure(3)
axis([0,MM,0,MM])
for i=1:MM
for j=1:MM
if G(i,j)==1
x1=j-1;y1=MM-i;
x2=j;y2=MM-i;
x3=j;y3=MM-i+1;
x4=j-1;y4=MM-i+1;
fill([x1,x2,x3,x4],[y1,y2,y3,y4],[0.2,0.2,0.2]);
hold on
else
x1=j-1;y1=MM-i;
x2=j;y2=MM-i;
x3=j;y3=MM-i+1;
x4=j-1;y4=MM-i+1;
fill([x1,x2,x3,x4],[y1,y2,y3,y4],[1,1,1]);
hold on
end
end
end
for k=1:K
PLK=PL(k,:);
minPLK=min(PLK);
pos=find(PLK==minPLK);
m=pos(1);
ROUT=ROUTES{k,m};
LENROUT=length(ROUT);
Rx=ROUT;
Ry=ROUT;
for ii=1:LENROUT
Rx(ii)=a*(mod(ROUT(ii),MM)-0.5);
if Rx(ii)==-0.5
Rx(ii)=MM-0.5;
end
Ry(ii)=a*(MM+0.5-ceil(ROUT(ii)/MM));
end
plot(Rx,Ry)
hold on
end
end
将上述算法应用于机器人路径规划，优化效果如下图所示
④ SGM半全局匹配算法求助
SGM半全局匹配算法用c++实现的代码如下，仅供参考（vs中编译通过）。
#include"stdafx.h"
#include<cstdio>
#include<cstring>
#include<iostream>
#include<cv.h>
#include<highgui.h>
#include<cmath>
usingnamespacestd;
constintWidth=1024;
constintHeight=1024;
intLvalue[Width][Width];
ucharC(intx,inty,intd,IplImage*matchImage,IplImage*baseImage)
{
uchar*pMatchPixel=NULL;
uchar*pBasePixel=NULL;
ucharmatchPixel=0;
ucharbasePixel=0;
ucharmatchMax=0;
ucharmatchMin=0;
uchartempMatch1=0;
uchartempMatch2=0;
uchardifPixel=0;
if(x+d<=matchImage->width&&x+d>=0)
{
pMatchPixel=(uchar*)matchImage->imageData+y*matchImage->widthStep+(x+d);
matchPixel=*pMatchPixel;
pBasePixel=(uchar*)baseImage->imageData+y*baseImage->widthStep+x;
basePixel=*pBasePixel;
//匹配影像上的像素点与其左、右像素点的平均值,线性插值的方法
tempMatch1=(*pMatchPixel+(*(pMatchPixel-1)))/2;
tempMatch2=(*pMatchPixel+(*(pMatchPixel+1)))/2;
matchMax=max(max(tempMatch1,tempMatch2),matchPixel);
matchMin=min(min(tempMatch1,tempMatch2),matchPixel);
//赋予视差d时的匹配代价C
//BT法
difPixel=max(max(basePixel-matchMax,matchMin-basePixel),0);
//AD法
//difPixel=abs(basePixel-matchPixel);
returndifPixel;
}
else
return255;
}
intmain()
{
IplImage*leftImage=cvLoadImage("headL.png",0);
IplImage*rightImage=cvLoadImage("headR.png",0);
intimageWidth=leftImage->width;
intimageHeight=leftImage->height;
intminLvalue=1000;
intminL1=1000;
intminL2=1000;
intP1=2;
intP2=5;
intdisparity=0;
intminDis=0;
intmaxDis=21;
intscale=12;
unsignedchar*pPixel=NULL;
#pragmaregion实现横向DP（以右影像为准,0度方向）
IplImage*MyDPImage_0=cvCreateImage(cvGetSize(leftImage),leftImage->depth,1);
cvZero(MyDPImage_0);
intt3=clock();
for(inti=0;i<imageHeight;i++)
{
for(intj=0;j<imageWidth;j++)
{
disparity=0;
minL1=1000;
minL2=1000;
for(intd=minDis;d<=maxDis;d++)
{
//初始化为代价函数的值
Lvalue[j][d]=C(j,i,d,leftImage,rightImage);
if(j>0)
{
minL1=min(minL1,Lvalue[j-1][d]);
}
}
for(intd=minDis;d<=maxDis;d++)
{
if(j>0)
{
minL2=min(minL2,Lvalue[j-1][d]);
minL2=min(minL2,(Lvalue[j-1][d+1]+P1));
minL2=min(minL2,(Lvalue[j-1][d-1]+P1));
minL2=min(minL2,(minL1+P2));
Lvalue[j][d]=Lvalue[j][d]+(minL2-minL1);
}
}
intmax=Lvalue[j][0];
for(intd=minDis;d<=maxDis;++d)
{
if(Lvalue[j][d]<max)
{
disparity=d;
max=Lvalue[j][d];
}
}
disparity=disparity*scale;
//生成视差图
pPixel=(uchar*)MyDPImage_0->imageData+i*MyDPImage_0->widthStep+j;
*pPixel=disparity;
}
}
intt4=clock();
cout<<"横向DP共用时:"<<t4-t3<<"ms"<<endl;
cvNamedWindow("MyDPImage_0",1);
cvShowImage("MyDPImage_0",MyDPImage_0);
cvSaveImage("MyDPImage_0.jpg",MyDPImage_0);
#pragmaendregion
#pragmaregion实现横向DP（以左影像为准，0度方向）
IplImage*MyDPImage_0_L=cvCreateImage(cvGetSize(leftImage),leftImage->depth,1);
cvZero(MyDPImage_0_L);
for(inti=0;i<imageHeight;i++)
{
for(intj=0;j<imageWidth;j++)
{
disparity=0;
minL1=1000;
minL2=1000;
for(intd=minDis;d<=maxDis;d++)
{
//初始化为代价函数的值
Lvalue[j][d]=C(j,i,-d,rightImage,leftImage);
if(j>0)
{
minL1=min(minL1,Lvalue[j-1][d]);
}
}
for(intd=minDis;d<=maxDis;d++)
{
if(j>0)
{
minL2=min(minL2,Lvalue[j-1][d]);
minL2=min(minL2,(Lvalue[j-1][d+1]+P1));
minL2=min(minL2,(Lvalue[j-1][d-1]+P1));
minL2=min(minL2,(minL1+P2));
Lvalue[j][d]=Lvalue[j][d]+minL2-minL1;
}
}
intmax=Lvalue[j][0];
for(intd=0;d<=maxDis;++d)
{
if(Lvalue[j][d]<max)
{
disparity=d;
max=Lvalue[j][d];
}
}
disparity=disparity*scale;
//生成视差图
pPixel=(uchar*)MyDPImage_0_L->imageData+i*MyDPImage_0_L->widthStep+j;
*pPixel=disparity;
}
}
cvNamedWindow("MyDPImage_0_L",1);
cvShowImage("MyDPImage_0_L",MyDPImage_0_L);
cvSaveImage("MyDPImage_0_L.jpg",MyDPImage_0_L);
#pragmaendregion
#pragmaregion一致性检查
uchar*pLeftDepthPixel=NULL;
uchar*pRightDepthPixel=NULL;
ucharleftDepthPixel=0;
ucharrightDepthPixel=0;
uchardifDepthPixel=0;
IplImage*CheckImage_0=cvCloneImage(MyDPImage_0);
cvZero(CheckImage_0);
for(inti=0;i<imageHeight;i++)
{
for(intj=0;j<imageWidth;j++)
{
pRightDepthPixel=(uchar*)MyDPImage_0->imageData+i*MyDPImage_0->widthStep+j;
rightDepthPixel=*pRightDepthPixel;
if(j+rightDepthPixel/scale<imageWidth)
{
pLeftDepthPixel=(uchar*)MyDPImage_0_L->imageData+i*MyDPImage_0_L->widthStep+j+rightDepthPixel/scale;
leftDepthPixel=*pLeftDepthPixel;
}
else
leftDepthPixel=0;
difDepthPixel=abs((leftDepthPixel-rightDepthPixel)/scale);
pPixel=(uchar*)CheckImage_0->imageData+i*CheckImage_0->widthStep+j;
if(difDepthPixel<=1)
{
*pPixel=rightDepthPixel;
}
else
{
//否则，当前像素的视差值为零
*pPixel=0;
}
}
}
intlp,rp;
intlx,rx;
for(inti=0;i<imageHeight;i++)
{
for(intj=0;j<imageWidth;j++)
{
inttem=*((uchar*)CheckImage_0->imageData+i*CheckImage_0->widthStep+j);
if(tem==0)
{
lp=rp=0;
lx=j;
rx=j;
if(lx-1<0)
lp=*((uchar*)CheckImage_0->imageData+i*CheckImage_0->widthStep+lx);
while((lp==0)&&(lx-1>=0))
lp=*((uchar*)CheckImage_0->imageData+i*CheckImage_0->widthStep+(--lx));
if(rx+1>=imageWidth)
rp=*((uchar*)CheckImage_0->imageData+i*CheckImage_0->widthStep+rx);
while((rp==0)&&(rx+1<imageWidth))
rp=*((uchar*)CheckImage_0->imageData+i*CheckImage_0->widthStep+(++rx));
if(lp>rp)
{
*((uchar*)CheckImage_0->imageData+i*CheckImage_0->widthStep+j)=rp;
}
else
{
*((uchar*)CheckImage_0->imageData+i*CheckImage_0->widthStep+j)=lp;
}
}
}
}
cvSmooth(CheckImage_0,CheckImage_0,CV_MEDIAN,3,0,0);
cvNamedWindow("CheckImage_0",1);
cvShowImage("CheckImage_0",CheckImage_0);
cvSaveImage("CheckImage_0.jpg",CheckImage_0);
#pragmaendregion
cout<<"完成！"<<endl;
cvWaitKey(0);
return0;
}
⑤ 画椭圆的算法
#include<stdio.h>
#include<graphics.h>
#include<math.h>
#include<time.h>
#include<conio.h>
void ellipsepoint(int x,int y,int value,int rx,int ry)
{
putpixel((int)rx+x,(int)ry+y,value);
putpixel((int)rx-x,(int)ry+y,value);
putpixel((int)rx+x,(int)ry-y,value);
putpixel((int)rx-x,(int)ry-y,value);
}
void MidPointEllipse(int a,int b,int value,int rx,int ry)
{
long x=0;
long y=b;
long sa=a*a,sb=b*b;
long xp=(long)((float)sa/(float)sqrt((float)(sa+sb)));
long yp=(long)((float)sb/(float)sqrt((float)(sa+sb)));
long d=sb-sa*(b-0.25);
ellipsepoint(x,y,value,0,0);
while(x<xp)
{
if(d<0)
{ d+=sb*(2*x+3);
x++; }
else
{
d+=sb*(2*x+3)+sa*(-2*y+2);
x++;
y--;
}
ellipsepoint(x,y,value,rx,ry);
}
x=a;y=0;d=4*sa+sb-4*a*sb;
while(y<yp)
{
if(d<0)
{
d+=sa*(2*y+3);
y++;
}
else
{
d+=sa*(2*y+3)+sb*(2-2*x);
y++;
x--;
}
ellipsepoint(x,y,value,rx,ry);
}
}
int main()
{
int gdriver,gmode;
gdriver=VGA;
gmode=VGAHI;
registerbgidriver(EGAVGA_driver);
initgraph(&gdriver,&gmode,"");
MidPointEllipse(50,30,5,100,100);
getch();
closegraph();
return 0;
}
⑥ 求大神给出基于粒子群算法的多目标搜索算法的完整程序。。。从目标函数到最后。。
%% 该函数演示多目标perota优化问题
%清空环境
clc
clear
load data
%% 初始参数
objnum=size(P,1); %类中物品个数
weight=92; %总重量限制
%初始化程序
Dim=5; %粒子维数
xSize=50; %种群个数
MaxIt=200; %迭代次数
c1=0.8; %算法参数
c2=0.8; %算法参数
wmax=1.2; %惯性因子
wmin=0.1; %惯性因子
x=unidrnd(4,xSize,Dim); %粒子初始化
v=zeros(xSize,Dim); %速度初始化
xbest=x; %个体最佳值
gbest=x(1,:); %粒子群最佳位置
% 粒子适应度值
px=zeros(1,xSize); %粒子价值目标
rx=zeros(1,xSize); %粒子体积目标
cx=zeros(1,xSize); %重量约束
% 最优值初始化
pxbest=zeros(1,xSize); %粒子最优价值目标
rxbest=zeros(1,xSize); %粒子最优体积目标
cxbest=zeros(1,xSize); %记录重量，以求约束
% 上一次的值
pxPrior=zeros(1,xSize);%粒子价值目标
rxPrior=zeros(1,xSize);%粒子体积目标
cxPrior=zeros(1,xSize);%记录重量，以求约束
%计算初始目标向量
for i=1:xSize
for j=1:Dim %控制类别
px(i) = px(i)+P(x(i,j),j); %粒子价值
rx(i) = rx(i)+R(x(i,j),j); %粒子体积
cx(i) = cx(i)+C(x(i,j),j); %粒子重量
end
end
% 粒子最优位置
pxbest=px;rxbest=rx;cxbest=cx;
%% 初始筛选非劣解
flj=[];
fljx=[];
fljNum=0;
%两个实数相等精度
tol=1e-7;
for i=1:xSize
flag=0; %支配标志
for j=1:xSize
if j~=i
if ((px(i)<px(j)) && (rx(i)>rx(j))) ||((abs(px(i)-px(j))<tol)...
&& (rx(i)>rx(j)))||((px(i)<px(j)) && (abs(rx(i)-rx(j))<tol)) || (cx(i)>weight)
flag=1;
break;
end
end
end
%判断有无被支配
if flag==0
fljNum=fljNum+1;
% 记录非劣解
flj(fljNum,1)=px(i);flj(fljNum,2)=rx(i);flj(fljNum,3)=cx(i);
% 非劣解位置
fljx(fljNum,:)=x(i,:);
end
end
%% 循环迭代
for iter=1:MaxIt
% 权值更新
w=wmax-(wmax-wmin)*iter/MaxIt;
%从非劣解中选择粒子作为全局最优解
s=size(fljx,1);
index=randi(s,1,1);
gbest=fljx(index,:);
%% 群体更新
for i=1:xSize
%速度更新
v(i,:)=w*v(i,:)+c1*rand(1,1)*(xbest(i,:)-x(i,:))+c2*rand(1,1)*(gbest-x(i,:));
%位置更新
x(i,:)=x(i,:)+v(i,:);
x(i,:) = rem(x(i,:),objnum)/double(objnum);
index1=find(x(i,:)<=0);
if ~isempty(index1)
x(i,index1)=rand(size(index1));
end
x(i,:)=ceil(4*x(i,:));
end
%% 计算个体适应度
pxPrior(:)=0;
rxPrior(:)=0;
cxPrior(:)=0;
for i=1:xSize
for j=1:Dim %控制类别
pxPrior(i) = pxPrior(i)+P(x(i,j),j); %计算粒子i 价值
rxPrior(i) = rxPrior(i)+R(x(i,j),j); %计算粒子i 体积
cxPrior(i) = cxPrior(i)+C(x(i,j),j); %计算粒子i 重量
end
end
%% 更新粒子历史最佳
for i=1:xSize
%现在的支配原有的，替代原有的
if ((px(i)<pxPrior(i)) && (rx(i)>rxPrior(i))) ||((abs(px(i)-pxPrior(i))<tol)...
&& (rx(i)>rxPrior(i)))||((px(i)<pxPrior(i)) && (abs(rx(i)-rxPrior(i))<tol)) || (cx(i)>weight)
xbest(i,:)=x(i,:);%没有记录目标值
pxbest(i)=pxPrior(i);rxbest(i)=rxPrior(i);cxbest(i)=cxPrior(i);
end
%彼此不受支配，随机决定
if ~( ((px(i)<pxPrior(i)) && (rx(i)>rxPrior(i))) ||((abs(px(i)-pxPrior(i))<tol)...
&& (rx(i)>rxPrior(i)))||((px(i)<pxPrior(i)) && (abs(rx(i)-rxPrior(i))<tol)) || (cx(i)>weight) )...
&& ~( ((pxPrior(i)<px(i)) && (rxPrior(i)>rx(i))) ||((abs(pxPrior(i)-px(i))<tol) && (rxPrior(i)>rx(i)))...
||((pxPrior(i)<px(i)) && (abs(rxPrior(i)-rx(i))<tol)) || (cxPrior(i)>weight) )
if rand(1,1)<0.5
xbest(i,:)=x(i,:);
pxbest(i)=pxPrior(i);rxbest(i)=rxPrior(i);cxbest(i)=cxPrior(i);
end
end
end
%% 更新非劣解集合
px=pxPrior;
rx=rxPrior;
cx=cxPrior;
%更新升级非劣解集合
s=size(flj,1);%目前非劣解集合中元素个数
%先将非劣解集合和xbest合并
pppx=zeros(1,s+xSize);
rrrx=zeros(1,s+xSize);
cccx=zeros(1,s+xSize);
pppx(1:xSize)=pxbest;pppx(xSize+1:end)=flj(:,1)';
rrrx(1:xSize)=rxbest;rrrx(xSize+1:end)=flj(:,2)';
cccx(1:xSize)=cxbest;cccx(xSize+1:end)=flj(:,3)';
xxbest=zeros(s+xSize,Dim);
xxbest(1:xSize,:)=xbest;
xxbest(xSize+1:end,:)=fljx;
%筛选非劣解
flj=[];
fljx=[];
k=0;
tol=1e-7;
for i=1:xSize+s
flag=0;%没有被支配
%判断该点是否非劣
for j=1:xSize+s
if j~=i
if ((pppx(i)<pppx(j)) && (rrrx(i)>rrrx(j))) ||((abs(pppx(i)-pppx(j))<tol) ...
&& (rrrx(i)>rrrx(j)))||((pppx(i)<pppx(j)) && (abs(rrrx(i)-rrrx(j))<tol)) ...
|| (cccx(i)>weight) %有一次被支配
flag=1;
break;
end
end
end
%判断有无被支配
if flag==0
k=k+1;
flj(k,1)=pppx(i);flj(k,2)=rrrx(i);flj(k,3)=cccx(i);%记录非劣解
fljx(k,:)=xxbest(i,:);%非劣解位置
end
end
%去掉重复粒子
repflag=0; %重复标志
k=1; %不同非劣解粒子数
flj2=[]; %存储不同非劣解
fljx2=[]; %存储不同非劣解粒子位置
flj2(k,:)=flj(1,:);
fljx2(k,:)=fljx(1,:);
for j=2:size(flj,1)
repflag=0; %重复标志
for i=1:size(flj2,1)
result=(fljx(j,:)==fljx2(i,:));
if length(find(result==1))==Dim
repflag=1;%有重复
end
end
%粒子不同，存储
if repflag==0
k=k+1;
flj2(k,:)=flj(j,:);
fljx2(k,:)=fljx(j,:);
end
end
%非劣解更新
flj=flj2;
fljx=fljx2;
end
%绘制非劣解分布
plot(flj(:,1),flj(:,2),'o')
xlabel('P')
ylabel('R')
title('最终非劣解在目标空间分布')
disp('非劣解flj中三列依次为P，R，C')
⑦ RxRx0.617=公斤是不是钢筋重量的标准算法
是标准算法。正确的说。每米
⑧ rx global mask什么意思
rx-global-mask
接收全球面膜
双语例句
1. There's no Rx for unemployment.
失业问题无法解决。
来自《权威词典》
2. Secondly, RX, GMRF and SEM, three representational algorithms of anomaly detection are studied.
其次, 研究了具有代表性的三种奇异检测算法:RX算法 、 基于高斯马尔可夫随机场模型(GMRF)的检测算法和基于随机最大期望(SEM)分类的检测算法.
⑨ 用什么算法确定保险用户是欺诈用户
1、 数据类型转换：因”gender”实际代表男性、女性用户，故需要将“gender”变量转换成因子变量，且因子水平用“F”替换1，“M”替换2；将“fraudRisk”变量也转换成因子变量。
答：首先将数据导入到R中，并查看数据维度：
> # 导入数据> ccFraud <- read.csv("ccFraud.csv")> # 查看数据维度> str(ccFraud)'data.frame': 10000000 obs. of 9 variables:$ custID : int 1 2 3 4 5 6 7 8 9 10 ...$ gender : int 1 2 2 1 1 2 1 1 2 1 ...$ state : int 35 2 2 15 46 44 3 10 32 23 ...$ cardholder : int 1 1 1 1 1 2 1 1 1 1 ...$ balance : int 3000 0 0 0 0 5546 2000 6016 2428 0 ...$ numTrans : int 4 9 27 12 11 21 41 20 4 18 ...$ numIntlTrans: int 14 0 9 0 16 0 0 3 10 56 ...$ creditLine : int 2 18 16 5 7 13 1 6 22 5 ...$ fraudRisk : int 0 0 0 0 0 0 0 0 0 0 ...
ccFraud数据集一共有一千万行9列，各列均为整型变量。按照题目要求先将“gender”变量转变为因子型，且因子水平用“F”替换1，“M”替换2。代码如下：
> ccFraud$gender <- factor(ifelse(ccFraud$gender==1,'F','M'))> str(ccFraud)'data.frame': 10000000 obs. of 9 variables:$ custID : int 1 2 3 4 5 6 7 8 9 10 ...$ gender : Factor w/ 2 levels "F","M": 1 2 2 1 1 2 1 1 2 1 ...$ state : int 35 2 2 15 46 44 3 10 32 23 ...$ cardholder : int 1 1 1 1 1 2 1 1 1 1 ...$ balance : int 3000 0 0 0 0 5546 2000 6016 2428 0 ...$ numTrans : int 4 9 27 12 11 21 41 20 4 18 ...$ numIntlTrans: int 14 0 9 0 16 0 0 3 10 56 ...$ creditLine : int 2 18 16 5 7 13 1 6 22 5 ...$ fraudRisk : int 0 0 0 0 0 0 0 0 0 0 ...
将“fraudRisk”变量也转换成因子变量，代码如下：
>ccFraud$fraudRisk <- as.factor(ccFraud$fraudRisk)>str(ccFraud)'data.frame': 10000000 obs. of 9 variables: $ custID : int 1 2 3 4 5 6 7 8 9 10 ... $ gender : Factor w/ 2 levels "F","M": 1 2 2 1 1 2 1 1 2 1... $ state : int 35 2 2 15 46 44 3 10 32 23... $ cardholder : int 1 1 1 1 1 2 1 1 1 1 ... $ balance : int 3000 0 0 0 0 5546 2000 60162428 0 ... $ numTrans : int 4 9 27 12 11 21 41 20 4 18... $ numIntlTrans: int 14 0 9 0 16 0 0 3 10 56 ... $ creditLine : int 2 18 16 5 7 13 1 6 22 5 ... $ fraudRisk : Factor w/ 2 levels"0","1": 1 1 1 1 1 1 1 1 1 1 ...
2、 数据探索：查看“fraudRisk”变量中0、1的频数及占比情况。
答：此题是送分题，通过table函数、prop.table函数即可实现。代码如下：
> table(ccFraud$fraudRisk) 0 1 9403986 596014 > prop.table(table(ccFraud$fraudRisk)) 0 1 0.9403986 0.0596014
3、 数据分区：需要按照变量fraudRisk来进行等比例抽样，其中80%作为训练集train数据，20%作为测试集test数据。
答：由于要根据fraudRisk变量进行等比例抽样，我们用caret包的createDataPartition函数实现。代码如下：
> library(caret)载入需要的程辑包：lattice载入需要的程辑包：ggplot2> idx <- createDataPartition(ccFraud$fraudRisk,p=0.8,list=F)> train <- ccFraud[idx,]> test <- ccFraud[-idx,]> prop.table(table(train$fraudRisk)) 0 1 0.94039851 0.05960149 > prop.table(table(test$fraudRisk)) 0 1 0.94039897 0.05960103
4、 建立模型：利用至少三种常见的分类算法（如KNN近邻算法、决策树算法、随机森林等）对数据建立预测模型。
答：由于数据量较大，学员反馈运行慢，这边利用MicrosoftML包来跑模型。关于MRS的快速入门请查阅之前文章：https://ask.hellobi.com/blog/xiejiabiao/8559
> # 模型一：利用MicrosoftML包的rxFastTrees()函数构建快速决策树模型> (a <- Sys.time()) #模型运行前时间[1] "2017-09-03 23:32:04 CST"> treeModel <- rxFastTrees(fraudRisk ~ gender + cardholder + balance + numTrans+ + numIntlTrans + creditLine,data = train)Not adding a normalizer.Making per-feature arraysChanging data from row-wise to column-wiseBeginning processing data.Rows Read: 8000001, Read Time: 0, Transform Time: 0Beginning processing data.Processed 8000001 instancesBinning and forming Feature objectsReserved memory for tree learner: 79664 bytesStarting to train ...Not training a calibrator because it is not needed.Elapsed time: 00:01:04.6222538> (b <- Sys.time()) #模型运行后时间[1] "2017-09-03 23:33:09 CST"> b-a # 模型运行时长Time difference of 1.086313 mins> # 模型二：利用MicrosoftML包的rxFastForest()函数构建快速随机森林模型> (a <- Sys.time()) #模型运行前时间[1] "2017-09-03 23:33:31 CST"> forestModel <- rxFastForest(fraudRisk ~ gender + cardholder + balance + numTrans+ + numIntlTrans + creditLine,data = train)Not adding a normalizer.Making per-feature arraysChanging data from row-wise to column-wiseBeginning processing data.Rows Read: 8000001, Read Time: 0, Transform Time: 0Beginning processing data.Processed 8000001 instancesBinning and forming Feature objectsReserved memory for tree learner: 79664 bytesStarting to train ...Training calibrator.Beginning processing data.Rows Read: 8000001, Read Time: 0, Transform Time: 0Beginning processing data.Elapsed time: 00:01:25.4585776> (b <- Sys.time()) #FastTrees模型运行后时间[1] "2017-09-03 23:34:57 CST"> b-a # 模型运行时长Time difference of 1.433823 mins> # 模型三：利用MicrosoftML包的rxLogisticRegression()函数构建快速逻辑回归模型> (a <- Sys.time()) #模型运行前时间[1] "2017-09-03 23:34:57 CST"> logitModel <- rxLogisticRegression(fraudRisk ~ gender + cardholder + balance + numTrans+ + numIntlTrans + creditLine,data = train)Automatically adding a MinMax normalization transform, use 'norm=Warn' or 'norm=No' to turn this behavior off.Beginning processing data.Rows Read: 8000001, Read Time: 0, Transform Time: 0Beginning processing data.Beginning processing data.Rows Read: 8000001, Read Time: 0, Transform Time: 0Beginning processing data.Beginning processing data.Rows Read: 8000001, Read Time: 0, Transform Time: 0Beginning processing data.LBFGS multi-threading will attempt to load dataset into memory. In case of out-of-memory issues, turn off multi-threading by setting trainThreads to 1.Beginning optimizationnum vars: 8improvement criterion: Mean ImprovementL1 regularization selected 8 of 8 weights.Not training a calibrator because it is not needed.Elapsed time: 00:00:19.5887244Elapsed time: 00:00:00.0383181> (b <- Sys.time()) #模型运行后时间[1] "2017-09-03 23:35:17 CST"> b-a # 模型运行时长Time difference of 20.27396 secs>
逻辑回归模型运行时间最短，消耗20.3秒，其次是决策树，消耗1.08分钟，时间最长的是随机森林，时长为1.4分钟。
5、 模型评估：利用上面建立的预测模型（至少第四点建立的三个模型），对训练集和测试集数据进行预测，评估模型效果，最终选择最优模型作为以后的业务预测模型。（提示：构建混淆矩阵）
答：针对上面的三种模型，我们分别对train、test数据集进行预测并评估。
# 利用决策树模型对数据进行预测，并计算误差率> treePred_tr <- rxPredict(treeModel,data = train)Beginning processing data.Rows Read: 8000001, Read Time: 0, Transform Time: 0Beginning processing data.Elapsed time: 00:00:52.1015119Finished writing 8000001 rows.Writing completed.> t <- table(train$fraudRisk,treePred_tr$PredictedLabel)> t 0 1 0 7446742 76447 1 253008 223804> (paste0(round((sum(t)-sum(diag(t)))/sum(t),3)*100,"%")) #计算决策树对train数据集的预测误差率[1] "4.1%"> treePred_te <- rxPredict(treeModel,data = test)Beginning processing data.Rows Read: 1999999, Read Time: 0, Transform Time: 0Beginning processing data.Elapsed time: 00:00:13.4980323Finished writing 1999999 rows.Writing completed.> t1 <- table(test$fraudRisk,treePred_te$PredictedLabel)> t1 0 1 0 1861406 19391 1 63176 56026> (paste0(round((sum(t1)-sum(diag(t1)))/sum(t1),3)*100,"%")) #计算决策树对test数据集的预测误差率[1] "4.1%"> # 利用随机森林模型对数据进行预测，并计算误差率> forestPred_tr <- rxPredict(forestModel,data = train)Beginning processing data.Rows Read: 8000001, Read Time: 0.001, Transform Time: 0Beginning processing data.Elapsed time: 00:00:56.2862657Finished writing 8000001 rows.Writing completed.> t <- table(train$fraudRisk,forestPred_tr$PredictedLabel)> t 0 1 0 7508808 14381 1 373777 103035> (paste0(round((sum(t)-sum(diag(t)))/sum(t),3)*100,"%")) #计算随机森林对train数据集的预测误差率[1] "4.9%"> forestPred_te <- rxPredict(forestModel,data = test)Beginning processing data.Rows Read: 1999999, Read Time: 0.001, Transform Time: 0Beginning processing data.Elapsed time: 00:00:14.0430130Finished writing 1999999 rows.Writing completed.> t1 <- table(test$fraudRisk,forestPred_te$PredictedLabel)> t1 0 1 0 1877117 3680 1 93419 25783> (paste0(round((sum(t1)-sum(diag(t1)))/sum(t),3)*100,"%")) #计算随机森林对test数据集的预测误差率[1] "1.2%"> # 利用逻辑回归模型对数据进行预测，并计算误差率> logitPred_tr <- rxPredict(logitModel,data = train)Beginning processing data.Rows Read: 8000001, Read Time: 0.001, Transform Time: 0Beginning processing data.Elapsed time: 00:00:08.1674394Finished writing 8000001 rows.Writing completed.> t <- table(train$fraudRisk,logitPred_tr$PredictedLabel)> t 0 1 0 7444156 79033 1 250679 226133> (paste0(round((sum(t)-sum(diag(t)))/sum(t),3)*100,"%")) #计算逻辑回归对train数据集的预测误差率[1] "4.1%"> logitPred_te <- rxPredict(logitModel,data = test)Beginning processing data.Rows Read: 1999999, Read Time: 0, Transform Time: 0Beginning processing data.Elapsed time: 00:00:02.0736547Finished writing 1999999 rows.Writing completed.> t1 <- table(test$fraudRisk,logitPred_te$PredictedLabel)> t1 0 1 0 1860885 19912 1 62428 56774> (paste0(round((sum(t1)-sum(diag(t1)))/sum(t),3)*100,"%")) #计算逻辑回归对test数据集的预测误差率[1] "1%"
从训练集和测试集的预测误差率来看，对于此份数据，逻辑回归是最优的选择。