modules/contrib/src/chowliutree.cpp

   1 /*M///////////////////////////////////////////////////////////////////////////////////////
   2 //
   3 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
   4 //
   5 //  By downloading, copying, installing or using the software you agree to this license.
   6 //  If you do not agree to this license, do not download, install,
   7 //  copy or use the software.
   8 //
   9 // This file originates from the openFABMAP project:
  10 // [http://code.google.com/p/openfabmap/]
  11 //
  12 // For published work which uses all or part of OpenFABMAP, please cite:
  13 // [http://ieeexplore.ieee.org/xpl/articleDetails.jsp?arnumber=6224843]
  14 //
  15 // Original Algorithm by Mark Cummins and Paul Newman:
  16 // [http://ijr.sagepub.com/content/27/6/647.short]
  17 // [http://ieeexplore.ieee.org/xpl/articleDetails.jsp?arnumber=5613942]
  18 // [http://ijr.sagepub.com/content/30/9/1100.abstract]
  19 //
  20 //                           License Agreement
  21 //
  22 // Copyright (C) 2012 Arren Glover [aj.glover@qut.edu.au] and
  23 //                    Will Maddern [w.maddern@qut.edu.au], all rights reserved.
  24 //
  25 //
  26 // Redistribution and use in source and binary forms, with or without modification,
  27 // are permitted provided that the following conditions are met:
  28 //
  29 //   * Redistribution's of source code must retain the above copyright notice,
  30 //     this list of conditions and the following disclaimer.
  31 //
  32 //   * Redistribution's in binary form must reproduce the above copyright notice,
  33 //     this list of conditions and the following disclaimer in the documentation
  34 //     and/or other materials provided with the distribution.
  35 //
  36 //   * The name of the copyright holders may not be used to endorse or promote products
  37 //     derived from this software without specific prior written permission.
  38 //
  39 // This software is provided by the copyright holders and contributors "as is" and
  40 // any express or implied warranties, including, but not limited to, the implied
  41 // warranties of merchantability and fitness for a particular purpose are disclaimed.
  42 // In no event shall the Intel Corporation or contributors be liable for any direct,
  43 // indirect, incidental, special, exemplary, or consequential damages
  44 // (including, but not limited to, procurement of substitute goods or services;
  45 // loss of use, data, or profits; or business interruption) however caused
  46 // and on any theory of liability, whether in contract, strict liability,
  47 // or tort (including negligence or otherwise) arising in any way out of
  48 // the use of this software, even if advised of the possibility of such damage.
  49 //
  50 //M*/
  51
  52 #include "precomp.hpp"
  53 #include "opencv2/contrib/openfabmap.hpp"
  54
  55 namespace cv {
  56
  57 namespace of2 {
  58
  59 ChowLiuTree::ChowLiuTree() {
  60 }
  61
  62 ChowLiuTree::~ChowLiuTree() {
  63 }
  64
  65 void ChowLiuTree::add(const Mat& imgDescriptor) {
  66     CV_Assert(!imgDescriptor.empty());
  67     if (!imgDescriptors.empty()) {
  68         CV_Assert(imgDescriptors[0].cols == imgDescriptor.cols);
  69         CV_Assert(imgDescriptors[0].type() == imgDescriptor.type());
  70     }
  71
  72     imgDescriptors.push_back(imgDescriptor);
  73
  74 }
  75
  76 void ChowLiuTree::add(const std::vector<Mat>& _imgDescriptors) {
  77     for (size_t i = 0; i < _imgDescriptors.size(); i++) {
  78         add(_imgDescriptors[i]);
  79     }
  80 }
  81
  82 const std::vector<cv::Mat>& ChowLiuTree::getImgDescriptors() const {
  83     return imgDescriptors;
  84 }
  85
  86 Mat ChowLiuTree::make(double infoThreshold) {
  87     CV_Assert(!imgDescriptors.empty());
  88
  89     unsigned int descCount = 0;
  90     for (size_t i = 0; i < imgDescriptors.size(); i++)
  91         descCount += imgDescriptors[i].rows;
  92
  93     mergedImgDescriptors = cv::Mat(descCount, imgDescriptors[0].cols,
  94         imgDescriptors[0].type());
  95     for (size_t i = 0, start = 0; i < imgDescriptors.size(); i++)
  96     {
  97         Mat submut = mergedImgDescriptors.rowRange((int)start,
  98             (int)(start + imgDescriptors[i].rows));
  99         imgDescriptors[i].copyTo(submut);
 100         start += imgDescriptors[i].rows;
 101     }
 102
 103     std::list<info> edges;
 104     createBaseEdges(edges, infoThreshold);
 105
 106     // TODO: if it cv_asserts here they really won't know why.
 107
 108     CV_Assert(reduceEdgesToMinSpan(edges));
 109
 110     return buildTree(edges.front().word1, edges);
 111 }
 112
 113 double ChowLiuTree::P(int a, bool za) {
 114
 115     if(za) {
 116         return (0.98 * cv::countNonZero(mergedImgDescriptors.col(a)) /
 117             mergedImgDescriptors.rows) + 0.01;
 118     } else {
 119         return 1 - ((0.98 * cv::countNonZero(mergedImgDescriptors.col(a)) /
 120             mergedImgDescriptors.rows) + 0.01);
 121     }
 122
 123 }
 124 double ChowLiuTree::JP(int a, bool za, int b, bool zb) {
 125
 126     double count = 0;
 127     for(int i = 0; i < mergedImgDescriptors.rows; i++) {
 128         if((mergedImgDescriptors.at<float>(i,a) > 0) == za &&
 129             (mergedImgDescriptors.at<float>(i,b) > 0) == zb) {
 130                 count++;
 131         }
 132     }
 133     return count / mergedImgDescriptors.rows;
 134
 135 }
 136 double ChowLiuTree::CP(int a, bool za, int b, bool zb){
 137
 138     int count = 0, total = 0;
 139     for(int i = 0; i < mergedImgDescriptors.rows; i++) {
 140         if((mergedImgDescriptors.at<float>(i,b) > 0) == zb) {
 141             total++;
 142             if((mergedImgDescriptors.at<float>(i,a) > 0) == za) {
 143                 count++;
 144             }
 145         }
 146     }
 147     if(total) {
 148         return (double)(0.98 * count)/total + 0.01;
 149     } else {
 150         return (za) ? 0.01 : 0.99;
 151     }
 152 }
 153
 154 cv::Mat ChowLiuTree::buildTree(int root_word, std::list<info> &edges) {
 155
 156     int q = root_word;
 157     cv::Mat cltree(4, (int)edges.size()+1, CV_64F);
 158
 159     cltree.at<double>(0, q) = q;
 160     cltree.at<double>(1, q) = P(q, true);
 161     cltree.at<double>(2, q) = P(q, true);
 162     cltree.at<double>(3, q) = P(q, true);
 163     //setting P(zq|zpq) to P(zq) gives the root node of the chow-liu
 164     //independence from a parent node.
 165
 166     //find all children and do the same
 167     std::vector<int> nextqs = extractChildren(edges, q);
 168
 169     int pq = q;
 170     std::vector<int>::iterator nextq;
 171     for(nextq = nextqs.begin(); nextq != nextqs.end(); nextq++) {
 172         recAddToTree(cltree, *nextq, pq, edges);
 173     }
 174
 175     return cltree;
 176
 177
 178 }
 179
 180 void ChowLiuTree::recAddToTree(cv::Mat &cltree, int q, int pq,
 181                                std::list<info>& remaining_edges) {
 182
 183     cltree.at<double>(0, q) = pq;
 184     cltree.at<double>(1, q) = P(q, true);
 185     cltree.at<double>(2, q) = CP(q, true, pq, true);
 186     cltree.at<double>(3, q) = CP(q, true, pq, false);
 187
 188     //find all children and do the same
 189     std::vector<int> nextqs = extractChildren(remaining_edges, q);
 190
 191     pq = q;
 192     std::vector<int>::iterator nextq;
 193     for(nextq = nextqs.begin(); nextq != nextqs.end(); nextq++) {
 194         recAddToTree(cltree, *nextq, pq, remaining_edges);
 195     }
 196 }
 197
 198 std::vector<int> ChowLiuTree::extractChildren(std::list<info> &remaining_edges, int q) {
 199
 200     std::vector<int> children;
 201     std::list<info>::iterator edge = remaining_edges.begin();
 202
 203     while(edge != remaining_edges.end()) {
 204         if(edge->word1 == q) {
 205             children.push_back(edge->word2);
 206             edge = remaining_edges.erase(edge);
 207             continue;
 208         }
 209         if(edge->word2 == q) {
 210             children.push_back(edge->word1);
 211             edge = remaining_edges.erase(edge);
 212             continue;
 213         }
 214         edge++;
 215     }
 216
 217     return children;
 218 }
 219
 220 bool ChowLiuTree::sortInfoScores(const info& first, const info& second) {
 221     return first.score > second.score;
 222 }
 223
 224 double ChowLiuTree::calcMutInfo(int word1, int word2) {
 225     double accumulation = 0;
 226
 227     double P00 = JP(word1, false, word2, false);
 228     if(P00) accumulation += P00 * std::log(P00 / (P(word1, false)*P(word2, false)));
 229
 230     double P01 = JP(word1, false, word2, true);
 231     if(P01) accumulation += P01 * std::log(P01 / (P(word1, false)*P(word2, true)));
 232
 233     double P10 = JP(word1, true, word2, false);
 234     if(P10) accumulation += P10 * std::log(P10 / (P(word1, true)*P(word2, false)));
 235
 236     double P11 = JP(word1, true, word2, true);
 237     if(P11) accumulation += P11 * std::log(P11 / (P(word1, true)*P(word2, true)));
 238
 239     return accumulation;
 240 }
 241
 242 void ChowLiuTree::createBaseEdges(std::list<info>& edges, double infoThreshold) {
 243
 244     int nWords = imgDescriptors[0].cols;
 245     info mutInfo;
 246
 247     for(int word1 = 0; word1 < nWords; word1++) {
 248         for(int word2 = word1 + 1; word2 < nWords; word2++) {
 249             mutInfo.word1 = (short)word1;
 250             mutInfo.word2 = (short)word2;
 251             mutInfo.score = (float)calcMutInfo(word1, word2);
 252             if(mutInfo.score >= infoThreshold)
 253             edges.push_back(mutInfo);
 254         }
 255     }
 256     edges.sort(sortInfoScores);
 257 }
 258
 259 bool ChowLiuTree::reduceEdgesToMinSpan(std::list<info>& edges) {
 260
 261     std::map<int, int> groups;
 262     std::map<int, int>::iterator groupIt;
 263     for(int i = 0; i < imgDescriptors[0].cols; i++) groups[i] = i;
 264     int group1, group2;
 265
 266     std::list<info>::iterator edge = edges.begin();
 267     while(edge != edges.end()) {
 268         if(groups[edge->word1] != groups[edge->word2]) {
 269             group1 = groups[edge->word1];
 270             group2 = groups[edge->word2];
 271             for(groupIt = groups.begin(); groupIt != groups.end(); groupIt++)
 272             if(groupIt->second == group2) groupIt->second = group1;
 273             edge++;
 274         } else {
 275             edge = edges.erase(edge);
 276         }
 277     }
 278
 279     if(edges.size() != (unsigned int)imgDescriptors[0].cols - 1) {
 280         return false;
 281     } else {
 282         return true;
 283     }
 284
 285 }
 286
 287 }
 288
 289 }