MathGroup Archive 2007

[Date Index] [Thread Index] [Author Index]

Search the Archive

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


  • Prev by Date: Re: piecewise functions from logical relationships (ie. solving with constraints)
  • Next by Date: Re: Re: LegendreP error (bug?) in Mathematica
  • Previous by thread: Re: mutual information of red and green channel
  • Next by thread: Linux V6: slow input