% ODOT - Demonstrates odot and oslash operators on 1D signal
%
% Usage:  [smooth, energy] = odot(f, K)
%
% Arguments:    f - a 1D signal
%               K - optional 'Weiner' type factor to condition the results
%                   where division by 0 occurs in the 'oslash' operation.
%                   K defaults to 'eps', If oscillations appear in the
%                   plots try increasing the value of K
%
% Returns:     energy - the Local Energy of the signal.
%              smooth - the smooth component of the signal obtained by
%                       performing the 'oslash' operator between the
%                       signal and its Local Energy.
% Plots:
%           Signal            Hilbert Transform of Signal
%           Local Energy      Hilbert Transform of Energy
%           Smooth Component  Reconstruction
%
% Smooth         = signal 'oslash' energy
% Reconstruction = energy  'odot'  smooth
%
% Glitches in the results will be seen at points where the Local Energy
% peaks - these points cause numerical grief.  These problems can be
% alleviated by smoothing the signal slightly and/or increasing the
% parameter K.
%
% This code only works for 1D signals - I am not sure how you would
% implement it for 2D images...

% Reference:  Robyn Owens.  "Feature-Free Images", Pattern Recognition
% Letters. Vol 15. pp 35-44, 1994.

% Copyright (c) 2004 Peter Kovesi
% School of Computer Science & Software Engineering
% The University of Western Australia
% http://www.csse.uwa.edu.au/
% 
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
% 
% The above copyright notice and this permission notice shall be included in 
% all copies or substantial portions of the Software.
%
% The Software is provided "as is", without warranty of any kind.

% March 2004

function [smooth, energy] = odot(f, K)
    
    if nargin == 1
	K = eps;
    end
    
    N = length(f);
    
    if rem(N,2) > 0  % odd No of elements - trim the last one off to make the
		     % number of elements even for simplicity.
	N = N - 1;
	f = f(1:N);
    end
    
    F = fft(f);
    
    F(1) = 0;                        % Kill DC component
    f = real(ifft(F));               % Reconstruct signal minus DC component
    
    % Perform 90 degree phase shift on the signal by multiplying +ve
    % frequencies of the fft by i and the -ve frequencies by -i, and then
    % inverting. 
    
    phaseshift = [ ones(1,N/2)*i ones(1,N/2)*(-i) ];
    
    % Hilbert Transform of signal
    h = real(ifft(F.*phaseshift));        
    
    energy = sqrt(f.^2 + h.^2);           % Energy
    % Hilbert Transform of Energy
    energyh = real(ifft(fft(energy).*phaseshift)); 
					  
    % smooth = signal 'oslash' energy
    divisor = energy.^2 + energyh.^2;   
    % Where divisor << K,  weinercorrector -> 0/K
    % Where divisor >> K,  weinercorrector -> 1    
    weinercorrector = divisor.^2 ./ ((divisor.^2)+K);
    smooth = (f.*energy + energyh.*h)./divisor .* weinercorrector;
    
    % Hilbert transform of smooth component
    smoothh=real(ifft(fft(smooth).*phaseshift));
    
    % Reconstruction = energy odot smooth
    recon = (smooth.*energy - smoothh.*energyh);     
    
    subplot(3,2,1), plot(f),title('Signal'); 
    subplot(3,2,2), plot(h),title('Hilbert Transform of Signal'); 
    subplot(3,2,3), plot(energy),title('Local Energy'); 
    subplot(3,2,4), plot(energyh),title('Hilbert Transform of  Energy'); 
    subplot(3,2,5), plot(smooth),title('Smooth Component'); 
    subplot(3,2,6), plot(recon),title('Reconstruction');