From: riQui on
This is a small module where you can find some simple dip-related
functions I've put together taking a couple courses.
You can find some chromatic, erosion effects and a bit of signal-
related code in it.

Comments are often written in Italian but descriptions are in English.
Please let me know if you find errors!

-----------------------------------------------------------------------

(*:Title:Lame DIP*)

(* :Author: Riccardo V. Vincelli, Universita' di Milano-Bicocca *)
(* :Mathematica Version: 7.0 *)
(* :Package Version: 1.0 *)
(* :Context: LameDIP`*)
(* :Summary:
This package is intended to be a first helping hand to those who
want to delve into digital image processing using Mathematica rather
than another system.
It presents some basic functions and can be seen as a lame alternative
to the Digital Image Processing package by Wolfram.
Contents:
-Some eyecandy functions to generate impressive images
(decompositions, neonish effects etc...) by applying builtin functions
to images
-Basic statistics to compare images: focusing on the differences
between images, the contrast of an image.
-Color space transformations operating on triples; you can always try
something like applying them to image data and create a new image from
the result...
*)
(* :History:
I put together stuff I wrote
*)
(* :Updates:
lookup http://reim.webhop.net
*)
(* Botherings? Mail:
Riccardo V. Vincelli - r.vincelli(a)campus.unimib.it
*)
(*License: assume FreeBSD License (though I don't care that much)*)

(**Lame DIP**)

BeginPackage["LameDIP`"]

Print[Date[][[1]],"-",Date[][[2]],"-",Date[][[3]]]
Print["package by: Riccardo aka reim Vincelli"]
Print["Type ?LameDIP`* to visualize package functions"]
Print["Now loading ..."]

Unprotect[
CMYKColorSeparate,
CMYKDominantAdd,
FourierPhase,
FourierSpectrum,
GrayscaleClassicAnaglyph,
HSVToRGB,
ImagePlanes,
MeanSquareError,
PeakSignalToNoiseRatio,
RandomGradient,
RGBColorSeparate,
RGBDominantAdd,
RGBToHSV,
RGBToYCbCr,
RGBToYIQ,
RMSContrast,
RandomPhotoshop,
YCbCrToRGB,
YIQToRGB
]

(*CMYKColorSeparate*)
CMYKColorSeparate::usage = "CMYK version of RGBColorSeparate."

(*CMYKDominantAdd*)
CMYKDominantAdd::usage = "CMYK version of RGBDominantAdd."

(*FourierPhase*)
FourierPhase::usage = "Come implementazione praticamente uguale al
prossimo, ritorna la fase del punto dfp nel piano polare"

(*FourierSpectrum*)
FourierSpectrum::usage = "This function computes the spectrum of the
discrete 'Fourier' transform applied to an image; the result, which is
an array, is finally plotted. As usual, the plot
brings in its center the middlepoint. Recall that, according to
'Fourier' analysis, we're doing nothing but porting the original
discrete signal (the argument image) into the frequency domain;
the term 'Fourier' spectrum denotes the module of the polar form of
the transform."

(*GrayscaleClassicAnaglyph*)
GrayscaleClassicAnaglyph::usage =
"This algorithm produces one of those old fashioned redgreen 3d
pictures, by superimposing a red and a cyan image out of phase; best
viewed with red (left) and cyan (right) 3d glasses, go get some!
Warning: the coefficients for the overlay coordinates as well as the
one needed to cutoff the monochromatic part of the picture were
stimated with common sense from a sample picture: quality may lack..."

(*HSVToRGB*)
HSVToRGB::usage =
"HSVToRGB[img] converts a hue-saturation-value triple to a red-green-
blue one."

(*ImagePlanes*)
ImagePlanes::usage =
"ImagePlanes[img] gives eight bitplanes of img; img must be a byte-
valued image having one or three channels (generalizations are
trivial, check the code)."

(*MeanSquareError*)
MeanSquareError::usage =
"MeanSquareError[imgx,imgy] gives the MSE between images imgx and imgy
of the same kind obviously (useful to quantify the difference between
two differently quantized copies of an image)."

(*PeakSignalToNoiseRatio*)
PeakSignalToNoiseRatio::usage =
"PeakSignalToNoiseRatio[imgx,imgy] is an image statistic based on
MeanSquareError."

(*RandomGradient*)
RandomGradient::usage =
"Applies a randomly picked gradient among the ones available on the
platform."

(*RandomPhotoshop*)
RandomPhotoshop::usage =
"RandomPhotoshop[img] applies a random-pick method or function to img
returning somewhat of artistic version of it. Consider that only a
small group of applications is proposed, who knows you can find
exoteric maps I don't even know doing the job."

(*RGBColorSeparate*)
RGBColorSeparate::usage =
"RGBColorSeparate[img] assumes img to be an rgb one and it returns a
list of the three chromatic components."

(*RGBDominantAdd*)
RGBDominantAdd::usage = "Modifies the image in order to obtain an x-
ish or a x^-1-ish one, where x is obviously one of R, G, B; the
specified channel is scaled by the given factor and the other two by
the inverse of it. Zero is a nonsense value and negative values result
in a black image;
notice that y>1 actually shades in the specified channel whereas y<1
exalts the complementary tone (ie (R, 0.5)->cyan hue)."

(*RGBToHSV*)
RGBToHSV::usage =
"RGBToHSV[rgb] converts the triple from the rgb to the hsv space."

(*RGBToYCbCr*)
RGBToYCbCr::usage =
"RGBToYCbCr[rgb] converts the triple from the rgb to the YCbCr, a
space where Y, Cb and Cr are the luma and the two chrominance
components respectively."

(*RGBToYIQ*)
RGBToYIQ::usage =
"RGBToYIQ[rgb] converts the triple from the rgb to the YIQ, a space
where the Y, I and Q are the luma and the two chrominance components
respectively (as in the YCbCr space)."

(*RMSContrast*)
RMSContrast::usage =
"RMSContrast[img] gives the contrast of img computed through the root
mean square statistic; img must be a three-channel one
(generalizations are trivial, check the code)."

(*YCbCrToRGB*)
YCbCrToRGB::usage =
"YCbCrToRGB[ycbcr] converts a ycbcr triple to a rgb one."

(*YIQToRGB*)
YIQToRGB::usage =
"YIQToRGB[yiq] converts a yiq triple to a rgb one."

Begin["`Private`"]

CMYKColorSeparate[img_] :=
If[ImageChannels[img] == 4,
{
ImageApply[{#[[1]], #[[2]] 0, #[[3]] 0, #[[4]] 0} &, img],
ImageApply[{#[[1]] 0, #[[2]], #[[3]] 0, #[[4]] 0} &, img],
ImageApply[{#[[1]] 0, #[[2]] 0, #[[3]], #[[4]] 0} &, img],
ImageApply[{#[[1]] 0, #[[2]] 0, #[[3]] 0, #[[4]]} &, img]
},
Print["Error: incompatible picture; check the number of channels "]
]

CMYKDominantAdd[img_, chan_, val_] :=
If[
ImageChannels[img] ==
4 && (chan == "C " || chan == "M " || chan == "Y " || chan == "K
"),
Switch[chan,
"C ", ImageApply[{#[[1]] val, #[[2]] (1/val), #[[3]] (1/
val), #[[4]] (1/val)} &, img],
"M ", ImageApply[{#[[1]] (1/val), #[[2]] val, #[[3]] (1/
val), #[[4]] (1/val)} &, img],
"Y ", ImageApply[{#[[1]] (1/val), #[[2]] (1/
val), #[[3]] val, #[[4]] (1/val)} &, img],
"K ", ImageApply[{#[[1]] (1/val), #[[2]] (1/val), #[[3]] (1/
val), #[[4]] val} &, img]
],
Print["Error: incompatible picture or wrong params; check the \
number of channels or the image and the channel and value \
parameters"]
]

FourierPhase[img_]:=
(*Dal codice per dare l' array spettro cambia veramente poco: qua
abbiamo da applicare la fase al posto che l' assoluto, non utilizziamo
il padding: pena, immagine poco dettagliata) e trascuriamo anche la
log (non siamo nel regno dell' oscurita' col sole in mezzo come per lo
spettro)*)
If[
ImageQ[img],
Module[{

tutta=Arg[Fourier[ImageData[ColorConvert[img,"Grayscale"]],FourierParameters-
>{-1,1}]]
},
Module[{
m=Dimensions[tutta][[1]],(*righe*)
n=Dimensions[tutta][[2]] (*colonne*)
},
Module[{
altosx=Take[tutta,{1,Floor[m/2]},{1,Floor[n/2]}],
altodx=Take[tutta,{1,Floor[m/2]},{Floor[n/2],n}],
bassosx=Take[tutta,{Floor[m/2],m},{1,Floor[n/2]}],
bassodx=Take[tutta,{Floor[m/2],m},{Floor[n/2],n}]
},
ImageCrop[ArrayPlot[Log[Join[Join[bassodx,bassosx,
2],Join[altodx,altosx,2]]],Frame->False]]
]
]
],
Print["Error: not an image "]
]

FourierSpectrum[img_]:=
(*L' idea e' applicare la dft (o la fft? non so sotto cosa applichi
poi il sistema) ad un' immagine; i parametri sono specifici per il
data processing generico.
Lavoriamo su greyscale perche' l' idea di potenza del segnale
'colorato' non e' persa, ragionevolmente il grigio e' una media ed
anche perche' via ArrayPlot e' l' unico modo che ho trovato :L
L' immagine data e' prima di essere passata bordata di bianco
(windowing function) per uno spettro visivamente migliore; tali bordi
sono infine croppati*)

If[
ImageQ[img],
Module[{
tutta=Abs[Fourier[ImageData[ImagePad[ColorConvert[img,"Grayscale"],
256,White]],FourierParameters->{-1,1}]]
},
Module[{
m=Dimensions[tutta][[1]],(*righe*)
n=Dimensions[tutta][[2]] (*colonne*)
},
(*Siccome vogliamo avere il 'punto di luce' al centro, dobbiamo
ritagliare e ricomporre i quattro quadranti della matrice...*)
Module[{
altosx=Take[tutta,{1,Floor[m/2]},{1,Floor[n/2]}],
altodx=Take[tutta,{1,Floor[m/2]},{Floor[n/2],n}],
bassosx=Take[tutta,{Floor[m/2],m},{1,Floor[n/2]}],
bassodx=Take[tutta,{Floor[m/2],m},{Floor[n/2],n}]
},
ImageCrop[ArrayPlot[Log[Join[Join[bassodx,bassosx,
2],Join[altodx,altosx,2]]],Frame->False]]
]
]
],
Print["Error: not an image"]
]

GrayscaleClassicAnaglyph[img_]:=
If[
ImageChannels[img]==3,
Module[{
imgr=ImageApply[{#[[1]],#[[2]]0,#[[3]]0}
&,ColorConvert[ColorConvert[img,"Grayscale "],"RGB "]],
imgc=ImageApply[{#[[1]]0,#[[2]],#[[3]]}
&,ColorConvert[ColorConvert[img,"Grayscale "],"RGB "]]
},
Module[{
a=ImageCompose[imgr,{imgc,0.5},{Ceiling[ImageDimensions[img]
[[1]]*183/343]-12,Ceiling[ImageDimensions[img][[2]]*255/512]}]
},
ImageTake[a,ImageDimensions[a][[2]],{17*ImageDimensions[a][[1]]/
343,ImageDimensions[a][[1]]}]
]
],
Print["Error: incompatible picture; check the number of channels "]
]

HSVToRGB[hsv_] :=
Module[{
h = hsv[[1]],
s = hsv[[2]],
v = hsv[[3]]
},
If[
h < 0 || h >= 360 || s < 0 || s > 1 || v < 0 || v > 1,
Print["Error: out of range; correct intervals are: 0<=h<360, \
0<=s,v<=1"],
Module[{
hi = Mod[Floor[h/60], 6],
f = h/60 - Floor[h/60],
p = v (1 - s)},
Module[{
q = v (1 - f s),
t = v (1 - (1 - f) s)
},
Switch[hi,
0, {v, t, p},
1, {q, v, p},
2, {p, v, t},
3, {p, q, v},
4, {t, p, v},
5, {v, p, q}
]
]
]
]
]

ImagePlanes[img_] :=
If[
ImageChannels[img] == 3,
Module[{
a = ColorSeparate[img]
}, {
ImageApply[If[# <= 0.5,#*0,#*0 + 1] &, a[[1]]],
ImageApply[If[# <= 0.25,#*0,#*0 + 1] &, a[[1]]],
ImageApply[If[# <= 0.125,#*0,#*0 + 1] &, a[[1]]],
ImageApply[If[# <= 0.0625,#*0,#*0 + 1] &, a[[1]]],
ImageApply[If[# <= 0.03125,#*0,#*0 + 1] &, a[[1]]],
ImageApply[If[# <= 0.015625,#*0,#*0 + 1] &, a[[1]]],
ImageApply[If[# <= 0.0078125,#*0,#*0 + 1] &, a[[1]]],
ImageApply[If[# <= 0.00390625,#*0,#*0 + 1] &, a[[1]]],
ImageApply[If[# <= 0.5,#*0,#*0 + 1] &, a[[2]]],
ImageApply[If[# <= 0.25,#*0,#*0 + 1] &, a[[2]]],
ImageApply[If[# <= 0.125,#*0,#*0 + 1] &, a[[2]]],
ImageApply[If[# <= 0.0625,#*0,#*0 + 1] &, a[[2]]],
ImageApply[If[# <= 0.03125,#*0,#*0 + 1] &, a[[2]]],
ImageApply[If[# <= 0.015625,#*0,#*0 + 1] &, a[[2]]],
ImageApply[If[# <= 0.0078125,#*0,#*0 + 1] &, a[[2]]],
ImageApply[If[# <= 0.00390625,#*0,#*0 + 1] &, a[[2]]],
ImageApply[If[# <= 0.5,#*0,#*0 + 1] &, a[[3]]],
ImageApply[If[# <= 0.25,#*0,#*0 + 1] &, a[[3]]],
ImageApply[If[# <= 0.125,#*0,#*0 + 1] &, a[[3]]],
ImageApply[If[# <= 0.0625,#*0,#*0 + 1] &, a[[3]]],
ImageApply[If[# <= 0.03125,#*0,#*0 + 1] &, a[[3]]],
ImageApply[If[# <= 0.015625,#*0,#*0 + 1] &, a[[3]]],
ImageApply[If[# <= 0.0078125,#*0,#*0 + 1] &, a[[3]]],
ImageApply[If[# <= 0.00390625,#*0,#*0 + 1] &, a[[3]]]
}
],
If[
ImageChannels[img] == 1,
{
ImageApply[If[# <= 0.5,#*0,#*0 + 1] &, img],
ImageApply[If[# <= 0.25,#*0,#*0 + 1] &, img],
ImageApply[If[# <= 0.125,#*0,#*0 + 1] &, img],
ImageApply[If[# <= 0.0625,#*0,#*0 + 1] &, img],
ImageApply[If[# <= 0.03125,#*0,#*0 + 1] &, img],
ImageApply[If[# <= 0.015625,#*0,#*0 + 1] &, img],
ImageApply[If[# <= 0.0078125,#*0,#*0 + 1] &, img],
ImageApply[If[# <= 0.00390625,#*0,#*0 + 1] &, img]
},
Print["Error: incompatible picture; check the number of channels "]
]
]

MeanSquareError[imgx_,imgy_]:=
(*L'idea e' quella di fare l' mse sui tre canali separatamente e
sommare*)
If[ImageQ[imgx]&&ImageQ[imgy]&&ImageDimensions[imgx]==ImageDimensions[imgy]&&ImageType[imgx]==ImageType[imgy]&&ImageChannels[imgx]==ImageChannels[imgy],
(*applichiamo il prodotto tra le dimensioni dell' immagine:*)
Module[{n=Apply[(#1 #2)&,ImageDimensions[imgx]]},
(*sommiamo le tre componenti mse, una per ogni canale*)
Apply[(#1[[1]]+#2[[2]]+#3[[3]])&,
(*applichiamo la definizione: differenza tra valori pixel*)
(1/n) Apply[((#1[[1]]-#2[[1]])^2)&,{ImageData[imgx],ImageData[imgy]}];
(1/n) Apply[((#1[[2]]-#2[[2]])^2)&,{ImageData[imgx],ImageData[imgy]}];
(1/n) Apply[((#1[[3]]-#2[[3]])^2)&,
{ImageData[imgx],ImageData[imgy]}]]],
Print["Error: incompatible picture couple; check they share the same
dimensions, type, number of channels "]]

PeakSignalToNoiseRatio[imgx_, imgy_] := 10 Log[10, (1/
MeanSquareError[imgx, imgy])]

RGBColorSeparate[img_] :=
If[ImageChannels[img] == 3,
{
ImageApply[{#[[1]], #[[2]] 0, #[[3]] 0} &, img],
ImageApply[{#[[1]] 0, #[[2]], #[[3]] 0} &, img],
ImageApply[{#[[1]] 0, #[[2]] 0, #[[3]]} &, img]
},
Print["Error: incompatible picture; check the number of channels "]
]

RandomPhotoshop[img_] :=
If[
ImageQ[img],
Module[{
r = RandomInteger[{1, 15}]
},
Switch[r,
1, ImageApply[Round[#] &, img],
2, ImageApply[Floor[#] &, img],
3, ImageApply[SquareWave[#] &, img],
4, ImageApply[TriangleWave[#] &, img],
5, ImageApply[SawtoothWave[#] &, img],
6, ImageApply[Cos[#] &, img],
7, ImageApply[Tan[#] &, img],
8, ImageApply[Sinc[#] &, img],
9, ImageApply[Erfc[#] &, img],
10, ImageApply[InverseErfc[#] &, img],
11, ImageApply[AiryAi[#] &, img],
12, ImageApply[EllipticK[#] &, img],
13, ImageApply[HeavisidePi[#] &, img],
14, ImageApply[HeavisideLambda[#] &, img],
15, ImageApply[Mean, img]
]
],
Print["Error: not an image "]
]

RandomGradient[img_]:=
If[
ImageQ[img],
Module[{cdg=ColorData["Gradients "]},
Module[{g=cdg[[RandomInteger[{1,Length[cdg]}]]]},

ImageApply[List@@ColorData[ToString[g],#]&,ColorConvert[img,"Grayscale
"]]]],
Print["Error: not an image "]]

RGBDominantAdd[img_, chan_, val_] :=
If[
ImageChannels[img] ==
3 && (chan == "R " || chan == "G " || chan == "B "),
Switch[chan,
"R ", ImageApply[{#[[1]] val, #[[2]] (1/val), #[[3]] (1/val)} &,
img],
"G ", ImageApply[{#[[1]] val, #[[2]] (1/val), #[[3]] (1/val)} &,
img],
"B ", ImageApply[{#[[1]] (1/val), #[[2]] (1/val), #[[3]] val} &,
img]
],
Print["Error: incompatible picture or wrong params; check the \
number of channels or the image and the channel and value \
parameters"]
]

RGBToHSV[rgb_] :=
Module[{
r = rgb[[1]],
g = rgb[[2]],
b = rgb[[3]]
},
If[
r < 0 || r > 1 || g < 0 || g > 1 || b < 0 || b > 1,
Print["Error: out of range; correct intervals are: 0<=r,g,b<=1"],
Module[{
max = Max[rgb],
min = Min[rgb],
h,
s,
v
},
Switch[max,
min, h = 0,
r, h = Mod[60 ((g - b)/(max - min)) + 360, 360],
g, h = 60 ((b - r)/(max - min)) + 120,
b, h = 60 ((r - g)/(max - min)) + 240
];
If[
max == 0,
s = 0,
s = 1 - min/max
];
v = max;
Return[{h, s, v}]
]
]
]

RGBToYCbCr[rgb_] :=
Module[{
r = rgb[[1]],
g = rgb[[2]],
b = rgb[[3]]
},
If[
r < 0 || r > 1 || g < 0 || g > 1 || b < 0 || b > 1,
Print["Error: out of range; correct intervals are: 0<=r,g,b<=1"],
Return[
{0.2989 r + 0.5866 g + 0.1145 b,
-0.1687 r - 0.3312 g + 0.5000 b,
0.5000 r - 0.4183 g - 0.0816 b}
]
]
]

RGBToYIQ[rgb_] :=
Module[{
r = rgb[[1]],
g = rgb[[2]],
b = rgb[[3]]
},
If[
r < 0 || r > 1 || g < 0 || g > 1 || b < 0 || b > 1,
Print["Error: out of range; correct intervals are: 0<=r,g,b<=1"],
Return[
{0.299 r + 0.587 g + 0.114 b,
0.595716 r - 0.274453 g - 0.321263 b,
0.211456 r - 0.522591 g + 0.311135 b}
]
]
]

RMSContrast[img_] :=
If[
ImageChannels[img] == 3,
Module[{
x = Flatten[ImageData[ColorSeparate[img][[1]]]],
y = Flatten[ImageData[ColorSeparate[img][[2]]]],
z = Flatten[ImageData[ColorSeparate[img][[3]]]]
},
{StandardDeviation[x], StandardDeviation[y], StandardDeviation[z]}
],
If[ImageChannels[img] == 1,
StandardDeviation[Flatten[ImageData[img]]],
Print["Error: incompatible picture; check the number of channels
"]]]

YCbCrToRGB[ycbcr_] :=
Module[{
y = ycbcr[[1]],
cb = ycbcr[[2]],
cr = ycbcr[[3]]
},
If[
y < 0 || y > 1 || cb < 0 || cb > 1 || cr < 0 || cr > 1,
Print["Error: out of range; correct intervals are: \
0<=y,cb,cr<=1"],
Return[
{y + 0.0000 cb + 1.4022 cr,
y - 0.3456 cb - 0.7145 cr,
y + 1.7710 cb + 0.0000 cr}
]
]
]

YIQToRGB[yiq_] :=
Module[{
y = yiq[[1]],
i = yiq[[2]],
q = yiq[[3]]
},
If[
y < 0 || y > 1 || i < 0 || i > 1 || q < 0 || q > 1,
Print["Error: out of range; correct intervals are: 0<=y,i,q<=1"],
Return[
{y + 0.9563 i + 0.621 q,
y - 0.2721 i - 0.6474 q,
y - 1.107 i + 1.7046 q}
]
]
]

End[ ]
Print[" ready."]
EndPackage[ ]





-----------------------------------------------------------------------------------------------------------------------