imlook3d.m

CT等医学图像的分割及其三维重建的MATLAB源代码

img = getimage(gcbf);
if (htype==1)
    figure;hist(double(img(:)),64);
    title('Slice histogram');
    ylabel('Frequency');
    xlabel('Intensity');
elseif (htype == 2)
    figure;hist(log(double(img(:))+1),64);     
else
    htype = 0;
end
return;

% --- Executes during object creation, after setting all properties.
function popup_view_CreateFcn(hObject, eventdata, handles)
if ispc
    set(hObject,'BackgroundColor','white');
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
%-------------------------------------------------------------
% --- Executes on selection change in popup_view.
%-------------------------------------------------------------
function popup_view_Callback(hObject, eventdata, handles)
v = get(hObject,'value');
if v ~=1
    nimg = convert_view(handles.oimage.CData,v);
else
    nimg = handles.oimage.CData;
end
[r,c,z]=size(nimg);
handles.image.CData = nimg;
handles.image.Xdata = [1 c];
handles.image.Ydata = [1 r];
handles.viewtype = v;
handles.imSize = [r,c,z];
reset(handles.axes1); axis off;
set(handles.SliceNumSlider,'Min',1,'Max',z,'Value',...
    round(z/2),'SliderStep',[1.0/double(z-1) 1.0/double(z-1)]);
set(handles.ImgObject,'Xdata',[1 c],'Ydata',[1 r]);
% Set the properties of the axes

set(handles.axes1,'XLim',[1 c],'YLim',[1 r],'PlotBoxAspectRatio',[c r 1]);

if ~isempty(handles.imageInfo)
    Sthick = handles.imageInfo.SliceThickness;
    if v == 2
        psize = handles.imageInfo.PixelSpacing(1);
        set(handles.axes1,'DataAspectRatio',[1 psize/Sthick  1]);
    elseif v == 3
        psize = handles.imageInfo.PixelSpacing(2);
        set(handles.axes1,'DataAspectRatio',[1 psize/Sthick  1]);
    else
        psize = handles.imageInfo.PixelSpacing;
        set(handles.axes1,'DataAspectRatio',[psize(1) psize(2)  1]);
    end
end
set(handles.SliceNumEdit,'String',num2str(round(z/2)));
cimg = handles.image.CData(:,:,round(z/2));
handles.currentSliceNumber = round(z/2);
set(handles.ImgObject,'Cdata',cimg);
guidata(hObject, handles);
drawnow
%-------------------------------------------------------------
% function s= convert_view(img,desired_view)
% img: input image
% It assumes that the original image is in transverse view(1)
% 1: transverse, 2: coronal, 3: sagittal views
% s: output image
function s= convert_view(img,desired_view)
[r,c,z]=size(img);
switch desired_view
    case 2
        s = permute(img,[3 2 1]);
        s=flipdim(s,1);
    case 3
        s = permute(img,[3 1 2]);
        s=flipdim(s,1);
    otherwise;
end
% --------------------------------------------------------------------
function readOrdinaryImageSeries_Callback(hObject, eventdata, handles)
[filename, pathname] = uigetfile({'*.*'},'Image');
if isequal(filename,0)|isequal(pathname,0)
    setstatus('File not found');
else
    cd(pathname)
    k = findstr(filename,'.');
    fnames = dir(['*' filename(k+1:end)]);
    timg = imread(filename);
    [r,c,z] = size(timg);
    numSlices = size(fnames,1);
    nimg = zeros([r,c,numSlices]);
    if z == 3
        hw=warndlg('Color images will be converted to gray');
        uiwait;
        nimg(:,:,1) = rgb2gray(timg);
    else
        nimg(:,:,1) = timg;
    end
    hwb = waitbar(0,'Reading Images, please wait...');
    for i=2:numSlices
        waitbar(i/numSlices, hwb);
        timg = imread(fnames(i).name);
        [r,c,z] = size(timg);
        if z == 3
            nimg(:,:,i) = rgb2gray(timg);
        else
            nimg(:,:,i) = timg;
        end
    end
    close(hwb);
    clear timg;
    [r,c,z]=size(nimg);
    handles.image.CData = nimg;
    handles.image.Xdata = [1 c];
    handles.image.Ydata = [1 r];
    handles.imSize = [r,c,z];
    handles.SliceThickness = 1;
    if z==1
        sliderStep = 1;
    else
        sliderStep = 1.0/double(z-1);
    end
    set(handles.SliceNumSlider,'Min',1,'Max',z,'Value',...
        round(z/2),'SliderStep',[sliderStep sliderStep]);
    set(handles.ImgObject,'Xdata',[1 c],'Ydata',[1 r]);
    % Set the properties of the axes
    set(handles.axes1,'XLim',[1 c],'YLim',[1 r],'PlotBoxAspectRatio',[c r 1]);
    set(handles.SliceNumEdit,'String',num2str(round(z/2)));
    cimg = handles.image.CData(:,:,round(z/2));
    set(handles.ImgObject,'Cdata',cimg);
    handles.oimage = handles.image;
    guidata(hObject, handles);
    clear nimg;
    if z > 1
        turnDisplayButtons(handles,'on');
        set(handles.popupmenu_imageModality,'visible','off');
    end
    drawnow  
end
% --------------------------------------------------------------------
function saveImageIntoWorkspace_Callback(hObject, eventdata, handles)
assignin('base','IIMG',handles.image.CData);
% --------------------------------------------------------------------
function imfileMenu_Callback(hObject, eventdata, handles)
% --------------------------------------------------------------------
function ImMaps_Callback(hObject, eventdata, handles)
% --------------------------------------------------------------------
% --------------------------------------------------------------------
function readWorkspaceImage_Callback(hObject, eventdata, handles)
Vars = evalin('base','who;');
if ~isempty (Vars)
    [s,v] = listdlg('PromptString','Select an image:',...
        'SelectionMode','single',...
        'Listsize',[150,200],...
        'ListString',Vars);
    if (v)
        imlook3d(evalin('base',Vars{s}));
    end
else
    disp('No workspace variable');
end

% --------------------------------------------------------------------
% Reads series of dicom images into a 3D array
% --------------------------------------------------------------------
function readDicomPETSeries_Callback(hObject, eventdata, handles)
[filename, pathname] = uigetfile({'*.dcm';'*.*'},'Image');
if isequal(filename,0)|isequal(pathname,0)
    setstatus('File not found');
else
    cd(pathname)
    k = findstr(filename,'.');
    fnames = dir([filename(1:2) '*']);
    info = dicominfo(filename);
    r = info.Width;
    c = info.Height;
    numSlices = length(fnames);
    nimg = zeros([r,c,numSlices]);
    hw = waitbar(0,'Reading volume image, please wait...');
    for i=1:numSlices
        info = dicominfo(fnames(i).name);
        nimg(:,:,info.ImageIndex) = info.RescaleSlope*dicomread(info)+info.RescaleIntercept;
        waitbar(i/numSlices,hw);
    end
    close(hw);
    
    fields = {'Filename','FileModDate','FileSize','ImageType','Width','Height','BitDepth',...
        'NumberOfFrames','Rows','Columns','BitsAllocated','SmallestImagePixelValue',...
        'LargestImagePixelValue','PixelDataGroupLength','StudyDescription','SeriesDescription',...
        'StudyID','ImageID','PatientID','SliceThickness','StudyDate','AcquisitionDate',...
        'AcquisitionTime','PixelSpacing'};
    
    myinfo = keepfield(info,fields);
    
    handles.imageInfo = myinfo;
 
    [r,c,z]=size(nimg);
    handles.image.CData = nimg;
    handles.image.Xdata = [1 c];
    handles.image.Ydata = [1 r];
    handles.imSize = [r,c,z];
    set(handles.SliceNumSlider,'Min',1,'Max',z,'Value',...
        round(z/2),'SliderStep',[1.0/double(z-1) 1.0/double(z-1)]);
    set(handles.ImgObject,'Xdata',[1 c],'Ydata',[1 r]);
    % Set the properties of the axes
    set(handles.axes1,'XLim',[1 c],'YLim',[1 r],'PlotBoxAspectRatio',[c r 1]);
    set(handles.SliceNumEdit,'String',num2str(round(z/2)));
    cimg = handles.image.CData(:,:,round(z/2));
    set(handles.ImgObject,'Cdata',cimg);
    handles.oimage = handles.image;
    guidata(hObject, handles);
    clear nimg;
    if z>1
        turnDisplayButtons(handles,'on')
    end
    handles.imageModality = 1;
    set(handles.popupmenu_imageModality,'visible','on');
    set(handles.popupmenu_imageModality,'value',1);
    drawnow  
end
% --------------------------------------------------------------------
function saveCurrentSlice_Callback(hObject, eventdata, handles)
[filename, pathname,filterindex] = uiputfile({'*.png';'*.tif';'*.jpg'},'Image');
fformat = {'png','tif','jpg'}; 
if isequal(filename,0) | isequal(pathname,0)
    warndlg('No file name entered!');
else
    I= getframe(handles.axes1);
    I=frame2im(I);
    buffer=pwd;
    cd(pathname);
    imwrite(I,filename,fformat{filterindex});
    cd(buffer);
end
% --------------------------------------------------------------------
function adjust_intensity_window_Callback(hObject, eventdata, handles)
setappdata(0  ,  'hMainGui'    , gcf);
hMainGui = getappdata(0, 'hMainGui');
maxINT = max(max(max(handles.image.CData)));
minINT = min(min(min(handles.image.CData)));
cclim = get(handles.axes1,'CLim');
intensityParam = struct('lowerLimit',cclim(1),'upperLimit',cclim(2),'maxIntensity',maxINT, 'minIntensity',minINT);
setappdata(hMainGui,  'intensityParam', intensityParam);
adjustIntensityRangeGui;
uiwait;
hMainGui = getappdata(0,'hMainGui');
handles.intensityParam = getappdata(hMainGui,  'intensityParam');
guidata(hObject, handles);
% -----------------------------------------------------------------
function global_intensity_windowing_Callback(hObject, eventdata, handles)
str = get(gcbo,'Checked');
if strcmp(str,'off')
    handles.global_windowing =1;
    set(gcbo,'Checked','on');
else
    handles.global_windowing =0;
    set(gcbo,'Checked','off');
end
guidata(hObject,handles);
% --------------------------------------------------------------------
function renderVolume_Callback(hObject, eventdata, handles)
img= handles.image.CData;
if size(img,3) > 1
    prompt = {'Enter Lower Threshold:','Enter Upper Threshold:'};
    dlg_title = 'Enter Thresholds'; num_lines = 1; def = {'1','2'};
    answer = inputdlg(prompt,dlg_title,num_lines,def);
    T(1) = str2num(answer{1,1});
    T(2) = str2num(answer{2,1});
    nimg = img >= T(1) & img <= T(2);
    s = size(img);
    vlim=[1 s(2) 1 s(1) 1 s(3)];
    figure;
    % set(get(hndl,'CurrentAxes'),'YDir','rev');
    [x y z D] = subvolume(nimg,vlim);
    p = patch(isosurface(x,y,z,D,0),'FaceColor','r','EdgeColor','none');
    p2 = patch(isocaps(x,y,z,D,0),'FaceColor','interp','EdgeColor','none');
    set(gca,'XTick',[],'YTick',[],'ZTick',[],'Box','on','Color','c');
    isonormals(x,y,z,D,p);
    axis tight;
    daspect([1 1 1])
    lighting gouraud
else
    msgbox('This is not a volume image','Error','modal');
end

% --------------------------------------------------------------------
% READS CT DICOM IMAGES
% --------------------------------------------------------------------
function readDicomCTSeries_Callback(hObject, eventdata, handles)
[filename, pathname] = uigetfile({'*.dcm';'*.*'},'Image');
if isequal(filename,0)|isequal(pathname,0)
    setstatus('File not found');
else
    cd(pathname)
    k = findstr(filename,'.');
    fnames = dir(['*' filename(k(end)+1:end)]);
    info = dicominfo(filename);
    r = info.Width;
    c = info.Height;
    numSlices = length(fnames);
    nimg = zeros([r,c,numSlices]);
    hwb = waitbar(0,'Reading Images, please wait...');
    for i=1:numSlices
        waitbar(i/numSlices, hwb);
        info = dicominfo(fnames(i).name);
        nimg(:,:,i) = info.RescaleSlope*dicomread(info)+info.RescaleIntercept;
    end
    close(hwb);
    
    fields = {'Filename','FileModDate','FileSize','ImageType','Width','Height','BitDepth',...
        'NumberOfFrames','Rows','Columns','BitsAllocated','SmallestImagePixelValue',...
        'LargestImagePixelValue','PixelDataGroupLength','StudyDescription','SeriesDescription',...
        'StudyID','ImageID','PatientID','SliceThickness','StudyDate','AcquisitionDate','AcquisitionTime','PixelSpacing'};
    
    myinfo = keepfield(info,fields);
    handles.imageInfo = myinfo;
    
    [r,c,z]=size(nimg);
    handles.image.CData = nimg;
    handles.image.Xdata = [1 c];
    handles.image.Ydata = [1 r];
    handles.imSize = [r,c,z];
    set(handles.SliceNumSlider,'Min',1,'Max',z,'Value',...
        round(z/2),'SliderStep',[1.0/double(z-1) 1.0/double(z-1)]);
    set(handles.ImgObject,'Xdata',[1 c],'Ydata',[1 r]);
    % Set the properties of the axes
    set(handles.axes1,'XLim',[1 c],'YLim',[1 r],'PlotBoxAspectRatio',[c r 1]);
    set(handles.SliceNumEdit,'String',num2str(round(z/2)));
    cimg = handles.image.CData(:,:,round(z/2));
    set(handles.ImgObject,'Cdata',cimg);
    handles.oimage = handles.image;
    guidata(hObject, handles);
    clear nimg;
    handles.imageModality = 4;
    if z>1
        turnDisplayButtons(handles,'on')
    end
    handles.imageModality = 2;
    set(handles.popupmenu_imageModality,'visible','on');
    set(handles.popupmenu_imageModality,'value',2);