Threshold-based classification of single model: real-data example

This notebook  assumes that the definitions in the notebook "statistical_modeling.nb" have already been evaluated.

Global variables

File information

\(\(dir = "/Users/nechyba/www/eel6825/course_materials/t5.em_examples/two_objects_data";\)\)

segmentedFile1 = "kellogs1_segmented.ppm";
segmentedFile2 = "tide1_segmented.ppm";
imageType = "PPM";

SetDirectory[dir];

Replacement color

replacementColor = {0, 255, 0};

Plotting options

pointStyle1  =  {Orange,  
          PointSize[ .01]};
pointStyle2  =  {Purple,  PointSize[ .01]};

Classification region colors

classColor1  =  RGBColor[ .5, 1, 1];
classColor2  =  RGBColor[1, 1,  .5]; 
    classColors  =  {classColor1,  classColor2};

pr = {{0, 255}, {0, 255}};

Local definitions

Segment loaded image

SegmentObject[img_, thresh_:  245,  
          color_:  
            replacementColor] := 
Module[{img2}, 
img2  =  
            Map[If[#[[1]]  >  thresh  &&  #[[2]]  >  
                      thresh  &&  #[[3]]  >  thresh,  
                  replacementColor, #] & , 
              img[[1, 1]], {2}]; 
ReplacePart[img, 
            img2, {1, 1}]
];

Segment file image

\(\(ConvertImage[outFile_, inFile_, \ type_, \ thresh_:  245, 
          color_:  
            replacementColor] := 
Module[{img}, 
img = Import[inFile, type]; 
img2 = 
            SegmentObject[img, thresh, color]; 
Export[
            outFile, img2, "PPM"];
];\)\)

Extract RGB data from segmented images

ExtractRGB[img_, 
          color_:  
            replacementColor] := 
Module[{rgb}, 
rgb = 
            Partition[Flatten[img[[1, 1]]], 
              3]; 
DeleteCases[rgb, 
            color]
];

Get representative sample of data

RepresentativeSample[dat_,  
          n_] := 
Module[{d2, 
            m}, 
m  =  
            Ceiling[Length[dat]/n]; 
d2  =  
            First  /@  Partition[dat,  m]];

Data plotting routines

Plot representative sample data (columns a and b from dat ) from a larger data set in 2d

\(\(PlotRepresentativeSample2D[dat_, \ n_, {a_, b_}, \ 
          opt___] := 
Module[{d2, m, 
            rgblabel}, 
rgblabel\  = \ {"Red", \ "Green", \ "Blue"}; 
m\  = \ 
            Ceiling[Length[dat]/n]; 
d2\  = \ 
            First\  /@ \ Partition[dat, \ m]; 
Show[
            Graphics[
              Join[opt, \(Point[{#[\([a]\)], #[\([b]\)]}] &\) /@ \ d2]], \ 
            Frame -> True, AspectRatio -> Automatic, \ 
            FrameLabel -> {rgblabel[\([a]\)], \ rgblabel[\([b]\)]}, \ 
            DisplayFunction -> Identity, \ 
            PlotRange -> {{0, 255}, {0, 255}}]];\)\)

Plot representative sample data  from a larger data set in 3d

PlotRepresentativeSample3D[dat_,  n_,  
          opt___] := 
Module[{d2, 
            m}, 
m  =  
            Ceiling[Length[dat]/n]; 
d2  =  
            First  /@  Partition[dat,  m]; 
Show[
            Graphics3D[Join[opt, Point /@  d2]],  
            DisplayFunction -> Identity, 
            PlotRange -> {{0, 255}, {0, 255}, {0, 255}}]];

data2pixel[{r_, g_, b_, x_, y_}] := {RGBColor[r/255. , g/255. , 
          b/255. ],  Rectangle[{x, -y - 1}, {x + 1, -y}]}

Plot multiple results for mixture model

EMAllPlot[{μ_, Σ_, priors_},  
          gin_] := 
Module[{pr, g2, g4, 
            g5}, 
pr  =  {{0, 1}, {0, 
                1}}; 
g2  =  
            PlotMeans[{μ, Σ, priors},  pr,  Blue,  nS,  
              FrameTicks -> None]; 
g4  =  
            PlotContourPDF[{μ, Σ, priors},  pr, 50, nS,  
              FrameTicks -> None]; 
g5  =  
            PlotPDF[{μ, Σ, priors},  pr, 50, nS,  
              Axes -> None]; 
Show[
            GraphicsArray[{{g2, gin}, {g4, g5}}], yS,  
            imSize]
];

Classify image (2d)

OneClassClassifyImage2DAugmented[im_, pdf_, threshold_, 
          which_: {Red,  Green},  
          scale_:  1] := 
Module[{class, colors, pos1, 
            pos2, lpdf, gr, 
            bl}, 
gr = (List @@ (pointStyle1  //  
                      First))*255; 
bl = (List @@ (pointStyle2 // First))*255; 
colors = {Red, Green, 
              Blue}; 
pos1 = Position[colors, 
                  which[[1]]] // Flatten  // 
              First; 
pos2 = Position[colors, 
                  which[[2]]] // Flatten  // 
              First; 
class  =  
            Map[If[# == replacementColor, replacementColor, 
                  If[(pdf  /.  {x -> #[[pos1]]/scale  //  N, 
                            y -> #[[pos2]]/scale  //  N})  >  
                      threshold,  gr,  bl]] &,  
              im[[1, 1]], {2}]; 
ReplacePart[im, 
            class, {1, 1}]
];

Full classify image (2d)

OneClassClassify2[model_,  threshold_,  obj_,  
          obj2_] := 
Module[{pdf, class1, 
            class2}, 
pdf = 
            MixtureLikelihoodNative[{x, y}, 
              model]; 
class1 = 
            OneClassClassifyImage2DAugmented[obj, pdf, 
              threshold, {Red, Green}, 255. ];  
          class2 = 
            OneClassClassifyImage2DAugmented[obj2, pdf, 
              threshold, {Red, Green}, 255. ]; {class1, 
            class2}
];

OneClassError[class1_,  
          class2_] := 
Module[{cnt1, cnt2, 
            bar}, 
cnt1 = 
            CountPixels[class1]; 
cnt2 = 
            CountPixels[
              class2]; 
bar  =  {cnt1,  cnt2}  //  
              Transpose; 
BarChart[bar[[1]],  
            bar[[2]],  
            BarStyle -> {pointStyle1[[1]], 
                pointStyle2[[1]]},  
Frame -> 
              True, BarLabels -> {"Kellogs", "Tide"},  
            imSize]
];

Count pixels in classified images

CountPixels[classim_, 
          color1_: (List  @@  (pointStyle1  //  First))*255,  
          color2_: (List  @@  (pointStyle2  //  First))*255, 
          backgroundColor_:  
            replacementColor] := 
Module[{cnts, n1, n2, 
            n3}, 
cnts = 
            Frequencies[
              Partition[Flatten[classim[[1, 1]]], 
                3]]; 
n1 = 
            cnts[[Position[cnts, color1]  //  Flatten  //  First, 
                1]]; 
n2 = 
            cnts[[Position[cnts, color2]  //  Flatten  //  First, 
                1]]; 
n3 = 
            cnts[[Position[cnts, backgroundColor]  //  Flatten  //  
                  First, 1]]; 
{100*n1/(n1 + n2),  
              100*n2/(n1 + n2)}  //  
            N
];

3d EM initialization function

Scalar covariance initialization 3d

\(\(EMInitializationScalar3D[x_, \ k_] := 
        Module[{n, \ rl, \ μ, \ σ, \ priors, \ sm, \ del2}, 
	
          n\  = \ Length[x]; 
	
          rl\  = \ Table[{Random[Integer, {1, n}]}, \ {k}]; 
	μ\  = \ 
            Partition[Flatten[\(Take[x, \ #] &\)\  /@ \ rl], \ 3]; 
	
          sm\  = \ 
            Table[\(dist[#, \ μ[\([i]\)]] &\)\  /@ \ μ, \ {i, \ 1, \ 
                k}]; 
	
          del2\  = \ 
            Table[{i, i}\ \ , \ {i, \ 1, 
                k}]; 
	σ\  = \ \(\((#*#)\) &\)\  /@ \ \((Min\  /@ \ 
                  Delete[sm, \ del2])\); 
	
          priors\  = \ Table[1/k, \ {k}]\  // \ N; 
	{μ, \ σ, \ 
            priors}
];\)\)

Full covariance initialization 3d

EMInitializationFull3D[x_, 
          k_] := 
Module[{μ, σ, priors, 
            dim}, 
{μ, σ, priors}  =  
            EMInitializationScalar3D[x,  k]; 
dim  =  
            Length[First[μ]]; 
{μ, (# IdentityMatrix[dim]) &  /@  σ, 
            priors}
];

Load and visualize data

Load segmented images

obj1 = Import[segmentedFile1, "PPM"];
obj2 = Import[segmentedFile2, "PPM"];
Show[GraphicsArray[{{obj1}, {obj2}}],  
        ImageSize -> 288];

Extract RGB data

data1 = ExtractRGB[obj1];
data2 = ExtractRGB[obj2];

2d mixture-of-Gaussians one-class classification (RG plane)

Visualize distributions

which = {1, 2};
gobj1  =  
        PlotRepresentativeSample2D[data1,  2000, which, 
          pointStyle1];
gobj2  =  
        PlotRepresentativeSample2D[data2,  2000, which, 
          pointStyle2];
g1 = 
        Show[gobj1, FrameTicks -> False, FrameLabel -> None, 
          nS];
g2 = 
        Show[gobj2, FrameTicks -> False, FrameLabel -> None, 
          nS];
gbase = Show[g1, g2, nS];
Show[GraphicsArray[{g1, g2}],  imSize,  yS];

EM experiment #1: object 1 (Kellogs car)

Extract 2d data

d2d1 = ({#[[which[[1]]]], #[[which[[2]]]]} & /@  (data1/255. ));

Choose number of Gaussians in mixture model

k = 4;

Assumes full covariance matrices

initFull  =  EMInitializationFull[d2d1, k];

Compute EM algorithm

Timing[em  =  EM[d2d1, initFull, 0.001];]

Save solution

model1 =  em  //  Last;

Final model

EMAllPlot[model1, g1];

One-class image classification for different thresholds (0.5 to 5.0)

thresholdList = { .5, 1. , 1.5, 2. , 2.5, 3. , 3.5, 4. , 4.5, 
          5. };
results = OneClassClassify2[model1, #, obj1, obj2] &  /@  
          thresholdList;
Show[GraphicsArray[{{obj1, #[[1]]}, {obj2, #[[2]]}}],  
            ImageSize -> 400] &  /@  results;

Error plot for different threshold (0.5 to 5.0)

OneClassError[#[[1]], #[[2]]] & /@ results;