Unsharp Masking producing different results (imsharpen)

Hi, I tried to do unsharp masking for image sharpening with these two methods. Mathematically these 2 methods looking same but giving different results. Can anyone please tell me why these are producing different results.

Mathematically both B and B1 looks same to me but when I ran B and B1 are different .

Could anyone please help me where I am missing?

Answer 1

[DGM]

The error is largely due to order of operations and numeric class. You're working with a uint8 image. Tools like imfilter() will return an image of the same class as they're given, but internally, operations are done in floating point. That means you're dealing with the data being requantized on the output.

If you want to do a bunch of sequential operations on an image, it's best to cast the image as double before working on it, otherwise the rounding and truncation errors will keep stacking up. Cases where you're doing differences and sums between nominally similar images are a good example of where you can introduce errors. Your operation is I + a*(I - Igauss). Say a given pixel has the following theoretical result: 80 + 1*(80 - 100) = 60. If I and Igauss are both uint8, they can't have a negative difference, so the result is 80 instead.

It gets more complicated if the scaling factor is considered. If (I - Igauss) is negative and I and Igauss are uint8, the scaling factor will do nothing, since (in my example) the difference is truncated at zero.

Here's what you're starting with:

% Method1:
I = imread('pout.tif');
radius = 2;
amount = 1;

% Gaussian blurring filter
filtRadius = ceil(radius*2); 
filtSize = 2*filtRadius + 1; % this is a good general support size
gaussFilt = fspecial('gaussian',[filtSize filtSize],radius);

% High-pass filter
sharpFilt = zeros(filtSize,filtSize); % zeros (null filter)
sharpFilt(filtRadius+1,filtRadius+1) = 1; % 1 in the center (neutral filter)
sharpFilt = sharpFilt - gaussFilt; % difference with gaussian (I-Ig)
sharpFilt = amount*sharpFilt; % scaling (a*(I-Ig))
sharpFilt(filtRadius+1,filtRadius+1) = sharpFilt(filtRadius+1,filtRadius+1) + 1; % (I + a*(I-Ig))

% Method 1:
% all operations are applied in float to the filter
% we're never operating directly using integer-class arrays
B = imfilter(I,sharpFilt,'replicate','conv');

% Method 2: 
% the uint8 image is blurred
% but then we do a potentially lossy difference in uint8
% we attempt to scale integers by a potentially fractional float
% last, we do a potentially lossy sum in uint8
B1 = I + amount*(I - imfilter(I,gaussFilt,'replicate','conv')); % note what i changed here!

% Computing MSE
immse(B,B1)

ans = 8.1409

Now let's make one minor change

I = im2double(I); % make sure the image is double float first

% same thing as before
B = imfilter(I,sharpFilt,'replicate','conv');
B1 = I + amount*(I - imfilter(I,gaussFilt,'replicate','conv'));
immse(B,B1)

ans = 3.0421e-32

That's a tiny error, but bear in mind, the error scale also corresponds to the change in scale that's characteristic with the change in numeric class. Eventually at the end of our workflow, we're (probably) going to convert back to uint8. What is the actual realized error on the output?

B = im2uint8(B); % cast and rescale
B1 = im2uint8(B1); % cast and rescale
immse(B,B1) % that tiny error disappears when everything is quantized again

ans = 0

Note the change I made on the B1 line. Gotta keep your parameters in consistent use.