How to degrade red channel in underwater images to artificially increase color distortion?

I've got some good quality underwater images that have small negligible blurring and haze. Therefore, I want to use them as reference image and distort their red channel. By doing this, I can numerally (instead of visually) measure error after applying my restoring methods. For example, I want to numerally measure the effect of increasing all red intensities (by a small constant value) on error between the initial image and the restored image (since I have the reference image). On the other hand, error could be simply measured by mean square error or taking norm.

I've written the following code which uses medium transmission to degrade the red channel. Medium transmission is computed by a variant of dark channel prior proposed in " Automatic Red-Channel Underwater Image Restoration ".

Unfortunately, my code is too slow. Is there any better way to increase speed efficiency?

Use this script to run the mediumtransmissionMat() function. Assign the path to your image into fileNameDataSet.

Notice: the parameter gsdestruction is set to 3 but feel free to change it with any greater odd number.

im=imread(fileNameDataSet);
im2 = im;
%%%%%% destructing image
gsdestruction = 3;
[medtransMat , globalBackgLight]  =  mediumtransmissionMat ( im , gsdestruction ) ;
im2( : , : , 1) = im2uint8( im2double( im( : , : , 1) ) .* medtransMat - ( globalBackgLight(1) ) * ( 1 - medtransMat) );

Definition of mediumtransmissionMat:

function [medtransMat , globalBackgLight]  =  mediumtransmissionMat ( im , gs ) 
    %%%%%% gs must be an odd num
    im = im2double(im);
    imheight = size ( im , 1 ) ;
    imwidth = size ( im , 2 ) ;
    %% find brightest pixel in the dark channel- red
    maxvalred  =  max  ( max ( im (  :  ,  :  ,  1  )   )   ) ;
    indxmaxvalred = find ( im ( : , : , 1 )  == maxvalred ) ;
    %% green global background light- scalar
    greenglobalBackLight = im ( indxmaxvalred ( 1 ) +imheight*imwidth ) ;
    %% blue global background light- scalar
    blueglobalBackLight = im ( indxmaxvalred ( 1 ) +2*imheight*imwidth ) ;
    globalBackgLight = double ( [maxvalred , greenglobalBackLight , blueglobalBackLight] );
    %%  medium transmission size is as im
    mediumtransmission = zeros ( imheight , imwidth ) ;
    %%  gs is grid size
    for r = 1:imheight
        for c = 1:imwidth
    %%      four corners of rect
          rmin = r- ( gs-1 ) /2;
          rmax = r+ ( gs-1 ) /2;
          cmin = c- ( gs-1 ) /2;
          cmax = c+ ( gs-1 ) /2;
    %%      if rmin is out of boundary
          if  ( rmin<1 ) 
            rmin = 1;
            rmax = rmin + gs - 1;
          endif
    %%      if cmin is out of boundary
          if  ( cmin<1 ) 
            cmin = 1;
            cmax = cmin+gs - 1;
          endif
    %%      if rmax is out of boundary
          if  ( rmax>imheight ) 
            rmax = imheight;
            rmin = rmax - gs - 1;
          endif
    %%      if cmax is out of boundary
          if  ( cmax>imwidth ) 
            cmax = imwidth;
            cmin = cmax - gs - 1;
          endif    
    %%    get gs by gs patch from image
          patchImage =  im ( rmin:rmax , cmin:cmax , 1:3 );
    %%    find the min of each patch- scalar
          minpatchred = min ( min ( 1 - patchImage ( : , : , 1 ) ) );
          minpatchgreen = min ( min ( patchImage ( : , : , 2 )  )  ) ;
          minpatchblue = min ( min ( patchImage ( : , : , 3 )  )  ) ;
    %%    normalize green and blue by local back light
          patchImage ( : , : , 1 ) = patchImage ( : , : , 1 ) / ( 1 - maxvalred );
          patchImage ( : , : , 2 )  = patchImage ( : , : , 2 ) /greenglobalBackLight;
          patchImage ( : , : , 3 )  = patchImage ( : , : , 3 ) /blueglobalBackLight;
    %%    min of green and blue patches-scalar
          minpatch = min ( minpatchgreen , minpatchblue ) ;
          minpatch = min ( minpatch , minpatchred );      
    %%    medium transmission- scalar
          mediumtransmission ( r , c )  = 1 - minpatch;     
        endfor % for c    
    endfor  % for r
    medtransMat = mediumtransmission;
    end  % end of function

Answer 1

[vahid rowghanian]

I googled on how to increase the speed of scanning images and I found circshift() function which speeds up local neighbour processing. this function also shrank the code lines a lot.

so, here is the medium transmission calculation:

function medtransMat =  mediumtransmissionMat ( im , gs ) 
%%%%%% gs must be an odd num
%%% im must be normalize [0-1]


half=floor(gs*gs/2);
  %%%%% add padding to im
paddedim=padarray(im,[half,half],'both');
immin =ones(size(paddedim));    
imheight = size ( im , 1 ) ;
imwidth = size ( im , 2 ) ;

%% find brightest pixel in the dark channel- red
maxvalred=max(max(im(:,:,1))) ;
indxmaxvalred = find ( im ( : , : , 1 )  == maxvalred ) ;
%% green global background light- scalar
greenglobalBackLight = im ( indxmaxvalred ( 1 ) +imheight*imwidth ) ;
%% blue global background light- scalar
blueglobalBackLight = im ( indxmaxvalred ( 1 ) +2*imheight*imwidth ) ;
globalBackgLight = double ( [maxvalred , greenglobalBackLight , blueglobalBackLight] );

%%%% loop in green and blue channels to find min
for k=2:3  
  for i=-half:half
    for j=-half:half
      immin(:,:,k)=min(immin(:,:,k) , circshift(paddedim(:,:,k),[i,j]));
    endfor
  endfor
endfor

%%%% loop in red channel only due to different equation
for i=-half:half
  for j=-half:half
    immin(:,:,1)=min( immin(:,:,1) , ones(size(paddedim(:,:,1)))-circshift(paddedim(:,:,1),[i,j]));
  endfor
endfor
%%%% normalize 0-1
immin(:,:,1)=immin(:,:,1)/(1-globalBackgLight(1));
immin(:,:,2)=immin(:,:,2)/globalBackgLight(2);
immin(:,:,3)=immin(:,:,3)/globalBackgLight(3);

%%% remove the padding
imminnopad=immin(half+1:end-half,half+1:end-half,:);
%%% find min in 3 channels
medtransMat=min(imminnopad,[],3);
medtransMat=medtransMat/max(medtransMat(:));
medtransMat=1-medtransMat;


end % end of function