How to calculate the (affine or projective) distortion tensor for an RGB color image (without flickering) in frequency domain?

Matlab 2016b on Linux

Hello experts for Matlab,

I only deal with Matlab occasionally and image processing is my hobby. As a civil engineer I only have basic knowledge of linear algebra and need a tip from experts for my job.

The aim is to calculate the (affine or projective) distortion tensor for an RGB color image (three band image) with an interference signal on one channel. In the attached example, the interference is on the green channel (G) - here simply mirrored.

R-Band

G-Band (distorted channel)

B-Band

The hypothesis is:

1. The frequencies (energy) in the three channels to the RGB image are similar.
2. The frequencies (energy) in the three channels to the RGB image are similar.
3. In the distorted image, only the phase (position vector to the image pixel) changes compared to the other two channels.
4. In the distorted image, the frequencies (energy in the image) change only slightly compared to the other two channels.
5. The phase from the distorted image can be substituted by the phase of the two undistorted channels.
6. The cross multiplication of the frequency of the distorted band with the phase of the undistorted band leads to the back-transformation of the distorted image (matching).

In my example {code #1} the combination of the three channels (RGB bands) results in a color image with glittering signals.

Does anyone know how this flickering occurs and what the cause is?

Does anyone know the solution how to remove the flickering (flares)?

Can the bands (RGB) be normalized before a Fourier transformation so that the phases can be clearly mapped to the frequencies?

The second script {code #2} can be used to generate the three channel images (input to the first script).

I'm very grateful for clues. Also for hints on already asked questions in the forum. There's so much on this subject. I certainly haven't found all the important clues yet.

Chris

% MATLAB R2016B
% Test Tutorial - 11.05.2018 cn - Gaussian Phase-Shift Image-Matching -
% (c) Christian Naujoks, Eberbach, Germany  2018,
% c.naujoks[*]microanaglyph.com

%{
References:
1975 C. Kuglin, D. Hines - The phase correlation image alignment method

Matlab Example: openExample('images/UsePhaseCorrelationAsPreprocessingStepInRegistrationExample')
%}

% project location
% driveLocation = '/01.bench/20180904/fft_phase'; 
clear all; close all; imtool close all;

% 1 Read Image (RGB uint8)
ImageBand_R_undistorted = imread('ImageBand_R_undistorted.jpg');
ImageBand_G_distorted   = imread('ImageBand_G_distorted.jpg');
ImageBand_B_undistorted = imread('ImageBand_B_undistorted.jpg');

% 2 Transform Image-Domain (Time- or Space-Domain) to FFourier-Domain
FFkoef_R = fft2(ImageBand_R_undistorted);
FFkoef_x = fft2(ImageBand_G_distorted); 
FFkoef_B = fft2(ImageBand_B_undistorted);

      % >>>>>>> Reconstruct Image only on green Channel
      % 2.2 Manpulation Frequency- / Phase-Domain
      % compute magnitude and phase i.e. energy and position
      % Matching the distorted G-Band on Phase R or B
      % cross magnitudes and phases of the two images; convert to uint8 images
      % Reconstruct the Image direct invers from FFT / Frequency and Phase

      % >>>>>>> Reconstruct Image only on green Channel
      energy_R=abs(FFkoef_R); phase_R=atan2(imag(FFkoef_R),real(FFkoef_R)); position_R=exp(1i*phase_R);
      energy_x=abs(FFkoef_x); phase_x=atan2(imag(FFkoef_x),real(FFkoef_x)); position_x=exp(1i*phase_x); 
      energy_B=abs(FFkoef_B); phase_B=atan2(imag(FFkoef_B),real(FFkoef_R)); position_B=exp(1i*phase_B);
      % .*    Array multiplication. 
      %       A.*B is the element-by-element product of the arrays A and B. 
      %       A and B must have the same size, unless one of them is a scalar. 
      %ImageChannel_G_CrossMagPhas=(uint8((ifft2(time(energy_x,position_B)))));  
      ImageBand_G_CrossMagPhas2=(uint8((ifft2(energy_x.*position_B))));  

      % only a 1D FFT - not a 2D FFT:
      ImageBand_G_CrossMagPhas1=(uint8((ifft(abs(fft(ImageBand_G_distorted)).*(exp(1i*(atan2(imag(fft(ImageBand_B_undistorted)),real(fft(ImageBand_R_undistorted))))))))));  
      %ImageBand_G_CrossMagPhas2=(uint8((ifft2(abs(fft2(ImageBand_G_distorted)).*(exp(1i*(atan2(imag(fft2(ImageBand_B_undistorted)),real(fft2(ImageBand_R_undistorted))))))))));  
      ImageBand_G_CrossMagPhas3=(uint8((ifft(abs(fft(ImageBand_G_distorted)).*(exp(1i*(atan2(imag(fftn(ImageBand_B_undistorted)),real(fft(ImageBand_R_undistorted))))))))));  

% Summary RGB Band

figure('Position', get(0,'ScreenSize'), 'name','Summary Tutorial FFT Restore')
subplot(2,3,1);  imshow(ImageBand_R_undistorted); title('input: R-Band');
subplot(2,3,2);  imshow(ImageBand_G_distorted); title('input: G-Band (distorted)');
subplot(2,3,3);  imshow(ImageBand_B_undistorted); title('input: B-Band');
subplot(2,3,4);  imshow(cat(3, ImageBand_R_undistorted, ImageBand_G_CrossMagPhas1, ImageBand_B_undistorted)); title('1D-FFT: Output from Feq (G) & Phase (R or B)');
subplot(2,3,5);  imshow(cat(3, ImageBand_R_undistorted, ImageBand_G_CrossMagPhas2, ImageBand_B_undistorted)); title('2D-FFT: Output from Feq (G) & Phase (R or B)');
subplot(2,3,6);  imshow(cat(3, ImageBand_R_undistorted, ImageBand_G_CrossMagPhas3, ImageBand_B_undistorted)); title('nD-FFT: Output from Feq (G) & Phase (R or B)');

{code #2} - run first to produce the three input images (3 band) for the first script.

% MATLAB R2016B
% Test Tutorial - 11.05.2018 cn - Gaussian Phase-Shift Image-Matching -
% (c) Christian Naujoks, Eberbach, Germany  2018,
% c.naujoks[*]microanaglyph.com
% project location
% driveLocation = '/01.bench/20180904/fft_phase'; 
clear all; close all; imtool close all;

% 1 Read Image (RGB uint8)
%ImageIn = imread('lenna.jpg','jpg'); 
%ImageIn = imread('lighthouse.png');
%load mandrill;  Image = ind2rgb(X,map); 
Image = imread('kobi.png');

figure('name','Perform the crop operation per Mouse by double-clicking in the crop rectangle');
Image = imcrop( Image ); 
%Image = imcrop( Image, [0 0 8 8] ); % limited 8 x 8 pixel
delete(get(groot,'CurrentFigure'));
[Image_Rows, Image_Columns, numberOfColorBands] = size(Image); 

% 1.2 Split RGB Channel
% ---------------------------
ImageBand_R_undistorted = Image(:,:,1);
ImageBand_G_undistorted = Image(:,:,2);
    theta = 0; 
    % >>>>>>> Strech Image only green Band (distorted image)
    %ImageBand_G_distorted = ImageBand_G_undistorted
    %ImageBand_G_distorted = imresize(imtransform(ImageBand_G_undistorted, maketform('affine',[1 0 0; -.75 1 0; 0 0.2 1])), [Image_Rows Image_Columns]); 
    %ImageBand_G_distorted = imresize(imtransform(ImageBand_G_undistorted, maketform('affine',[-.97 0 0; .5 1 0; 0 0 1])), [Image_Rows Image_Columns]); 
    ImageBand_G_distorted = imresize(imtransform(ImageBand_G_undistorted, maketform('affine',[-1 0 0; 0 1 0; 0 0 1])), [Image_Rows Image_Columns]); 
    %ImageBand_G_distorted = imresize(imtransform(ImageBand_G_undistorted, maketform('affine',[-1 0 0; 0 -1 0; 0 0 1])), [Image_Rows Image_Columns]); % (rotation 180 Grad) this case  is nearly correct!
    %ImageBand_G_distorted = imresize(imwarp(ImageBand_G_undistorted, projective2d([cosd(theta) -sind(theta) 0.001; sind(theta) cosd(theta) 0.01; 0 0 1])), [Image_Rows Image_Columns]); 
    %ImageBand_G_distorted = imresize(imwarp(ImageBand_G_undistorted, projective2d([cosd(theta) -sind(theta) 0.01; sind(theta) cosd(theta) 0.0; 0 0 1])), [Image_Rows Image_Columns]);
  ImageBand_B_undistorted = Image(:,:,3);

      %%FFT Fourier Domain - Frequency & Phase
      % take FFT2 / IFFT2
      %freq     =  fft2(import);   phase     = angle(freq);   image = uint8(ifft2(angle(phase)));
      %freq_inv = ifft2(freq);     phase_inv = angle(phase);

      % 2.1 Transform Image-Domain (Time- or Space-Domain) to FFourier-Domain
      FFkoef_0 = fft2(rgb2gray(Image)); 
      FFkoef_R = fft2(ImageBand_R_undistorted);
      FFkoef_G = fft2(ImageBand_G_undistorted);
      FFkoef_x = fft2(ImageBand_G_distorted); 
      FFkoef_B = fft2(ImageBand_B_undistorted);
      FFkoef_1 = fft2(rgb2gray(cat(numberOfColorBands, ImageBand_R_undistorted, ImageBand_G_undistorted, ImageBand_B_undistorted)));
      FFkoef_2 = fft2(rgb2gray(cat(numberOfColorBands, ImageBand_G_distorted, ImageBand_G_undistorted, ImageBand_B_undistorted)));

2.2 Manpulation Frequency- / Phase-Domain compute magnitude and phase i.e. energy and position Matching the distorted G-Channel on Phase R or B cross magnitudes and phases of the two images; convert to uint8 images Reconstruct the Image direct invers from FFT / Frequency and Phase

      % >>>>>>> Reconstruct Image only on green Channel
      ImageBand_G_CrossMagPhas=(uint8(abs((ifft2((abs(FFkoef_x)).*(exp(1i*(atan2(imag(FFkoef_R),real(FFkoef_R))))))))));  %imshow(CrossMagPhax); pause; % here is something mystify !
   %   ImageBand_G_CrossMagPhas=(uint8(abs((ifft2((abs(FFkoef_x)).*(exp(1i*(atan2(imag(FFkoef_B),real(FFkoef_B))))))))));  %imshow(CrossMagPhax); pause; % here is something mystify !

      % comment next line with double %%to calculate in single step mode
      %{   
      energy_R=abs(FFkoef_R); phase_R=atan2(imag(FFkoef_R),real(FFkoef_R)); position_R=exp(1i*phase_R);
      energy_G=abs(FFkoef_G); phase_G=atan2(imag(FFkoef_G),real(FFkoef_G)); position_G=exp(1i*phase_G); 
      energy_x=abs(FFkoef_x); phase_x=atan2(imag(FFkoef_x),real(FFkoef_x)); position_x=exp(1i*phase_x); 
      energy_B=abs(FFkoef_B); phase_B=atan2(imag(FFkoef_B),real(FFkoef_R)); position_B=exp(1i*phase_B);
      % .*    Array multiplication. 
      %       A.*B is the element-by-element product of the arrays A and B. 
      %       A and B must have the same size, unless one of them is a scalar. 
      %ImageBand_G_CrossMagPhas=(uint8((ifft2(time(energy_x,position_B)))));  
      ImageBand_G_CrossMagPhas=(uint8((ifft2(energy_x.*position_B))));  
      %imshow(ImageBand_G_CrossMagPhas); pause; % here is something mystify !

      % output values
      %ImageBand_G_undistorted
      %ImageBand_G_distorted
      %FFkoef_x

      disp('>>>>>>>>>>>>:::::::::::  (Amplitude.*Position)');  (energy_x.*position_B)
  %    ifft2(xA.*BB)
      disp('Spatial image out undistorted Fequency and undistorted Phase:  uint8((ifft2(GA.*BB)))'); uint8((ifft2(energy_G.*position_B)))
      disp('Spatial image out distorted Fequency and undistorted Phase:    uint8((ifft2(xA.*BB)))'); uint8((ifft2(energy_x.*position_B)))

      disp('Energie (freq) out undistorted Fequency and undistorted Phase:  energy_GA=abs(FFkoef_G)'); energy_G
      disp('Energie (freq) out undistorted Fequency and undistorted Phase:  energy_BA=abs(FFkoef_B)'); energy_B
      disp('Energie (freq) out distorted Fequency and undistorted Phase:    energy_xA=abs(FFkoef_x)'); energy_x
      %}
      FFkoef_g = fft2(ImageBand_G_CrossMagPhas); %S_IG2 = abs(FFkoef_g);
      FFkoef_3 = fft2(rgb2gray(cat(numberOfColorBands, ImageBand_G_CrossMagPhas, ImageBand_G_undistorted, ImageBand_B_undistorted)));
% 3 .Summary RGB Channel
figure('Position', get(0,'ScreenSize'), 'name','Summary Tutorial FFT Transform')

subplot(5,9,1);  imshow(Image); title('Original Image (I0)'); 
subplot(5,9,2);  imshow(ImageBand_R_undistorted); title('R-Band'); 
subplot(5,9,3);  imshow(ImageBand_G_undistorted); title('G-Band original'); 
subplot(5,9,6);  imshow(ImageBand_B_undistorted); title('B-Band'); 
subplot(5,9,7);  imshow(cat(numberOfColorBands, ImageBand_R_undistorted, ImageBand_G_undistorted, ImageBand_B_undistorted)); title('RGB-Assemble');
subplot(5,9,8);  imagesc(abs(log(1+abs(fftshift(FFkoef_1))))); title('RGB-Frequency');
subplot(5,9,9);  imagesc(angle(fftshift(FFkoef_1))); title('RGB-Phase');

%subplot(5,9,10);  imshow(ImageInRGB); title('Gray Image (I0)');  
subplot(5,9,11); imshow(ImageBand_R_undistorted); 
subplot(5,9,13); imshow(ImageBand_G_distorted); title('G-Band distorted (+Noise)');
subplot(5,9,15); imshow(ImageBand_B_undistorted); 
subplot(5,9,16); imshow(cat(numberOfColorBands, ImageBand_R_undistorted, ImageBand_G_distorted, ImageBand_B_undistorted));
subplot(5,9,17); imagesc(abs(log(1+abs(fftshift(FFkoef_2)))));
subplot(5,9,18); imagesc(angle(fftshift(FFkoef_2)));

%subplot(5,9,15); imshow(ImOut0); title('Image Reconstruction from Phase');
subplot(5,9,20); imshow(ImageBand_R_undistorted); 
subplot(5,9,23); imshow(ImageBand_G_CrossMagPhas); title('G-Band matching on R- or B-Phase');
subplot(5,9,24); imshow(ImageBand_B_undistorted); 
subplot(5,9,25); imshow(cat(numberOfColorBands, ImageBand_R_undistorted, ImageBand_G_CrossMagPhas, ImageBand_B_undistorted)); title('Assemble with Gaussian-Matching - fail');
subplot(5,9,45); imshow(cat(numberOfColorBands, ImageBand_R_undistorted, ImageBand_G_CrossMagPhas, ImageBand_B_undistorted)); title('Output Gaussian-Matched-Image');
subplot(5,9,26); imagesc(abs(log(1+abs(fftshift(FFkoef_3)))));
subplot(5,9,27); imagesc(angle(fftshift(FFkoef_3)));

subplot(5,9,28); imagesc(abs(log(1+abs(fftshift(FFkoef_0))))); title('abs Powerspectrum / Freq'); 
subplot(5,9,29); imagesc(abs(log(1+abs(fftshift(FFkoef_R))))); title('- Band');
subplot(5,9,30); imagesc(abs(log(1+abs(fftshift(FFkoef_G))))); 
subplot(5,9,31); imagesc(abs(log(1+abs(fftshift(FFkoef_x))))); %imtool(abs(log(1+abs(fftshift(FFkoef_x)))));
subplot(5,9,32); imagesc(abs(log(1+abs(fftshift(FFkoef_g))))); 
subplot(5,9,33); imagesc(abs(log(1+abs(fftshift(FFkoef_B))))); 

subplot(5,9,37); imagesc(angle(fftshift(FFkoef_0))); title('Phase');
subplot(5,9,38); imagesc(angle(fftshift(FFkoef_R))); title('- Band');
subplot(5,9,39); imagesc(angle(fftshift(FFkoef_G))); 
subplot(5,9,40); imagesc(angle(fftshift(FFkoef_x))); %imtool(angle(fftshift(FFkoef_x))); 
subplot(5,9,41); imagesc(angle(fftshift(FFkoef_g))); 
subplot(5,9,42); imagesc(angle(fftshift(FFkoef_B))); 

imwrite(ImageBand_R_undistorted,'ImageBand_R_undistorted.jpg');
imwrite(ImageBand_G_distorted,'ImageBand_G_distorted.jpg');
imwrite(ImageBand_B_undistorted,'ImageBand_B_undistorted.jpg');
imwrite(cat(numberOfColorBands, ImageBand_R_undistorted, ImageBand_G_CrossMagPhas, ImageBand_B_undistorted),'ImageBand_catRGB.jpg');

Translated with www.DeepL.com/Translator