how to estimate bright channel for color image in Matlab??

how to estimate the bright channel prior for color image (convert the image of the spatial domain into Brightness Domain) using the following formula:

Does anyone know how to do this in Matlab??

[Rik]

Do you have a bit more context for this? The math notation is not clear enough on its own for me to understand it. It looks like you would do max(A,3), but I'm really sure.

[ghada sandoub]

this is the form of estimating bright channel. where Ic is a color channel of I and Ω(x) is a local patch centered at x. y is a pixel in the local path Ω(x). the bright channel is the outcome of two maximum operators. The operator "max c∈{r,g,b}” is performed on each pixel for the three channels. The operator “max y∈Ω(x)” is a maximum filter. the local patch size for both operation was set as 15 × 15. so, i asked about how to do this in matlab?

[Rik]

So you want a window average, but with a maximum, with a 15-by-15-by-size(I,3) window?

[ghada sandoub]

I think that..

Answer 1

[Image Analyst]

The inner operation is dilation. The outer one is max. So simply dilate each color channel (take the max in a 15x15 window) with imdilate(). Like for the red channel.

dilatedRed = imdilate(redChannel, ones(3));

Same for blue and green, then recombine into an rgb image and use max

brightnessImage = max(rgbImage, [], 3);

A full demo is attached.

But I don't recall seeing this formula you gave. Taking the max is what you do to get L in the RGB-to-HSL, or to get the V in RGB-to-HSV formula. But they don't do dilation first. Not sure what that gets you, do you know? Why not just do it without dilation and use the built-in rgb2hsv() function

hsvImage = rgb2hsv(rgbImage); % Convert to HSV color space.
brightnessImage = hsvImage(:, :, 3); % Extract the V channel.

[Rik]

It is quite useful in cases like this, but why would a dilation be equivalent to a moving max filter? I would have expected something similar to a convolution to happen instead. Is that just the limitation in my knowledge of the field, or is there a deeper reason why it works like this?

[Image Analyst]

On a grayscale image, a dilation is a moving max filter. That's what it does. A convolution can't do that - it's just a weighted sum of the pixel values in the window. No way to get the max out of that operation. You could do it with nlfilter() - take the local max - but then you'd essentially be doing the dilation, just on your own. Also, a morphological erosion is the same as a moving min.

[ghada sandoub]

Thanks alot for your response. so, in another case if i want to estimate the dark channel prior for the color image with the following formula:

should i use erosion instead of dilation?? and make the min operation first then make the erosion or what should i do in this case?

[Image Analyst]

Yes, except the order is switched. It is:

Imin = imerode(min(Ig, [], 3), ones(windowSize));

Where are you getting these formulas from? I can understand that the dark current signal might be estimated as the min of the red, green, and blue pixels in the demosaicing window but I really doubt any camera uses a demosaicing algorithm where the interpolated pixels are the uniformly weighted average inside a 15x15 window.

[ghada sandoub]

I'm working now in researching and I get these formulas from a paper that interested in low light image enhancement. this paper uses the bright channel formula to estimate the insufficient illumination in the low light image.

Answer 2

[Rik]

Below you find 2.5 strategies. The first should be a lot faster, especially for larger images. With the last it is easier to see that the method returns the expected outcome.

clc
A=randi(20000,300,300,3);%random image
window=15;

%strategy 1:
tic
GetSingleLayerMax=@(A,window) nlfilter(A,[window window],@(x) max(x(:)));
maxIM=GetSingleLayerMax(A(:,:,1),window);
for n=2:size(A,3)
    tmp=cat(3,maxIM,GetSingleLayerMax(A(:,:,n),window));
    maxIM = max(tmp,[],3);
end
toc

%strategy 1b:
tic
GetSingleLayerMax2=@(A,window) colfilt(A,[window window],'sliding',@max);
maxIM_1b=GetSingleLayerMax2(A(:,:,1),window);
for n=2:size(A,3)
    tmp=cat(3,maxIM_1b,GetSingleLayerMax2(A(:,:,n),window));
    maxIM_1b = max(tmp,[],3);
end
toc

%strategy 2:
tic
maxIM2=zeros(size(A,1),size(A,2));
tmp=(1:window);tmp=tmp-mean(tmp);
tmp=ceil(tmp);%in case of even window size
[w_x,w_y,w_z]=meshgrid(tmp,tmp,1:size(A,3));
w=[w_x(:) w_y(:) w_z(:)];
for r=1:size(A,1)
    for c=1:size(A,2)
        %get all indices
        tmp=w;
        tmp(:,1)=w(:,1)+r;
        tmp(:,2)=w(:,2)+c;
        
        %remove invalid indices
        tmp(any(tmp<1,2),:)=[];
        tmp(tmp(:,1)>size(A,1),:)=[];
        tmp(tmp(:,2)>size(A,2),:)=[];
        
        %convert to linear indices
        ind=sub2ind(size(A),tmp(:,1),tmp(:,2),tmp(:,3));
        
        %compute the max
        maxIM2(r,c)=max(A(ind));
    end
end
toc

if isequal(maxIM,maxIM2) && isequal(maxIM_1b,maxIM2)
    disp('everything finished succesfully')
else
    disp('something went wrong (check for float rounding errors)')
end

[ghada sandoub]

thanks alot