% dip05_5.m                gr09mar15  gr23oct23
% Laplacians and sharpening

%%%%%%% note
iptsetpref('ImshowBorder','tight');
iptsetpref('ImshowInitialMagnification','fit');
%%%%%%%%

clear all;

[x,y] = meshgrid(1:100,1:200);
z = (x).^2 + (y-100).^2;
halfdisk = double(z<6000)*255;
figure(1);  imshow(halfdisk,[]);
return;
%%
 f = halfdisk;
% f = double(imread('.\img\lena.tif')); f=f(200:350,200:350);
[nr,nc] = size(f);

% f = f + randn(size(f))*7;
                                  % Laplacians
lapl1a = [0 1 0; 1 -4 1; 0 1 0]; % correct
g1a = filter2(lapl1a,f);
lapl1b = - lapl1a;              % wrong sign
g1b = filter2(lapl1b,f);
lapl2a = [1 1 1; 1 -8 1; 1 1 1];  % augmented
g2a = filter2(lapl2a,f);
lapl2b = - lapl2a;           % augmented, wrong sign
g2b = filter2(lapl2b,f);

imshow([3*g1a+128, 3*g1b+128; 1.5*g2a+128, 1.5*g2b+128],[0 255]);

%%
lambda = 1;                        % Unsharp Masking
gum = f - lambda*g1a;   % try g1b too
imshow([f, gum],[0 255]);

%%
figure
g1dz = g1a .* (abs(g1a)>40);        % Unsharp Masking with dead zone
gumdz = f - lambda*g1dz;
imshow([f, gumdz],[0 255]);

%%
sobx = [1 2 1; 0 0 0; -1 -2 -1];   soby = sobx';    % with Sobel sub-operators
gsob = filter2(sobx,f) + filter2(soby,f);
gumsob = f + .3*gsob;
imshow([f, gumsob],[0 255]);

%%
mu = 30;   lambda = 1;         % Rational Unsharp Masking  
% mu = 20;   lambda = 2;        % (better results than dead-zone-UM
A = f;                          %  but computationally heavier)
a = [A(:,1),A(:,1:nc-1)];    b = A;    c = [A(:,2:nc),A(:,nc)];
d = [A(1,:);A(1:nr-1,:)];                   e = [A(2:nr,:);A(nr,:)];
zx = (a-c).^2 ./ mu.^2;     zy = (d-e).^2 ./ mu.^2;
zz = (2*b-a-c) .* (zx.^2 + 4*zx) ./ (3*zx.^2 + 2) + (2*b-d-e) .* (zy.^2 + 4*zy) ./ (3*zy.^2 + 2);
B = b + lambda * zz;
B(B<0) = 0;   B(B>255) = 255;
imshow([f B],[0 255]);
%%
%  plot RUM control function
zx = ([0:255]).^2 ./ mu.^2;    ctrl = (zx.^2 + 4*zx) ./ (3*zx.^2 + 2);
figure(2);  plot([0:255],ctrl,'LineWidth',3); axis tight;
xlabel('Absolute lum. difference'); ylabel('Control term');
return;
