# The correct way to linearize colorspace before resizing, blurring etc

There was an interesting discussion on MathGroups dealing with the fact that image-processing functions in Mathematica (and many other software, including Adobe Photoshop) work with RGB, Grayscale etc. intensity values as if they would be linear and additive while in fact these values are powers of the physical intensity values and consequently must be linearized before making additive operations on them. This topic is expanded in the linked article where examples of incorrect default image resizing and blurring are given along with the general explanation of the correct algorithm and images generated with it.

Matthias Odisio (Wolfram Research) replied:

This genuine issue will be properly addressed in a future release of Mathematica.

So, Mathematica 9 is released. But I cannot find any example in the Documentation on how to linearize a colorspace correctly. Obviously, the linearization algorithm must depend on the colorspace used.

The question is: is there efficient and straightforward way to linearize a colorspace in Mathematica correctly?

## UPDATE

Here is an excerpt from the documentation for Adobe After Effects which highlights some benefits of linear color space:

By performing operations in a linear color space, you can prevent certain edge and halo artifacts, such as the fringing that appears when high-contrast, saturated colors are blended together. Many color operations benefit from working in a linear color space, including those operations involved in image resampling, blending between layers with blending modes, motion blur, and anti-aliasing.

<...>

Note: A linearized working color space works best with higher color depths—16 bpc and 32 bpc—and is not recommended for 8-bpc color.

Also, good explanation of the difference between linear RGB and sRGB color spaces can be found here.

-
A very interesting question highlighting an issue of which I was only partially aware. –  Mr.Wizard Dec 3 '12 at 4:31
I've tested this with (industrial) cameras, grayvalues were practically linear. So this is probably not a property of the color space, but of the camera type and analog/digital preprocessing done in the camera. But if it isn't linear - why not just use ImageApply before resizing? –  nikie Dec 3 '12 at 7:42
@nikie ImageApply is an option but it is not straighforward because the used must define transformation functions for every colorspace by himself based on the colorspace specifications. It is not obvious and requires substantial work for anyone who is not deeply in image processing (like me). I also concerned about efficiency. –  Alexey Popkov Dec 3 '12 at 9:06
I don't have Mathematica 9, but the online documentation shows that ColorConvert has been updated to include CIE XYZ, and there is also a new symbol ColorProfileData for representing ICC color profiles. So I think you should be able to ColorConvert from whatever color space to XYZ, do the resizing and blurring, then ColorConvert back to the original color space (or any other you choose). –  Simon Woods Dec 3 '12 at 15:23
@Simon Woods Unfortunately, ImageResize gives wrong result after converting to the "XYZ" colospace: img=Import["http://www.4p8.com/eric.brasseur/gamma_dalai_lama_gray.jpg"];ImageR‌​esize[ColorConvert[img,"XYZ"],Scaled[1/2]]. From the other side, linearized sRGB colorspace produces expected result when ImageResize is applied (see code in the linked discussion). –  Alexey Popkov Dec 3 '12 at 20:23

Here is my attempt to figure out how the correct colorspace linearization should be made. I used specially designed test images by Eric Brasseur for comparison of two colorspace linearization algorithms. The first algorithm is just an implementation of the corresponding formulae from the Specification of sRGB made by Jari Paljakka who started the discussion on MathGroups. This algorithm does not take into acount the alpha channel (and probably will work incorrectly with it). The second algorithm utilizes new image processing functionality of Mathematica 9: the support of colorprofiles.

Both algorithms assume that the input image is in the sRGB colorspace which is the most commonly-used color space and also is the standard default color space for the Internet. More than 90% of all images in the Internet are sRGB-encoded.

The problem with the sRGB colorspace is that it does not have pure gamma curve and hence cannot be correctly linearized just by applying gamma to RGB values. But there is a pure gamma based colorspace: Adobe RGB 1998. So I decided to convert an image to Adobe RGB 1998, then linearize the colorspace by applying gamma using ImageAdjust, resize it with ImageResize (which operates under the assumption of linearity for the pixel data), then apply ImageAdjust with inverse gamma of Adobe RGB 1998 and finally convert from Adobe RGB 1998 to sRGB.

Here is a comparison of the results (these screenshots should be seen in original size with 100% resolution; the code follows):

The NASA image "Earth's City Lights" is a very extreme case where non-linear colorspace effects have a big impact on the results of resizing the image (reference from here):

(*Color profile based approach*)
(*The AdobeRGB1998.icc profile is from \
(*http://scarse.sourceforge.net/goodies/profiles/*)
sRGB = Import[
"http://scarse.sourceforge.net/goodies/profiles/spaces/sRGB.icm"];
sRGB2Linear[sRGBimg_Image] :=
ColorConvert[
Image[sRGBimg, "Real",(*ColorSpace->sRGB,*)Interleaving -> True],
sRGB -> adobeRGB1998], {0, 0, 563/256}]
linear2sRGB[linearRGB_Image] :=
ColorConvert[
Interleaving -> True], adobeRGB1998 -> sRGB]
linearResize[sRGBimg_Image, scaling_] :=
linear2sRGB[ImageResize[sRGB2Linear[sRGBimg], scaling]]
linearBlur[sRGBimg_Image, r_] :=
linear2sRGB[Blur[sRGB2Linear[sRGBimg], r]]
linear2Grayscale[sRGBimg_Image] :=
ColorConvert[sRGB2Linear[sRGBimg], "Grayscale"], {0, 0, 256/563}]

(*functional approach*)
srgb2linear =
Compile[{{Csrgb, _Real, 1}},
With[{α = 0.055},
Table[Piecewise[{{C/12.92,
C <= 0.04045}, {((C + α)/(1 + α))^2.4,
C > 0.04045}}], {C, Csrgb}]], RuntimeAttributes -> {Listable}];
linear2srgb =
Compile[{{Clinear, _Real, 1}},
With[{α = 0.055},
Table[Piecewise[{{12.92*C,
C <= 0.0031308}, {(1 + α)*C^(1/2.4) - α,
C > 0.0031308}}], {C, Clinear}]],
RuntimeAttributes -> {Listable}];
linearresize[sRGBimg_Image, scaling_] :=
Image[linear2srgb[
ImageData[
ImageResize[
Image[srgb2linear[
ImageData[ColorConvert[Image[sRGBimg, "Real"], "RGB"],
Interleaving -> True]], ColorSpace -> "RGB"], scaling],
Interleaving -> True]], ColorSpace -> "RGB"]
linearblur[sRGBimg_Image, r_] :=
Image[linear2srgb[
ImageData[
Blur[Image[
srgb2linear[
ImageData[ColorConvert[Image[sRGBimg, "Real"], "RGB"],
Interleaving -> True]], ColorSpace -> "RGB"], r],
Interleaving -> True]], ColorSpace -> "RGB"]
linear2grayscale[sRGBimg_Image] :=
Image[linear2srgb[
ImageData[
ColorConvert[
Image[srgb2linear[
ImageData[ColorConvert[Image[sRGBimg, "Real"], "RGB"],
Interleaving -> True]], ColorSpace -> "RGB"], "Grayscale"],
Interleaving -> True]]]

(*comparison*)
testImages =
Import /@ \
{"http://www.4p8.com/eric.brasseur/gamma_dalai_lama_gray.jpg",
"http://www.4p8.com/eric.brasseur/gamma_colors.jpg",
"http://www.4p8.com/eric.brasseur/gamma_2.2.jpg"};
testImages[[3]] = ImageRotate[testImages[[3]]];
correctResults =
Import /@ \
{"http://www.4p8.com/eric.brasseur/gamma_dalai_lama_gray_good.jpg",
"http://www.4p8.com/eric.brasseur/gamma_colors_good.jpg",
"http://www.4p8.com/eric.brasseur/gamma_good_2.2.jpg"};
correctResults[[3]] = ImageRotate[correctResults[[3]]];
fullSizeFix =
Style[Image[#, Magnification -> 1], Magnification -> 1] &;
correctResults = fullSizeFix /@ correctResults;
Grid[Prepend[
Join[Map[fullSizeFix, {ImageResize[#, Scaled[1/2]],
linearResize[#, Scaled[1/2]], linearresize[#, Scaled[1/2]]} & /@
testImages, {2}], List /@ correctResults, 2],
Style[#, Bold, FontFamily -> "Times",
TextAlignment -> Center] & /@ {"Just\nImageResize",
"ColorProfile-based\nlinearized resizing",
"Functional\nlinearized resizing",
"Expected result\n(Eric Brasseur)"}], Frame -> All]

calliphora =
Calliphora_sp_Portrait.jpg"];
calliphoraCorrect =
Import["http://www.4p8.com/eric.brasseur/gamma_calliphora_sp_\
portrait.1.jpg"];
calliphoraSize = ImageDimensions[calliphoraCorrect];
saturn = Import[
"http://photojournal.jpl.nasa.gov/jpegMod/PIA11667_modest.jpg"];
saturnCorrect =
Import["http://www.4p8.com/eric.brasseur/gamma_388135main_pia11667.\
1.jpg"];
saturnSize = ImageDimensions[saturnCorrect];
Grid[{Style[#, Bold, FontFamily -> "Times",
TextAlignment -> Center] & /@ {"Just\nImageResize",
"ColorProfile-based\nlinearized resizing",
"Functional\nlinearized resizing",
"Expected result\n(Eric Brasseur)"},
fullSizeFix /@ {ImageResize[calliphora, calliphoraSize],
linearResize[calliphora, calliphoraSize],
linearresize[calliphora, calliphoraSize], calliphoraCorrect},
fullSizeFix /@ {ImageResize[saturn, saturnSize],
linearResize[saturn, saturnSize],
linearresize[saturn, saturnSize], saturnCorrect}}]

Grid[{Style[#, Bold, FontFamily -> "Times",
TextAlignment -> Center] & /@ {"Just\nBlur[#,2]&",
"ColorProfile-based\nlinearized blur",
"Functional\nlinearized blur",
"Expected result\n(Eric Brasseur)"},
fullSizeFix /@ {Blur[testImages[[1]], 2],
linearBlur[testImages[[1]], 2], linearblur[testImages[[1]], 2],
Import["http://www.4p8.com/eric.brasseur/gamma_dalai_lama_gray_\
blurred.1.jpg"]}}]

gamma4 = Import[
"http://www.4p8.com/eric.brasseur/gamma_color_to_grayscale_start.\
jpg"];
Grid[{Style[#, Bold, FontFamily -> "Times",
TextAlignment ->
Center] & /@ {"Just\nColorConvert[#,\"Grayscale\"]&",
"ColorProfile-based\nlinearization", "Functional\nlinearization",
"Expected result\n(Eric Brasseur)"},
fullSizeFix /@ {ColorConvert[gamma4, "Grayscale"],
linear2Grayscale[gamma4], linear2grayscale[gamma4],
Import["http://www.4p8.com/eric.brasseur/gamma_color_to_grayscale_\
good.jpg"]}}]

earthLights =
Import["http://eoimages.gsfc.nasa.gov/ve//1438/earth_lights_4800.\
tif"];
fullSizeFix =
Style[Image[#, Magnification -> 1], Magnification -> 1] &;
Grid[{fullSizeFix /@ {ImageResize[earthLights, 500],
linearResize[earthLights, 500], linearresize[earthLights, 500]},
Style[#, Bold, FontFamily -> "Times",
TextAlignment -> Center] & /@ {"Just\nImageResize",
"ColorProfile-based\nlinearized resizing",
"Functional\nlinearized resizing"}
} // Transpose]

-
Everything looks good except the halftone patterns at the bottom of the first image. The two middle results do not match the target on the right. (+1) –  Mr.Wizard Dec 10 '12 at 13:57
@Mr.Wizard I suspect that the reason is that the Eric Brasseur's test images were lossy JPEG-encoded by him in order to get reasonable size of the web-page (he noticed this in the text). The "expected" resulting images probably were generated from the original lossless images. I cannot be completely sure but I think that lossy compression may cause the mentioned differences. –  Alexey Popkov Dec 10 '12 at 17:54
@Mr.Wizard But of course this does not explain differences in other images, especially for the Saturn image. Eric did not mention what algorithm he used for image resizing. It is possible that selecting other resizing method could give closer match... –  Alexey Popkov Dec 10 '12 at 18:01
The solution based on color profile is very fine indeed. The differences with the so-called "expected result" have no significance. –  Matthias Odisio Dec 18 '12 at 21:26