Generate a nonuniform illumination (bias field) 2D image

Question

I am creating a nonuniform illumination in MATLAB R2015. The nonuniform illumination (known as bias field) is defined as

Nonuniform Illumination (intensity inhomogeneity) which manifests itself as slow intensity variations in the same class over the image domain.

The non-linear degree of intensity inhomogeneity is indicated by the range of values of the bias field in the interval [1 − α, 1 + α] with α > 0

I used the below MATLAB code to generate the intensity inhomogeneity image in case of α = 0.2:

cols=256;rows=256;
alpha=0.2;
x = linspace(1-alpha,1+alpha,cols);
y = linspace(1-alpha,1+alpha,row2);
% Transform into a grid
[X,Y] = meshgrid(x,y);
bias=X;
imshow(bias,[]) 

However, I am confused about the term 'non-linear'. It means that the value in the range [1 − α, 1 + α] needs to be slowly changed in a non-linear manner. In the above code, it seems that it is linearly changing. Is this correct? Could you suggest me a method to create an image based on the definition?

This is an expected example which I need to obtain

or

Answer 1

I provide an example for the generation of a nonuniform illumination using polynomials.

First, I take a vertical lineout (see first image) out of your last image, X being the pixel position and Y being the image intensity. For the horizontal direction I am using a simple parabola. Using kroneker product, the image is generated.

In order to ensure that the image intensity is between 1-alpha and 1+alpha, a streching is implemented as proposed by @Carl_Witthoft.

See code and resulting image below:

N=3; % degree of polynomial fit
alpha=0.2;

c = polyfit(X,Y,N);

y_f = zeros(size(Y));

for i=0:N
    y_f = y_f + c(N-i+1)*X.^i;
end

figure;
plot(X,Y, 'DisplayName', 'lineout'); hold on;
plot(X,y_f, 'DisplayName', 'fit'); legend show;

y_f = y_f./100;

xr = linspace(-4,4,150);
yr = 2-0.1*xr.^2;

img = kron(yr,y_f);

%% ensure that all values are between 1-\alpha and 1+\alpha

minV = min(img(:));
maxV = max(img(:));
mV = maxV-minV;

img = img-minV;                 % now from 0 to maxV-minV
img = img/(maxV-minV)*2*alpha;  % now from 0 to 2*alpha
img = img + 1-alpha;            % now from 1-alpha to 1+alpha

figure; imshow(img,[])

Answer 2

An usual approach is to use two-dimensional polynomials of low degrees (lets say up to 5). For your images, probably a polynomial of 3rd degree would be sufficient to simulate/estimate the bias field.

If your images happen to have more structure (i.e. magnitude images), it is feasable to smooth them strongly (e.g. Gaussian blur with a large FWHM), since those field variations live on a larger scale than anatomical variations.