imgfilterexpofn.m

非常好的数字处理教程

function imgfilterexpofn(action,datastruct)
    if nargin < 1
        action='init';
    end

    % http://www-mmdb.iai.uni-bonn.de/lehre/BIT/ss03_DSP_Vorlesung/matlab_demos/
	% Include lowpass, etc
    
	name = mfilename;
    figname = [name(1:end-2) '_fig'];
    f=findobj('Tag',figname);
    handles = get(f,'UserData');

    switch action
		case 'help'
			display_help(figname);
        case 'init'
            reset(handles.spectrum);
			set(handles.filtertag,'Visible','off');
			set(handles.filterval,'Visible','off');
			setdefaults;
		case 'load'
			handles.imagedata = load_imagedata;
            if ~isfield(handles.imagedata, 'filenamepath')
                return;
			elseif isfield(handles,'figwin');
                try,
                    close(handles.figwin);
                catch
                end
            end
			handles.original = handles.imagedata;
            make2DSpecPlot(handles);
			handles.figwin = showimage(handles.imagedata, ...
				get(handles.scale,'Value'), get(handles.inverse,'Value'));
			place_header(handles.figwin, handles);
            movegui(handles.figwin,'onscreen');
		case 'readinput'
            reset(handles.spectrum);
			set(handles.filtertag,'Visible','off');
			set(handles.filterval,'Visible','off');
			setdefaults;
			handles.imagedata = datastruct;
            clear datastruct;
            make2DSpecPlot(handles);
			handles.original = handles.imagedata;
			handles.figwin = showimage(handles.imagedata, ...
				get(handles.scale,'Value'), get(handles.inverse,'Value'));
			place_header(handles.figwin, handles);
            movegui(handles.figwin,'onscreen');
		case {'scale','inverse'}
			if isfield(handles, 'imagedata')
				figure(handles.figwin);
				X = get(gca, 'XLim');
				Y = get(gca, 'YLim');
				showimage(handles.imagedata, get(handles.scale,'Value'), get(handles.inverse,'Value'), handles.figwin);
				set(gca,'XLim',X,'YLim',Y);
			end
        case {'colormap','inversecolor'}
			if isfield(handles, 'imagedata')
				figure(handles.imgfilterexpo_fig);
				axes(handles.spectrum);
				contents = get(handles.colormap,'String');
				cmap = colormap(lower(contents{get(handles.colormap,'Value')}));
				if (get(handles.inversecolor,'Value'))
					colormap(flipud(cmap));
				else
					colormap(cmap);
				end
			end
		case 'zoom_reset'
			if isfield(handles, 'imagedata')
				showimage(handles.imagedata, get(handles.scale,'Value'), get(handles.inverse,'Value'), handles.figwin);
			end
        case 'plot'
            contents = get(handles.filtermenu,'String');
	        kernel_size = str2double(get(handles.kernel_size,'String'));
            switch (lower(contents{get(handles.filtermenu,'Value')}))
		        case 'average'
					mask = ones(kernel_size,kernel_size)./kernel_size^2;
				case 'gaussian'
					var = str2double(get(handles.filterval,'String'));
					G = gauss2d(kernel_size,kernel_size,var,var);
					mask = G/sum(sum(G));
		        otherwise
                    warndlg('Not able to perform this.');
					return;
            end
			figure;
			freqz2(mask);
		case 'apply_filter'
			if isfield(handles, 'imagedata')
				handles.imagedata2 = handles.imagedata;
				figure(handles.figwin);
				X = get(gca, 'XLim');
				Y = get(gca, 'YLim');
				handles.imagedata = apply_filter(handles);
                make2DSpecPlot(handles);
				showimage(handles.imagedata, get(handles.scale,'Value'), get(handles.inverse,'Value'), handles.figwin);
				set(gca,'XLim',X,'YLim',Y);
			end
		case 'apply_noise'
			if isfield(handles, 'imagedata')
				handles.imagedata2 = handles.imagedata;
				figure(handles.figwin);
				X = get(gca, 'XLim');
				Y = get(gca, 'YLim');
				handles.imagedata = apply_noise(handles);
                make2DSpecPlot(handles);
				showimage(handles.imagedata, get(handles.scale,'Value'), get(handles.inverse,'Value'), handles.figwin);
				set(gca,'XLim',X,'YLim',Y);
			end
		case {'undo_filter', 'undo_noise'}
			if isfield(handles, 'imagedata')
				figure(handles.figwin);
				X = get(gca, 'XLim');
				Y = get(gca, 'YLim');
				switch action
					case 'undo_filter'
						handles = undo('filter',handles);
					case 'undo_noise'
						handles = undo('noise',handles);
				end
                make2DSpecPlot(handles);
				showimage(handles.imagedata, get(handles.scale,'Value'), get(handles.inverse,'Value'), handles.figwin);
				set(gca,'XLim',X,'YLim',Y);
			end
		case {'imgaliasexpo','imgspecexpo'}
			if isfield(handles,'imagedata')
				imagedata = handles.imagedata;
				switch action
					case 'imgaliasexpo'
						imgaliasexpogui(imagedata);
					case 'imgspecexpo'
						imgspectrumexpogui(imagedata);
				end
			end
		case 'save'
			save_imagedata(handles.imagedata);
		case 'print'
			print_figure(f);
		case 'close'
			close_figure(f,figname(1:end-4));
			if isfield(handles, 'figwin'),
				close_figure(handles.figwin);
			end
            return;
	end
	set(f,'UserData',handles);

%-------------------------------------------------------------
function make2DSpecPlot(handles)
	figure(handles.imgfilterexpo_fig);
	axes(handles.spectrum);
	spectrum = log(abs(fftshift(fft2(handles.imagedata.data))));
	% Get rid of inf
	points = isinf(spectrum);
	spectrum(points) = 0;
	imagesc(spectrum);
	contents = get(handles.colormap,'String');
	cmap = colormap(lower(contents{get(handles.colormap,'Value')}));
	if (get(handles.inversecolor,'Value'))
		colormap(flipud(cmap));
	else
		colormap(cmap);
	end
	axis image;
    
%-------------------------------------------------------------
function imagedata = apply_filter(handles)
	imagedata = handles.imagedata;
	imagedata.data = double(imagedata.data);
	A = imagedata.data;
	[rows,columns] = size(A);

	kernel_size = str2double(get(handles.kernel_size,'String'));
	contents = get(handles.filtermenu,'String');
	h = waitbar(0,'Filtering');
	% Extend borders to accomodate kernels
	k = floor(str2double(get(handles.kernel_size,'String'))/2);
	% Create place to store values
	A = zeros(rows+2*k,columns+2*k);
	% Put the original image in the center
	A(k+1:rows+k,k+1:columns+k) = imagedata.data;

	% Four corners
	A(1:k,1:k) = imagedata.data(1,1);
	A(1:k,columns+k+1:columns+2*k) = imagedata.data(1,columns);
	A(rows+k+1:rows+2*k, columns+k+1:columns+2*k) = ...
		imagedata.data(rows,columns);
	A(rows+k+1:rows+2*k,1:k) = imagedata.data(rows,1);

	% Handle the four sides
	for i = 1:k
		A(i,k+1:columns+k) = imagedata.data(1,1:columns);
		A(k+1:rows+k,columns+k+i) = ...
			imagedata.data(1:rows,columns);
		A(rows+k+i,k+1:columns+k) = ...
			imagedata.data(rows,1:columns);
		A(k+1:rows+k,i) = imagedata.data(1:rows,1);
	end

	switch (lower(contents{get(handles.filtermenu,'Value')}))
		case 'median'
			if get(handles.horizontal,'Value')
				for i = 1:rows
					for j = 1:columns
						% Determine pixel block.
						pixel_block = reshape(A(i+k,j:j+2*k),kernel_size,1);
						% Find median without using median function.
						sorted_block = sort(pixel_block);
						median_position = ceil(kernel_size / 2); 
						block_median = sorted_block (median_position);
						B(i,j) = block_median;
					end
					waitbar(i/rows,h);
				end
			elseif get(handles.vertical,'Value')
				for i = 1:rows
					for j = 1:columns
						% Determine pixel block.
						pixel_block = reshape(A(i:i+2*k,j+k),kernel_size,1);
						% Find median without using median function.
						sorted_block = sort(pixel_block);
						median_position = ceil(kernel_size / 2); 
						block_median = sorted_block (median_position);
						B(i,j) = block_median;
					end
					waitbar(i/rows,h);
				end
			elseif get(handles.kernel,'Value')
				for i = 1:rows
					for j = 1:columns
						% Determine pixel block.
						pixel_block = reshape(A(i:i+2*k,j:j+2*k),kernel_size^2,1);
						% Find median without using median function.
						sorted_block = sort(pixel_block);
						median_position = ceil(kernel_size^2 / 2); 
						block_median = sorted_block (median_position);
						B(i,j) = block_median;
					end
					waitbar(i/rows,h);
				end
			end
		case 'average'
			b_k = ones(1,kernel_size)./kernel_size;
			if get(handles.horizontal,'Value')
				for i = 1:rows
					C = filter(b_k,1,A(i+k,:));
					B(i,:) = C(1,2*k:columns+2*k);
				end
			elseif get(handles.vertical,'Value')
				for i = 1:columns
					C = filter(b_k,1,A(:,i+k));
					B(:,i) = C(2*k:rows+2*k,1);
				end
			elseif get(handles.kernel,'Value')
				mask = ones(kernel_size,kernel_size)./kernel_size^2;
				B = zeros(rows, columns);
				for i=1:rows,
					for j=1:columns,
						current_block = A(i:i+2*k, j:j+2*k);
						B(i,j) = sum(sum(current_block.*mask));
					end
					waitbar(i/rows,h);
				end
				% Handle by convolution, but how to take care of the sides?
				%mask = ones(kernel_size,kernel_size)./kernel_size^2;
				%B = conv2(A,mask,'same');
			end
		case 'gaussian'
			var = str2double(get(handles.filterval,'String'));
			b_k = gauss2d(1,kernel_size,var,var)./kernel_size;
			if get(handles.horizontal,'Value')
				for i = 1:rows
					C = filter(b_k,1,A(i+k,:));
					B(i,:) = C(1,2*k:columns+2*k);
				end
			elseif get(handles.vertical,'Value')
				for i = 1:columns
					C = filter(b_k,1,A(:,i+k));
					B(:,i) = C(2*k:rows+2*k,1);
				end
			elseif get(handles.kernel,'Value')
				% Make Gaussian mask
				G = gauss2d(kernel_size,kernel_size,var,var);
				% Normalize
				G = G/sum(sum(G));
				for i=1:rows,
					for j=1:columns,
						current_block = A(i:i+2*k, j:j+2*k);
						B(i,j) = sum(sum(current_block.*G));
					end
					waitbar(i/rows,h);
				end
				%B = conv2(A,G,'same');
			end
		case 'del gaussian'
			var = str2double(get(handles.filterval,'String'));
			if get(handles.horizontal,'Value')
				G = d2dgauss(kernel_size,var,kernel_size,var,0);
				for i=1:rows,
					for j=1:columns,
						current_block = A(i:i+2*k, j:j+2*k);
						B(i,j) = sum(sum(current_block.*G));
					end
					waitbar(i/rows,h);
				end
			elseif get(handles.vertical,'Value')
				G = d2dgauss(kernel_size,var,kernel_size,var,pi/2);
				for i=1:rows,
					for j=1:columns,
						current_block = A(i:i+2*k, j:j+2*k);
						B(i,j) = sum(sum(current_block.*G));
					end
					waitbar(i/rows,h);
				end
			elseif get(handles.kernel,'Value')
				G = d2dgauss(kernel_size,var,kernel_size,var,pi/2);
				for i=1:rows,
					for j=1:columns,
						current_block = A(i:i+2*k, j:j+2*k);
						B1(i,j) = sum(sum(current_block.*G));
					end
					waitbar(i/(2*rows),h);
				end
				G = d2dgauss(kernel_size,var,kernel_size,var,0);
				for i=1:rows,
					for j=1:columns,
						current_block = A(i:i+2*k, j:j+2*k);
						B2(i,j) = sum(sum(current_block.*G));
					end
					waitbar((i+rows)/(2*rows),h);
				end
				B = sqrt(B1.*B1 + B2.*B2);
	
			end
				
		case 'laplacian'
			b_k = laplace2d(1,kernel_size)./kernel_size;
			if get(handles.horizontal,'Value')
				for i = 1:rows
					C = filter(b_k,1,A(i+k,:));
					B(i,:) = C(1,2*k:columns+2*k);
				end
			elseif get(handles.vertical,'Value')
				for i = 1:columns
					C = filter(b_k,1,A(:,i+k));
					B(:,i) = C(2*k:rows+2*k,1);
				end
			elseif get(handles.kernel,'Value')
				var = str2double(get(handles.filterval,'String'));
				L = laplace2d(kernel_size,kernel_size,var,var);
				L = L/sum(sum(L));
				for i=1:rows,
					for j=1:columns,
						current_block = A(i:i+2*k, j:j+2*k);
						B(i,j) = sum(sum(current_block.*L));
					end
					waitbar(i/rows,h);
				end
				%B = conv2(A,L,'same');
			end
	end
	close(h);

	imagedata.data = B;
	set(handles.undo,'String','Undo');


%-------------------------------------------------------------
function imagedata = apply_noise(handles)
	imagedata = handles.imagedata;

	contents = get(handles.noisemenu,'String');
	switch (lower(contents{get(handles.noisemenu,'Value')}))
		case 'salt and pepper'
			imagedata.data = noise(imagedata.data, 'sp', str2double(get(handles.noiseval,'String')));
		case 'gaussian'
			imagedata.data = noise(imagedata.data, 'ag', str2double(get(handles.noiseval,'String')));
	end
	set(handles.undo2,'String','Undo');


%-------------------------------------------------------------
function handles = undo(type,handles)
	switch type
		case 'filter'
			hObject = handles.undo;
		case 'noise'
			hObject = handles.undo2;
	end
			
	if (strcmp(get(hObject,'String'),'Undo'))
		% Store old data for redo
		handles.imagedata3 = handles.imagedata;
		% Reload old data
		handles.imagedata = handles.imagedata2;
		set(hObject,'String','Redo');
	else
		handles.imagedata = handles.imagedata3;
		set(hObject,'String','Undo');
	end