MATLAB绘图基础
 
2 MATLAB的立体绘图
2.1 生成网格坐标矩阵的函数
2.2 画三维曲面的函数
2.3 三维旋转体的绘制
柱面图
球面图
    
2.4 三维图形的处理
视点处理
色彩处理
消隐处理
裁剪处理
   
   
❤隐函数作图
 

 
2 MATLAB的立体绘图
 
2.1 生成网格坐标矩阵的函数
 
[X,Y]=meshgrid(x,y)
生成X-Y平面的网格坐标矩阵
[X,Y,Z]=sphere(n)
生成球面的三维坐标矩阵
[X,Y,Z]=cylinder(r,n)
生成柱面的三维坐标矩阵
[X,Y]=meshgrid(x,y);Z=peaks(X,Y)
peaks多峰函数，生成多峰曲面的坐标矩阵
例子：
 
%例一
theta=0:pi/50:6*pi;
x=cos(theta); y=sin(theta); z=0:300;
plot3(x,y,z);
%例二
x=-3:0.1:3; y=-3:0.1:3;
[X,Y]=meshgrid(x,y); Z=X.^2+Y.^2;
surf(X,Y,Z);
 
2.2 画三维曲面的函数
 
plot3(X,Y,Z)
三维曲面
mesh(X,Y,Z)
三维网格曲面
meshz(X,Y,Z)
可将曲面加上围裙
meshc(X,Y,Z)
同时画出网状图与等高线
surf(X,Y,Z)
三维填充曲面（更精细）
surfc(X,Y,Z)
同时画出曲面图与等高线
bar3(A)
三维柱状图
stem3(A)
三维棒状图
pie3(A)
三维饼状图
fill3(A)
三维填充图
waterfall(X,Y,Z)
三维瀑布图，可在x方向或y方向产生水流效果
contour3(X,Y,Z)
三维等高线图
contour(X,Y,Z)
画出曲面等高线在XY平面的投影
例子：
 
%mesh和plot是三度空间立体绘图的基本命令，mesh可画出立体网状图，plot则可画出立体曲面图，两者产生的图形都会依高度而有不同颜色。
%plot3可画出三度空间中的曲线。
t=linspace(0,20*pi, 501);   
subplot(3,4,1),plot3(t.*sin(t), t.*cos(t), t);title('plot3');  
%也可同时画出两条三度空间中的曲线 
t=linspace(0, 10*pi, 501);  
subplot(3,4,2),plot3(t.*sin(t), t.*cos(t), t, t.*sin(t), t.*cos(t), -t);title('plot3');
%mesh可画出立体网状图。
x=linspace(-2, 2, 25); %在x轴上取25点  
y=linspace(-2, 2, 25); %在y轴上取25点  
[xx,yy]=meshgrid(x, y); %xx和yy都是21x21的矩阵  
zz=xx.*exp(-xx.^2-yy.^2); %计算函数值，zz也是21x21的矩阵  
subplot(3,4,3),mesh(xx, yy, zz);title('mesh'); %画出立体网状图
%peaks函数可产生一个凹凸有致的曲面，包含了三个局部极大点及三个局部极小点。  
subplot(3,4,4),peaks;
%meshz可将曲面加上围裙。
[x,y,z]=peaks;  
subplot(3,4,5),meshz(x,y,z);title('meshz');
axis([-inf inf -inf inf -inf inf]);  
%meshc同时画出网状图与等高线。 
[x,y,z]=peaks;  
subplot(3,4,6),meshc(x,y,z);title('meshc');  
axis([-inf inf -inf inf -inf inf]);   
%surf画出立体曲面图
x=linspace(-2, 2, 25); %在x轴上取25点  
y=linspace(-2, 2, 25); %在y轴上取25点  
[xx,yy]=meshgrid(x, y); %xx和yy都是21x21的矩阵   
zz=xx.*exp(-xx.^2-yy.^2); %计算函数值，zz也是21x21的矩阵   
subplot(3,4,7),surf(xx, yy, zz);title('surf'); %画出立体曲面图
%surfc同时画出曲面图与等高线。  
[x,y,z]=peaks;  
subplot(3,4,8),surfc(x,y,z);title('surfc');
axis([-inf inf -inf inf -inf inf]); 
%waterfall可在x方向或y方向产生水流效果。 
[x,y,z]=peaks;  
subplot(3,4,9),waterfall(x,y,z);title('waterfall-x');
axis([-inf inf -inf inf -inf inf]); 
%在y方向产生水流效果 
[x,y,z]=peaks;  
subplot(3,4,10),waterfall(x',y',z');title('waterfall-y');
axis([-inf inf -inf inf -inf inf]);  
%contour3画出曲面在三度空间中的等高线。
subplot(3,4,11),contour3(peaks, 20);title('contour3');
axis([-inf inf -inf inf -inf inf]);
%contour画出曲面等高线在XY平面的投影。
subplot(3,4,12),contour(peaks, 20);title('contour');
 
 
2.3 三维旋转体的绘制
 
柱面图
 
由cylinder函数实现，调用格式为：
 
[x,y,z]=cylinder(R,n); %其中R是一个向量，存放柱面各个等间隔高度上的半径，n表示在圆柱圆周上有n个间隔点，默认有20个间隔点。
如：cylinder(3)生成一个圆柱，cylinder([10,1])生成一个圆锥，而t=0:pi/100:4*pi; R=sin(t); cylinder(R,30);生成一个正弦圆柱面。
[X,Y,Z]=cylinder(R)或[X,Y,Z]=cylinder此形式为默认N=20且R=[1 1]
 
例子：
 
t=0:pi/20:2*pi;
[x,y,z]=cylinder(2+sin(t),30);
subplot(1,3,1);surf(x,y,z);
axis('equal');axis('square'); %控制坐标轴的大小相同
subplot(1,3,2);
x=0:pi/20:pi*3;
r=5+cos(x);
[a,b,c]=cylinder(r,30);
mesh(a,b,c);
axis('equal');axis('square');
subplot(1,3,3);
r=abs(exp(-0.25*t).*sin(t));
t=0:pi/12:3*pi;
r=abs(exp(-0.25*t).*sin(t));
[X,Y,Z]=cylinder(r,30);
mesh(X,Y,Z);
colormap([1 0 0]);
axis('equal');axis('square');
 
 
球面图
 
由sphere函数实现，调用格式为：
 
[X,Y,Z]=sphere(N); %此函数生成3个(N+1)*(N+1)的矩阵,利用函数surf(X,Y,Z) 可绘制出圆心位于原点、半径为1的单位球体。
n决定了球面的圆滑程度，其默认值为20。若n值取的比较小，则绘制出多面体的表面图。
[X,Y,Z]=sphere; %此形式使用了默认值N=20。
sphere(N); %只是绘制了球面图而不返回任何值。
 
例子：
 
subplot(1,2,1);
[x,y,z]=sphere;
surf(x,y,z);
axis('equal');axis('square'); %控制坐标轴的大小相同
subplot(1,2,2);
[a,b,c]=sphere(40);
t=abs(c);
surf(a,b,c,t);
axis('equal');axis('square');
colormap('hot');
 
 
2.4 三维图形的处理
 
视点处理
 
设置视点的函数view，调用格式为：
 
view(az,el); %其中az为方位角，el为仰角，它们均以度为单位。系统默认的视点定义为方位角为-37.5度，仰角30度。
 
例子：
 
subplot(2,2,1);mesh(peaks);view(-37.5,30);title('-37.5°,30°');
subplot(2,2,2);mesh(peaks);view(0,90);title('0°,90°');
subplot(2,2,3);mesh(peaks);view(90,0);title('90°,0°');
subplot(2,2,4);mesh(peaks);view(-7,-10);title('-7°,-10°');
 
 
色彩处理
 
内建矩阵
 
colormap hot
 
三维图形表面的着色
 
shading faceted; %将每个网格片用其高度对应的颜色进行着色，网格线是黑色
shading flat; %将每个网格片用同一个颜色进行着色，且网格线也用相应的颜色
shading interp; %在网格片内采用颜色插值处理
 
消隐处理
 
例子：
 
z=peaks(50);
subplot(2,1,1);
mesh(z);title('消隐前的网图')
hidden off
subplot(2,1,2)
mesh(z);title('消隐后的网图')
hidden on
colormap([0 0 1]);
 
 
裁剪处理
 
将图形中需要裁剪部分对应的函数值设置成NaN，使函数值为NaN的部分将不显示出来，从而达到对图形进行裁剪的目的。
 例子：
 
%削掉正弦波顶部或底部大于0.5的部分。
subplot(1,2,1);
x=0:pi/10:4*pi;
y=sin(x);
i=find(abs(y)>0.5);
x(i)=NaN;
plot(x,y);
axis([0 15 -0.4 0.4]);axis square
%绘制两个球面，其中一个在另一个里面，将外面的球裁掉一部分，以便能看到里面的球。
subplot(1,2,2);
[x,y,z]=sphere(25);%生成外面的大球
z1=z;
z1(:,1:4)=NaN;%将大球裁去一部分
c1=ones(size(z1));
surf(3*x,3*y,3*z1,c1);       %生成里面的小球
hold on
z2=z;
c2=2*ones(size(z2));
c2(:,1:4)=3*ones(size(c2(:,1:4)));
surf(1.5*x,1.5*y,1.5*z2,c2);
colormap([0 1 0;0.5 0 0;1 0 0]);
grid on;hold off;axis equal
 
 
❤隐函数作图
 
Matlab提供了ezplot函数绘制隐函数图形，用法如下：
 ①对于函数f=f(x)，ezplot的调用格式为：
 
ezplot(f); %在默认区间(-2pi,2pi)绘制图形
ezplot(f,[a,b]); %在区间(a,b)绘制图形
 
②对于隐函数f=f(x,y)，ezplot的调用格式为：
 
ezplot(f); %在默认区间(-2pi,2pi)绘制f(x,y)=0的图形
ezplot(f,[xmin,xmax,ymin,ymax]); %在区间[xmin,xmax,ymin,ymax]绘制图形
ezplot(f,[a,b]); %在区间(a,b)绘制图形
 
③对于参数方程x=x(t)，y=y(t)，ezplot函数的调用格式为：
 
ezplot(x,y); %在默认区间绘制x=x(t),y=y(t)图形
ezplot(x,y,[tmin,tmax]); %在区间(tmin,tmax)绘制x=x(t),y=y(t)图形
 
其他隐函数绘图还有ezpolar，ezcontour，ezplot3，ezmesh，ezmeshc，ezsurf，ezsurfc。
 例子：
 
subplot(2,2,1);
ezplot('x^2+y^2-9');axis equal;
subplot(2,2,2);
ezplot('x^3+y^3-5*x*y+1/5')
subplot(2,2,3);
ezplot('cos(tan(pi*x))',[0,1]);
subplot(2,2,4);
ezplot('8*cos(t)','4*sqrt(2)*sin(t)',[0,2*pi]);
 

参考材料 | 文献:
 代码参考：Heart Model for Valentine’s Day
 创意来源：Matlab 可以做一些什么有意思的事？
 参考文献：Wolfram MathWorld Heart Surface
 
 
代码运行效果:
 
 完整代码:
 
function heartSudoku2
sizeTheta=30;
sizeFai=40;
theta=linspace(0,pi,sizeTheta)';
nudge=0.0001; % used to avoid the operlapping
fai=linspace(0+nudge,pi/2-nudge,round(sizeFai/4));
a=9/4;
b=9/80;
A=1+(a-1)*(sin(theta).^2)*(sin(fai).^2);
B=(sin(theta).^2.*cos(theta).^3) * (1+(b-1)*(sin(fai).^2));
rou = zeros(size(A));
for iLoop=1:numel(A)
    curA=A(iLoop);
    curB=B(iLoop);
    % this is the polar coordinates version of the sextic equation found on
    % http://mathworld.wolfram.com/HeartSurface.html
    polyFactors=[curA^3,-curB,-3*curA^2,0,3*curA,0,-1];
    solutions=roots(polyFactors);    
    realRou=real(solutions(abs(imag(solutions))<1e-9));
    rou(iLoop)=realRou(realRou>0);    
end

% x,y,z for the quater of whole heart
x=rou.*(sin(theta)*cos(fai));
y=rou.*(sin(theta)*sin(fai));
z=rou.*(cos(theta)*ones(size(fai)));

% x,y,z for the whole heart
cordX=[x,-fliplr(x),-x,fliplr(x)];
cordY=[y,fliplr(y),-y,-fliplr(y)];
cordZ=[z,fliplr(z),z,fliplr(z)];

heartSudokuFig=figure('units','pixels',...
    'position',[200 50 630 630],...
    'Numbertitle','off',...
    'menubar','none',...
    'resize','off',...
    'name','heart sudoku',...
    'color',[0.99 0.99 0.99]);
for i=1:3
    for j=1:3
      heartSudokuAxes(i,j)=axes('Units','pixels',...
      'parent',heartSudokuFig,...
      'Position',[(i-1)*210-25,(j-1)*210-25-10,210+50,210+50],...
      'Color',[0.99 0.99 0.99],...
      'Box','on', ...
      'XLim',[-1.28,1.28],'YLim',[-1.28,1.28],'ZLim',[-1.28,1.28],...
      'XColor','none','YColor','none','ZColor','none');
      hold(heartSudokuAxes(i,j),'on')
      camtarget([0, 0, 0]);
      view(30,20);
      
    end
end
markerColor=[248 150 170]./255;
lineColor=[248 150 170]./255;
lineWidth=0.8;
%[1 1]=====================================================================
for iLoop=1:3:size(cordX,1)
    plot3(heartSudokuAxes(1,1),cordX(iLoop,:),cordY(iLoop,:),cordZ(iLoop,:),...
        'Marker','.','MarkerSize',3.5,'MarkerEdgeColor',markerColor,'LineStyle','none');
    plot3(heartSudokuAxes(1,1),cordX(iLoop,:),cordY(iLoop,:),cordZ(iLoop,:),...
        'color',lineColor,'LineStyle','-','LineWidth',lineWidth);
end
for iLoop=1:4:size(cordX,2)
    plot3(heartSudokuAxes(1,1),cordX(:,iLoop),cordY(:,iLoop),cordZ(:,iLoop),...
        'color',lineColor,'LineStyle','-','LineWidth',lineWidth);
end
%[1 2]=====================================================================
for iLoop=1:3:size(cordX,1)
    plot3(heartSudokuAxes(1,2),cordX(iLoop,:),cordY(iLoop,:),cordZ(iLoop,:),...
        'Marker','.','MarkerSize',3.5,'MarkerEdgeColor',markerColor,'LineStyle','none');
    plot3(heartSudokuAxes(1,2),cordX(iLoop,:),cordY(iLoop,:),cordZ(iLoop,:),'color',...
        lineColor,'LineStyle','-','LineWidth',lineWidth);
end
for iLoop=1:size(cordX,2)
    plot3(heartSudokuAxes(1,2),cordX(:,iLoop),cordY(:,iLoop),cordZ(:,iLoop),'color',...
        lineColor,'LineStyle','-','LineWidth',lineWidth);
end
%[1 3]=====================================================================
for iLoop=1:size(cordX,1)
    plot3(heartSudokuAxes(1,3),cordX(iLoop,:),cordY(iLoop,:),cordZ(iLoop,:),...
        'Marker','.','MarkerSize',3.5,'MarkerEdgeColor',markerColor,'LineStyle','none');
    plot3(heartSudokuAxes(1,3),cordX(iLoop,:),cordY(iLoop,:),cordZ(iLoop,:),...
        'color',lineColor,'LineStyle','-','LineWidth',lineWidth);
end
for iLoop=1:size(cordX,2)
    plot3(heartSudokuAxes(1,3),cordX(:,iLoop),cordY(:,iLoop),cordZ(:,iLoop),...
        'color',lineColor,'LineStyle','-','LineWidth',lineWidth);
end
%[2 1]=====================================================================
for iLoop=[2 8 15 25]
    plot3(heartSudokuAxes(2,1),cordX(iLoop,:),cordY(iLoop,:),cordZ(iLoop,:),...
        'Marker','.','MarkerSize',3.5,'MarkerEdgeColor',markerColor,'LineStyle','none');
    plot3(heartSudokuAxes(2,1),cordX(iLoop,:),cordY(iLoop,:),cordZ(iLoop,:),...
        'color',lineColor,'LineStyle','-','LineWidth',lineWidth);
end
for iLoop=1:2:size(cordX,2)
    plot3(heartSudokuAxes(2,1),cordX(:,iLoop),cordY(:,iLoop),cordZ(:,iLoop),...
        'color',lineColor,'LineStyle','-','LineWidth',lineWidth);
end
%[2 2]=====================================================================
for iLoop=1:2:size(cordX,1)
    plot3(heartSudokuAxes(2,2),cordX(iLoop,:),cordY(iLoop,:),cordZ(iLoop,:),...
        'Marker','.','MarkerSize',3.5,'MarkerEdgeColor',markerColor,'LineStyle','none');
    plot3(heartSudokuAxes(2,2),cordX(iLoop,:),cordY(iLoop,:),cordZ(iLoop,:),...
        'color',lineColor,'LineStyle','-','LineWidth',lineWidth);
end
for iLoop=1:2:size(cordX,2)
    plot3(heartSudokuAxes(2,2),cordX(:,iLoop),cordY(:,iLoop),cordZ(:,iLoop),...
        'color',lineColor,'LineStyle','-','LineWidth',lineWidth);
end
%[2 3]=====================================================================
for iLoop=[1 3 5 9 20 25 27]
    plot3(heartSudokuAxes(2,3),cordX(iLoop,:),cordY(iLoop,:),cordZ(iLoop,:),...
        'Marker','.','MarkerSize',3.5,'MarkerEdgeColor',markerColor,'LineStyle','none');
    plot3(heartSudokuAxes(2,3),cordX(iLoop,:),cordY(iLoop,:),cordZ(iLoop,:),...
        'color',lineColor,'LineStyle','-','LineWidth',lineWidth);
end
for iLoop=1:size(cordX,2)
    plot3(heartSudokuAxes(2,3),cordX(:,iLoop),cordY(:,iLoop),cordZ(:,iLoop),...
        'color',lineColor,'LineStyle','-','LineWidth',lineWidth);
end
%[3 1]=====================================================================
for iLoop=[2 6 8 15 ]
    plot3(heartSudokuAxes(3,1),cordX(iLoop,:),cordY(iLoop,:),cordZ(iLoop,:),...
        'Marker','.','MarkerSize',3.5,'MarkerEdgeColor',markerColor,'LineStyle','none');
    plot3(heartSudokuAxes(3,1),cordX(iLoop,:),cordY(iLoop,:),cordZ(iLoop,:),...
        'color',lineColor,'LineStyle','-','LineWidth',lineWidth);
end
for iLoop=[8 9 18 21 30 32 35 40]
    plot3(heartSudokuAxes(3,1),cordX(:,iLoop),cordY(:,iLoop),cordZ(:,iLoop),...
        'color',lineColor,'LineStyle','-','LineWidth',lineWidth);
end
%[3 2]=====================================================================
for iLoop=[8 12 15 18 23 25]
    plot3(heartSudokuAxes(3,2),cordX(iLoop,:),cordY(iLoop,:),cordZ(iLoop,:),...
        'Marker','.','MarkerSize',3.5,'MarkerEdgeColor',markerColor,'LineStyle','none');
    plot3(heartSudokuAxes(3,2),cordX(iLoop,:),cordY(iLoop,:),cordZ(iLoop,:),...
        'color',lineColor,'LineStyle','-','LineWidth',lineWidth);
end
for iLoop=[8 9 18 21 25 28 30 32 35 40]
    plot3(heartSudokuAxes(3,2),cordX(:,iLoop),cordY(:,iLoop),cordZ(:,iLoop),...
        'color',lineColor,'LineStyle','-','LineWidth',lineWidth);
end
%[3 3]=====================================================================
for iLoop=[8 18 25]
    plot3(heartSudokuAxes(3,3),cordX(iLoop,:),cordY(iLoop,:),cordZ(iLoop,:),...
        'Marker','.','MarkerSize',3.5,'MarkerEdgeColor',markerColor,'LineStyle','none');
    plot3(heartSudokuAxes(3,3),cordX(iLoop,:),cordY(iLoop,:),cordZ(iLoop,:),...
        'color',lineColor,'LineStyle','-','LineWidth',lineWidth);
end
for iLoop=1:3:40
    plot3(heartSudokuAxes(3,3),cordX(:,iLoop),cordY(:,iLoop),cordZ(:,iLoop),...
        'color',lineColor,'LineStyle','-','LineWidth',lineWidth);
end
end
 
 
注：
 由于各心形的形状要体现随机性，故各个心形的数据是单独编辑的，如果想要获得规整的心形，不妨采用如下方法：
 
 效果：
 
 另：
 对于参考代码，源文件中附带以下版权声明，依据声明中要求在此进行展示，同时提醒代码在进行再次发布和利用时须附带此声明：
 
 
 
Copyright © 2009, Xin Zhao
 All rights reserved.
 Redistribution and use in source and binary forms, with or without
 modification, are permitted provided that the following conditions are
 met:
 
 
* Redistributions of source code must retain the above copyright
  notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
  notice, this list of conditions and the following disclaimer in
  the documentation and/or other materials provided with the distribution
* Neither the name of the The MathWorks, Inc. nor the names
  of its contributors may be used to endorse or promote products derived
  from this software without specific prior written permission.
 
 
 
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS”
 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 POSSIBILITY OF SUCH DAMAGE.
 

 最近在做基于双目视觉的三维重建。比较opencv和matlab工具箱的立体标定结果精度时，发现貌似如果手工选取角点不那么离谱的话，matlab标定结果精度更高也更鲁棒。就想先用matlab标定好相机，再把结果供opencv函数加载使用。如何将Matlab标定结果的.mat文件转成需要的CvMat矩阵，就是本篇博客所要讲的。
 
 主要参考：http://www.jianshu.com/p/ad6a2f8a3fc8
 
 这里分以下三步来讲：
 
 1、先在matlab中加载.mat文件成matlab中的矩阵;
 
 (1)生成MAT文件
 eg1:
 >>A=[1 2;3 4];%matlab中矩阵按列保存而非opencv中按行保存
 >>save mydata A;%mydata表示文件名
 
 eg2:
 >>A=[1 2;3 4];
 >>B=[2 3;5 6];
 >>save mydata A B;%保存多个矩阵直接将矩阵名都写上,简单方法是save mydata;回车即可则将运行空间中的所有变量都保存到文件中.
 
 (2)从.mat文件中读取所有数据
 首先将想打开的mat文件所在的目录设置为当前工作目录(cd 路径);
 然后执行：load mydata;即可,想看哪个变量的值直接输入变量名回车即可.
 
 若只想读取mat文件指定数据,eg:load mydata A;则运行空间中将不出现B
 
 2、再将该矩阵按yml格式保存成yml文件；
 matlab中int类型对应i,float类型对应f,double类型对应d.
 定义了函数的.m脚本文件,若想在命令行窗口中直接调用,即使当前运行路径与脚本文件相同也需要将脚本文件所在目录添加到设置路径中才能够直接调用，否则会说该函数未定义.(在脚本中调用某个脚本函数,只需要将两个脚本文件放在同一目录下即可)
 
 *************************创建脚本matlab2opencv.m**************************************
 
 %该函数将matlab中的mat矩阵保存成yml文件
 
 function matlab2opencv( variable, fileName, flag)
 
 [rows cols] = size(variable);
 
 % Beware of Matlab's linear indexing
 variable = variable';
 
 % Write mode as default
 if ( ~exist('flag','var') )
     flag = 'w';
 end
 
 if ( ~exist(fileName,'file') || flag == 'w' )
     % New file or write mode specified
     file = fopen( fileName, 'w');%不存在则创建写入模式
     fprintf( file, '%%YAML:1.0\n');
 else
     % Append mode
     file = fopen( fileName, 'a');%存在则追加模式
 end
 
 % Write variable header
 fprintf( file, '%s: !!opencv-matrix\n', inputname(1));
 fprintf( file, '    rows: %d\n', rows);
 fprintf( file, '    cols: %d\n', cols);
 fprintf( file, '    dt: d\n');%double类型
 fprintf( file, '    data: [ ');
 
 % Write variable data
 for i=1:rows*cols
     fprintf( file, '%.16d', variable(i));%16表示小数点后有16位
     if (i == rows*cols), break, end
     if mod(i+1,4) == 0
         fprintf( file, ',\n        ');
     else
         fprintf( file, ', ');
     end
 end
 
 fprintf( file, ']\n');
 
 fclose(file);
 
 ************************************************************************
 
 接着就可以从.mat立体标定结果文件中提取需要的信息保存到CvMat中：
 
 *******************************创建脚本savemat2yml.m****************************
 %利用matlab2opencv函数从matlab标定结果矩阵中提取所需的数值保存到对应的yml文件中(按opencv函数所需矩阵的顺序设置值),以供opencv直接利用cvLoad等加载调用.
 
 cd D:\MATLAB\R2014a\toolbox\TOOLBOX_calib\calib_images;%切换到立体标定结果的目录下
 load Calib_Results_stereo.mat;%加载标定结果矩阵
 
 %左摄像机内参数矩阵3x3
 %fx 0 cx
 %0 fy cy
 %0 0 1
 M1(1,1)=fc_right(1); M1(1,2)=0;           M1(1,3)=cc_right(1);
 M1(2,1)=0;           M1(2,2)=fc_right(2); M1(2,3)=cc_right(2);
 M1(3,1)=0;           M1(3,2)=0;           M1(3,3)=1;
 
 %左畸变参数向量1x5(貌似也可以是5x1则直接用),记录的顺序为k1,k2,p1,p2,k3
 D1(1,1)=kc_right(1); 
 D1(1,2)=kc_right(2); 
 D1(1,3)=kc_right(3); 
 D1(1,4)=kc_right(4); 
 D1(1,5)=kc_right(5); %一般为0
 
 %右摄像机内参数矩阵3x3
 %fx 0 cx
 %0 fy cy
 %0 0 1
 M2(1,1)=fc_right(1); M2(1,2)=0;           M2(1,3)=cc_right(1);
 M2(2,1)=0;          M2(2,2)=fc_right(2);  M2(2,3)=cc_right(2);
 M2(3,1)=0;          M2(3,2)=0;           M2(3,3)=1;
 
 %右畸变参数向量1x5(貌似也可以是5x1则直接用),记录的顺序为k1,k2,p1,p2,k3
 D2(1,1)=kc_right(1); 
 D2(1,2)=kc_right(2); 
 D2(1,3)=kc_right(3); 
 D2(1,4)=kc_right(4); 
 D2(1,5)=kc_right(5);
 
 %摄像机间的旋转矩阵3x3,可直接用
 %摄像机间的平移矩阵3x1,也可直接用
 
 %调用matlab2opencv函数保存矩阵到yml文件
 matlab2opencv(M1,'M1.yml');%注意该函数脚本要与本脚本在同一目录下或者该函数脚本已设置路径了
 matlab2opencv(D1,'D1.yml');
 matlab2opencv(M2,'M2.yml');
 matlab2opencv(D2,'D2.yml');
 matlab2opencv(R,'R.yml');
 matlab2opencv(T,'T.yml');
 ***********************************************************************
 
 3、最后opencv加载该yml文件，并以Mat、CvMat的格式保存到一个变量中.
 若yml文件中只有一个矩阵可以直接调用cvSave和cvLoad保存和加载yml文件.eg:
 
 
 

  
 
 
 
 1 #include <cv.h>
 2 #include <highgui.h>
 3 #include<iostream>
 4 using namespace std;  5  6 ostream& operator<<(ostream& out,CvMat* mat)  7 {  8 cout<<"C++方式输出矩阵:"<<endl;  9 cout<<"mat.rows="<<mat->rows<<",mat.cols="<<mat->cols<<"\n["<<endl; 10 for(int i=0;i<mat->rows;i++) 11  { 12 for(int j=0;j<mat->cols;j++) 13  { 14 //设置输出宽度为25,按左对齐,数值的有效数字位数为20(整数+小数) 15 cout<<setw(25)<<setiosflags(ios::left)<<setprecision(20)<<CV_MAT_ELEM(*mat,double,i,j); 16  } 17 if(i==mat->rows-1) 18 cout<<"\n]\n\n"; 19 else 20 cout<<endl; 21  } 22 return out; 23 } 24 25 26 void main() 27 { 28 CvMat* M1=(CvMat*)cvLoad("C:/Users/shark/Desktop/M1.yml"); 29 cout<<M1<<endl; 30  getchar(); 31 }
 
 

  
 
 
 
 博主matlab入门级,解决该问题还是花费了一定时间,望本博客对可能同是菜鸟的你有所帮助。

 
论文写作指导：请加QQ229366758
 
[摘 要] 本文是对一个静态的情况下货品分配货位的探讨；利用权重赋值的方法把多目标问题转化为求解一个目标的情况进行求解；并使用matlab软件进行迭代求解，得到一个很好的收敛优化解。
 
下载论文网 http://www.xzlunwen.com
 
[关键词] matlab 立体仓库 货位优化 仿真
 
前言
 
自动化立体仓库又称为自动化仓储AS/RS(Automatic Storage and Retrieval System)系统，它可以按照指令自动完成货物的存取，并能对库存货物进行自动管理，实现自动化作业，具有作业效率高，占地面积小，存储容量大，计算机自动控制等优势，受到企业越来越多的重视。本文主要是就货品静态货位分配进行优化，以达到提高仓库存储效率的目的。
 
1.基本假设：
 
1)货架是矩形，货架的I/O处位于货架的左下角。货架共有N列，M层，货位长度为w，高度为h。
 
2)系统中有多种货物，且每个托盘只存取有一种类型的货物，并且所有的托盘是相同的。货架上所有货位在尺寸上都是相同的，并且任意一个托盘都可以被存放于货架上的任意一个货位。
 
3)堆垛机可以在水平方向上和垂直方向上同时运动。为了简化处理，设堆垛机只有一种速度，即其最大速度，最大速度分别为Vh、Ah。
 
2.立体仓库基本参数：
 
堆垛机运行水平速度：120m/min；托盘货位高：1050mm；堆垛机运行垂直速度：30m/min；货架层数：10；托盘货位长：1300mm；货架列数：20；托盘货位宽：1100mm；托盘型号：1000*1200mm。
 
1 建立目标函数
 
1.1 货架货位编号
 
根据堆垛机运行速度原则为货位进行编号，堆垛机水平运行速度120m/min、垂直运行速度为30m/min。假定条件为堆垛机水平和垂直加速度一样，则根据货架货位到达出入库口的堆垛机运行时间tij=[L*i/vx, H*(j-1)/vy]为在不考虑堆垛机启动和制动的情况下,将第i列第j层货位上的货品搬运到出入库台所用的时间；vx，vy为堆垛机的水平和垂直运行速度；L，H为单个货位的长度和高度)由短到长进行货架的货位编号，如果堆垛机垂直运行和水平运行时间一样则以水平运行时间优先编号，得到编号如下：
 
其中左下侧为货物I/O侧。
 
1.2 货品编号
 
根据货品的COI值为Ii(COI为立方体索引号，Heskett给出的COI的计算公式为：Ii=Ci/fi，其中Ci为某种货品储存总量所需的库存容量；fi为某种货品的出库频率)货品COI值越小，出入库频率越高，应该离出入口越近，所给编号也就越小。
 
使用COI值是因为在货品进行存储时，货架上由一个初始的存货量，虽然某货品出入库频率低，但如果其初始存储量比较大，则把其放在O/I近的货位会使后续进出库货物所用时间更长。
 
出入库货品频率和容量：
 
出入库频率――10类货物相对比的频率，他们的加权和为1；
 
所占货格数――表示当该类货物存储量最大时需要占据的货格数量。
 
1.3 优化函数
 
根据出入库堆垛机运行时间最少原则，当堆垛机运行时间t和货品COI的乘积之和Q最小时(出入库频率高的货品对应近的货位)，我们得到一个优化目标条件：
 
minQ= tij*COIij
 
tij―堆垛机运行到第i层第j列货位时运行时间；
 
COIij―第i层第j列货位对应的货品COI值。
 
根据上轻下重的原则，仓库每排货架共有10层20列，设地面层为第一层，离出入口最近的列为第一列，处在第i层第j列的货品重量记为Wij，我们货位分配优化目标是每个托盘上货品的质量与其所在层的乘积之和S最小。得到第二个优化目标条件：
 
minS= Wij*(i-1)
 
由上述优化函数可以知道，货位分配需要同时考虑货架的稳定性和存取效率，这是个组合多目标优化问题。
 
2 货位分配的计算过程
 
1)把多目标优化问题简化成单目标优化
 
(1)和(2)两目标函数分别是求堆垛机利用率最高和货架重心最低的函数，我们通过对量目标函数赋给权重，不同的仓库类型对两目标的重要性具有不同的要求，这里为了研究的需要，赋给权重相等，分别为0.5.那么多目标函数就变成：
 
minH=S+Q= Wij*(i-1)+ tij*COIij
 
2)进行初始分配，给货品随机赋给一个货位，根据表1-2的内容，假设货架满货时(也就是没有空货位存在)。不同的货品所占货位格数已给出。
 
3)货品赋给的货位正确与否的判断条件，任意选取第i行第j列货品其对应的重量为Wij，COI值为COIij，使用该货品与第m行第n列货品其进行比较对应的重量为Wmn，COI值为COImn，把第i行第j列货品的H值与第m行第n列货品的H值进行求和得到式(5)：
 
he1=Wij*(i-1)+tij*COIij+Wmn*(m-1)+tmn*COImn
 
然后再把第i行第j列货品与第m行第n列货品互换，同时计算其和值得到式(6)：
 
he2=Wij*(m-1)+tmn*COIij+Wmn*(i-1)+tij*COImn
 
把式4和式5进行比较，如果式4大于式5，则第i行第j列货品与第m行第n列货品不进行互换，反之，进行互换。
 
4)重复第三步骤，直到每个货位上的货品与其他每个货位上货品进行比较完成，这样得到一个最优的货位分配方案，得出最终解。
 
3 基于matlab的静态货位分配优化仿真
 
3.1 求解
 
通过上面求解过程，在Matlab中进行编码计算，计算代码如下：
 
%初始化
 
for i=1:10
 
for j=1:20
 
for m=i:10
 
for n=1:20
 
chuizhi=i*1.05/0.3;
 
shuiping=j*1.1/1.20;
 
if chuizhi>shuiping
 
p=chuizhi;
 
else p=shuiping;
 
end
 
chuizhi=m*1.05/0.30;
 
shuiping=n*1.1/1.20;
 
if chuizhi>shuiping
 
q=chuizhi;
 
else q=shuiping;
 
end
 
heij=p*COI(i,j)+WEIGHT(i,j)*(i-1)*1.05+q*COI(m,n)+WEIGHT(m,n)*(m-1)*1.05;
 
hemn=q*COI(i,j)+WEIGHT(i,j)*(m-1)*1.05+p*COI(m,n)+WEIGHT(m,n)*(i-1)*1.05;
 
…………
 
end
 
end
 
end
 
k=(i-1)*20+j;
 
Q(k)=ones(1,10)*(t.*COI)*ones(20,1);
 
S(k)=(0:9)*(s2.*WEIGHT)*ones(20,1);
 
H(k)=ones(1,10)*(t.*COI)*ones(20,1)+(0:9)*WEIGHT*ones(20,1)
 
%计算目标函数最终解的和
 
end
 
end
 
plot(Q);
 
plot(S);
 
plot(H);
 
%图表输出
 
3.2 结果分析
 
计算过程中目标函数H值的变化
 
从图中可以看出：通过函数的迭代，目标函数总值H随着迭代次数的增加而曲线下降，在次数180到200此之间，函数值基本上没有太大的变化，可以看出迭代到了这个之间已经基本达到收敛；目标函数总值H开始迭代时为1.2554*106，迭代完成后H值为8.2887*105；迭代开始时目标函数Q值为1.1798*106，迭代结束值变为7.7184*105。迭代开始时目标函数S值为7.271.3*105，迭代结束后变为5.312*105。验证H=Q+S,迭代开始时 1.2554*106=1.1798*106+7.271.3*105，迭代结束后8.2887*105=7.7184*105+5.312*105。函数最终达到很好的收敛。
 
4 总结
 
本文是在一个静态的情况下对货品进行分配货位的求解；利用权重赋值的方法把多目标问题转化为求解一个目标的情况进行求解；使用matlab软件进行迭代求解，并得到一个很好的收敛优化解。
 
在实际使用过程中存在很多不确定因素，包括权重的取值、货架的多层次、堆垛机的实际速度等等都是需要考虑的问题，本文求解仅作一个探讨。继续深入研究应使用动态仿真方法来实现。
 
参 考 文 献
 
[1]柳赛男，柯映林，李江雄，等.基于调度策略的自动化仓库系统优化问题研究.计算机集成制造系统，2006
 
[2]贾煜亮.自动化立体仓库中货位实时分配优化问题研究. 清华大学工学硕士学位论文.2007
 
[3]邵琳.自动化立体仓库堆垛机单一作业循环时间简析.China storage & transport magazine.2007
 
[4]陈玉金.立体仓库堆垛机设计.现代物流.2006