snibgo's ImageMagick pages

Chromatic aberration

... and correcting it.

• References
• The problem
• Fixing the problem
• Transverse correction by scaling
• Purple fringing
• Scripts

References

• Correcting Barrel Distortion, Helmut Dersch, 1999.
• Wikipedia: Chromatic aberration
• Cambridge in colour: camera lens corrections
• PanoTools: Chromatic aberration
• Hugin: Correcting Chromatic Aberrations
• Wikipedia: Purple fringing
• Lensfun at GitHub
• Hugin at Sourceforge

The problem

A single-element lens bends light (refraction), with the intention that rays from a distant object hitting different parts of the lens are focused on a single point on the sensor. The amount of bending depends on the colour of the light, and on the glass used (its index of refraction, IoR, which varies according to wavelength). Blue light bends more than red light. So if white light passes through a single-element lens the red, green and blue components will bend by different amounts, and hence will focus at different places. This is chromatic aberration (CA).

Crop from the left side of photo taken with single-element lens (185mm f/4.3), showing chromatic, comatic and other aberrations.
%IM%convert ^ %PICTLIB%\20180509\AGA_3631_sRGB.tiff ^ -strip ^ -crop 400x315+830+1738 +repage ^ cab_single.jpg

A different colour may focus at a different distance from the lens (axial CA, aka longitudinal CA) or a different distance from the centre of the sensor (transverse CA, aka lateral CA), or both. The distances depend on the lens aperture but also on the distance of the light source from the plane that is in focus.

Simulating chromatic aberration
Simulate axial CA: blur the red and blue channels %IMG7%magick ^ -size 400x300 ^ xc:Black ^ -fill White ^ -draw "circle 150,150 250,150" ^ -draw "circle 350,50 360,50" ^ -draw "circle 350,250 352,250" ^ -channel RGB -separate +channel ^ +write cab_sim.miff ^ ( -clone 0 ^ -blur 0x10 ^ ) ^ -swap 0,3 +delete ^ ( +clone ^ -blur 0x10 ^ ) ^ +swap +delete ^ -combine ^ cab_simaca.png
Simulate transverse CA: scale the red and blue channels %IMG7%magick ^ cab_sim.miff ^ ( -clone 0 ^ -distort SRT 2000,0,1.005,0 ^ ) ^ -swap 0,3 +delete ^ ( +clone ^ -distort SRT 2000,0,0.995,0 ^ ) ^ +swap +delete ^ -combine ^ cab_simtca.png
Simulate both axial and transverse CA %IMG7%magick ^ cab_sim.miff ^ ( -clone 0 ^ -blur 0x10 ^ -distort SRT 2000,0,1.005,0 ^ ) ^ -swap 0,3 +delete ^ ( +clone ^ -blur 0x10 ^ -distort SRT 2000,0,0.995,0 ^ ) ^ +swap +delete ^ -combine ^ cab_simatca.png

In film photography, CA is difficult to correct in post. In digital photography, transverse CA can be reduced by a geometrical distortion of the red, green and blue components of the image.

Camera lenses reduce CA by using multiple elements with different IoRs. But the problem can't be entirely removed. It is one of the many trade-offs that a lens designer has to make.

The problem shows as false colours at high-contrast edges. It is most noticable at edges furthest from the image centre.

Axial CA:

• is worse at long focal lengths;
• is reduced at small apertures;
• is difficult to correct;

Transverse CA:

• is worse at short focal lengths;
• is not affected by aperture;
• shows as red and green fringing on opposite sides of light patches on dark backgrounds (or dark patches on light background);
• is quite easy to correct;

set SRCNEF=%PICTLIB%20120918\DSC_0314.NEF

%DCRAW% -6 -T -w -O cab_1.tiff %SRCNEF%

set CaW=9
set CaH=9
set CaX=54
set CaY=4881

set sCROP=-crop %CaW%x%CaH%+%CaX%+%CaY% +repage

set sPROC=-strip %sCROP% -scale 4000%%

%IM%convert cab_1.tiff %sPROC% cab_1.png

Fixing the problem

But does it need fixing? On many images, the effect of CA is not noticable, and fixing it doesn't improve the image. Sometimes CA adds interest to an image. It can be treated like any photographic quality just like depth of field or slow shutter speeds or anything else. So leaving CA alone is always an option.

But sometimes we do want to fix it.

Possible solutions include:

1. Proprietary software from the lens manufacturer or elsewhere.
2. The Lensfun database or software or both, or software that uses Lensfun.
3. General-purpose software such as ImageMagick.

The manufacturer should be in the best position to know the characteristics of a lens, including tolerances, and is in the best position to know how to correct aberrations. But manufacturers don't publish these details, so we have a blunt choice of using their software or not.

Modelling transverse CA

As a first approximation, we model chromatic aberration as a geometric distortion of two channels with respect to the third. In theory this will fix transverse CA, but not axial CA. By convention we distort the red and blue channels to align with the green channel.

We assume the distortion that corrects the aberration is a radial movement towards or away from the centre. (More strictly, towards or away the optical axis, the location on the sensor where the lens axis intersects.) We assume circular symmetry, so the distortion is independent of the polar angle.

rd = f(ru)

... where r_u is the undistorted normalised radius, r_d is the distorted radius, and f() is a function that is zero at zero radius, and otherwise returns a value for r_d that is close to r_u. "Normalised" means the radius is one at the centre of the long sides. (Beware: other definitions of "normalised" can be used.) A polynomial is one such function:

rd = a*ru4 + b*ru3 + c*ru2 + d*ru

This is the same formula that is sometimes used to correct barrel/pincushion distortion. Effectively, we are correcting barrel/pincushion distortion in two channels with respect to the third.

If d is close to one and the other coefficients are close to zero, r_d will be close to r_u.

If we require that at r_u=1, there is no resizing so r_d=1 (the radius is unchanged at r==1) then:

a + b + c + d = 1

For overall barrel/pincushion distortion, some software will default d = 1-(a+b+c), and some software will always use this value without allowing the user to change it.

There will be one set of coefficients for the red channel, and a different set of coefficients for the green channel.

Using ImageMagick

We assume the only non-zero coefficient is d, so the transformation is:

rd = d*ru

This is a simple resizing.

Blah get the multipliers. Resize image with those parameters, or feed back to dcraw.

Use linear input?

The script chromAberr.bat finds scaling factors for the red and blue channels that make them most closely match the green channel within a defined crop area. It is fast and gives a good result for that area and other areas at a similar radius from the centre. But it might make CA worse at other radii.

call %PICTBAT%chromAberr ^
  cab_1.tiff ^
  %CaW% %CaH% %CaX% %CaY% ^
  cab_env

set cab_env 

cab_env_CaRad=0.98421
cab_env_ScaleB=1.000116101932579
cab_env_ScaleR=0.9998098897033573

The script caResize.bat applies these resizing factors to the red and green channels with ImageMagick's -distort SRT. We crop and enlarge to see the effect:

call %PICTBAT%caResize ^ cab_1.tiff ^ cab_corrim.tiff ^ %cab_env_ScaleR% %cab_env_ScaleB% %IMG7%magick ^ cab_corrim.tiff ^ %sPROC% ^ cab_corrim_sm.png

As expected, the result for this cropped area is good. We could extend the analysis to sample multiple areas to obtain a resizing that was optimised for the entire image. If we sampled four or more areas, we could use simultaneous equations to obtain the abcd parameters for a barrel/pincushion distortion of the red and blue channels.

We can directly use the two resizing factors as parameters to the dcraw -C option:

%DCRAW% ^ -6 -T -w ^ -C %cab_env_ScaleR% %cab_env_ScaleB% ^ -O cab_1d.tiff ^ %SRCNEF% %IM%convert cab_1d.tiff %sPROC% cab_1d.png

The result is similar to the caResize.bat result, but they are not identical. In dcraw.c, the function scale_colors() uses simple bilinear interpolation between four input pixels. IM's -distort SRT uses higher-quality EWA resampling.

Using tca_correct

The Hugin toolset includes tca_correct.exe that samples locations of an image to find the barrel/pincushion distortion to make the red and blue channels align with the green channel.

The Hugin equation for barrel/pincushion distortion works in the opposite direction to IM. [Is it???] The Hugin equation is:

rdest = a*rsrc4 + b*rsrc3 + c*rsrc2 + d*rsrc

Questions: where is r=1? Which way round are r parameters?

According to Lens correction model, r=1 for the largest circle that fits in the image, so r=1 at a distance half the smallest side from the centre.

In tca_correct.exe, the four parameters are named a, b, c and v. We can optimize all four of these...

for /F "usebackq tokens=*" %%L in (`%HUG%tca_correct ^
  -o abcv ^
  cab_1.tiff`) do set TCA_CORR1=%%L

echo TCA_CORR1=%TCA_CORR1% 

TCA_CORR1=-r 0.0003954:-0.0015797:0.0019806:0.9993112 -b -0.0003378:0.0011113:-0.0012434:1.0004323  

%HUG%fulla ^ %TCA_CORR1% ^ --dont-rescale ^ -o cab_corrtc1.tiff ^ cab_1.tiff %IMG7%magick ^ cab_corrtc1.tiff ^ %sPROC% ^ cab_corrtc1_sm.png

... or just the v parameter:

for /F "usebackq tokens=*" %%L in (`%HUG%tca_correct ^
  -o v ^
  cab_1.tiff`) do set TCA_CORR2=%%L

echo TCA_CORR2=%TCA_CORR2% 

TCA_CORR2=-r 0.0000000:0.0000000:0.0000000:1.0000817 -b 0.0000000:0.0000000:0.0000000:0.9999544  

%HUG%fulla ^ %TCA_CORR2% ^ --dont-rescale ^ -o cab_corrtc2.tiff ^ cab_1.tiff %IMG7%magick ^ cab_corrtc2.tiff ^ %sPROC% ^ cab_corrtc2_sm.png

Or just get the multipliers and feed back to dcraw. Or use in IM.

Transverse correction by scaling

Shift lens: guestimating the lens axis.

Purple fringing

Scripts

For convenience, .bat scripts are also available in a single zip file. See Zipped BAT files.

chromAberr.bat

rem From image %1
rem crop parameters %2 to %5
rem calculates chromatic aberration parameters by trial and error
rem setting environment variables prefixed with %6.

@if "%1"=="" findstr /B "rem @rem" %~f0 & exit /B 1

@setlocal enabledelayedexpansion

rem @call echoOffSave


set INFILE=%1

set CaW=%2
if "%CaW%"=="." set CaW=
if "%CaW%"=="" set CaW=w

set CaH=%3
if "%CaH%"=="." set CaH=
if "%CaH%"=="" set CaH=h

set CaX=%4
if "%CaX%"=="." set CaX=
if "%CaX%"=="" set CaX=0

set CaY=%5
if "%CaY%"=="." set CaY=
if "%CaY%"=="" set CaY=0

set ENVPREF=%6
if "%ENVPREF%"=="." set ENVPREF=
if "%ENVPREF%"=="" set ENVPREF=cab__

for /F "usebackq" %%L in (`%IM%identify ^
  -ping ^
  -format "WW=%%w\nHH=%%h\nHalfSht=%%[fx:min(w,h)/2]\nCaDX=%%[fx:(w/2-%CaX%-%CaW%/2)]\nCaDY=%%[fx:(h/2-%CaY%-%CaH%/2)]\n" ^
  %INFILE%`) do set %%L

for /F "usebackq" %%L in (`%IM%identify ^
  -ping ^
  -format "CaRad=%%[fx:hypot(2*%CaDX%/%WW%,2*%CaDY%/%HH%)]\n" ^
  xc:`) do set %%L

for /F "usebackq" %%L in (`%IM%identify ^
  -ping ^
  -format "CaRad=%%[fx:hypot(%CaDX%/%WW%,%CaDY%/%HH%)/%HalfSht%]\n" ^
  xc:`) do set %%L

for /F "usebackq" %%L in (`%IM%identify ^
  -ping ^
  -format "CaRad=%%[fx:2*hypot(%CaDX%,%CaDY%)/hypot(%WW%,%HH%)]\n" ^
  xc:`) do set %%L

for /F "usebackq" %%L in (`%IM%identify ^
  -ping ^
  -format "PorM=%%[fx:hypot(%CaW%,%CaH%)/hypot(%CaDX%,%CaDY%)/4]\n" ^
  xc:`) do set %%L

for /F "usebackq" %%L in (`%IM%identify ^
  -ping ^
  -format "F1=%%[fx:1-%PorM%]\nF2=%%[fx:1+%PorM%]\n" ^
  xc:`) do set %%L

set sCROP=-crop %CaW%x%CaH%+%CaX%+%CaY% +repage

set sPROC=-strip %sCROP% -scale 5000%%

set TMPPREF=\temp\ch_
set TMPEXT=miff

%IM%convert ^
  %INFILE% ^
  -strip %sCROP% -separate ^
  +adjoin ^
  %TMPPREF%_%%d.%TMPEXT%

echo WW=%WW% HH=%HH% CaDX=%CaDX% CaDY=%CaDY% CaRad=%CaRad% >cab_var.lis

rem set F1=0.994
rem set F2=1.010

call %PICTBAT%whatScaleT2 ^
  %TMPPREF%_0.%TMPEXT% ^
  %TMPPREF%_1.%TMPEXT% ^
  %F1% %F2% %CaDX% %CaDY%

if ERRORLEVEL 1 exit /B 1

set ScaleR=%wstSCALE%

call %PICTBAT%whatScaleT2 ^
  %TMPPREF%_2.%TMPEXT% ^
  %TMPPREF%_1.%TMPEXT% ^
  %F1% %F2% %CaDX% %CaDY%

if ERRORLEVEL 1 exit /B 1

set ScaleB=%wstSCALE%


call echoRestore

endlocal & set %ENVPREF%_CaRad=%CaRad%& set %ENVPREF%_ScaleR=%ScaleR%& set %ENVPREF%_ScaleB=%ScaleB%

caResize.bat

rem From image %1
rem makes output %2
rem %3 scale factor for red channel
rem %4 scale factor for blue channel


@if "%1"=="" findstr /B "rem @rem" %~f0 & exit /B 1

@setlocal enabledelayedexpansion

@call echoOffSave

call %PICTBAT%setInOut %1 car

if not "%2"=="" if not "%2"=="." set OUTFILE=%2

set ScaleR=%3
if "%ScaleR%"=="." set ScaleR=
if "%ScaleR%"=="" set ScaleR=1

set ScaleB=%4
if "%ScaleB%"=="." set ScaleB=
if "%ScaleB%"=="" set ScaleB=1

%IMG7%magick ^
  %INFILE% ^
  -channel RGB ^
  -separate ^
  +channel ^
  ( -clone 0 ^
    -distort SRT %ScaleR%,0 ) ^
  -swap 0,3 +delete ^
  ( -clone 2 ^
    -distort SRT %ScaleB%,0 ) ^
  +swap +delete ^
  -combine ^
  %OUTFILE%

call echoRestore

endlocal & set carOUTFILE=%OUTFILE%

whatScaleT2.bat

@rem Given same-size images %1 and %2,
@rem finds scale S in "-distort SRT S,0" for %1 to best match %2, and a score.
@rem
@rem Parameters:
@rem   %3 min scale (default 0.5)
@rem   %4 max scale (default 1.5)
@rem   %5 CX } Centre for scaling
@rem   %6 CY }

@rem Unlike whatScale.bat,
@rem this works by brute force trial-and-error of the entire images.
@rem So this may be more suitable when not much of either image is in focus.


@if "%2"=="" findstr /B "rem @rem" %~f0 & exit /B 1

@setlocal enabledelayedexpansion

@call echoOffSave


if "%wstMETRIC%"=="" set wstMETRIC=RMSE

set wstMIN_SCALE=%3
if "%wstMIN_SCALE%"=="." set wstMIN_SCALE=
if "%wstMIN_SCALE%"=="" set wstMIN_SCALE=0.5

set wstMAX_SCALE=%4
if "%wstMAX_SCALE%"=="." set wstMAX_SCALE=
if "%wstMAX_SCALE%"=="" set wstMAX_SCALE=1.5

if "%wstMIN_SCALE%"=="" exit /B 1

set CX=%5
if "%CX%"=="." set CX=
if "%CX%"=="" set CX=0

set CY=%6
if "%CY%"=="." set CY=
if "%CY%"=="" set CY=0

echo CX=%CX%  CY=%CY%

set PRECISION=%wstPRECISION%
if "%PRECISION%"=="" set PRECISION=-precision 6

set TMPEXT=.miff
set TMPDIR=%TEMP%

call %PICTBAT%setInOut %1 wst
set IN_A=%INFILE%
set IN_A_TMP=%TMPDIR%\%~n1_wst_inA%TMPEXT%
set IN_A_ONE=%TMPDIR%\%~n1_wst_oneA%TMPEXT%

call %PICTBAT%setInOut %2 wst
set IN_B=%INFILE%
set IN_B_TMP=%TMPDIR%\%~n2_wst_inB%TMPEXT%
set IN_B_ONE=%TMPDIR%\%~n2_wst_oneB%TMPEXT%

set DEBUG_FILE=wst_%~n1_%~n2_dbg%EXT%

echo IN_A=%IN_A%  IN_B=%IN_B%

if "%~x1"=="%TMPEXT%" (
  set IN_A_TMP=%IN_A%
) else (
  %IM%convert %IN_A% %IN_A_TMP%

  if ERRORLEVEL 1 exit /B 1
)

if "%~x2"=="%TMPEXT%" (
  set IN_B_TMP=%IN_B%
) else (
  %IM%convert %IN_B% %IN_B_TMP%

  if ERRORLEVEL 1 exit /B 1
)

for /F "usebackq" %%L ^
in (`%IM%identify ^
  -precision 16 ^
  -format "WW1=%%w\nHH1=%%h\nInvCX=%%[fx:-(%CX%)]\nInvCY=%%[fx:-(%CY%)]" ^
  %IN_A_TMP%`) ^
do set %%L

if "%WW1%"=="" exit /B 1

echo InvCX=%InvCX%  InvCY=%InvCY%

%IM%convert ^
  %IN_A_TMP% ^
  -distort SRT %CX%,%CY%,1,0 ^
  %IN_A_ONE%

%IM%convert ^
  %IN_B_TMP% ^
  -distort SRT %CX%,%CY%,1,0 ^
  %IN_B_ONE%

set SCALE0=%wstMIN_SCALE%
set SCALE1=0
set SCALE2=0
set SCALE3=%wstMAX_SCALE%

set A_SCALE=%TMPDIR%\%~n1_wst_A_S%TMPEXT%
set B_SCALE=%TMPDIR%\%~n2_wst_B_S%TMPEXT%


call :TryScale %SCALE0%
set COMP0=!COMP!

call :TryScale %SCALE3%
set COMP3=!COMP!

set nIter=0
set FINISHED=0


:loop

if %SCALE0%==%SCALE2% set FINISHED=1

for /F "usebackq" %%L ^
in (`%IM%identify %PRECISION% ^
  -precision 16 ^
  -format "SCALE1=%%[fx:%SCALE0%*pow(%SCALE3%/%SCALE0%,1/3)]\nSCALE2=%%[fx:%SCALE0%*pow(%SCALE3%/%SCALE0%,2/3)]" ^
  xc:`) ^
do set %%L

if %SCALE0%==%SCALE1% set FINISHED=1
if %SCALE1%==%SCALE2% set FINISHED=1
if %SCALE2%==%SCALE3% set FINISHED=1

call :TryScale %SCALE1%
set COMP1=!COMP!

call :TryScale %SCALE2%
set COMP2=!COMP!

echo Scale: %SCALE0% %SCALE1% %SCALE2% %SCALE3%
echo Comp: %COMP0% %COMP1% %COMP2% %COMP3%

for /F "usebackq" %%L ^
in (`%IM%identify ^
  -format "High0=%%[fx:%COMP0%>%COMP3%?1:0]\nHigh3=%%[fx:%COMP0%<%COMP3%?1:0]" ^
  xc:`) ^
do set %%L

if %High0%==1 (
  set SCALE0=%SCALE1%
  set COMP0=%COMP1%
) else if %High3%==1 (
  set SCALE3=%SCALE2%
  set COMP3=%COMP2%
) else (
  set FINISHED=1
)

set /A nIter+=1

if %FINISHED%==0 goto loop

echo nIter=%nIter%

call echoRestore

endlocal & set wstSCALE=%SCALE1%& set wstCOMP=%COMP1%

exit /B 0


rem -------------------------------------
rem Subroutines

:TryScale
set SC=%1

for /F "usebackq" %%L ^
in (`%IM%identify ^
  -precision 16 ^
  -format "Gtr1=%%[fx:%SC%>1?1:0]\nInvSc=%%[fx:1.0/%SC%]" ^
  xc:`) ^
do set %%L

if %Gtr1%==1 (

  %IM%convert ^
    %IN_A_TMP% ^
    -distort SRT %CX%,%CY%,%SC%,0 ^
    %A_SCALE%

  for /F "tokens=2 usebackq delims=() " %%L ^
in (`%IM%compare ^
  %A_SCALE% ^
  %IN_B_ONE% ^
  -metric %wstMETRIC% ^
  NULL: 2^>^&1`) ^
do set COMP=%%L

) else (
  %IM%convert ^
    %IN_B_TMP% ^
    -distort SRT %CX%,%CY%,%InvSc%,0 ^
    %B_SCALE%

  for /F "tokens=2 usebackq delims=() " %%L ^
in (`%IM%compare ^
  %IN_A_ONE% ^
  %B_SCALE% ^
  -metric %wstMETRIC% ^
  NULL: 2^>^&1`) ^
do set COMP=%%L

)
exit /B 0

tcaCorr.bat

rem From image %1,
rem makes output %2
rem using transverse chromatic aberration (TCA) parameters %3
rem as created by Hugin tca_correct.exe, probably quoted.

@if "%1"=="" findstr /B "rem @rem" %~f0 & exit /B 1

@setlocal enabledelayedexpansion

rem @call echoOffSave

call %PICTBAT%setInOut %1 tcac

if not "%2"=="" if not "%2"=="." set OUTFILE=%2

set PARAMS=%~3

for /F "tokens=1-10 delims=: " %%A in ("%PARAMS%") do (
  if not "%%A"=="-r" exit /B 1
  if not "%%F"=="-b" exit /B 1
  set ar=%%B
  set br=%%C
  set cr=%%D
  set dr=%%E
  set ab=%%G
  set bb=%%H
  set cb=%%I
  set db=%%J
)

for /F "usebackq" %%L in (`%IM%identify ^
  -format "sumr=%%[fx:%ar%+(%br%)+(%cr%)+(%dr%)]\nsumb=%%[fx:%ab%+(%bb%)+(%cb%)+(%db%)]\n" ^
  xc:`) do set %%L

echo sumr=%sumr% sumb=%sumb%

%IMG7%magick ^
  %INFILE% ^
  -channel RGB ^
  -separate ^
  +channel ^
  -verbose ^
  ( -clone 0 ^
    -distort BarrelInverse %ar%,%br%,%cr%,%dr% ) ^
  -swap 0,3 +delete ^
  ( -clone 2 ^
    -distort BarrelInverse %ab%,%bb%,%cb%,%db% ) ^
  +verbose ^
  +swap +delete ^
  -combine ^
  %OUTFILE%

call echoRestore

endlocal & set tcacOUTFILE=%OUTFILE%

%IM%identify -version

Version: ImageMagick 6.9.5-3 Q16 x86 2016-07-22 http://www.imagemagick.org
Copyright: Copyright (C) 1999-2015 ImageMagick Studio LLC
License: http://www.imagemagick.org/script/license.php
Visual C++: 180040629
Features: Cipher DPC Modules OpenMP 
Delegates (built-in): bzlib cairo flif freetype jng jp2 jpeg lcms lqr openexr pangocairo png ps rsvg tiff webp xml zlib

%IMG7%magick identify -version

Version: ImageMagick 7.0.7-28 Q16 x64 2018-03-25 http://www.imagemagick.org
Copyright: Copyright (C) 1999-2015 ImageMagick Studio LLC
License: http://www.imagemagick.org/script/license.php
Visual C++: 180040629
Features: Cipher DPC HDRI Modules OpenMP 
Delegates (built-in): bzlib cairo flif freetype gslib heic jng jp2 jpeg lcms lqr openexr pangocairo png ps raw rsvg tiff webp xml zlib

%HUG%tca_correct -h | findstr version 

tca_correct version 2018.0.0.5abfb4de7961 built by Thomas

Page version v1.0 14-May-2018.

Page created 18-May-2018 10:08:53.

Copyright © 2018 Alan Gibson.