A module for some digital image processing-related tasks
- To: mathgroup at smc.vnet.net
- Subject: [mg108361] A module for some digital image processing-related tasks
- From: riQui <reimone at gmail.com>
- Date: Mon, 15 Mar 2010 00:06:10 -0500 (EST)
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 at 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[ ] -----------------------------------------------------------------------------------------------------------------------