DICOMファイルを開く

SPMにはDICOMファイルのheader情報を得る関数が用意されています。 これを使ってDICOM fileを開いてみます。 MatlabのImage Processing Toolboxの関数dicominfoに相当するのが spm_dicom_headers.mですが、 Image Processing Toolboxのdicomreadに相当する関数が見当たらないので以下のようにしました。
[filename,pathname]=uigetfile('*'); % DICOM ファイルを選択
hdr=spm_dicom_headers(strcat(pathname,filename));
fid=fopen(hdr{1}.Filename,'rb');
fseek(fid,0,'eof'); % ファイル終端から0 byte移動
Vsize=[hdr{1}.Rows,hdr{1}.Columns,hdr{1}.BitsAllocated];
Vsize(3)=Vsize(3)/8;
fseek(fid,-prod(Vsize),'eof');
V=fread(fid,inf,'int16'); % Vsize(3)をみて2byte整数と判明したと仮定
fclose(fid);
V=reshape(V,Vsize(1:2));
figure;imagesc(V');colormap(gray);daspect([1,1,1]);

Analyse形式と異なり、反転180度回転したものになっています。
Analyse形式にするには
rot90(Z',2);
とする必要があります。ああ面倒くさい!
Analyse形式の*.imgファイルを作ってみます。 尚、読み込むファイル名はxxx1.dcm〜xxx100.dcmなどとなっていると仮定します。
[filename,pathname]=uigetfile('*.dcm');
loadname=strcat(pathname,filename);
extension=loadname(length(loadname)-3:end);
loadname=loadname(1:end-5);
for n=1:1000;...
newname=strcat(pathname,sprintf('%d',n),extension);...
if exist(newname);hdr{n}=newname;else break;end;...
end;
header=spm_dicom_headers(hdr{1});
Vsize=[header{1}.Rows,header{1}.Columns,header{1}.BitsAllocated];
Vsize(3)=Vsize(3)/8;
datasize=prod(Vsize);
V=zeros(Vsize(1),Vsize(2),n-1);
for nn=1:n-1;...
fid=fopen(hdr{nn},'rb');...
fseek(fid,-datasize,'eof');...
W=fread(fid,inf,'int16');...% 2byte整数と仮定
V(:,:,nn)=rot90(reshape(W',Vsize(1:2)),2);...
end;
V=int16(V);
[filename,pathname]=uiputfile('*.img');
fid=fopen(strcat(pathname,filename),'wb');
fwrite(fid,V,'int16');
fclose(fid);
figure;imagesc(rot90(squeeze(V(:,128,:))));colormap(gray);daspect([1,1,1]);