Why the error?

As a result of executing the following code, there is an error between a and b.

a and b have different values.

I used R2023a.

Please tell me the cause.

img=rand(3248,2048);
img2=img(1613:3248,:);
filter=ones(25,25);

img=imfilter(img,filter,0);
img2=imfilter(img2,filter,0);
a=img(1625:3248,:);
b=img2(13:end,:);

disp(max(abs(a-b),[],'all'));

     0

[Anjaneyulu Bairi]

Can you post the screenshot of the error here

[Dyuman Joshi]

"... there is an error between a and b."

What error?

[Azusa Tsukahara]

@Anjaneyulu Bairi

@Dyuman Joshi

Thank you for your answer.

I added the information.

[Dyuman Joshi]

Just to be sure, How is tmp defined?

Is it

tmp = max(abs(a-b),[],'all')

[Azusa Tsukahara]

tmp.m is the code above.

disp(max(abs(a-b),[],'all'));

is

1.1369e-13

Answer 1

[DGM]

The error is (almost) entirely a consequence of the difference in the way you're handling the region boundaries. When the image is filtered by imfilter(), it will be padded by half the filter width prior to processing. The areas of the output image that are within this same radius are influenced by whatever the padding content is.

With imfilter(), the default padding is zero. This is usually a bad choice, as it causes noticeable vignetting. While using replication or mirroring for the padding will produce less visually-noticeable effects in the output, the result will not generally be the same as if the original image content were used to begin with.

rng(3)
img = im2uint8(rand(5)); % a 5x5 test image
fk = ones(3)/9; % a 3x3 box average filter

% crop first, then filter
a = img(2:4,:); % middle 3 rows
a = imfilter(a,fk); % filter with default padding (zero)
imshow(a,'border','tight') % top and bottom rows are dark due to the influence of the padding

% crop first, but filter with different padding
b = img(2:4,:); % middle 3 rows
b = imfilter(b,fk,'replicate'); % filter with better padding
imshow(b,'border','tight') % top and bottom rows come from replicated image data

% filter first (with better padding), then crop
c = imfilter(img,fk,'replicate'); % filter with better padding
c = c(2:4,:); % middle 3 rows
imshow(c,'border','tight') % top and bottom rows come from original image data

% the interior of the middle rows are all equal, 
% because they're not influenced by the edge content
range([a(2,:); b(2,:); c(2,:)],1)

ans = 1×5

42 0 0 0 37

% but _none_ of the three methods are equal
% because the edge handling is different
% edge replication _is not_ equivalent to edge preservation
[immse(a,b) immse(b,c) immse(a,c)]

ans = 1×3

1.0e+03 * 1.4274 0.3520 2.4783

So that's the source of your error. I'd avoid using default padding with imfilter(), but if 'replicate' or 'symmetric' aren't good enough, and you want to crop prior to filtering, you'll have to over-crop by at least half the filter width and then crop again after filtering. That way you can be assured that the internal padding method doesn't influence the output.

Note that in both your example and mine, the ROI is the entire width of the image. In this case, the region can't be over-cropped on that axis, so the left and right edges will always be influenced by the internal padding method.

[Azusa Tsukahara]

Thank you for your answer!!

By using 'replicate', the difference between a and b became 0.

```you want to crop prior to filtering, you'll have to over-crop by at least half the filter width and then crop again after filtering```

In the example above, I think I'm doing what you said.

Also shown in the example below.

param.norishiro_V =  param.wV; %over-crop size

is not enough?

set(0,'DefaultFigureWindowStyle','docked');
set(0, 'defaultFigureColormap', jet(256));
clear;

% Param
param.sH = 2048; %horizontal image size
param.sV = 3248; %vertical image size
param.wH = 12;   % horizontal window size (windows size = 2*wH+1)
param.wV = 12;   % vertical window size (windows size = 2*wV+1)
param.division_calc_number = 4;% division number
param.division_height = param.sV/param.division_calc_number;
param.norishiro_V =  param.wV; %over-crop size


rng(1);
img=rand(param.sV,param.sH);
filter=ones(2*param.wV+1,2*param.wH+1);

% Calc at once
result_once=imfilter(img,filter,0);

% divide the image
% Calc by dividing
result_div=zeros(param.sV,param.sH);
for i=1:param.division_calc_number
    vst = max(1, (i-1).*param.division_height+1 - param.norishiro_V);
    ven = min(param.sV, i.*param.division_height + param.norishiro_V);
    %Extract divided image
    %over-crop half the filter width
    img_tmp=img(vst:ven,:);
    result_tmp=imfilter(img_tmp,filter,0);

    %Crop
    vst_out = (i-1) * param.division_height+1;
    ven_out = i * param.division_height;
    result_div(vst_out:ven_out,:)=result_tmp(1+vst_out-vst:param.division_height+vst_out-vst,:);
end

% Check if A and B match
% 0 : A = B
% else : have errors in the processing of imfilter() ???
disp(max(abs(result_once-result_div),[],'all'));
%figure;imagesc(result_once-result_div);axis image;

[DGM]

Oh good gravy you're right. I didn't notice that was the difference in the index ranges.

What I don't understand is where the problem is occurring if it's showing zero error. It shows zero error in your original example (the part executed in the editor, not the screenshot). This new example also seems to match when I try it.

% Param
param.sH = 512; %horizontal image size
param.sV = 384; %vertical image size
param.wH = 12;   % horizontal window size (windows size = 2*wH+1)
param.wV = 12;   % vertical window size (windows size = 2*wV+1)
param.division_calc_number = 4;% division number
param.division_height = param.sV/param.division_calc_number;
param.norishiro_V =  param.wV; %over-crop size

rng(1);
img=rand(param.sV,param.sH);
filter=ones(2*param.wV+1,2*param.wH+1);

% Calc at once
result_once=imfilter(img,filter,0);

% divide the image
% Calc by dividing
result_div=zeros(param.sV,param.sH);
for i=1:param.division_calc_number
    vst = max(1, (i-1).*param.division_height+1 - param.norishiro_V);
    ven = min(param.sV, i.*param.division_height + param.norishiro_V);
    
    %Extract divided image
    %over-crop half the filter width
    img_tmp=img(vst:ven,:);
    result_tmp=imfilter(img_tmp,filter,0);
    
    %Crop
    vst_out = (i-1) * param.division_height+1;
    ven_out = i * param.division_height;
    result_div(vst_out:ven_out,:)=result_tmp(1+vst_out-vst:param.division_height+vst_out-vst,:);
end

% Check if A and B match
% 0 : A = B
% else : have errors in the processing of imfilter() ???
err = max(abs(result_once-result_div),[],'all');
bothmatch = 0 == err

bothmatch = logical
   1

So is the error (the screenshot) something that happens intermittently or am I misinterpreting the question? I mean, Star Strider had a point regarding float rounding error, though I'd have to double check to make sure that there isn't a way that overcropping could let that happen.

[Azusa Tsukahara]

The error (the screenshot) happens all the time in my environment.

If I run your code, it will be below.

My environment is 9.14.0.2206163 (R2023a).

I created another environment(23.2.0.2380103 (R2023b)) and ran it with similar results.

[DGM]

Maybe I'm missing something obvious, but I don't see where the problem could be. If it's down to float rounding, maybe it's hardware dependent.

FWIW, I'm testing on R2019b in Linux on a ~12 year old AMD board. Don't know if that's even relevant.

[Azusa Tsukahara]

Thank you for your response.

I appreciate your answer.

I'll check float rounding and hardware.

P.S.:

The problem did not occur in another environment(9.6.0.1072779 (R2019a)).

Your example matches on R2019a.

[Azusa Tsukahara]

The problem

R2019a : don’t occur

R2022a : don’t occur

R2022b : occur

R2023a : occur

R2023b : occur

According to the following URL, imfilter() has changed in version R2022b.

多次元イメージの N 次元フィルター処理 - MATLAB imfilter - MathWorks 日本

I think the change contains a problem.

[DGM]

Ah. Maybe it's handling things differently internally. I'm going to have to find a way to dig into that if I want to be able to talk about its behavior.

Answer 2

[Hassaan]

img = rand(3248,2048);
filter = ones(25,25);

% Apply the filter to the whole image with 'same' to ensure output size matches input
img_filtered = imfilter(img, filter, 'same', 'replicate');

% Extract the parts of the image to compare
a = img_filtered(1625:3248,:);
img2 = img_filtered(1613:3248,:);
b = img2(13:end,:);

% Ensure that 'a' and 'b' are the same size
assert(isequal(size(a), size(b)), 'Arrays a and b must be the same size.');

% Compute the maximum absolute difference
disp(max(abs(a - b), [], 'all'));

In this corrected code, the 'same' option in the imfilter function call ensures the output image is the same size as the input. The 'replicate' option is used to handle the boundaries by replicating the edge values of the image, which should be done consistently for both parts of the image.

By extracting a and b after applying the filter to the entire image, we ensure that they have the same boundary conditions. Finally, we assert that a and b are the same size before performing the element-wise subtraction and computing the maximum absolute difference. If they are not the same size, MATLAB will throw an error message.

---------------------------------------------------------------------------------------------------------------------------------------------------

If you find the solution helpful and it resolves your issue, it would be greatly appreciated if you could accept the answer. Also, leaving an upvote and a comment are also wonderful ways to provide feedback.

Professional Interests

• Technical Services and Consulting
• Embedded Systems | Firmware Developement | Simulations
• Electrical and Electronics Engineering

Feel free to contact me.

[DGM]

The example is obfuscating the problem.

The output variables a and b are not identical because of the difference in the specified edge handling. They're identical because they're simply the exact same region cropped out of one array.

% Extract the parts of the image to compare
a = img_filtered(1625:3248,:);
% img2 = img_filtered(1613:3248,:); % this is just obfuscation
% b = img2(13:end,:);
b = img_filtered(1625:3248,:); % the addressing is identical to a

Explicitly specifying 'same' isn't necessary or relevant, as it's already the default behavior of imfilter().