%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % MatLab convolution demonstration by JV Stone, January 2010. % Copyright: This software is licensed under the Creative Commons License % (CC-GNU GPL version 2.0 or later). Acknowledgement of use should cite % the book associated with this software, which is: % Seeing: The Computational Approach to Biological Vision % by JP Frisby and JV Stone, MIT Press, 2010. % Contact: j.v.stone at shef.ac.uk %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % In order to see all figures neatly laid out, see comment about the function % tilefigs at the end of this block of code. % HTML FOR WEB % % %
% This software is licensed under the CC-GNU GPL version 2.0 or later. % % Clear all variables. clear all; %--------------------------- % Set various parameters for figures. jfontsize = 15; jlinewidth = 2; fignum = 1; az = 36; % azimuth for 3D plots. el = 30; % elevation for 3D plots. %--------------------------- % Set standard deviation of larger standard deviation (sd) used in filter. sd0 = 6; % make filter 3 times as long as sd. a = 3*sd0; inc = 1; x = -a:inc:a; % x contains n numbers from -a to a. y = x; [one n] = size(x); n2 = ceil(n/2); % need this later. % Put x and y values into a grid. [X Y] = meshgrid(x); %--------------------------- %%%%%%%%%%%%%% % Make small 2D gaussian filter. %%%%%%%%%%%%%% % Set standard deviation of smaller Gaussian. % Ratio of large sd/small sd = 1.6.(the DoG then appproximates the 2nd % derivative of a Gaussian). sd = sd0/1.6; sdsmall = sd; Z = exp(-(X.^2 + Y.^2)./(2*sd^2)); % calculate weights in filter % Normalise filter Z so that all weights in Z sum to one. Z = Z./sum(abs(Z(:))); Z1 = Z; %--------------------------- % Plot 2D filter Z1. figure(fignum); fignum=fignum+1; clf; hold on; % plot 2D Gaussian filter as a surface. surfl(X,Y,Z1); shading interp; colormap(gray); % plot cross-section of 2D Gaussian cross_section1=Z1(n2,:); plot3(a*ones(size(Z1(n2,:))),-a:inc:a,cross_section1,'r','LineWidth',jlinewidth); view(az,el); % set viewpoint, type 'rotate3d on' to change view with mouse title('Small 2D Gaussian filter Z1 (showing cross-section in red)','FontSize',jfontsize); xlabel('x-position','FontSize',jfontsize); ylabel('y-position','FontSize',jfontsize); zlabel('filter weight','FontSize',jfontsize); grid on; set(gca,'FontSize',jfontsize); set(gca,'Linewidth',jlinewidth); %%%%%%%%%%%%%% %%%%%%%%%%%%%% % Make large gaussian filter. %%%%%%%%%%%%%% % Set sd of larger Gaussian. sd = sd0; sdbig = sd; Z = exp(-(X.^2+Y.^2)./(2*sd^2)); Z = Z./sum(abs(Z(:))); Z2 = Z; %--------------------------- % Plot 2D filter Z2 upside down to make the point that we will subtract Z2 from Z1 figure(fignum); fignum=fignum+1; clf ; hold on; Z2 = -Z2; surfl(X,Y,Z2); shading interp;colormap(gray); cross_section2=Z2(n2,:); plot3(a*ones(size(Z2(n2,:))),-a:inc:a,cross_section2,'r','LineWidth',jlinewidth); g2 = Z2(n2,:); view(az,el); Z2 = -Z2; title('Large 2D filter Z2, plotted upside down','FontSize',jfontsize); xlabel('x-position','FontSize',jfontsize); ylabel('y-position','FontSize',jfontsize); zlabel('filter weight','FontSize',jfontsize); grid on; set(gca,'FontSize',jfontsize); set(gca,'Linewidth',jlinewidth); %%%%%%%%%%%%%% % Show cross-section of small and large guassians and their difference. figure(fignum); fignum=fignum+1; clf ; cross_section1=Z1(n2,:); cross_section2=Z2(n2,:); cross_section_diff=cross_section1-cross_section2; plot(cross_section_diff,'k--','LineWidth',jlinewidth); hold on; plot(cross_section1,'r','LineWidth',jlinewidth); plot(cross_section2,'b','LineWidth',jlinewidth); title('Red=small Guass, Blue=big Gauss, Black=difference (DOG)','FontSize',jfontsize); axis off; %%%%%%%%%%%%%% % Difference large and small 2D Gaussian filters to get difference of gaussian (DOG) filter. %%%%%%%%%%%%%% %--------------------------- % Make DOG filter as DIFFERENCE of two Gaussian filters. filter = Z1-Z2; filter = real(filter); % remove any imaginery components. filter = filter./sum(abs(filter(:))); % ensure filter weights sum to one. filter = filter-mean(filter(:)); % ensure filter has zero mean, %(so that the excitatory and inhibitory parts cancel each other). Z3=filter; % Plot 2D DOG filter Z3. figure(fignum); fignum=fignum+1; clf ; hold on; surfl(X,Y,Z3); shading interp;colormap(gray); plot3(a*ones(size(Z3(n2,:))),-a:inc:a,Z3(n2,:),'r','LineWidth',jlinewidth); g2=Z2(n2,:); view(az,el); title('2D DOG filter, with cross-section shown in red','FontSize',jfontsize); xlabel('x-position','FontSize',jfontsize); ylabel('y-position','FontSize',jfontsize); zlabel('filter weight','FontSize',jfontsize); grid on; set(gca,'FontSize',jfontsize); set(gca,'Linewidth',jlinewidth); zmin=min(cross_section_diff)*1.5; zmax=max(cross_section_diff)*1.5; set(gca,'ZLim',[zmin zmax]); axis off; %--------------------------- % Get a standard matlab image to filter. load clown; figure(fignum); fignum=fignum+1; clf; imagesc(X); colormap(gray);axis equal;axis off; title('Unfiltered image','FontSize',jfontsize); % CHECK %op = jconv_dog2d(X,sdsmall); %f(76);imagesc(op);axis equal; axis off;colormap(gray); %title('from jconv\_dog2d','FontSize',jfontsize); %%%%%%%%%%%%%% % Convolve image with DOG filter. %%%%%%%%%%%%%% result = conv2(X,filter,'same'); figure(fignum); fignum=fignum+1; imagesc(result); colormap(gray); axis equal; axis off; title('Filtered image','FontSize',jfontsize); %%%%%%%%%%%%%% % Convolve image with DOG filter the SLOW way. %%%%%%%%%%%%%% doslow=0; if doslow % filter the slow way! figure(fignum); fignum=fignum+1; colormap(gray);axis equal;drawnow;axis equal; axis off; im=X; im1=zeros(size(im)); [nrf ncf]=size(filter); nrf2=ceil(nrf/2); nrf3=nrf2-1; ncf2=ceil(ncf/2); ncf3=ncf2-1; [nrim ncim]=size(im); b=round(nrf*1.2); % unfiltered border around image fprintf('Now convolving the slow way (so be patient) ...\n'); % for each pixel in image, im(rim0,cim0) ... for rim0=b:nrim-b for cim0=b:ncim-b % reset tot sum of products to zero tot=0; % now step through each element of filter % and multiply each pixel around rim0,cim0 by corresponding % value in filter, keeping record of the sum. for rf=-nrf3:nrf3 for cf=-ncf3:ncf3 subtot = im(rim0+rf,cim0+cf) * filter(nrf2+rf,ncf2+cf); tot=tot+subtot; end end % store sum of products in pixel of convolution array at im1(rim0,cim0) im1(rim0,cim0)=tot; f(6);imagesc(im1); drawnow; end end f(fignum);imagesc(im1); end % doslow %%%%%%%%%%%%%% % Get edge map by finding all points where the convolution image pass % through zero. % This is actually a tricky operation, and I have devised a quick and dirty % (and fairly opaque!) way to do it below. But bear in mind that the % underlying objective is to find each pixel which is located between two pixels % with a positive and a negative value (eg adjacent pixels with values 0.9 0.04 -0.6 % imply that the middle pixel is roughly at a zero-crossing). % Uncomment commented code below to see intermediate steps to finding zero-crossings. %%%%%%%%%%%%%% % binarise filtered image - the edges (zero-crossings) lie at the border of % each black-white contour. result_binary = double(result>0); figure(fignum); fignum=fignum+1; imagesc(result_binary); colormap(gray); axis equal; axis equal; axis off; title('Binary convolved image, with edges at black-white borders','FontSize',jfontsize); % make a small 3x3 filter to begin finding edges. filter = ones(3,3); filter(2,2) = 0; filter = filter./sum(abs(filter(:))); filter(2,2) = -1; filter = -filter; % Find places in result_binary where images goes from 0 to 1. result1 = conv2(result_binary,filter,'same'); % figure(fignum); fignum=fignum+1; % imagesc(result1);colormap(gray); axis equal;drawnow;axis equal; axis off; % title('result1','FontSize',jfontsize); % Find pixels in result1 images that are in range -t to +t, where t is a % threshold, set to 0.1 below. These are the zero-corssings in the % convolved image. figure(fignum); fignum=fignum+1; t=0.1; % set threshold for range around zero. result2 = ~(result1>-t & result1 tile figure windows usage: tilefigs ([nrows ncols],border_in pixels) % % Restriction: maximum of 100 figure windows % % Without arguments, tilefigs will determine the closest N x N grid % % % %Charles Plum Nichols Research Corp. % % 70 Westview Street % %Tel: (781) 862-9400 Kilnbrook IV % %Fax: (781) 862-9485 Lexington, MA 02173 % % % maxpos = get (0,'screensize'); % determine terminal size in pixels % maxpos(4) = maxpos(4) - 25; % hands = get (0,'Children'); % locate fall open figure handles % hands = sort(hands); % sort figure handles % numfigs = size(hands,1); % number of open figures % % % maxfigs = 100; % % % if (numfigs>maxfigs) % figure limit check % disp([' More than ' num2str(maxfigs) ' figures ... get serious pal']) % return % end % % % if nargin == 0 % maxfactor = sqrt(maxfigs); % max number of figures per row or column % sq = [2:maxfactor].^2; % vector of integer squares % sq = sq(find(sq>=numfigs)); % determine square grid size % gridsize = sq(1); % best grid size % nrows = sqrt(gridsize); % figure size screen scale factor % ncols = nrows; % figure size screen scale factor % elseif nargin > 0 % nrows = tile(1); % ncols = tile(2); % if numfigs > nrows*ncols % disp ([' requested tile size too small for ' ... % num2str(numfigs) ' open figures ']) % return % end % end % if nargin < 2 % border = 0; % else % maxpos(3) = maxpos(3) - 2*border; % maxpos(4) = maxpos(4) - 2*border; % end % xlen = fix(maxpos(3)/ncols) - 30; % new tiled figure width % ylen = fix(maxpos(4)/nrows) - 45; % new tiled figure height % % % % tile figures by postiion % % Location (1,1) is at bottom left corner % pnum=0; % for iy = 1:nrows % ypos = maxpos(4) - fix((iy)*maxpos(4)/nrows) + border +25; % figure location (row) % for ix = 1:ncols % xpos = fix((ix-1)*maxpos(3)/ncols + 1) + border+7; % figure location (column) % pnum = pnum+1; % if (pnum>numfigs) % break % else % figure(hands(pnum)) % set(hands(pnum),'Position',[ xpos ypos xlen ylen ]); % move figure % end % end % end % return % END tilefigs.m % % % ------------------------------------------------------------------------- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % END OF FILE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%