当前位置:文档之家› CT图像三维重建(附源码)

CT图像三维重建(附源码)

CT图像三维重建(附源码)
CT图像三维重建(附源码)

程序流图:

MATLAB 源码: clc;

clear all;

close all;

% load mri %载入mri 数据,是matlab 自带库

% ph = squeeze(D); %压缩载入的数据D ,并赋值给ph ph = phantom3d(128);

prompt={'Enter the Piece num(1 to 128):'}; %提示信息“输入1到27的片的数字”

name='Input number'; %弹出框名称

defaultanswer={'1'}; %默认数字

numInput=inputdlg(prompt,name,1,defaultanswer) %弹出框,并得到用户的输入信息 P= squeeze(ph(:,:,str2num(cell2mat(numInput))));%将用户的输入信息转换成数字,并从ph 中得到相应的片信息P

imshow(P) %展示图片P

D = 250; %将D 赋值为250,是从扇束顶点到旋转中心的像素距离。 dsensor1 = 2; %正实数指定扇束传感器的间距2

F1 = fanbeam(P,D,'FanSensorSpacing',dsensor1); %通过P ,D 等计算扇束的数据值 dsensor2 = 1; %正实数指定扇束传感器的间距1

F2 = fanbeam(P,D,'FanSensorSpacing',dsensor2); %通过P ,D 等计算扇束的数据值

dsensor3 = 0.25 %正实数指定扇束传感器的间距0.25 [F3, sensor_pos3, fan_rot_angles3] = fanbeam(P,D,...

'FanSensorSpacing',dsensor3); %通过P ,D 等计算扇束的数据值,并得到扇束传感器的位置sensor_pos3和旋转角度fan_rot_angles3

figure, %创建窗口

imagesc(fan_rot_angles3, sensor_pos3, F3) %根据计算出的位置和角度展示F3的图片

colormap(hot); %设置色图为hot colorbar; %显示色栏

xlabel('Fan Rotation Angle (degrees)') %定义x 坐标轴

ylabel('Fan Sensor Position (degrees)') %定义y 坐标轴

output_size = max(size(P)); %得到P 维数的最大值,并赋值给output_size Ifan1 = ifanbeam(F1,D, ... 'FanSensorSpacing',dsensor1,'OutputSize',output_size);

%根据扇束投影数据F1及D 等数据重建图像

figure, imshow(Ifan1) %创建窗口,并展示图片Ifan1

title('图一');

disp('图一和原图的性噪比为:');

result=psnr1(Ifan1,P);

Ifan2 = ifanbeam(F2,D, ...

'FanSensorSpacing',dsensor2,'OutputSize',output_size);

%根据扇束投影数据F2及D 等数据重建图像

figure, imshow(Ifan2) %创建窗口,并展示图片Ifan2

disp('图二和原图的性噪比为:');

result=psnr1(Ifan2,P);

title('图二');

Ifan3 = ifanbeam(F3,D, ...

生成128的输入图片数字对图片信息进行预处

用函数fanbeam 进行映射,得到扇束的数据,并用函数ifanbeam 根据扇束投影数据重建图像,并计算重建图像和原图的

结束

'FanSensorSpacing',dsensor3,'OutputSize',output_size);

%根据扇束投影数据F3及D等数据重建图像figure, imshow(Ifan3) %创建窗口,并展示图片Ifan3

title('图三');

disp('图三和原图的性噪比为:');

result=psnr1(Ifan3,P);

function [p,ellipse]=phantom3d(varargin)

%PHANTOM3D Three-dimensional analogue of MATLAB Shepp-Logan phantom

% P = PHANTOM3D(DEF,N) generates a 3D head phantom that can

% be used to test 3-D reconstruction algorithms.

%

% DEF is a string that specifies the type of head phantom to generate.

% Valid values are:

%

% 'Shepp-Logan' A test image used widely by researchers in

% tomography

% 'Modified Shepp-Logan' (default) A variant of the Shepp-Logan phantom

% in which the contrast is improved for better

% visual perception.

%

% N is a scalar that specifies the grid size of P.

% If you omit the argument, N defaults to 64.

%

% P = PHANTOM3D(E,N) generates a user-defined phantom, where each row

% of the matrix E specifies an ellipsoid in the image. E has ten columns,

% with each column containing a different parameter for the ellipsoids:

%

% Column 1: A the additive intensity value of the ellipsoid

% Column 2: a the length of the x semi-axis of the ellipsoid

% Column 3: b the length of the y semi-axis of the ellipsoid

% Column 4: c the length of the z semi-axis of the ellipsoid

% Column 5: x0 the x-coordinate of the center of the ellipsoid

% Column 6: y0 the y-coordinate of the center of the ellipsoid

% Column 7: z0 the z-coordinate of the center of the ellipsoid

% Column 8: phi phi Euler angle (in degrees) (rotation about z-axis)

% Column 9: theta theta Euler angle (in degrees) (rotation about x-axis) % Column 10: psi psi Euler angle (in degrees) (rotation about z-axis)

%

% For purposes of generating the phantom, the domains for the x-, y-, and % z-axes span [-1,1]. Columns 2 through 7 must be specified in terms

% of this range.

%

% [P,E] = PHANTOM3D(...) returns the matrix E used to generate the phantom. %

% Class Support

% -------------

% All inputs must be of class double. All outputs are of class double.

%

% Remarks

% -------

% For any given voxel in the output image, the voxel's value is equal to the

% sum of the additive intensity values of all ellipsoids that the voxel is a

% part of. If a voxel is not part of any ellipsoid, its value is 0.

%

% The additive intensity value A for an ellipsoid can be positive or negative;

% if it is negative, the ellipsoid will be darker than the surrounding pixels.

% Note that, depending on the values of A, some voxels may have values outside % the range [0,1].

%

% Example

% -------

% ph = phantom3d(128);

% figure, imshow(squeeze(ph(64,:,:)))

%

% Copyright 2005 Matthias Christian Schabel (matthias @ stanfordalumni . org) % University of Utah Department of Radiology

% Utah Center for Advanced Imaging Research

% 729 Arapeen Drive

% Salt Lake City, UT 84108-1218

%

% This code is released under the Gnu Public License (GPL). For more information, % see : https://www.doczj.com/doc/1717250795.html,/copyleft/gpl.html

%

% Portions of this code are based on phantom.m, copyrighted by the Mathworks %

[ellipse,n] = parse_inputs(varargin{:});

p = zeros([n n n]);

rng = ( (0:n-1)-(n-1)/2 ) / ((n-1)/2);

[x,y,z] = meshgrid(rng,rng,rng);

coord = [flatten(x); flatten(y); flatten(z)];

p = flatten(p);

for k = 1:size(ellipse,1)

A = ellipse(k,1); % Amplitude change for this ellipsoid

asq = ellipse(k,2)^2; % a^2

bsq = ellipse(k,3)^2; % b^2

csq = ellipse(k,4)^2; % c^2

x0 = ellipse(k,5); % x offset

y0 = ellipse(k,6); % y offset

z0 = ellipse(k,7); % z offset

phi = ellipse(k,8)*pi/180; % first Euler angle in radians

theta = ellipse(k,9)*pi/180; % second Euler angle in radians

psi = ellipse(k,10)*pi/180; % third Euler angle in radians

cphi = cos(phi);

sphi = sin(phi);

ctheta = cos(theta);

stheta = sin(theta);

cpsi = cos(psi);

spsi = sin(psi);

% Euler rotation matrix

alpha = [cpsi*cphi-ctheta*sphi*spsi cpsi*sphi+ctheta*cphi*spsi spsi*stheta;

-spsi*cphi-ctheta*sphi*cpsi -spsi*sphi+ctheta*cphi*cpsi cpsi*stheta;

stheta*sphi -stheta*cphi ctheta];

% rotated ellipsoid coordinates

coordp = alpha*coord;

idx = find((coordp(1,:)-x0).^2./asq + (coordp(2,:)-y0).^2./bsq + (coordp(3,:)-z0).^2./csq <= 1);

p(idx) = p(idx) + A;

end

p = reshape(p,[n n n]);

return;

function out = flatten(in)

out = reshape(in,[1 prod(size(in))]);

return;

function [e,n] = parse_inputs(varargin)

% e is the m-by-10 array which defines ellipsoids

% n is the size of the phantom brain image

n = 128; % The default size

e = [];

defaults = {'shepp-logan', 'modified shepp-logan', 'yu-ye-wang'};

for i=1:nargin

if ischar(varargin{i}) % Look for a default phantom

def = lower(varargin{i});

idx = strmatch(def, defaults);

if isempty(idx)

eid = sprintf('Images:%s:unknownPhantom',mfilename);

msg = 'Unknown default phantom selected.';

error(eid,'%s',msg);

end

switch defaults{idx}

case 'shepp-logan'

e = shepp_logan;

case 'modified shepp-logan'

e = modified_shepp_logan;

case 'yu-ye-wang'

e = yu_ye_wang;

end

elseif numel(varargin{i})==1

n = varargin{i}; % a scalar is the image size elseif ndims(varargin{i})==2 && size(varargin{i},2)==10

e = varargin{i}; % user specified phantom

else

eid = sprintf('Images:%s:invalidInputArgs',mfilename);

msg = 'Invalid input arguments.';

error(eid,'%s',msg);

end

end

% ellipse is not yet defined

if isempty(e)

e = modified_shepp_logan;

end

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Default head phantoms: % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function e = shepp_logan

e = modified_shepp_logan;

e(:,1) = [1 -.98 -.02 -.02 .01 .01 .01 .01 .01 .01];

return;

function e = modified_shepp_logan

%

% This head phantom is the same as the Shepp-Logan except

% the intensities are changed to yield higher contrast in

% the image. Taken from Toft, 199-200.

%

% A a b c x0 y0 z0 phi theta psi

% -----------------------------------------------------------------

e = [ 1 .6900 .920 .810 0 0 0 0 0 0

-.8 .6624 .874 .780 0 -.0184 0 0 0 0

-.2 .1100 .310 .220 .22 0 0 -18 0 10

-.2 .1600 .410 .280 -.22 0 0 18 0 10

.1 .2100 .250 .410 0 .35 -.15 0 0 0

.1 .0460 .046 .050 0 .1 .25 0 0 0

.1 .0460 .046 .050 0 -.1 .25 0 0 0

.1 .0460 .023 .050 -.08 -.605 0 0 0 0

.1 .0230 .023 .020 0 -.606 0 0 0 0

.1 .0230 .046 .020 .06 -.605 0 0 0 0 ];

return;

function e = yu_ye_wang

%

% Yu H, Ye Y, Wang G, Katsevich-Type Algorithms for Variable Radius Spiral Cone-Beam CT %

% A a b c x0 y0 z0 phi theta psi

% -----------------------------------------------------------------

e = [ 1 .6900 .920 .900 0 0 0 0 0 0

-.8 .6624 .874 .880 0 0 0 0 0 0

-.2 .4100 .160 .210 -.22 0 -.25 108 0 0

-.2 .3100 .110 .220 .22 0 -.25 72 0 0

.2 .2100 .250 .500 0 .35 -.25 0 0 0

.2 .0460 .046 .046 0 .1 -.25 0 0 0

.1 .0460 .023 .020 -.08 -.65 -.25 0 0 0

.1 .0460 .023 .020 .06 -.65 -.25 90 0 0

.2 .0560 .040 .100 .06 -.105 .625 90 0 0

-.2 .0560 .056 .100 0 .100 .625 0 0 0 ]; return;

% func——计算两幅图像的psnr值

function result=psnr1(in1,in2)

z=mse(in1,in2);

result=10*log10(255.^2/z)

% plot(result)

function z=mse(x,y)

if ndims(x)==3

x=rgb2gray(x);

end

if ndims(y)==3

y=rgb2gray(y);

end

x=double(x);

y=double(y);

[m1,n1]=size(x);

[m2,n2]=size(y);

m=min(m1,m2);

n=min(n1,n2);

z=0;

for i=1:m

for j=1:n

z=z+(x(i,j)-y(i,j)).^2;

end

end

z=z/(m*n);

图像三维重建技术

1概述 随着计算机软硬件技术的快速发展,大规模复杂场景的实时绘制已经成为可能,这也加快了虚拟现实技术的发展,又对模型的复杂度和真实感提出了新的要求。虚拟场景是虚拟现实系统的重要组成部分,它的逼真度将直接影响整个虚拟现实系统的沉浸感。客观世界在空间上是三维的,而现有的图像采集装置所获取的图像是二维的。尽管图像中含有某些形式的三维空间信息,但要真正在计算机中使用这些信息进行进一步的应用处理,就必须采用三维重建技术从二维图像中合理地提取并表达这些 三维信息。 三维建模工具虽然日益改进,但构建稍显复杂的三维模型依旧是一件非常耗时费力的工作。而很多要构建的三维模型都存在于现实世界中,因此三维扫描技术和基于图像建模技术就成了人们心目中理想的建模方式;又由于前者一般只能获取景物的几何信息,而后者为生成具有照片级真实感的合成图像提供了一种自然的方式,因此它迅速成为目前计算机图形学领域中的研究热点。 2三维建模技术 三维重建技术能够从二维图像出发构造具有真实感的三维图形,为进一步的场景变化和组合运算奠定基础,从而促进图像和三维图形技术在航天、造船、司法、考古、 工业测量、 电子商务等领域的深入广泛的应用。3基于图像的三维重建技术 基于图像的建模最近几年兴起的一门新技术,它使用直接拍摄到的图像,采用尽量少的交互操作,重建场 景。 它克服了传统的基于几何的建模技术的许多不足,有无比的优越性。传统的三维建模工具虽然日益改进,但构建稍显复杂的三维模型依旧是一件非常耗时费力的工作。考虑到我们要构建的很多三维模型都能在现实世界中找到或加以塑造,因此三维扫描技术和基于图像建模技术就成了人们心目中理想的建模方式;又由于前者一般只能获取景物的几何信息,而后者为生成具有照片级真实感的合成图像提供了一种自然的方式,因此它迅速成为目前计算机图形学领域中的研究热点。 4 基于图像重建几何模型的方法 4.1 基于侧影轮廓线重建几何模型 物体在图像上的侧影轮廓线是理解物体几何形状的 一条重要线索1当以透视投影的方式从多个视角观察某一空间物体时,在每个视角的画面上都会得到一条该物体的侧影轮廓线,这条侧影轮廓线和对应的透视投影中心共同确定了三维空间中一个一般形状的锥体1显然,该物体必将位于这个锥体之内;而所有这些空间锥体的交则构成了一个包含该物体的空间包络1这个空间包络被称为物体的可见外壳,当观察视角足够多时,可见外壳就可以被认为是该物体的一个合理的逼近。鉴于此类算法一般需要大量的多视角图像,因此图像的定标工作就变得非常复杂。 4.2采用立体视觉方法重建几何模型 基于立体视觉重建三维几何是计算机视觉领域中的经典问题,被广泛应用于自动导航装置。近年来,立体视觉 图像三维重建技术 康皓,王明倩,王莹莹 (装甲兵技术学院电子工程系,吉林长春130117) 摘要:基于图像的三维重建属于计算机视觉中的一个重要的研究方向,从提出到现在已有十多年的历史。文章首先对三维重建技术做了详细阐述,并着重从计算机图形学的研究角度对基于图像建模技术进行了综述,介绍了 具有代表性的基于图像建模的方法及其最新研究进展,给出了这些方法的基本原理, 并对这些方法进行分析比较,最后对基于图像建模技术的未来研究给出了一些建议和应解决的问题。关键词:三维建模技术;图像建模技术;计算机图形学;虚拟现实中图分类号:TP271文献标识码:A 文章编号1006-8937(2009)11-0042-02 Three-dimensional image reconstruction technique KANG Hao,WANG Ming-qian,WANG Ying-ying (DepartmentofElectronicEngineering,ArmoredInstituteofTechnology,Changchun,Jilin130117,China) Abstract:Image-based Three-dimensional reconstruction is an important research direction in computer vision ,from now more than ten years'history.This article first describes three-dimensional reconstruction technique in detail and review image-based modeling techniques from the perspective of computer graphics research,introduce a representative of the method of image-based modeling and the latest research progress,give the basic principles of these methods,analysis and compare these methods,finally,give a number of recommendations and problems which should be solved on image-based modeling technology for future research. Keywords:three-dimensional modeling techniques;image modeling techniques;computer graphics;virtual reality 收稿日期:2009-03-19 作者简介:康皓(1978-),女,吉林长春人,硕士研究生,讲师,研 究方向:计算机辅助设计与编程。 TECHNOLOGICAL DEVELOPMENT OF ENTERPRISE 2009年6月Jun.2009 企业技术开发 第28卷

matlab相关图形实现代码

根据数据点绘制饼图和针状图: x=[1 2 3 4 5 6]; >> subplot(2,2,1);pie(x); >> subplot(2,2,2);pie3(x); >> subplot(2,2,3);stem(x); >>subplot(2,2,4);stem3(x); 5% 10% 14% 19% 24% 29% 24% 29% 19% 5%14% 10%0 2 4 6 2 4 6 5 10 01 2 05 10

根据数据点绘制向量场图、羽状图和罗盘图: x=[1 2 3 4 5 6];y=[1 2 3 4 5 6]; u=[1 2 3 4 5 6];v=[1 2 3 4 5 6]; subplot(2,2,1);quiver(x,y,u,v); subplot(2,2,2);quiver(x,y,u,v,'r'); subplot(2,2,3);feather(u,v); subplot(2,2,4);compass(u,v); 024680 246 802468 246 80 5 10 15 2 4 6 5 10 30 210 60240 90270 120 300 150330 180

rand(m,n)产生m ×n 均匀分布的随机矩阵,元素取值在0.0~1.0。 randn 函数:产生标准正态分布的随机数或矩阵的函数。 Y = randn(m,n) 或 Y = randn([m n])返回一个m*n 的随机项矩阵。 > theta=10*rand(1,50); %确定50个随机数theta >> Z=peaks; %确定Z 为峰值函数peaks >> x=0:0.01:2*pi;y=sin(x); %确定正弦函数数据点x.y >> t=randn(1000,1); %确定1000个随机数t >> subplot(2,2,1);rose(theta); %关于(theta )的玫瑰花图 >> subplot(2,2,2);area(x,y); %关于(x,y)的面积图 >> subplot(2,2,3);contour(Z); %关于Z 的等值线图(未填充) >> subplot(2,2,4);hist(t); %关于t 的柱状图 5 10 30 210 60 240 90270 120300150330 18000246 -1 -0.500.5 110 20 30 40 10 2030 40-4 -2 2 4 100 200 300

数字图像处理_旋转与幅度谱(含MATLAB代码)

数字图像处理实验一 15生医 一、实验内容 产生右图所示图像 f1(m,n),其中图像大小为256 ×256,中间亮条为128×32,暗处=0,亮处=100。 对其进行FFT: ①同屏显示原图f1(m,n)和FFT(f1)的幅度谱图; ②若令f2(m,n)=(-1)^(m+n)f1(m,n),重复 以上过程,比较二者幅度谱的异同,简述理由; ③若将f2(m,n)顺时针旋转90度得到f3(m,n),试显示FFT(f3)的 幅度谱,并与FFT(f2)的幅度谱进行比较; ④若将f1(m,n) 顺时针旋转90度得到f4(m,n),令f5(m,n) = f1(m,n) + f4(m,n),试显示FFT(f5)的幅度谱,指出其与 FFT(f1)和FFT(f4)的关系; ⑤若令f6(m,n)=f2(m,n)+f3(m,n),试显示FFT(f6)的幅度谱,并指出其与 FFT(f2)和FFT(f3)的关系,比较FFT(f6)和FFT(f5)的幅度谱。 二、运行环境 MATLAB R2014a 三、运行结果及分析 1.同屏显示原图f1(m,n)和FFT(f1)的幅度谱图:

50100150200250 100150200250 50100150200250 100150200250 2.令f2(m,n)=(-1)^(m+n )f1(m,n),对其进行FFT ,比较f2与f1幅度谱的异同,简述理由: 50100150200250 100150200250 50100150200250 100150200250 异同及理由:①空域:f2由于前边乘了系数(-1)^(m+n ),导致灰度值有正有负,而在MATLAB 的imshow 函数中默认把负值变为0(有些情况是取反),所以形成了如左图所示的黑白花纹。②频域:FFT(2)

三维重建综述

三维重建综述 三维重建方法大致分为两个部分1、基于结构光的(如杨宇师兄做的)2、基于图片的。这里主要对基于图片的三维重建的发展做一下总结。 基于图片的三维重建方法: 基于图片的三维重建方法又分为双目立体视觉;单目立体视觉。 A双目立体视觉: 这种方法使用两台摄像机从两个(通常是左右平行对齐的,也可以是上下竖直对齐的)视点观测同一物体,获取在物体不同视角下的感知图像,通过三角测量的方法将匹配点的视差信息转换为深度,一般的双目视觉方法都是利用对极几何将问题变换到欧式几何条件下,然后再使用三角测量的方法估计深度信息这种方法可以大致分为图像获取、摄像机标定、特征提取与匹配、摄像机校正、立体匹配和三维建模六个步骤。王涛的毕业论文就是做的这方面的工作。双目立体视觉法的优点是方法成熟,能够稳定地获得较好的重建效果,实际应用情况优于其他基于视觉的三维重建方法,也逐渐出现在一部分商业化产品上;不足的是运算量仍然偏大,而且在基线距离较大的情况下重建效果明显降低。 代表文章:AKIMOIO T Automatic creation of3D facial models1993 CHEN C L Visual binocular vison systems to solid model reconstruction 2007 B基于单目视觉的三维重建方法: 单目视觉方法是指使用一台摄像机进行三维重建的方法所使用的图像可以是单视点的单幅或多幅图像,也可以是多视点的多幅图像前者主要通过图像的二维特征推导出深度信息,这些二维特征包括明暗度、纹理、焦点、轮廓等,因此也被统称为恢复形状法(shape from X) 1、明暗度(shape from shading SFS) 通过分析图像中的明暗度信息,运用反射光照模型,恢复出物体表面法向量信息进行三维重建。SFS方法还要基于三个假设a、反射模型为朗伯特模型,即从各个角度观察,同一点的明暗度都相同的;b、光源为无限远处点光源;c、成像关系为正交投影。 提出:Horn shape from shading:a method for obtaining the shape of a smooth opaque object from one view1970(该篇文章被引用了376次) 发展:Vogel2008年提出了非朗伯特的SFS模型。 优势:可以从单幅图片中恢复出较精确的三维模型。 缺点:重建单纯依赖数学运算,由于对光照条件要求比较苛刻,需要精确知道光源的位置及方向等信息,使得明暗度法很难应用在室外场景等光线情况复杂的三维重建上。 2、光度立体视觉(photometric stereo) 该方法通过多个不共线的光源获得物体的多幅图像,再将不同图像的亮度方程联立,求解出物体表面法向量的方向,最终实现物体形状的恢复。 提出:Woodham对SFS进行改进(1980年):photometric method for determining surface orientation from multiple images(该文章被引用了891次) 发展:Noakes:非线性与噪声减除2003年; Horocitz:梯度场合控制点2004年; Tang:可信度传递与马尔科夫随机场2005年; Basri:光源条件未知情况下的三维重建2007年; Sun:非朗伯特2007年; Hernandez:彩色光线进行重建方法2007年;

matlab图像处理代码

附录 MATLAB图像处理命令  1.applylut  功能: 在二进制图像中利用lookup表进行边沿操作。 语法: A = applylut(BW,lut) 举例 lut = makelut('sum(x(:)) == 4',2); BW1 = imread('text.tif'); BW2 = applylut(BW1,lut); imshow(BW1) figure, imshow(BW2) 相关命令: makelut 2.bestblk  功能: 确定进行块操作的块大小。 语法: siz = bestblk([m n],k) [mb,nb] = bestblk([m n],k) 举例 siz = bestblk([640 800],72) siz = 64 50 相关命令: blkproc 3.blkproc  功能:

MATLAB 高级应用——图形及影像处理 320 实现图像的显式块操作。 语法: B = blkproc(A,[m n],fun) B = blkproc(A,[m n],fun,P1,P2,...) B = blkproc(A,[m n],[mborder nborder],fun,...) B = blkproc(A,'indexed',...) 举例 I = imread('alumgrns.tif'); I2 = blkproc(I,[8 8],'std2(x)*ones(size(x))'); imshow(I) figure, imshow(I2,[]); 相关命令: colfilt, nlfilter,inline 4.brighten  功能: 增加或降低颜色映像表的亮度。 语法: brighten(beta) newmap = brighten(beta) newmap = brighten(map,beta) brighten(fig,beta) 相关命令: imadjust, rgbplot 5.bwarea  功能: 计算二进制图像对象的面积。 语法: total = bwarea(BW) 举例 BW = imread('circles.tif'); imshow(BW);

SAR三维立体重建实验报告要点

SAR立体三维重建 姓名: ******* 学号: ********* 班级: ************* 指导教师: ******

1实验目的 1、理解基于合成孔径雷达立体像对的灰度信息进行三维重建的基本原理与方法; 2、了解ERDAS IMAGINE的基本功能,熟练掌握StereoSAR模块的使用方法; 3、理解SAR传感器几何模型及基于地面控制点(Ground Control Points, GCPs)几何模型精化的原理与方法; 4、通过真实SAR像对的数据处理,掌握SAR立体三维重建的基本流程。 2实验数据说明 本实验采用ERDAS IMAGINE软件的示例数据,RADASAT影像StereoSAR_Ref.img和StereoSAR_Match.img,这两景影像分别拍摄于1996年9月24日和1996年9月17日。 3实验原理 经过试验九的操作,使我们对InSAR提取测区DEM有了一定的掌握。而摄影测量中我们也学习了基于立体像对制作测区三维景观图,因此在此次实验中我们利用摄影测量的原理基于SAR影像进行三维重建。 3.1 SAR立体图像的获取 立体图像在摄影测量中称为立体相对。所谓立体相对是由不同摄站摄取的具有一定重叠的两张相片。因此雷达立体图像也可以定义为:由天线位置探测获取的具有一定影像重叠的两幅雷达图像[1]。 雷达立体图像的获取方式有两种:同侧立体观测和异侧立体观测。前者是指飞行器沿着不同的航线飞行(两次飞行方向可以相同或者相反),雷达从地物的

一侧对同一地区成像,同侧立体观测有可分为同一高度和不同高度两类;异侧立体观测是指雷达从地物的两侧分别对同一地区成像。 图 3.1-1 雷达立体图像获取方式 异侧立体观测获取的雷达立体图像视差明显,基高比(摄影基线与航高之比)大,有利于提高地物目标点高程的测量精度。但是地形起伏较大的地区,目标地物在立体像对的两幅图像上的相应影像不仅颜色差异很大,而且由于地形起伏引起的几何变形差异也很大。因为高出地面地物目标的左侧向着左航线的雷达天线,有效反射面积大,影像为浅色调;而对于右航线,该地物目标左侧反射信号弱或者为盲区,所以影像为深色调。同理,地物目标右侧在两幅图像上的色调与地物目标左侧色调刚好相反,由于异侧飞行所获取的雷达立体图像,使一幅图像上的阴影位于地物目标的一侧,而另一幅图像上的阴影在地物目标的另一侧,这就给立体观察与测量带来了极大的困难。另外,采用斜距显示的雷达图像由于地形起伏的影响,同一地物目标在立体像对的两幅图像上的变形差异也很大。因此,异侧(对侧)获取的雷达立体图像,只适于平坦或丘陵而不适合山地的立体观察与测量[2]。 同侧同高度或者不同高度获取的雷达立体图像,视差和基高比虽然比异侧获取的雷达立体图像要小,但两幅图像上相应影像的色调和图像变形差异较小,只要对雷达工作参数进行适当选择,还是能获得较好的立体效应的,故在丘陵地、山地一般都采用同侧获取雷达立体图像进行地物目标的三维定位和立体测图[2]。 3.2 SAR立体图像的视差 SAR立体图像的视差指高出某一基准面的地物目标在两幅图像上的位移差。摄影测量中称之为立体像对的左右视差,是地物目标点高差的反应,由左右视差

数字图像处理matlab代码

一、编写程序完成不同滤波器的图像频域降噪和边缘增强的算法并进行比较,得出结论。 1、不同滤波器的频域降噪 1.1 理想低通滤波器(ILPF) I1=imread('eight.tif'); %读取图像 I2=im2double(I1); I3=imnoise(I2,'gaussian',0.01); I4=imnoise(I3,'salt & pepper',0.01); figure,subplot(1,3,1); imshow(I2) %显示灰度图像 title('原始图像'); %为图像添加标题 subplot(1,3,2); imshow(I4) %加入混合躁声后显示图像 title('加噪后的图像'); s=fftshift(fft2(I4)); %将灰度图像的二维不连续Fourier 变换的零频率成分 移到频谱的中心 [M,N]=size(s); %分别返回s的行数到M中,列数到N中n1=floor(M/2); %对M/2进行取整 n2=floor(N/2); %对N/2进行取整 d0=40; %初始化d0 for i=1:M for j=1:N d=sqrt((i-n1)^2+(j-n2)^2); %点(i,j)到傅立叶变换中心的距离 if d<=d0 %点(i,j)在通带内的情况 h=1; %通带变换函数 else %点(i,j)在阻带内的情况 h=0; %阻带变换函数 end s(i,j)=h*s(i,j); %ILPF滤波后的频域表示

end end s=ifftshift(s); %对s进行反FFT移动 s=im2uint8(real(ifft2(s))); %对s进行二维反离散的Fourier变换后,取复 数的实部转化为无符号8位整数 subplot(1,3,3); %创建图形图像对象 imshow(s); %显示ILPF滤波后的图像 title('ILPF滤波后的图像(d=40)'); 运行结果: 1.2 二阶巴特沃斯低通滤波器(BLPF) I1=imread('eight.tif'); %读取图像 I2=im2double(I1); I3=imnoise(I2,'gaussian',0.01); I4=imnoise(I3,'salt & pepper',0.01); figure,subplot(1,3,1); imshow(I2) %显示灰度图像 title('原始图像'); %为图像添加标题 subplot(1,3,2); imshow(I4) %加入混合躁声后显示图像 title('加噪后的图像'); s=fftshift(fft2(I4));%将灰度图像的二维不连续Fourier 变换的零频率成分 移到频谱的中心 [M,N]=size(s); %分别返回s的行数到M中,列数到N中n=2; %对n赋初值

三维血管的重建

血管的三维重建 摘要 对于血管的三维重建,本文研究了血管这一类特殊管道的中轴线及其半径的算法,绘制中轴线在XY 、YZ 、ZX 平面的投影图这些问题,问题分为三部分。 针对第一部分,先将100张切片图片在MATLAB 中导出生成0-1矩阵数据,在计算100张切片的最大内切圆半径及对应圆心坐标,为减小误差求100张切片最大内切圆的平均半径41666.29 d 。中轴线的曲线方程可在MATLAB 中拟合得到。 针对第二部分,得到中轴线曲线方程在MATLAB 中绘制出中轴线方程的空间曲线,之后将其投影在XY 、YZ 、ZX 平面上。 针对第三部分,对100张切片进行叠加重合,得到血管的三维立体图,再通过MATLAB 对血管的三维立体图进行优化完成血管的三维重建。 关键词:MATLAB 软件 管道半径中轴线曲线方程

一、问题重述 1.1基本情况 断面可用于了解生物组织、器官等的形态。如果用切片机连续不断地将样本切成数十、成百的平行切片,可依次逐片观察。根据拍照并采样得到的平行切片数字图象,运用计算机可重建组织、器官等准确的三维形态。 1.2相关信息 假设某些血管可视为一类特殊的管道,该管道的表面是由球心沿着某一曲线(称为中轴线)的球滚动包络而成。 现有某管道的相继100张平行切片图象,记录了管道与切片的交。图象文件名依次为0.bmp、1.bmp、…、99.bmp,格式均为BMP,宽、高均为512个象素(pixel)。 取坐标系的Z轴垂直于切片,第1张切片为平面Z=0,第100张切片为平面Z=99。Z=z切片图象中象素的坐标依它们在文件中出现的前后次序为(-256,-256,z),(-256,-255,z),…(-256,255,z), (-255,-256,z),(-255,-255,z),…(-255,255,z), …… (255,-256,z),(255,-255,z),…(255,255,z)。 1.3提出的问题 问题一:计算出管道的中轴线与半径,给出具体的算法。 问题二:绘制中轴线在XY、YZ、ZX平面的投影图。 问题三:绘制血管的三维重建立体图。

数字图像处理MATLAB相关代码

1.图像反转 MATLAB程序实现如下: I=imread('xian.bmp'); J=double(I); J=-J+(256-1); %图像反转线性变换 H=uint8(J); subplot(1,2,1),imshow(I); subplot(1,2,2),imshow(H); 2.灰度线性变换 MATLAB程序实现如下: I=imread('xian.bmp'); subplot(2,2,1),imshow(I); title('原始图像'); axis([50,250,50,200]); axis on; %显示坐标系 I1=rgb2gray(I); subplot(2,2,2),imshow(I1); title('灰度图像'); axis([50,250,50,200]); axis on; %显示坐标系 J=imadjust(I1,[0.1 0.5],[]); %局部拉伸,把[0.1 0.5]内的灰度拉伸为[0 1] subplot(2,2,3),imshow(J); title('线性变换图像[0.1 0.5]'); axis([50,250,50,200]); grid on; %显示网格线 axis on; %显示坐标系 K=imadjust(I1,[0.3 0.7],[]); %局部拉伸,把[0.3 0.7]内的灰度拉伸为[0 1] subplot(2,2,4),imshow(K); title('线性变换图像[0.3 0.7]'); axis([50,250,50,200]); grid on; %显示网格线 axis on; %显示坐标系 3.非线性变换 MATLAB程序实现如下: I=imread('xian.bmp'); I1=rgb2gray(I); subplot(1,2,1),imshow(I1); title('灰度图像'); axis([50,250,50,200]); grid on; %显示网格线 axis on; %显示坐标系 J=double(I1); J=40*(log(J+1)); H=uint8(J);

图像处理实例(含Matlab代码)

信号与系统实验报告——图像处理 学院:信息科学与工程学院 专业:2014级通信工程 组长:** 组员:** 2017.01.02

目录 目录 (2) 实验一图像一的细胞计数 (3) 一、实验内容及步骤 (3) 二、Matlab程序代码 (3) 三、数据及结果 (4) 实验二图像二的图形结构提取 (5) 一、实验内容及步骤 (5) 二、Matlab程序代码 (5) 三、数据及结果 (6) 实验三图像三的图形结构提取 (7) 一、实验内容及步骤 (7) 二、Matlab程序代码 (7) 三、数据及结果 (8) 实验四图像四的傅里叶变化及巴特沃斯低通滤波 (9) 一、实验内容及步骤 (9) 二、Matlab程序代码 (9) 三、数据及结果 (10) 实验五图像五的空间域滤波与频域滤波 (11) 一、实验内容及步骤 (11) 二、Matlab程序代码 (11) 三、数据及结果 (12)

实验一图像一的细胞计数 一、实验内容及步骤 将该图形进行一系列处理,计算得到途中清晰可见细胞的个数。 首先,由于原图为RGB三色图像处理起来较为麻烦,所以转为灰度图,再进行二值化化为黑白图像,得到二值化图像之后进行中值滤波得到细胞分布的初步图像,为了方便计数对图像取反,这时进行一次计数,发现得到的个数远远多于实际个数,这时在进行一次中值滤波,去掉一些不清晰的像素点,剩下的应该为较为清晰的细胞个数,再次计数得到大致结果。 二、Matlab程序代码 clear;close all; Image = imread('1.jpg'); figure,imshow(Image),title('原图'); Image=rgb2gray(Image); figure,imshow(Image),title('灰度图'); Theshold = graythresh(Image); Image_BW = im2bw(Image,Theshold); Reverse_Image_BW22=~Image_BW; figure,imshow(Image_BW),title('二值化图像'); Image_BW_medfilt= medfilt2(Image_BW,[3 3]); figure,imshow(Image_BW_medfilt),title('中值滤波后的二值化图像'); Reverse_Image_BW = ~Image_BW_medfilt; figure,imshow(Reverse_Image_BW),title('图象取反'); Image_BW_medfilt2= medfilt2(Reverse_Image_BW,[20 20]); figure,imshow(Image_BW_medfilt2),title('第二次中值滤波的二值化图像'); [Label, Number]=bwlabel(Image_BW_medfilt,8);Number [Label, Number]=bwlabel(Image_BW_medfilt2,8);Number

三维重建方法综述

三维重建方法综述 三维重建方法大致分为两个部分1、基于结构光的2、基于图片的。这里主要对基于图片的三维重建的发展做一下总结。基于图片的三维重建方法: 基于图片的三维重建方法又分为双目立体视觉;单目立体视觉。 A双目立体视觉: 这种方法使用两台摄像机从两个(通常是左右平行对齐的,也可以是上下竖直对齐的)视点观测同一物体,获取在物体不同视角下的感知图像,通过三角测量的方法将匹配点的视差信息转换为深度,一般的双目视觉方法都是利用对极几何将问题变换到欧式几何条件下,然后再使用三角测量的方法估计深度信息这种方法可以大致分为图像获取、摄像机标定、特征提取与匹配、摄像机校正、立体匹配和三维建模六个步骤。王涛的毕业论文就是做的这方面的工作。双目立体视觉法的优点是方法成熟,能够稳定地获得较好的重建效果,实际应用情况优于其他基于视觉的三维重建方法,也逐渐出现在一部分商业化产品上;不足的是运算量仍然偏大,而且在基线距离较大的情况下重建效果明显降低。 代表文章:AKIMOIOT Automatic creation of 3D facial models 1993 CHENCL Visual binocular vison systems to solid model reconstruction 2007 B基于单目视觉的三维重建方法: 单目视觉方法是指使用一台摄像机进行三维重建的方法所使用的图像可以是单视点的单幅或多幅图像,也可以是多视点的多幅图像前者主要通过图像的二维特征推导出深度信息,这些二维特征包括明暗度、纹理、焦点、轮廓等,因此也被统称为恢复形状法(shape from X) 1、明暗度(shape from shading SFS) 通过分析图像中的明暗度信息,运用反射光照模型,恢复出物体表面法向量信息进行三维重建。SFS方法还要基于三个假设a、反射模型为朗伯特模型,即从各个角度观察,同一点的明暗度都相同的;b、光源为无限远处点光源;c、成像关系为正交投影。 提出:Horn shape from shading:a method for obtaining the shape of a smooth opaque object from one view 1970(该篇文章被引用了376次) 发展:V ogel2008年提出了非朗伯特的SFS模型。优势:可以从单幅图片中恢复出较精确的三维模型。 缺点:重建单纯依赖数学运算,由于对光照条件要求比较苛刻,需要精确知道光源的位置及方向等信息,使得明暗度法很难应用在室外场景等光线情况复杂的三维重建上。 2、光度立体视觉(photometric stereo) 该方法通过多个不共线的光源获得物体的多幅图像,再将不同图像的亮度方程联立,求解出物体表面法向量的方向,最终实现物体形状的恢复。 提出:Woodham对SFS进行改进(1980年):photometric method for determining surface orientation from multiple images(该文章被引用了891次) 发展:Noakes:非线性与噪声减除2003年; Horocitz:梯度场合控制点2004年; Tang:可信度传递与马尔科夫随机场2005年;Basri:光源条件未知情况下的三维重建2007年;Sun:非朗伯特2007年; Hernandez:彩色光线进行重建方法2007年; Shi:自标定的光度立体视觉法2010年。 3、纹理法(shape from texture SFT) 通过分析图像中物体表面重复纹理单元的大小形状,恢复出物体法向深度等信息,得到物体的三维几何模型。

基于MATLAB图像处理报告

基于M A T L A B图像处理报告一、设计题目 图片叠加。 二、设计要求 将一幅礼花图片和一幅夜景图片做叠加运算,使达到烟花夜景的美图效果。 三、设计方案 、设计思路 利用matlab强大的图像处理功能,通过编写程序,实现对两幅图片的像素进行线性运算,利用灰度变换的算法使图片达到预期的效果。 、软件介绍 MATLAB是matrix&laboratory两个词的组合,意为矩阵工厂(矩阵实验室)。是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中,为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案,并在很大程度上摆脱了传统非交互式程序设计语言(如C、Fortran)的编辑模式,代表了当今国际科学计算软件的先进水平。 MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等,主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。 MATLAB的基本数据单位是矩阵,它的指令表达式与数学、工程中常用的形式十分相似,故用MATLAB来解算问题要比用C,FORTRAN等语言完成相同的事情简捷得多,并且MATLAB 也吸收了像Maple等软件的优点,使MATLAB成为一个强大的数学软件。在新的版本中也加入了对C,FORTRAN,C++,JAVA的支持。可以直接调用,用户也可以将自己编写的实用程序导入到MATLAB函数库中方便自己以后调用,此外许多的MATLAB爱好者都编写了一些经典的程序,用户直接进行下载就可以用。

SAR三维立体重建实验报告

SAR立体三维重建 : ******* 学号: *********

班级: ************* 指导教师: ****** 1实验目的 1、理解基于合成孔径雷达立体像对的灰度信息进行三维重建的基本原理与方法; 2、了解ERDAS IMAGINE的基本功能,熟练掌握StereoSAR模块的使用方法; 3、理解SAR传感器几何模型及基于地面控制点(Ground Control Points, GCPs)几何模型精化的原理与方法; 4、通过真实SAR像对的数据处理,掌握SAR立体三维重建的基本流程。 2实验数据说明 本实验采用ERDAS IMAGINE软件的示例数据,RADASAT影像StereoSAR_Ref.img和StereoSAR_Match.img,这两景影像分别拍摄于1996年9月24日和1996年9月17日。 3实验原理 经过试验九的操作,使我们对InSAR提取测区DEM有了一定的掌握。而摄影测量中我们也学习了基于立体像对制作测区三维景观图,因此在此次实验中我们利用摄影测量的原理基于SAR影像进行三维重建。 3.1 SAR立体图像的获取 立体图像在摄影测量中称为立体相对。所谓立体相对是由不同摄站摄取的具有一定重叠的两相片。因此雷达立体图像也可以定义为:由天线位置探测获取的具有一定影像重叠的两幅雷达图像[1]。

雷达立体图像的获取方式有两种:同侧立体观测和异侧立体观测。前者是指飞行器沿着不同的航线飞行(两次飞行方向可以相同或者相反),雷达从地物的一侧对同一地区成像,同侧立体观测有可分为同一高度和不同高度两类;异侧立体观测是指雷达从地物的两侧分别对同一地区成像。 图 3.1-1 雷达立体图像获取方式 异侧立体观测获取的雷达立体图像视差明显,基高比(摄影基线与航高之比)大,有利于提高地物目标点高程的测量精度。但是地形起伏较大的地区,目标地物在立体像对的两幅图像上的相应影像不仅颜色差异很大,而且由于地形起伏引起的几何变形差异也很大。因为高出地面地物目标的左侧向着左航线的雷达天线,有效反射面积大,影像为浅色调;而对于右航线,该地物目标左侧反射信号弱或者为盲区,所以影像为深色调。同理,地物目标右侧在两幅图像上的色调与地物目标左侧色调刚好相反,由于异侧飞行所获取的雷达立体图像,使一幅图像上的阴影位于地物目标的一侧,而另一幅图像上的阴影在地物目标的另一侧,这就给立体观察与测量带来了极大的困难。另外,采用斜距显示的雷达图像由于地形起伏的影响,同一地物目标在立体像对的两幅图像上的变形差异也很大。因此,异侧(对侧)获取的雷达立体图像,只适于平坦或丘陵而不适合山地的立体观察与测量[2]。 同侧同高度或者不同高度获取的雷达立体图像,视差和基高比虽然比异侧获取的雷达立体图像要小,但两幅图像上相应影像的色调和图像变形差异较小,只要对雷达工作参数进行适当选择,还是能获得较好的立体效应的,故在丘陵地、山地一般都采用同侧获取雷达立体图像进行地物目标的三维定位和立体测图[2]。 3.2 SAR立体图像的视差

matlab数字图像处理源代码

数字图像去噪典型算法及matlab实现 希望得到大家的指点和帮助 图像去噪是数字图像处理中的重要环节和步骤。去噪效果的好坏直接影响 到后续的图像处理工作如图像分割、边缘检测等。图像信号在产生、传输过程中都可能会受到噪声的污染,一般数字图像系统中的常见噪声主要有:高斯噪声(主要由阻性元器件内部产生)、椒盐噪声(主要是图像切割引起的黑图像上的白点噪声或光电转换过程中产生的泊松噪声)等; 目前比较经典的图像去噪算法主要有以下三种: 均值滤波算法:也称线性滤波,主要思想为邻域平均法,即用几个像素灰度 的平均值来代替每个像素的灰度。有效抑制加性噪声,但容易引起图像模糊, 可以对其进行改进,主要避开对景物边缘的平滑处理。 中值滤波:基于排序统计理论的一种能有效抑制噪声的非线性平滑滤波信号处理技术。中值滤波的特点即是首先确定一个以某个像素为中心点的邻域,一般为方形邻域,也可以为圆形、十字形等等,然后将邻域中各像素的灰度值排序,取其中间值作为中心像素灰度的新值,这里领域被称为窗口,当窗口移动时,利用中值滤波可以对图像进行平滑处理。其算法简单,时间复杂度低,但其对点、线和尖顶多的图像不宜采用中值滤波。很容易自适应化。 Wiener维纳滤波:使原始图像和其恢复图像之间的均方误差最小的复原方法,是一种自适应滤波器,根据局部方差来调整滤波器效果。对于去除高斯噪声效果明显。 实验一:均值滤波对高斯噪声的效果 l=imread('C:\Documents and 桌面\1.gif');% 读取图像

J=imnoise(l,'gaussian',0,0.005);% 加入均值为0 ,方差为 0.005 的高斯噪声subplot(2,3,1);imshow(l); title(' 原始图像'); subplot(2,3,2); imshow(J); ti tle('加入高斯噪声之后的图像’); %采用MATLAB 中的函数filter2 对受噪声干扰的图像进行均值滤波 K1=filter2(fspecial('average',3),J)/255; % 模板尺寸为3 K2=filter2(fspecial('average',5),J)/255;% 模板尺寸为5 K3=filter2(fspecial('average',7),J)/255; % 模板尺寸为7 K4= filter2(fspecial('average',9),J)/255; % 模板尺寸为9 subplot(2,3,3);imshow(K1); ti tle(' 改进后的图像1'); subplot(2,3,4); imshow(K2); title(' 改进后的图像2'); subplot(2,3,5);imshow(K3); title(' 改进后的图像3'); subplot(2,3,6);imshow(K4); title(' 改进后的图像4');

基于matlab的图像去雾算法详细讲解与实现-附matlab实现源代码

本文主要介绍基于Retinex理论的雾霭天气图像增强及其实现。并通过编写两个程序来实现图像的去雾功能。 1 Rentinex理论 Retinex(视网膜“Retina”和大脑皮层“Cortex”的缩写)理论是一种建立在科学实验和科学分析基础上的基于人类视觉系统(Human Visual System)的图像增强理论。该算法的基本原理模型最早是由Edwin Land(埃德温?兰德)于1971年提出的一种被称为的色彩的理论,并在颜色恒常性的基础上提出的一种图像增强方法。Retinex 理论的基本内容是物体的颜色是由物体对长波(红)、中波(绿)和短波(蓝)光线的反射能力决定的,而不是由反射光强度的绝对值决定的;物体的色彩不受光照非均性的影响,具有一致性,即Retinex理论是以色感一致性(颜色恒常性)为基础的。 根据Edwin Land提出的理论,一幅给定的图像S(x,y)分解成两幅不同的图像:反射物体图像R(x,y)和入射光图像L(x,y),其原理示意图如图8.3-1所示。 图-1 Retinex理论示意图 对于观察图像S中的每个点(x,y),用公式可以表示为: S(x,y)=R(x,y)×L(x,y) (1.3.1)实际上,Retinex理论就是通过图像S来得到物体的反射性质R,也就是去除了入射光L的性质从而得到物体原本该有的样子。 2 基于Retinex理论的图像增强的基本步骤 步骤一: 利用取对数的方法将照射光分量和反射光分量分离,即: S'(x, y)=r(x, y)+l(x, y)=log(R(x, y))+log(L(x, y)); 步骤二:用高斯模板对原图像做卷积,即相当于对原图像做低通滤波,得到低通滤波后的图像D(x,y),F(x, y)表示高斯滤波函数: D(x, y)=S(x, y) *F(x, y); 步骤三:在对数域中,用原图像减去低通滤波后的图像,得到高频增强的图像G (x, y): G(x,y)=S'(x, y)-log(D(x, y)) ;

相关主题
相关文档 最新文档