0.简单说明:
matlab读取DICOM医学图像文件的方式很简单直接,matlab自带dicomread函数。
>> help dicomread
dicomread - Read DICOM image
This MATLAB function reads the image data from the compliant Digital Imaging and
Communications in Medicine (DICOM) file filename.
X = dicomread(filename)
X = dicomread(info)
[X,map] = dicomread(...)
[X,map,alpha] = dicomread(...)
[X,map,alpha,overlays] = dicomread(...)
[...] = dicomread(filename, 'frames', v)
[...] = dicomread(___,'UseVRHeuristic',TF)
但是,如果想要另外一种医学图像类型.mha或者.mhd的话,matlab就没有自带对应函数了。如果想继续用dicomread()函数,就需要用itk自己转换成dicom了,反而麻烦了。所以直接access MetaImage为宜,经自己测试,附上代码如下:
1. 代码正文:
1. 测试脚本1:
mhaAccessDemo.m
close all
clear
clc
%% mha or mhd file access demo
hdInfo = mha_read_header('./liverMaskDemo.mha');
imgData = mha_read_volume(hdInfo);
imshow((imgData(:,:,round(end/2))));% 显示Axial(Z轴中间层数据)
测试结果图:
2. 图像头文件信息获取函数:
mha_read_header.m
function info =mha_read_header(filename)
% Function for reading the header of a Insight Meta-Image (.mha,.mhd) file
%
% info = mha_read_header(filename);
%
% examples:
% 1, info=mha_read_header()
% 2, info=mha_read_header('volume.mha');
if(exist('filename','var')==0)
[filename, pathname] = uigetfile('*.mha', 'Read mha-file');
filename = [pathname filename];
end
fid=fopen(filename,'rb');
if(fid<0)
fprintf('could not open file %s\n',filename);
return
end
info.Filename=filename;
info.Format='MHA';
info.CompressedData='false';
readelementdatafile=false;
while(~readelementdatafile)
str=fgetl(fid);
s=find(str=='=',1,'first');
if(~isempty(s))
type=str(1:s-1);
data=str(s+1:end);
while(type(end)==' '); type=type(1:end-1); end
while(data(1)==' '); data=data(2:end); end
else
type=''; data=str;
end
switch(lower(type))
case 'ndims'
info.NumberOfDimensions=sscanf(data, '%d')';
case 'dimsize'
info.Dimensions=sscanf(data, '%d')';
case 'elementspacing'
info.PixelDimensions=sscanf(data, '%lf')';
case 'elementsize'
info.ElementSize=sscanf(data, '%lf')';
if(~isfield(info,'PixelDimensions'))
info.PixelDimensions=info.ElementSize;
end
case 'elementbyteordermsb'
info.ByteOrder=lower(data);
case 'anatomicalorientation'
info.AnatomicalOrientation=data;
case 'centerofrotation'
info.CenterOfRotation=sscanf(data, '%lf')';
case 'offset'
info.Offset=sscanf(data, '%lf')';
case 'binarydata'
info.BinaryData=lower(data);
case 'compresseddatasize'
info.CompressedDataSize=sscanf(data, '%d')';
case 'objecttype',
info.ObjectType=lower(data);
case 'transformmatrix'
info.TransformMatrix=sscanf(data, '%lf')';
case 'compresseddata';
info.CompressedData=lower(data);
case 'binarydatabyteordermsb'
info.ByteOrder=lower(data);
case 'elementdatafile'
info.DataFile=data;
readelementdatafile=true;
case 'elementtype'
info.DataType=lower(data(5:end));
case 'headersize'
val=sscanf(data, '%d')';
if(val(1)>0), info.HeaderSize=val(1); end
otherwise
info.(type)=data;
end
end
switch(info.DataType)
case 'char', info.BitDepth=8;
case 'uchar', info.BitDepth=8;
case 'short', info.BitDepth=16;
case 'ushort', info.BitDepth=16;
case 'int', info.BitDepth=32;
case 'uint', info.BitDepth=32;
case 'float', info.BitDepth=32;
case 'double', info.BitDepth=64;
otherwise, info.BitDepth=0;
end
if(~isfield(info,'HeaderSize'))
info.HeaderSize=ftell(fid);
end
fclose(fid);
3. 图像数据(坐标、灰度值)函数:
mha_read_volume.m
function V = mha_read_volume(info)
% Function for reading the volume of a Insight Meta-Image (.mha, .mhd) file
%
% volume = tk_read_volume(file-header)
%
% examples:
% 1: info = mha_read_header()
% V = mha_read_volume(info);
% imshow(squeeze(V(:,:,round(end/2))),[]);
%
% 2: V = mha_read_volume('test.mha');
if(~isstruct(info)), info=mha_read_header(info); end
switch(lower(info.DataFile))
case 'local'
otherwise
% Seperate file
info.Filename=fullfile(fileparts(info.Filename),info.DataFile);
end
% Open file
switch(info.ByteOrder(1))
case ('true')
fid=fopen(info.Filename,'rb','ieee-be'); %去掉fopen第一个参数的转置'
otherwise
fid=fopen(info.Filename,'rb','ieee-le'); %ditto
end
switch(lower(info.DataFile))
case 'local'
% Skip header
fseek(fid,info.HeaderSize,'bof');
otherwise
fseek(fid,0,'bof');
end
datasize=prod(info.Dimensions)*info.BitDepth/8;
switch(info.CompressedData(1))
case 'f'
% Read the Data
switch(info.DataType)
case 'char'
V = int8(fread(fid,datasize,'char'));
case 'uchar'
V = uint8(fread(fid,datasize,'uchar'));
case 'short'
V = int16(fread(fid,datasize,'short'));
case 'ushort'
V = uint16(fread(fid,datasize,'ushort'));
case 'int'
V = int32(fread(fid,datasize,'int'));
case 'uint'
V = uint32(fread(fid,datasize,'uint'));
case 'float'
V = single(fread(fid,datasize,'float'));
case 'double'
V = double(fread(fid,datasize,'double'));
end
case 't'
switch(info.DataType)
case 'char', DataType='int8';
case 'uchar', DataType='uint8';
case 'short', DataType='int16';
case 'ushort', DataType='uint16';
case 'int', DataType='int32';
case 'uint', DataType='uint32';
case 'float', DataType='single';
case 'double', DataType='double';
end
Z = fread(fid,inf,'uchar=>uint8');
V = zlib_decompress(Z,DataType);
end
fclose(fid);
V = reshape(V,info.Dimensions);
function M = zlib_decompress(Z,DataType)
import com.mathworks.mlwidgets.io.InterruptibleStreamCopier
a=java.io.ByteArrayInputStream(Z);
b=java.util.zip.InflaterInputStream(a);
isc = InterruptibleStreamCopier.getInterruptibleStreamCopier;
c = java.io.ByteArrayOutputStream;
isc.copyStream(b,c);
M=typecast(c.toByteArray,DataType);
2. 写在最后:
周末偷闲临时小纪,更多资料大家可以访问mathworks官网去获取。