1 // This is brcv/seg/sdet/algo/sdet_EMD.cxx
2
3 //:
4 //\file
5
6 #include "sdet_EMD.h"
7
8 // ---------------------------------------------------------------------------------//
9 // GRAYSCALE DISTANCE FUNCTIONS
10 // ---------------------------------------------------------------------------------//
11
12 //: Compute Chi^2 distance between two histograms
sdet_chi_sq_dist(const sdet_bin hist1[],const sdet_bin hist2[])13 double sdet_chi_sq_dist(const sdet_bin hist1[], const sdet_bin hist2[])
14 {
15 double dist = 0.0;
16 for (int i=0; i<NBINS; i++)
17 dist += (hist1[i].weight - hist2[i].weight)*(hist1[i].weight - hist2[i].weight)/(hist1[i].weight + hist2[i].weight + 2e-16);
18
19 return dist;
20 }
21
22 //: Compute Bhattacharya distance between two histograms
sdet_bhat_dist(const sdet_bin hist1[],const sdet_bin hist2[])23 double sdet_bhat_dist(const sdet_bin hist1[], const sdet_bin hist2[])
24 {
25 double dist = 0.0;
26 for (int i=0; i<NBINS; i++)
27 dist += hist1[i].weight*hist2[i].weight;
28
29 return dist;
30 }
31
32 //: Compute EMD between two signatures (signature 1=dirt, signature 2 = hole)
33 //
34 // This implementation of the Earth Mover's Distance (EMD) is for one-
35 // dimensional datasets where the distance between clusters can be easily
36 // computed on the fly.
37 //
38 // implemented from the pseudocode in Ruzon's Thesis
sdet_gray_EMD(const sdet_bin dirt[],const sdet_bin hole[])39 double sdet_gray_EMD(const sdet_bin dirt[], const sdet_bin hole[])
40 {
41 int i = -1, j = -1;
42 double leftoverdirt = 0.0, leftoverhole = 0.0, work = 0.0;
43 double dirt_amt, hole_amt;
44
45 while (true) {
46
47 // Compute the amount of mass in the lowest numbered bin that hasn't
48 // been moved yet from the piles of dirt
49 if (leftoverdirt == 0.0)
50 {
51 //advance i to the next non-empty interval
52 i++;
53 while(dirt[i].weight == 0.0 && i < NBINS) i++;
54
55 //if no more intervals
56 if (i == NBINS)
57 return (work); //we're done
58 else
59 dirt_amt = dirt[i].weight;
60 }
61 else //use the amount that was left over from the last move
62 dirt_amt = leftoverdirt;
63
64 // Do the same for the holes
65 if (leftoverhole == 0.0)
66 {
67 //advance j to the next non-empty interval
68 j++;
69 while(hole[j].weight == 0.0 && j < NBINS) j++;
70
71 //if no more intervals
72 if (j == NBINS)
73 return (work); //we're done
74 else
75 hole_amt = hole[j].weight;
76 }
77 else //use the amount that was left over from the last move
78 hole_amt = leftoverhole;
79
80 // Compute the work done moving the smaller amount of mass and decide
81 // how much is left over in each bin.
82 double massmoved = std::min(dirt_amt, hole_amt);
83
84 work += massmoved * (1 - std::exp(-std::fabs(hole[j].value - dirt[i].value) / GAMMA));
85
86 leftoverdirt = dirt_amt - massmoved;
87 leftoverhole = hole_amt - massmoved;
88 }
89 }
90
91 // ---------------------------------------------------------------------------------//
92 // COLOR SIGNATURE DISTANCE FUNCTIONS
93 // ---------------------------------------------------------------------------------//
94
95 //: Compute Chi^2 distance between two histograms
sdet_color_chi_sq_dist(const sdet_color_sig &,const sdet_color_sig &)96 double sdet_color_chi_sq_dist(const sdet_color_sig & /*sig1*/, const sdet_color_sig & /*sig2*/)
97 {
98 double dist = 0.0;
99 //for (int i=0; i<NBINS; i++)
100 // dist += (hist1[i].weight - hist2[i].weight)*(hist1[i].weight - hist2[i].weight)/(hist1[i].weight + hist2[i].weight + 2e-16);
101
102 return dist;
103 }
104
105 //: Compute Bhattacharya distance between two histograms
sdet_color_bhat_dist(const sdet_color_sig &,const sdet_color_sig &)106 double sdet_color_bhat_dist(const sdet_color_sig & /*sig1*/, const sdet_color_sig & /*sig2*/)
107 {
108 double dist = 0.0;
109 //for (int i=0; i<NBINS; i++)
110 // dist += hist1[i].weight*hist2[i].weight;
111
112 return dist;
113 }
114
sdet_color_EMD(sdet_color_sig *,sdet_color_sig *,sdet_flow *,int *)115 double sdet_color_EMD(sdet_color_sig * /*sig1*/, sdet_color_sig * /*sig2*/, sdet_flow* /*flow*/, int* /*flowsize*/)
116 {
117 //This seems to be hard
118
119 return 0;
120 }
121