Computational Vision
U. Minn. Psy 5036
Daniel Kersten
Lecture 10: Image processing

Initialize

Read in Add-in packages:

Off[General::"spell1"];
SetOptions[ArrayPlot,ColorFunction->"GrayTones",DataReversed->True,Frame->False,AspectRatio->Automatic,Mesh->False,PixelConstrained->True,ImageSize->Small];

The input 64x64 image: face

face = ImageData[ImageReflect[ColorConvert[,"Grayscale"]]];
{width,height}=Dimensions[face];
gface=ArrayPlot[face]

Upcoming dates:

Final project outlines due: Nov. 13th.

Final projects

Format

Should be written like a scientific paper. 
Might require most of the code to be put in appendices. 
Can use modules you find elsewhere, but preserve copyrights, and reference
Your "audience" will be your class peers.

Review: Multiresolution, spatial filter models of early visual processing

Single-channel spatial filtering

One assumes a convolution filter of a fixed shape (e.g. the weights of a DOG filter) that gets applied across the whole image to predict an output "neural image".

images can be represented in the fourier domain in terms of: amplitude and phase spectra.

convolution of an image with a filter in the space domain is mathematically equivalent to: 

	image⊗kernel = (F^(-1))( F (image)× F(kernel) )
or using Mathematica functions:
	InverseFourier[Fourier[filter]*Fourier[image]]

convolution is used to model: blur, neural “images”, the linear stage in neural networks, deep convolutional neural networks

Multiple spatial frequency channels

Most neurons in primary visual cortex, in particular “simple cells” and “complex cells”,  have receptive fields whose sizes are limited, and are selective for particular spatial frequencies and orientations. This suggests that neurons in primary visual cortex can be grouped into “channels” consisting of population of neurons across retinotopic locations, but with similar frequency and orientation tuning. (See classic experiment of Campbell and Robson, 1968)

To model how these simple cell population transform a pattern of incoming image intensities, one assumes a set of convolution filters, each set has a fixed shape (e.g. the weights of gabor filters describe a fixed template) that get applied across the  image to predict a set of output "neural images". Each neural image represents information at a particular spatial scale and orientation. The collection of filters, identical except for position, is called a channel.

Global vs. local: 
	Sinusoidal basis functions are global filters. Global filtering is not a good model because they give up spatial localization, at the expense of spatial frequency information (the filters would have to extend across the whole space implying large receptive fields).

Multiresolution analysis with local filters

Instead of (global) sinusoidal basis functions, the convolutions use localized "gabor function" filters:

Clear[Grating,GratingPatch,kern];

Grating[x_,y_,fx_,fy_,phase_] := Cos[(2.0 Pi (fx x + fy y) + phase)];
GratingPatch[x_,y_,fx_,fy_,sig_,phase_] := Exp[-((x)^2 + (y)^2)/(2*sig^2)]*Grating[x,y,fx,fy,phase];
kern[fx_, fy_, sig_,phase_] := 
  Table[GratingPatch[x, y, fx, fy, sig,phase], {x, -1, 1, .1}, {y, -1, 1, .1}];

The following demo illustrates a single neural image for various choices of spatial frequency and orientation. In addition, the envelope width manipulates "bandwidth"--i.e. more cycles under the gabor means narrow spatial frequency tuning and thus narrow bandwidth. The phase slider moves one from "bar" to "edge" type receptive fields.

Manipulate[
GraphicsRow[{
ArrayPlot[kern[fr*Cos[theta], fr*Sin[theta],sig,phase]],ArrayPlot[ListConvolve[kern[fr*Cos[theta], fr*Sin[theta],sig,phase],face]]}],{{fr,1,"radial frequency"},.1,2},{{theta,.4,"orientation"},0,Pi},{{sig,.4,"envelope width"},.001,1},{{phase,0,"phase"},.0,Pi/2}]

zebra=ImageData[ImageReflect[ColorConvert[RemoveAlphaChannel[],"Grayscale"]]];
Dimensions[zebra]

Manipulate[
GraphicsRow[{
ArrayPlot[zebra],ArrayPlot[ListConvolve[kern[fr*Cos[theta], fr*Sin[theta],sig,phase],zebra]]}],{{fr,1,"radial frequency"},.1,2},{{theta,.4,"orientation"},0,Pi},{{sig,.4,"envelope width"},.001,1},{{phase,0,"phase"},.0,Pi/2}]

Self-similarity

When the filters have the same shape except for a change of scale (x→ α x), where α is a constant, they are called self-similar. The self - similar idea is important to vision because of the importance of scale sensitive and scale invariant processing. Further, the self - similar aspect of  neural filter models bears a close resemblance to the emerging mathematical field of wavelet analysis.

Human statistical efficiency for detecting gabor patches

Burgess, Wagner, Jennings and Barlow (1981) combined the SKE observer and spatial frequency analysis of human vision to find out which patterns human observers had the best constrast discrimination efficiency. They showed in a 1981 Science article that narrowly windowed sinusoids were detected with high efficiency (>70%) when added to static visual noise. These targets were detected more efficiently than simple disks of light.

You have all the tools to replicate the experiment of Burgess et al. You can compute d' for the ideal observer for signal-known-exactly patterns. And you can generate Gaussian-windowed sinusoids and add them to gaussian white noise. If you measure the percent correct, and convert that to d' for the human observer, you can calculate the absolute efficiency for human detection--and contribute to answer the question of what the eye sees best. Watson, Barlow & Robson (1983) found that that a 7 c/deg grating drifting at 4 Hz, (with a narrow gaussian envelope in space and time) was detected more efficiently than other patterns. They did not add artificial noise. They used photopic light levels, and taking into account the photon noise, they calculated the quantum efficiency at less than <0.05%.

Bottom-line: image coding in terms of scale and orientation:
	A model for human spatial image representation

At each spatial location, project the image onto a collection of basis vectors (i.e. compute the dot product) that span a range of spatial scales and orientations:

In general, these neural models of basis functions may be over-complete, and non-orthogonal. And there may be a range of phases. Above we show only the "sine-phase" or "edge-detectors" of Hubel and Wiesel.

How linear are V1 neurons?

There are kinds of non-linearities that have been introduced to refine, or model other types of neural populations: 1) rectification--neuronal firing rate is by definition non-negative. The combination of on-center, off-center responses can be treated as a theoretical unit to represent negative and positive signal values. 2) adding squared outputs of sine-phase and cosine-phase filter to produce  “contrast energy” filters; 3) Contrast normalization, where the output of an otherwise linear neuron is normalized by the “energy” outputs of ones nearby in space. We’ll discuss some of these in later lectures

What is a multiresolution scale/orientation representation good for?

What is the computational significance of a wavelet-like decomposition?

Analysis of what vision needs to recognize objects, etc..
	-> First stages of local feature extraction
		Edges, bars, ...over various spatial scales
	-> “front end” for subsequent processing 
		textures, curvature, corners, ...

Efficient coding?
	-> savings in neurons, or metabolic requirements?
	-> representations for efficient learning, subsequent coding?
	-> analysis of natural image statistics

Image processing tools

Mathematica  has a library of built-in functions for doing image manipulations

See the Basic Image Manipulation, Image Processing & Analysis, and the Image Filtering & Neighborhood Processing  guides. 
Here are some examples:

Image[Reverse[face]]

We’ve seen a  simple blur using BoxMatrix:

ImageConvolve[,BoxMatrix[2]/9]// ImageAdjust ;

BoxMatrix[2]/9

And the built-in Laplacian of a Gaussian filter which can be used to model center-surround, lateral inhibitory filtering:

ifiltered = ImageAdjust[LaplacianGaussianFilter[, 2] ]

ImageAdjust[ifiltered]

We can quickly generate the “neural images” for on- and off-center responses:

GraphicsRow[{ifiltered,ImageMultiply[Binarize[ifiltered ],ifiltered],ImageMultiply[ColorNegate[Binarize[ifiltered ]],ifiltered]}]

Point operations

Various contrast definitions

We’ve seen that for simple stimuli, contrast can be defined as:
	(Imax - Imin)/(Imax + Imin): Called "Michelson contrast". Particularly appropriate for gratings, or stimuli with primary luminance peak and valleys in the image.
	ΔI/(Subscript[I, background]): Often used in psychophysics of small points/disks against a larger background, or temporal increments/decrements (ΔI) relative to a base level ((Subscript[I, background])).
	ΔI/(Subscript[I, mean]): Gives same number as Michelson for gratings, when ΔI corresponds with the amplitude of the grating.

For complex stimuli, we could represent contrast in terms of standard statistical measures on luminance. Let the mean luminance: (Subscript[I, mean]) = (∑_(x,y)) I(x,y)/N, where N is the number of pixels. The variance of the luminance is:  (((∑_(x,y)) (((I(x,y)-(Subscript[I, mean])))^2))/N). Then the standard deviation could be used to provide an overall measure of contrast:

	(Sqrt[ (((∑_(x,y)) (((I(x,y)-(Subscript[I, mean])))^2))/N)])
(The formula for the unbiased variance estimate could also used, where N is replaced by N-1).
However, the above standard deviation definition depends on the units (e.g. luminance is measured in candelas/meter^2). Further, it is useful to have contrast definitions that come closer to capturing the perceptual aspects of contrast in an image.

This is the rationale for defining contrast at a point--called local contrast. The visual system is relatively insensitive to the mean luminance, suggesting that a useful measure of local contrast is luminance divided by the mean, whose ratio is dimensionless. In addition, the subjective, qualitative difference between “bright” and “dark” suggests positive and negative contrast, respectively. Local contrast at a point (x,y) as:

c(x,y) =  ((I(x,y) - (Subscript[I, mean]))/(Subscript[I, mean]))

The mean of c(x,y) is zero by definition, with positive and negative values corresponding to bright and dark relative to the mean. (Note that local contrast can be defined as function of time too, c(x,y,t).) 

The definition of local contrast at a point raises the question: over what range should (Subscript[I, mean]) be calculated? This question is relevant to the problems of local adaption to light level, and to tone mapping. The default for us will be to take the mean over the whole image.

Given c(x,y), we can calculate a summary measure of contrast. The variance of c(x,y) is the same as the variance of I(x,y)/(Subscript[I, mean]) and is given by: (((∑_(x,y)) (c^2)(x,y))/N). The root mean square (r.m.s.) of set of measurements (Subscript[m, i]^2) is by definition: (Sqrt[(((∑_i) (Subscript[m, i]^2))/N)]). Thus one can define r.m.s. contrast as the square root of the variance of c:

	(Sqrt[(((∑_(x,y)) (c^2)(x,y))/N)])

This definition of r.m.s. contrast provides us with a useful summary measure of overall "contrastiness" for a complex image. 

By analogy with physics  (substituting space for time) the square of r.m.s contrast is sometimes called "contrast power". And then "contrast energy" is defined as: contrast power x area. In psychophysics, "area" is often measured in squared degrees of visual angle. Contrast energy is a useful measure when one is concerned about how the size (or duration) of an image patch affects its visibility. If x and y are measured in degrees of visual angle, then, ∫(c^2)(x,y)dxdy is the contrast energy, which in the discrete case can be estimated as:

	(∑_(x,y)) (c^2)(x,y) x ((area/N)) = (∑_(x,y)) (c^2)(x,y) ΔxΔy, where Δx and Δy are the dimensions of a pixel in degrees.

Recall from Lecture 5: d'=((r.m.s. contrast x (Sqrt[N]))/σ)

If we extend the above contrast definitions to time, in terms of general appearance, the apparent contrast of an image doesn't increase with time, so contrast energy (contrast power x duration) is irrelevant. However, over short durations (less than 100 msec or so), human vision approximately integrates contrast power, so contrast energy is a better predictor of threshold. But for longer durations, contrast power is the better predictor.

Side note on the relation of contrast energy of an image to its fourier representation: Parseval's theorem says that the contrast energy of an image is equal to the integral of the square of the amplitude spectrum (called the power spectrum) over spatial frequency. Contrast energy doesn't depend on the phase relationships between the frequency components.

Contrast manipulations

Adjusting contrast (gain=1 leaves image unchanged, gain=0 reduces it to a uniform field):

μ=Mean[Flatten[face]];gain=0.045;
Manipulate[
ArrayPlot[gain (face-μ)+μ,Mesh->False,Frame->False,PlotRange->{Min[face],Max[face]}],
{{gain,0.045},0,1}]

Psychophysics and contrast

When measuring human visual sensitivity, it is important to carefully measure and calibrate the image stimuli. Because standard 8-bit computer displays resolve 256 graylevels, it can be useful to convert the stimuli into a range going from 0 to 255.  Scale so values are represented as graylevels between 0 and 255:

α = 255/(Max[face] - Min[face]);
β = -α  Min[face];

face256 = α face + β;

Gamma correction for displays

Computer operating systems allow one to adjust for various non-linearities between displays. A typical screen has a non-linear relationship between measured screen intensity and the voltage supplied. A traditional way of summarizing the non-linear relationship is in terms of "gamma": intensity = a x voltage^gamma. Let's assume that we are using intensity units that range from 0 to 255 and voltage units also going from 0 to 255.  intensity = 255^(1-gamma) x voltage^gamma:

output[input_,gamma_]:=(255^(1-gamma)) (input^gamma);
Plot[output[x,2],{x,0,255},Frame->True,ImageSize->Small]

The computer’s display card has a look-up-table (LUT) that can be loaded with the inverse gamma function to linearize the display. The LUT remaps input values to output voltages  so that there is a linear relationship between internal values and screen luminance.

Solve[output == 255^(1-gamma) input^gamma,input]

inverse[output_,gamma_]:=((((255^(-1+gamma)) output))^(1/gamma));
Plot[inverse[x,2],{x,0,255},Frame->True,ImageSize->Small]

One can do silly things too, like a many-to-one mapping:

(output2[input_]:=127. Sin[(input/8)]+127;);
Plot[output2[x],{x,0,255},Frame->True,ImageSize->Small]

ArrayPlot[output2[face256],PlotRange->{0,255}]

..but maybe this isn’t so silly. What does the fact that you can still see a recognizable face-like form in the above picture tell you about human vision?

Using gamma to do  point operations

You can also use the gamma transformation to do non-linear point operations on an image to concentrate some intensity values more than others.

Manipulate[
GraphicsRow[{Plot[output[x,gamma],{x,0,255},Frame->True,ImageSize->Small,PlotRange->{0,255}],ArrayPlot[output[face256,gamma],PlotRange->{0,255}]}],
{gamma,0.1,4}]

Sigmoidal contrast manipulation

Here is a gain function (called the "logistic function") that manipulates contrast smoothly--a "soft" threshold:

squash[x_,μ_,γ_]:=N[(1/(1+(E^(-γ (x-μ)))))];Plot[squash[x,128,1],{x,0,255},Frame->True,ImageSize->Small]

gain=0.045`;μ0=Mean[Flatten[face256]];
Manipulate[GraphicsRow[{Plot[squash[x,μ,γ],{x,0,255},Frame->True,PlotRange->{{0,255},{0,1}}],ArrayPlot[squash[(face256-μ0)+μ0,μ,γ],PlotRange->{Min[face],Max[face]}]}],{{μ,128},0,255},{{γ,.5},0,1}]

As γ grows, the sigmoid approaches a hard-threshold. Images that are forced to have only two values are called "Mooney images". 
We can make Mooney images more directly using  a function that takes an image and sets pixels bigger than  τ to 255, and if less than (or equal to) τ,  to 0:

Mooney[image_,τ_]:=Map[If[#>τ,255,0]&,image,{2}];
ImageReflect[Image[Mooney[face256,32]]]

See also: Binarize[], ColorQuantize[].  And the famous “Dalmation dog illusion” in perception books, or via a web search.

See Moore and Engel, 2001, and  Hegdé J &  Kersten D. (2007)  for applications of Mooney images to studying human vision. 

Gray levels are effectively quantized at low light levels. Mooney images mimic this effect but at high light levels.

Simple statistics

First-order, i.e. don't take into account relations between pixels. We’ve already seen a summary measure of overall “contrastiness”.

μ=Mean[Flatten[face]];
σ=Sqrt[Variance[Flatten[face]]];
rmscontrast = Sqrt[Variance[Flatten[face]]]/Mean[Flatten[face]]

Histograms

For images or for any data:

{ImageHistogram[],Histogram[Flatten[face]]}

You can tell that the image is quantized at a coarse level (less than 4 bits).  Alternatively, you could calculate the histogram with more basic functions. To do the pattern match below, the floating point numbers are first converted to integers using Round[ ]. If we normalize the histogram so that the sum is one, then we have a probability:

domain=Range[0,255];
Freq=Map[Count[Round[Flatten[face256]],#]&,domain];

ListPlot[(Freq/(Plus@@Freq)),PlotStyle->PointSize[0.02],PlotRange->{0,.3},ImageSize->Tiny]

Getting regions of images

Three ways to pull out a region:

ArrayPlot[Take[face256,{1,32},{1,64}]]

Or  to get rows 1 to 32 and columns 1 to 64 you can use:

ArrayPlot[face256[[1;;32,1;;64]]]

The above method is particularly useful to learn because it generalizes to general lists of lists, and has parallels in both python/numpy and matlab (see “slicing” in python).

We can also use the built-in image function ImageTake[] to get rows 32 through 64:

ImageTake[,{32,64}]

Getting coordinates by hand

ArrayPlot[face256,PixelConstrained->{1,1}]

Click on the image above to select it. Now bring up the Drawing Tools, under the Graphics menu. Now use the “cross hairs” tool to select the first the lower left point{x0,y0}, and then the upper right point {x1,y1} as the corners of the rectangular patch that you want. Copy and paste in the argument of the Round[] cell below. Here are coordinates for diagonal points for the eyes.

Round[ToExpression@{{"10.09","26.71"},{"36.1","48.99"}}];
Reverse[Transpose[%]];
ArrayPlot[Take[face256,%[[1]],%[[2]]]]

Or using Mathematica’s built-in graphics utility:

Click out the output of Image[]. This will bring up a whole set of image processing tools

Image[]

Geometric image manipulations using function interpolation

You can turn a discrete representation into a continuous one using ListInterpolation[]. Then one can do symbolic operations such as take a derivative, or do morphing.

faceFunction=ListInterpolation[Transpose[face],{{-1,1},{-1,1}},InterpolationOrder->1];
DensityPlot[faceFunction[x,y],{x,-1,1},{y,-1,1},PlotPoints->256,Mesh->False,AspectRatio->Automatic,Frame->None, ColorFunction->"GrayTones"]

Now you can calculate intensities “between pixels”:

faceFunction[.1,.23]

Morphing

faceFunction=ListInterpolation[Transpose[face],{{-1,1},{-1,1}}];

DensityPlot[faceFunction[Sign[x] (x^2),Sign[y] (y^2)],{x,-1,1},{y,-1,1},PlotPoints->100,Mesh->False,AspectRatio->Automatic,Frame->None, ColorFunction->"GrayTones",ImageSize->Small]

FINAL project idea: How does our ability to recognize objects degrade with geometrical deformations?

FINAL project idea: Use Mathematica’s Interpolation functions to model the log polar model of the retinotopy of primary visual cortex to show how the foveal representation gets expanded on cortex.

More filtering: Calculating the spatial gradient of an image using function interpolation

Both first and second spatial derivatives have been used to model edge detection/amplification by neural processes in primary visual cortex. 

Hubel and Wiesel’s “edge detector” can be interpreted as an approximation to a first order derivative (“gradient” in 2D, see ∇f), and their “bar detector” as a second derivative (see (∇^2)G below).

Im this section, we show how to use symbolic derivatives to find the gradient.

The gradient (or gradient vector field) of an image intensity function f,  is given by the vector ∇f.

(OverVector[∇f])=(((∂f)/(∂x)),((∂f)/(∂y)))

The vector length is |∇f|. At a given point, it has a maximum length in the direction of greatest change in intensity.

|∇f|=|(((∂f)/(∂x)),((∂f)/(∂y)))|=(Sqrt[(((((∂f)/(∂x))))^2),(((((∂f)/(∂y))))^2)])

Let filterface =f, where we've blurred out face a little to reduce quantization artifacts:

kernel={{1,1,1},{1,1,1},{1,1,1}}; (*Alternatively, you could use BoxFilter[]*)
filterface=ListConvolve[kernel,face];

faceFunction=ListInterpolation[Transpose[filterface],{{-1,1},{-1,1}}];

Clear[x,y]

nx[x_,y_]:= Evaluate[D[faceFunction[x,y],x]];
ny[x_,y_]:= Evaluate[D[faceFunction[x,y],y]];
ImageGradient[x_,y_]:=Evaluate[Sqrt[D[nx[x,y],x]^2+D[ny[x,y],y]^2] ];

Plot the rate of change in the x-direction, where intensity in the output image is proportional to the rate of change.

temp=Table[nx[x,y],{x,-1,1,.005},{y,-1,1,.005}];
ArrayPlot[Transpose[temp]]

You can use VectorPlot to visualize the gradient field:

VectorPlot[{nx[x,y],ny[x,y]},{x,-1,1},{y,-1,1},ImageSize->Small]

Plot the magnitude of the gradient to highlight regions of the image where contrast is changing the most rapidly:

DensityPlot[ImageGradient[x,y],{x,-1,1},{y,-1,1},PlotPoints->width,Mesh->False,Frame->False, ColorFunction->"Rainbow",ImageSize->Small]

Discrete spatial derivative filtering with built-in image functions

Use can also use image specific filters. GradientFilter[image,r] gives an image corresponding to the magnitude of the gradient of image, using discrete derivatives of a Gaussian of pixel radius r.

GradientFilter[,6]// ImageAdjust

The (∇^2)G operator (one of the forms used to represent a “mexican hat” filter) in the previous lecture first convolves the image with a Gaussian blur kernel, and then takes the Laplacian (∇^2). This filter effectively blurs an image before taking the derivatives.

Later we’ll see why when we study edge detection.

LaplacianGaussianFilter[, 2] // ImageAdjust

Manipulating color images

image=ExampleData[{"TestImage","Mandrill"}]

RGBvalues=ImageData[image];
Dimensions[RGBvalues]

{512,512,3}

{512,512,3}

{512,512,3}

{512,512,3}

{512,512,3}

reds=Map[#[[1]]&,N[RGBvalues],{2}];
greens=Map[#[[2]]&,N[RGBvalues],{2}];
blues=Map[#[[3]]&,N[RGBvalues],{2}];

GraphicsRow[{Image[reds],Image[greens],Image[blues]}]

You can also use the function:

ColorSeparate[];

A weighted average of RGB values to produce a luminance image:

grayscale=Map[((0.3#[[1]]+0.59#[[2]]+0.11#[[3]])/255)&,N[RGBvalues],{2}];
Image[grayscale]//ImageAdjust ;

Note the weights above are arbitrary, and the chosen values will depend on the color calibration.

Putting the R, G, B images back together:

r=reds;
g=greens;
b=blues;
temp2=Partition[Transpose[{Flatten[r], Flatten[g], Flatten[b]}], Dimensions[r][[2]]];
Image[temp2]

Next time

Efficient coding

References

Adelson, E. H., Simoncelli, E., & Hingorani, R. (1987). Orthogonal Pyramid Transforms for Image Coding. Paper presented at the Proc. SPIE - Visual Communication & Image Proc. II, Cambridge, MA.
Barlow, H. B., & Olshausen, B. A. (2004). Convergent evidence for the visual analysis of optic flow through anisotropic attenuation of high spatial frequencies. J Vis, 4(6), 415-426.
Daugman, J. G. (1988). An information-theoretic view of analog representation in striate cortex, Computational Neuroscience. Cambridge, Massachusetts: M.I.T. Press.
Engel, S. A., Glover, G. H., & Wandell, B. A. (1997). Retinotopic organization in human visual cortex and the spatial precision of functional MRI. Cereb Cortex, 7(2), 181-192.
Gold, J. M., Murray, R. F., Bennett, P. J., & Sekuler, A. B. (2000). Deriving behavioural receptive fields for visually completed contours. Curr Biol, 10(11), 663-666.
Hegdé, J., & Kersten, D. (2010). A Link between Visual Disambiguation and Visual Memory. Journal of Neuroscience, 30(45), 15124–15133. doi:10.1523/JNEUROSCI.4415-09.2010
Konishi, S. M., Yuille, A. L., Coughlan, J. M., & Zhu, S. C. (2003). Statistical edge detection: Learning and evaluating edge cues. IEEE Transactions on Pattern Analysis and Machine Intelligence, 25(1), 57-74.
Olman, C. A., & Kersten, D. (2004). Classification objects, ideal observers & generative models. Cognitive Science, 28, 227-239.
Moore, C., & Engel, S. A. (2001). Neural response to perception of volume in the lateral occipital complex. Neuron, 29(1), 277-286.
Schwartz, E. L. (1980). A quantitative model of the functional architecture of human striate cortex with application to visual illusion and cortical texture analysis. Biol Cybern, 37(2), 63-76.
http://library.wolfram.com/howtos/images/#histograms