Code covered by the BSD License  

Highlights from
Gaussian Filter, Determination of integer parameters

image thumbnail

Gaussian Filter, Determination of integer parameters

by

 

Program calculates integer parameters of n x n gaussian filters.

GaussFilter(varargin)
function varargout = GaussFilter(varargin)
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @GaussFilter_OpeningFcn, ...
                   'gui_OutputFcn',  @GaussFilter_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end


% --- Executes just before GaussFilter is made visible.
function GaussFilter_OpeningFcn(hObject, eventdata, handles, varargin)
global k kold;
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes GaussFilter wait for user response (see UIRESUME)
% uiwait(handles.figure1);
set(handles.dropSize,'String',{'3','5','7','9','11','13','15','17','19','21','23','25','27'});
set(handles.dropSize,'Value',6);
%k=6;
kold=1;
 contents = get(handles.dropSize,'String'); 
 k=(str2num(contents{get(handles.dropSize,'Value')})-1)/2;
SigIdeal=num2str(k/2.5); 

set(handles.txtSigmaIdeal,'String',{SigIdeal});
%set(handles.txtSigmaIdeal,'String',{'0.8'});
FilterCalculate(handles,'ProposeSeparable');
% --- Outputs from this function are returned to the command line.
function varargout = GaussFilter_OutputFcn(hObject, eventdata, handles) 
varargout{1} = handles.output;


% --- Executes on selection change in dropSize.
function dropSize_Callback(hObject, eventdata, handles)
global k ;
 contents = get(hObject,'String'); 
 k=(str2num(contents{get(hObject,'Value')})-1)/2;
SigIdeal=num2str(k/2.5); 

set(handles.txtSigmaIdeal,'String',{SigIdeal});

% --- Executes on button press in pbSigmaMin.
function pbSigmaMin_Callback(hObject, eventdata, handles)
FilterCalculate(handles,'SigmaMin');

% --- Executes on button press in pbSigmaPlus.
function pbSigmaPlus_Callback(hObject, eventdata, handles)
FilterCalculate(handles,'SigmaPlus');

% --- Executes on button press in pbSigmaIdeal.
function pbSigmaIdeal_Callback(hObject, eventdata, handles)
%global k kold Exponent2  Exponent2Min  Exponent2Max  SigmaU  Goalsum;
FilterCalculate(handles,'SigmaIdeal');

% --- Executes on button press in pbSigmaMinus.
function pbSigmaMinus_Callback(hObject, eventdata, handles)
FilterCalculate(handles,'SigmaMinus');

% --- Executes on button press in pbSigmaMax.
function pbSigmaMax_Callback(hObject, eventdata, handles)
FilterCalculate(handles,'SigmaMax');

% --- Executes on button press in pbProposeSeparable.
function pbProposeSeparable_Callback(hObject, eventdata, handles)
FilterCalculate(handles,'ProposeSeparable');

% --- Executes on button press in pbShowSeparable.
function pbShowSeparable_Callback(hObject, eventdata, handles)
FilterCalculate(handles,'ShowSeparable');

%Not used:================================
% --- Executes during object creation, after setting all properties.
function dropSize_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
    set(hObject,'BackgroundColor','white');
end

% --- Executes during object creation, after setting all properties.
function plot3D_CreateFcn(hObject, eventdata, handles)
% --- Executes during object creation, after setting all properties.
function tbl2Dfilter_CreateFcn(hObject, eventdata, handles)
function txtSigmaIdeal_Callback(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function txtSigmaIdeal_CreateFcn(hObject, eventdata, handles)

% --- Executes during object creation, after setting all properties.
function txtSigmaProposed_CreateFcn(hObject, eventdata, handles)
%End not used =======================
%End GUI ==================================================================

%Calculation===============================================================
function FilterCalculate(handles,Propose)
global  handlesGUI;
handlesGUI=handles;
    switch Propose
        case 'SigmaIdeal'
            ProposeSigmaIdeal;
        case 'SigmaMin'
            ProposeSigmaMin;
        case 'SigmaPlus'
            ProposeBigger;
        case 'SigmaMinus'
            ProposeSmaller;
        case 'SigmaMax'
            ProposeSigmaMax;
        case 'ProposeSeparable'
            ProposeSeparable;
        case 'ShowSeparable'
            ShowSeparable;
    end
%code for first 2 left buttons on GUI =========================   
    function ProposeSeparable
    global k F1D Sigma2Nr SeparableNr iSep countSeparable handlesGUI;
    set(handlesGUI.pbShowSeparable,'Enable','off');
    set(handlesGUI.txtChoice,'String','WAIT');
    set(handlesGUI.pbShowSeparable,'String','WAIT');
    set(handlesGUI.txtMessage,'String','WAIT');
    surf(handlesGUI.plot3D,[0 0;0 0]);
    set(handlesGUI.tbl2Dfilter,'Data',[0 0;0 0]);
    drawnow;
    Sigma1 = k / 3.2;
    SigmaOpt=k/2.5;
    SigmaMax=k/1.8;
    F1D=uint16(ones(600,k-1));%set all possible parameters to 1
    F1Dtemp=ones(1,k-1);
    Sigma2Nr=ones(600,1);%look up table sigma vers. row number
    SeparableNr=uint16(ones(40,3));%filter x, filter y, rating
    i=0;
    iSep=0;
    countSeparable=0;
    iOpt=0;%i for optimal sigma
    while Sigma1<SigmaMax
        i=i+1;
        %calculate next sigma by using old sigma and 
        %increasing center parameter of filter by one:
        Sigma2=k*sqrt(1/(2*log(exp(k*k/(2*Sigma1*Sigma1))-1)));
        Sigma2Nr(i)=Sigma2;
        Sigma2Nr(i+1)=Sigma2;
        Sigma2Nr(i+2)=Sigma2;        
        if iOpt==0 && Sigma2>=k/2.5
            iOpt=i;%get row number for optimal sigma
        end
        TwoSigmaSquared = 2 * Sigma2 * Sigma2;
        divisor = exp(-k * k / TwoSigmaSquared);%divisor makes corners of filter to value of 1        
        %F1D(1)=1; F1D(k+k+1)=1; already done

        F1D(i,1)=round(1/divisor);
        F1D(i+1,1)=F1D(i,1);
        F1D(i+2,1)=F1D(i,1);
 %make filter with ceil        
        boolNotGood=0;
        for j=1:k-1
            F1Dtemp(j)=exp(-(j * j) / TwoSigmaSquared) / divisor;
            F1D(i,j+1)=ceil(F1Dtemp(j));
            if F1D(i,j+1)>F1D(i,j)%Check if Filter is monoton descending,except the center may be equal
                boolNotGood=1;
            end
        end
        if boolNotGood==1
                i=i-1;% eliminate last values
        end
%make filter with round         
        countDuplicate=0;
        for j=1:k-1   
            F1D(i+1,j+1)=round(F1Dtemp(j));
            if F1D(i+1,j+1)==F1D(i,j+1)%Check for duplicate filter
                countDuplicate=countDuplicate+1;
            end
         end
        if countDuplicate==k-1 && F1D(i+1,1)==F1D(i,1)
            i=i-1;
        end
%make filter with floor         
        countDuplicate=0;
        for j=1:k-1      
            F1D(i+2,j+1)=floor(F1Dtemp(j));
            if F1D(i+2,j+1)==F1D(i+1,j+1)%Check for duplicate filter
                countDuplicate=countDuplicate+1;
            end            
        end
        if countDuplicate==k-1 && F1D(i+2,1)==F1D(i+1,1)
            i=i-1;
        end        

        i=i+2;
        Sigma1=Sigma2;
    end
     imax=i;
   %Now many posible 1D filters are in F1D
   
   %create 2d filters with the combination of all 1D filters:
     j=1;
     Fh=ones(1,k+k+1);%horizontal filter
     Fv=ones(1,k+k+1);%vertical filter
     
     while j<=imax%loop all filters
        for i2=2:k%stored filter F1D has only few values, create complete symetricalkernel
            Fh(i2)=F1D(j,k+2-i2);
            Fh(k+k+2-i2)=Fh(i2);
        end
        Fh(k+1)=F1D(j,1);%center value
        for i=1:j%loop all filter below or equal j-filter
            for i2=2:k
                Fv(i2)=F1D(i,k+2-i2);
                Fv(k+k+2-i2)=Fv(i2);
            end
            Fv(k+1)=F1D(i,1);%center value  
            summe=sum(sum(Fv'*Fh));
            if log2(summe)==round(log2(summe))
                countSeparable=countSeparable+1;
                SeparableNr(countSeparable,1)=j;
                SeparableNr(countSeparable,2)=i;
                if i==j
                  %very good rating for symetric filter (5000)
                  %substract distance from optimal sigma
                  SeparableNr(countSeparable,3)=5000-abs((1-Sigma2Nr(i)/SigmaOpt)*600); 
                else
                  %bad rating for distance to optimal sigma
                  %very bad rating for distance of sigma_x to sigma_y
                  SeparableNr(countSeparable,3)=4000-abs((1-Sigma2Nr(i)/SigmaOpt)*600)-abs((1-Sigma2Nr(j)/SigmaOpt)*600)-2000*abs((Sigma2Nr(i)-Sigma2Nr(j))/SigmaOpt);                                        
                end
            end
        end
        j=j+1;
     end%while

     if countSeparable>0
         SeparableNr=flipud(sortrows(SeparableNr(1:countSeparable,:),3));%sort rows descending
         iSep=1;
         ShowSeparable;
         set(handlesGUI.pbShowSeparable,'String',['show next of ' int2str(countSeparable) ' separable filter']);
         set(handlesGUI.pbShowSeparable,'Enable','on');
     else
         set(handlesGUI.pbShowSeparable,'String','no sep. filter found');
         set(handlesGUI.txtMessage,'String',['no sep. filter found' char(10) 'with divisor=2^x']);
     end

function ShowSeparable
global k  F1D Sigma2Nr SeparableNr iSep countSeparable Goalsum handlesGUI;
     Fh=ones(1,k+k+1);%horizontal filter
     Fv=ones(1,k+k+1);%vertical filter
     i=SeparableNr(iSep,1);
     j=SeparableNr(iSep,2);
     for i2=2:k
        Fh(i2)=F1D(j,k+2-i2);
        Fh(k+k+2-i2)=Fh(i2);
     end
    Fh(k+1)=F1D(j,1);%center value
    for i2=2:k
        Fv(i2)=F1D(i,k+2-i2);
        Fv(k+k+2-i2)=Fv(i2);
    end
    Fv(k+1)=F1D(i,1);%center value 
    Filter=Fv'*Fh;
    summe=sum(sum(Filter));
    Goalsum=summe;
    Exponent=round(log2(summe));

    set(handlesGUI.txtGoalSumExponent,'String',['Goalsum 2^' int2str(Exponent)]); 
     fullscreen = get(0,'ScreenSize');%[left, bottom, width, height]
    set(handlesGUI.figure1,'Position',[1 35 fullscreen(3) fullscreen(4)-60]); %maximze window
    surf(handlesGUI.plot3D,Filter);
    set(handlesGUI.tbl2Dfilter,'Data',Filter);%insert data into table
    %column width
    temp=floor(log10(Fv(k+1)*Fh(k+1)))+1;
    if temp<2;temp=2;end
    set(handlesGUI.tbl2Dfilter,'ColumnWidth',{5*temp+10});
    drawnow; % wait for: set(handlesGUI.tbl2Dfilter,'ColumnWidth',{5*temp+10})
    %table position, size
    %x from left, y from bottom, width, height
    Pos=get(handlesGUI.tbl2Dfilter,'Extent');
    Pos=[0.14 0.95-(k+k+2)*18/fullscreen(4) Pos(3) Pos(4)];
    set(handlesGUI.tbl2Dfilter,'Position',Pos);
    
    %3D plot position
    Pos=[0.14 0.95-(k+k+2)*18/fullscreen(4)-0.4 0.3 0.4];
    set(handlesGUI.plot3D,'Position',Pos);
    set(handlesGUI.txtTotalSum,'String',num2str(summe));
    set(handlesGUI.txtGoalSum,'String',num2str(Goalsum));
    if i==j
        set(handlesGUI.txtSigmaProposed,'String',num2str(Sigma2Nr(i)));
        set(handlesGUI.txtMessage,'String',['filter is symetric' char(10) 'take first row as kernel']);
        set(handlesGUI.txtMessage,'ForegroundColor',[0.0 0.7 0.0]);
    else
        set(handlesGUI.txtSigmaProposed,'String',['~' num2str((Sigma2Nr(i)+Sigma2Nr(j))/2)]);
        set(handlesGUI.txtMessage,'String',['take first row and' char(10) 'first column as kernels']);
        set(handlesGUI.txtMessage,'ForegroundColor',[0.0 0.0 0.7]);
    end
    set(handlesGUI.txtChoice,'String',[num2str(iSep) '. choice']);
    iSep=iSep+1;
    if iSep>countSeparable
        iSep=1;
    end
         set(handlesGUI.pbShowSeparable,'String',...
            ['show next of ' int2str(countSeparable) ' separable filter']);  
        
% code for other buttons  =======================================================      
function ProposeBigger
    global Exponent2;
    Exponent2 = Exponent2 - 1;
    iniPropose;
    Propose;
    ShowValues;
function ProposeSmaller
    global Exponent2;
    Exponent2 = Exponent2 + 1;
    iniPropose;
    Propose;
    ShowValues;

% 
function ProposeSigmaIdeal
global k Exponent2 SigmaU;
 iniPropose;
    SigmaU = k / 2.5;
    Exponent2 = 0; % calculate new Exponent (GoalSum)
    Propose;
    ShowValues;
 
function ProposeSigmaMin
    global k  Exponent2 SigmaU;
    iniPropose;
    SigmaU = k / 3;
    Exponent2 = 0; % calculate new Exponent (GoalSum)
    Propose;
    ShowValues;
% End Sub
 function ProposeSigmaMax
    global k Exponent2 SigmaU;
    iniPropose;
    SigmaU = k / 2;
    Exponent2 = 0; % calculate new Exponent (GoalSum)
    Propose;
    ShowValues;

    
function iniPropose()
global k kold Exponent2  Exponent2Min  Exponent2Max  SigmaU handlesGUI;
        set(handlesGUI.pbShowSeparable,'Enable','off'); 
        if k ~= kold % run only for new k
            Exponent2 = 0; % calculate new Exponent (GoalSum)
            SigmaU = k / 1.8;
            Propose;
            Exponent2Min = Exponent2;
            SigmaU = k / 3.2;
            Exponent2 = 0;% calculate new Exponent (GoalSum)
            Propose;
            Exponent2Max = Exponent2;
            kold=k;
        end

function Propose
global k Exponent2 SigmaU  Goalsum;
    SigmaMin = k / 4.3;%3.3
    SigmaMax = k / 1.0;%1.7
    summe = Summe1;
    if Exponent2 == 0 ; Exponent2 = round(log(summe) / log(2)); end
    Goalsum = 2 ^ Exponent2;
    for iteration = 1 : 100
        if summe > Goalsum 
            SigmaMin = SigmaU;
        else
            SigmaMax = SigmaU;
        end
        SigmaU = (SigmaMax + SigmaMin) / 2;
        Summeold = summe;
        summe = Summe1;
        if (summe == Goalsum) 
            break;
        end
        if (Summeold == summe) && iteration < 90 
            iteration = 90; % Try 10 more times
        end
    end% iteration
    
% End Sub
function ShowValues
global k Exponent2 SigmaU  Goalsum handlesGUI;
set(handlesGUI.txtGoalSumExponent,'String',['Goalsum 2^' int2str(Exponent2)]); 

    TwoSigmaSquared = 2 * SigmaU * SigmaU;
    divisor = exp(-k * k / TwoSigmaSquared);

    divisor = exp(-k * k * 2 / TwoSigmaSquared);
    summe=0;
    summeS=0;
    Filter1D=zeros(1,k+k+1);
    for i=-k:k
        Filter1D(k+1+i)=round(exp(-(i * i+k*k ) / TwoSigmaSquared) / divisor);
    end
    Filter=zeros(k+k+1);
    %FilterS=zeros(k+k+1);%separable filter
    for i = -k:k
        i1=k+1+i;
        
        for j = -k:k
            j1=k+1+j;
            Filter(i1,j1) = round(exp(-(i * i + j * j) / TwoSigmaSquared) / divisor);
            summe=summe+Filter(i1,j1);
            summeS=summeS+Filter1D(i1)* Filter1D(j1);
        end
    end% i
  
     fullscreen = get(0,'ScreenSize');%[left, bottom, width, height]
 %maximze window
     set(handlesGUI.figure1,'Position',[1 35 fullscreen(3) fullscreen(4)-60]); 
   surf(handlesGUI.plot3D,Filter);
    set(handlesGUI.tbl2Dfilter,'Data',Filter);
    %cell width
    temp=floor(log10(1/divisor))+1;
    if temp<2;temp=2;end
    cw1=ones(1,k+k+1)*5*temp+10;
    cw = num2cell(cw1);
    set(handlesGUI.tbl2Dfilter,'ColumnWidth',cw);
    %table position, size
    %x from left, y from bottom, width, height
    Pos=[0.14 0.95-(k+k+2)*18/fullscreen(4) ((k+k+1)*10*temp+40)/fullscreen(3) (k+k+2)*18/fullscreen(4)];
    set(handlesGUI.tbl2Dfilter,'Position',Pos);
    %3D plot position
    Pos=[0.14 0.95-(k+k+2)*18/fullscreen(4)-0.4 0.3 0.4];
    set(handlesGUI.plot3D,'Position',Pos);
    set(handlesGUI.txtTotalSum,'String',num2str(summe));
    set(handlesGUI.txtGoalSum,'String',num2str(Goalsum));
    set(handlesGUI.txtSigmaProposed,'String',num2str(SigmaU));
    if summeS==summe
        set(handlesGUI.txtMessage,'String',['filter is separable' char(10) 'take first row as kernel']);
        set(handlesGUI.txtMessage,'ForegroundColor',[0.0 0.7 0.0]);
    else
        set(handlesGUI.txtMessage,'String','filter is NOT separable');
        set(handlesGUI.txtMessage,'ForegroundColor',[1.0 0.0 0.0]);
    end

function summe=Summe1
global k SigmaU;
    TwoSigmaSquared = 2 * SigmaU * SigmaU;
    summe = 0;
 
    divisor = exp(-k * k * 2 / TwoSigmaSquared);
    for i = 1 : k
        for j = 0 : k
            summe = summe + round(exp(-(i * i + j * j) / TwoSigmaSquared) / divisor);
        end
    end
    summe = 4 * summe + round(1 / divisor);
    


Contact us