As far as I understand it, there isn't one equation; there are various techniques of differing complexity. I don't see any that have already been implemented on the File Exchange, but maybe I just missed something.
There are a number of papers describing compensation methods. You're free to search and implement one. Here's one I just grabbed that seems relatively simple:
At the very least, it has some amount of literature review that should help.
EDIT:
This is more or less what that paper does. I omitted the median+guided filter attempt at noise reduction. I doubt it's really necessary.
The entire core of the gamma correction scheme is far simpler than it's made out to be. I decided to be a little less heavy-handed with the contrast stretching and gamma than the paper was. The output is the weighted combination of the stretched image and a CLAHE image. Play with alph to vary the influence of the two. When alph = 0, the CLAHE process dominates, and the image will have relatively uniform local contrast. When alph = 1, shadows and highlights will be emphasized at the expense of local contrast in those regions.
% based on:
% Underwater image enhancement based on red channel
% weighted compensation and gamma correction model
% Xiang, et al
% DOI: 10.29026/oea.2018.180024
RGB = imread('bbbb.png');
RGB = im2double(RGB);
imshow(RGB)
% output adjustment parameters
contrastclip = 0.01; % [0 0.5]
gammaadj = 1; % [0 Inf], null condition is 1
alph = 0.5; % [0 1]
lam = [620 540 450]; % nanometers
% b(lam_ref) is unneeded, since it cancels anyway
blam = (-0.00113*lam + 1.62517);
% B_lam is top 0.1% of each channel
Blam = permute(quantile(RGB,1-0.001,[1 2]),[1 3 2]);
% attenuation coefficient ratios
cgcr = (blam(2)*Blam(1))/(blam(1)*Blam(2));
cbcr = (blam(3)*Blam(1))/(blam(1)*Blam(3));
% weights
wrgb = [1 cgcr cbcr]/(1 + cgcr + cbcr);
% construct new image
Rnew = imapplymatrix(wrgb,RGB);
RGBnew = RGB;
RGBnew(:,:,1) = Rnew;
imshow(RGBnew)
% adjust contrast/gamma
% adjust input levels and gamma
inlim = stretchlim(RGBnew,contrastclip);
RGBcont = imadjust(RGBnew,inlim,[0 1],gammaadj);
% CLAHE to reduce some haze and help flatten
RGBahq = zeros(size(RGBnew));
for c = 1:3
RGBahq(:,:,c) = adapthisteq(RGBnew(:,:,c),'distribution','rayleigh');
end
% lincomb of contrast/gamma adjustment & CLAHE adjustment
RGBnew = RGBcont*alph + RGBahq*(1-alph);
imshow(RGBnew)
