Feature test for sky/ground classification

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

Local definitions

Turn off integration warning message

Off[NIntegrate::ncvb];
Off[NIntegrate::slwcon];

Read whole image

\(\(ReadWholeImage[file_, \ 
          dir_: "."] := 
Module[{im, 
            d}, 
im = 
            Import[dir <> "/" <> 
                file <> ".ppm", \ "PPM"]; 
im
];\)\)

Read in split images

ReadSplitImages[file_,  
          dir_: "."] := 
Module[{im1, im2, raw1, raw2, 
            d1, d2}, 
im1 = 
            Import[dir <> "/" <> 
                file <> "_A.ppm", "PPM"]; 
im2 = 
            Import[dir <> "/" <> 
                file <> "_B.ppm", "PPM"]; 
raw1 = 
            im1[[1, 1]]; 
raw2 = 
            im2[[1, 1]]; 
d1 = 
            DeleteCases[
              Partition[Flatten[raw1], 3], {255, 255, 
                255}]; 
d2 = 
            DeleteCases[
              Partition[Flatten[raw2], 3], {255, 255, 
                255}]; 
{d1, d2, im1, 
            im2}
];

Compute 1D pdfs

ComputePDF1D[data_, 
          which_:  Red] := 
Module[{colors, pos, 
            d2}, 
colors = {Red, Green, 
              Blue}; 
pos = Position[colors, which] // 
                Flatten  // First; 
d2 = 
            N[(#[[pos]]) & /@ 
                data]; 
stats = {Mean[d2], 
              StandardDeviation[d2]}; 
pdf1D[stats, 
            x]
];

One-dimensional decision boundary

ComputeDecisionBoundary1D[pdf1_, pdf2_, 
          x_:  x] := 
Module[{sol, lpdf1, 
            lpdf2}, 
lpdf1 = 
            Log[pdf1] // PowerExpand; 
lpdf2 = 
            Log[pdf2] // 
              PowerExpand; 
sol = (#[[2]]) & /@ (Solve[lpdf1  ==  lpdf2, x]  //  
                  Flatten); 
If[
            Length[sol] == 1,  
sol  //  
              First, 
If[
              sol[[1]]  ≥ 0  &&  
                sol[[1]]  ≤  255,  sol  //  First,  
              sol  //  Last]]
];

Expected error (1d)

ExpectedError1D[pdf1_,  pdf2_,  
          x_:  x] := 
NIntegrate[
            Min[pdf1, pdf2], {x, 0, 255}]*100;

Expected error (2d)

ExpectedError2D[pdf1_,  pdf2_,  x_:  x,  
          y_:  y] := 
NIntegrate[
            Min[pdf1, pdf2], {x, 0, 255}, {y, 0, 255}]*100;

Classify image (1d)

ClassifyImage1D[im_, db_, 
          which_:  Red] := 
Module[{class, colors, pos, gr, 
            bl}, 
colors = {Red, Green, 
              Blue}; 
pos = Position[colors, which] // 
                Flatten  // 
              First; 
gr = (List @@ 
                  Orange)*255; 
bl = (List @@ 
                  Blue)*255; 
class = 
            Map[If[#[[pos]] < db, gr, bl] &, 
              im[[1, 1]], {2}]; 
ReplacePart[im, 
            class, {1, 1}]
];

Classify image (2d)

ClassifyImage2D[im_, pdf1_, pdf2_,  
          which_: {Red,  Green}] := 
Module[{class, 
            colors, pos1, pos2, lpdf, gr, 
            bl}, 
gr = (List @@ 
                  Orange)*255; 
bl = (List @@ 
                  Blue)*255; 
colors = {Red, Green, 
              Blue}; 
pos1 = Position[colors, 
                  which[[1]]] // Flatten  // 
              First; 
pos2 = Position[colors, 
                  which[[2]]] // Flatten  // 
              First; 
lpdf  =  (Log[pdf1]  // 
                  PowerExpand)  -  (Log[pdf2]  //  
                  PowerExpand); 
class  =  
            Map[If[(lpdf  /.  {x -> #[[pos1]], 
                          y -> #[[pos2]]})  >  0,  gr,  
                  bl] &,  
              im[[1, 1]], {2}]; 
ReplacePart[im, 
            class, {1, 1}]
];

Colors representing two classes

Data colors

data1Color = Green;
data2Color = Blue;

style1  =  {data1Color, 
          PointSize[ .02]};
style2  =  {data2Color, PointSize[ .02]};

Classification region colors

class1Color  =  RGBColor[1, 1,  .5];
class2Color  =  RGBColor[ .5, 1, 1];

classColors  =  {class1Color,  class2Color};

Default plot range

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

Image file names

Image file names

TrainingFile = "blue3";
TestFile = "blue4";

Set file names/paths

TrainingDirectory = "small_split_images";
WholeImageDirectory = "small_whole_images";

Set directory

SetDirectory["/Users/nechyba/Desktop/lecture5/sky_feature_selection/\
"];

One-channel classification experiments

Read in training and test images

{d1, d2, s1, s2} = 
        ReadSplitImages[TrainingFile,  
          TrainingDirectory];
trn = 
        ReadWholeImage[TrainingFile, 
          WholeImageDirectory];
tst = ReadWholeImage[TestFile, WholeImageDirectory];

Show training and test images

Show[GraphicsArray[{{trn, tst}, {s2, s1}}],  
        ImageSize -> 400];

Compute one-dimensional pdfs and corresponding error on 1d color channels

p1r = ComputePDF1D[d1, Red];
p1g = ComputePDF1D[d1, Green];
p1b = ComputePDF1D[d1, Blue];

    
p2r = ComputePDF1D[d2, Red];
p2g = ComputePDF1D[d2, Green];
p2b = ComputePDF1D[d2, Blue];

pl1 = 
        ShowTwoGaussianPDFs1D[p1r,  p2r, {0, 255}, x, 
          FrameTicks -> None, nS];
pl2 = 
        ShowTwoGaussianPDFs1D[p1g,  p2g, {0, 255}, x, 
          FrameTicks -> None, nS];
pl3 = 
        ShowTwoGaussianPDFs1D[p1b,  p2b, {0, 255}, x, 
          FrameTicks -> None, nS];

errors = {ExpectedError1D[p1r, p2r, x], ExpectedError1D[p1g, p2g, x], 
          ExpectedError1D[p1b, p2b, x]};
pl4 = 
        BarChart[errors,  Frame -> True,  
          BarStyle -> {Red, Green, Blue}, 
          FrameTicks -> {None, Automatic, None, Automatic}, 
          GridLines -> {None, Automatic}, nS];

Show one-dimensional pdfs

Show[GraphicsArray[{{pl1, pl2}, {pl3, pl4}}], 
        ImageSize -> 400];

Classify training and test image

cr1 = 
        ClassifyImage1D[trn, ComputeDecisionBoundary1D[p1r, p2r, x], 
          Red];
cg1 = 
        ClassifyImage1D[trn, ComputeDecisionBoundary1D[p1g, p2g, x], 
          Green];
cb1 = 
        ClassifyImage1D[trn, ComputeDecisionBoundary1D[p1b, p2b, x], 
          Blue];

cr2 = 
        ClassifyImage1D[tst, ComputeDecisionBoundary1D[p1r, p2r, x], 
          Red];
cg2 = 
        ClassifyImage1D[tst, ComputeDecisionBoundary1D[p1g, p2g, x], 
          Green];
cb2 = 
        ClassifyImage1D[tst, ComputeDecisionBoundary1D[p1b, p2b, x], 
          Blue];

Show classification results

Show[GraphicsArray[{{trn, cr1}, {cg1, cb1}}],  
        ImageSize -> 400];

Show[GraphicsArray[{{tst, cr2}, {cg2, cb2}}],  
        ImageSize -> 400];

Two-channel classification experiments

Show training and test images

Show[GraphicsArray[{{trn, tst}, {s2, s1}}],  
        ImageSize -> 400];

Visualize two-channel color distributios of training image

RG plane

gd1  =  
        PlotData2D[{#[[1]], #[[2]]} &  /@  d1  // N, pr, 
          style1, nS];
gd2  =  
        PlotData2D[{#[[1]], #[[2]]} &  /@  d2  // N, pr, 
          style2, nS];
distRG  =  
        Show[gd1, gd2, FrameTicks -> None,  yS];

RB plane

gd1  =  
        PlotData2D[{#[[1]], #[[3]]} &  /@  d1  // N, pr, 
          style1, nS];
gd2  =  
        PlotData2D[{#[[1]], #[[3]]} &  /@  d2  // N, pr, 
          style2, nS];
distRB  =  
        Show[gd1, gd2,  FrameTicks -> None, yS];

GB plane

gd1  =  
        PlotData2D[{#[[2]], #[[3]]} &  /@  d1  // N, pr, 
          style1, nS];
gd2  =  
        PlotData2D[{#[[2]], #[[3]]} &  /@  d2  // N, pr, 
          style2, nS];
distGB  =  
        Show[gd1, gd2,  FrameTicks -> None, yS];

Compute two-dimensional pdfs and corresponding error on 2d color channel combinations

p1rg = 
        gaussian[{#[[1]], #[[2]]} &  /@  d1  // N, {x, 
            y}];
p1rb = 
        gaussian[{#[[1]], #[[3]]} &  /@  d1  // N, {x, 
            y}];
p1gb = 
        gaussian[{#[[2]], #[[3]]} &  /@  d1  // N, {x, 
            y}];

    
p2rg = 
        gaussian[{#[[1]], #[[2]]} &  /@  d2  // N, {x, 
            y}];
p2rb = 
        gaussian[{#[[1]], #[[3]]} &  /@  d2  // N, {x, 
            y}];
p2gb = 
        gaussian[{#[[2]], #[[3]]} &  /@  d2  // N, {x, 
            y}];

Show two-dimensional pdfs (skip)

ShowTwoGaussianPDFs[p1rg, p2rg, {{0, 255}, {0, 255}}, x, y,  
        ImageSize -> 400];

ShowTwoGaussianPDFs[p1rb, p2rb, {{0, 255}, {0, 255}}, x, y,  
        ImageSize -> 400];

ShowTwoGaussianPDFs[p1gb, p2gb, {{0, 255}, {0, 255}}, x, y,  
        ImageSize -> 400];

Show decision regions

RG plane

decRG = 
        ShowTwoClassificationRegions[p1rg, p2rg, pr, x, y, classColors, 
          False, Axes -> False, Frame -> True,  
          FrameTicks -> None, nS];
Show[GraphicsArray[{{decRG, distRG}}], 
        ImageSize -> 400];

RB plane

decRB = 
        ShowTwoClassificationRegions[p1rb, p2rb, pr, x, y, classColors, 
          False, Axes -> False, Frame -> True,  
          FrameTicks -> None, nS];
Show[GraphicsArray[{{decRB, distRB}}], 
        ImageSize -> 400];

RB plane

decGB = 
        ShowTwoClassificationRegions[p1gb, p2gb, pr, x, y, classColors, 
          False, Axes -> False, Frame -> True,  
          FrameTicks -> None, nS];
Show[GraphicsArray[{{decGB, distGB}}], 
        ImageSize -> 400];

Compute expected error

errors = {ExpectedError2D[p1rg, p2rg, x, y], 
          ExpectedError2D[p1rb, p2rb, x, y], 
          ExpectedError2D[p1gb, p2gb, x, y]};
BarChart[errors,  Frame -> True,  
        BarStyle -> {Red, Green, Blue}, 
        FrameTicks -> {None, Automatic, None, Automatic}, 
        GridLines -> {None, Automatic}];

Classify training and test image

crg1 = 
        ClassifyImage2D[trn, p1rg, 
          p2rg, {Red, Green}];
crb1 = ClassifyImage2D[trn, p1rb, p2rb, {Red, Blue}]; 
    cgb1 = ClassifyImage2D[trn, p1gb, p2gb, {Green,  Blue}];

crg2 = 
        ClassifyImage2D[tst, p1rg, 
          p2rg, {Red, Green}];
crb2 = ClassifyImage2D[tst, p1rb, p2rb, {Red, Blue}]; 
    cgb2 = ClassifyImage2D[tst, p1gb, p2gb, {Green,  Blue}];

Show classification results

Show[GraphicsArray[{{trn, crg1}, {crb1, cgb1}}],  
        ImageSize -> 400];

Show[GraphicsArray[{{tst, crg2}, {crb2, cgb2}}],  
        ImageSize -> 400];