Re: mutual information of red and green channel
- To: mathgroup at smc.vnet.net
- Subject: [mg81332] Re: [mg81274] mutual information of red and green channel
- From: DrMajorBob <drmajorbob at bigfoot.com>
- Date: Wed, 19 Sep 2007 05:35:33 -0400 (EDT)
- References: <29233217.1190160008456.JavaMail.root@m35>
- Reply-to: drmajorbob at bigfoot.com
You divided by zero and/or calculated Log[0]. The Sum can't include cases
where either of the marginals is zero... even though the joint pdf will be
zero there as well.
Here's an alternative solution for the whole thing, in case it gives you
any ideas.
First I get a file name:
<<Histograms`
filename="C:\\Documents and Settings\\bobby\\My Documents\\My Pictures\\";
{FileNameSetter[Dynamic[filename]],Dynamic[filename]}
Click on the button and navigate to a picture.
This pulls out the red-green channels:
picture = Import[filename];
redGreen = Flatten[picture[[1, 1, All, All, {1, 2}]], 1];
There was a LOT of white in my example, so I used DeleteCases to prune the
data a bit:
redGreen = DeleteCases[redGreen, {255, _} | {_, 255}];
Histogram3D is slow for a picture like mine, but perhaps it's worth seeing
the result anyway:
Histogram3D[redGreen]
This creates the joint PDF in terms of an InterpolatingFunction:
joint = Tally@redGreen;
joint = Sort@
Join[joint, {#, 0} & /@
Complement[Flatten[Table[{i, j}, {i, 0, 255}, {j, 0, 255}], 1],
joint[[All, 1]]]];
joint[[All, 2]] = #/Total@# &[N@joint[[All, 2]]];
joint = Interpolation@joint;
I reused the variable name "joint" to save memory.
This plots the PDF:
Plot3D[joint[x, y], {x, 0, 255}, {y, 0, 255}, PlotRange -> All]
(InterpolatingOrder defaulted to 3; this slightly affects this plot, but
it won't affect the calculations below.)
info[x,y] are the marginal PDFs and the joint information term at {x,y},
based on your formula:
Clear[red, green, info]
red[x_] := red[x] = Sum[joint[x, y], {y, 0, 255}]
green[y_] := green[y] = Sum[joint[x, y], {x, 0, 255}]
info[x_Integer, y_Integer] := Module[{j = joint[x, y]},
If[j == 0, 0, j Log[2, j/(red[x] green[y])]]]
Notice how I avoided dividing by zero or calculating Log[0]. If either
marginal is zero, so is the joint.
And finally, here's the sum:
Sum[info[x, y], {x, 0, 255}, {y, 0, 255}]
3.20989
I don't calculate information every day, so excuse it if I've arrived at
something off-the-wall. Some of the steps may help you anyway.
Bobby
On Mon, 17 Sep 2007 23:40:54 -0500, <vickyisai at gmail.com> wrote:
> Hello
>
> I am using the following code to calculate mutual information of red
> and green channel of an image.
>
> I am calculating Ixy= Sum(i)Sum(j) P(REDi,GREENj) Log[P(REDi,GREENj)/
> P(REDi)/P(GREENj) ]
>
> (that is this formula
> http://upload.wikimedia.org/math/c/b/b/cbb518be041b820958181a932a5cd4ff.png
> )
>
> I am getting an error of "Indeterminate expression 0(-=E2=88=9E=
) =
> encountered.
> "
>
> I believe I am getting this because I am taking Probability of Redi
> from a histogram that I have computed (value of REDi divided by total
> no of pixels ....as my probability of REDi. ame for GREENj
> But that histogram has 0 at some places which is leading to divide by
> zero.
>
> (i plotted that histogram to check that)
>
> How should I calculated mutual information of red and green channel of
> an image.
>
> ( I use Mathematica 5.2 and dont have image processing package)
> Thank You
>
> Code is below:
> ---------------------------------------------------------------------- ----
> image=Import["c:/fish.ppm"];
> (* Show[image]; *)
>
>
>
> part=image[[1,1]];
>
> red=part[[All,All,1]];
>
> histRed=Table[0,{256}];
>
> For[i=1,i<Dimensions[red][[1]],++i,
> For[j=1,j<Dimensions[red][[2]],++j,
> indxr=Floor[red[[i,j]]]+1;
> histRed[[indxr]]=histRed[[indxr]]+1;
> ]
> ]
> (* histR=ListPlot[histRed,PlotJoined=EF=94=A2True]; *)
>
> green=part[[All,All,2]];
> histGreen=Table[0,{256}];
>
> For[i=1,i<Dimensions[green][[1]],++i,
> For[j=1,j<Dimensions[green][[2]],++j,
> indxg=Floor[green[[i,j]]]+1;
> histGreen[[indxg]]=histGreen[[indxg]]+1;
> ]
> ]
> (* histG=ListPlot[histGreen,PlotJoined=EF=94=A2True]; *)
>
>
> (* Joint Prob. Matrix *)
> joinRG=Table[0,{256},{256}];
>
> For[i=1,i<Dimensions[red][[1]],++i,
> For[j=1,j<Dimensions[green][[2]], ++j,
> idx1=Floor[red[[i,j]]]+1;
> idx2=Floor[green[[i,j]]]+1;
> joinRG[[idx1,idx2]]=joinRG[[idx1,idx2]]+1;
> ]
> ]
>
> (* Ixy= Sum(i)Sum(j) P(xi,yj) Log[P(xi,yj)/P(xi)/P(yj) *)
>
> For[i=1,i<Dimensions[histRed][[1]],++i,
> For[j=1,j<Dimensions[histGreen][[1]],++j,
> Ixy = Ixy + joinRG[[i,j]] / 65536 * Log[2 , ,joinRG[[i,j]] /
> histRed[[i]] / histGreen[[j]] ]
> ]
> ]
>
>
>
-- =
DrMajorBob at bigfoot.com
