center_surround.cpp
Go to the documentation of this file.
00001 /*********************************************************************
00002  * Software License Agreement (BSD License)
00003  *
00004  *  Copyright (c) 2010, Georgia Institute of Technology
00005  *  All rights reserved.
00006  *
00007  *  Redistribution and use in source and binary forms, with or without
00008  *  modification, are permitted provided that the following conditions
00009  *  are met:
00010  *
00011  *   * Redistributions of source code must retain the above copyright
00012  *     notice, this list of conditions and the following disclaimer.
00013  *   * Redistributions in binary form must reproduce the above
00014  *     copyright notice, this list of conditions and the following
00015  *     disclaimer in the documentation and/or other materials provided
00016  *     with the distribution.
00017  *   * Neither the name of the Georgia Institute of Technology nor the names of
00018  *     its contributors may be used to endorse or promote products derived
00019  *     from this software without specific prior written permission.
00020  *
00021  *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
00022  *  "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
00023  *  LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
00024  *  FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
00025  *  COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
00026  *  INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
00027  *  BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
00028  *  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
00029  *  CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
00030  *  LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
00031  *  ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
00032  *  POSSIBILITY OF SUCH DAMAGE.
00033  *********************************************************************/
00034 
00035 #include <cpl_visual_features/saliency/center_surround.h>
00036 #include <opencv2/highgui/highgui.hpp>
00037 
00038 #include <iostream>
00039 #include <sstream>
00040 
00041 // #define DISPLAY_SALIENCY_MAPS
00042 // #define DISPLAY_COLOR_MAPS
00043 // #define DISPLAY_GABOR_FILTERS
00044 
00045 #define R_INDEX 2
00046 #define G_INDEX 1
00047 #define B_INDEX 0
00048 
00049 using cv::Mat;
00050 using std::vector;
00051 namespace cpl_visual_features
00052 {
00053 
00054 //
00055 // Constructors and Initialization functions
00056 //
00057 CenterSurroundMapper::CenterSurroundMapper(int min_c, int max_c, int min_delta,
00058                                            int max_delta) :
00059     min_c_(min_c), max_c_(max_c), min_delta_(min_delta), max_delta_(max_delta),
00060     N_(4), gabor_size_(7)
00061 {
00062   num_scales_ = max_c_ + max_delta_;
00063   generateGaborFilters();
00064 }
00065 
00066 void CenterSurroundMapper::generateGaborFilters()
00067 {
00068   Mat base_x(gabor_size_, gabor_size_, CV_64FC1);
00069   Mat base_y(gabor_size_, gabor_size_, CV_64FC1);
00070 
00071   // Populate the base x and y matrices
00072   for (int i = 0; i < base_x.cols; ++i)
00073   {
00074     for (int j = 0; j < base_x.rows; ++j)
00075     {
00076       base_x.at<double>(i,j) = i - gabor_size_/2;
00077       base_y.at<double>(i,j) = j - gabor_size_/2;
00078     }
00079   }
00080 
00081   for (int alpha = 0; alpha < N_; ++alpha)
00082   {
00083     float theta = M_PI / N_*alpha;
00084     Mat x_theta = base_x *  cos(theta) + base_y * sin(theta);
00085     Mat y_theta = base_x * -sin(theta) + base_y * cos(theta);
00086 
00087     Mat gabor_alpha(gabor_size_, gabor_size_, CV_64FC1, 1.0);
00088     gabor_alpha *= 1/(2*M_PI);
00089 
00090     Mat to_exp(gabor_size_, gabor_size_, CV_64FC1);
00091     to_exp = (x_theta.mul(x_theta) + y_theta.mul(y_theta))*-0.5;
00092 
00093     for(int i = 0; i < gabor_size_; i++)
00094     {
00095       for(int j = 0; j < gabor_size_; j++)
00096       {
00097         gabor_alpha.at<double>(i,j) *= exp(to_exp.at<double>(i,j))*
00098             cos(x_theta.at<double>(i,j)*M_PI*2);
00099       }
00100     }
00101 
00102     gabor_filters_.push_back(gabor_alpha);
00103   }
00104 }
00105 
00106 //
00107 // High level interaction functions
00108 //
00109 
00110 Mat CenterSurroundMapper::operator()(Mat& frame, bool use_gradient)
00111 {
00112   Mat saliency_map = getSaliencyMap(frame);
00113   Mat gradient_map(saliency_map.rows, saliency_map.cols, CV_32FC1);
00114   if (use_gradient) // TODO: Make nicer functionality for specifying this map
00115   {
00116     saliency_map.copyTo(gradient_map);
00117 
00118     for (int i = 0; i < gradient_map.rows; ++i)
00119     {
00120       for (int j = 0; j < gradient_map.cols; ++j)
00121       {
00122         gradient_map.at<float>(i,j) = i*1.0/gradient_map.rows;
00123         // gradient_map.at<float>(i,j) = (0.5*i/gradient_map.rows +
00124         //                                0.5 - std::abs(gradient_map.cols/2 - j)
00125         //                                *1.0/gradient_map.cols);
00126       }
00127     }
00128 
00129     saliency_map = saliency_map*0.75 + gradient_map*0.25;
00130   }
00131 
00132 
00133   Mat scaled = scaleMap(saliency_map);
00134   Mat display_scaled = upSampleResponse(scaled, min_delta_, frame.size());
00135 #ifdef DISPLAY_SALIENCY_MAPS
00136   if (use_gradient)
00137     cv::imshow("Top Down map", gradient_map);
00138   cv::imshow("Saliency", scaled);
00139   cv::imshow("Map", saliency_map);
00140   cv::waitKey();
00141 #endif // DISPLAY_SALIENCY_MAPS
00142 
00143   return display_scaled;
00144 }
00145 
00146 
00147 Mat CenterSurroundMapper::operator()(Mat& frame, Mat& depth_map)
00148 {
00149   // TODO: Split frame_saliency into the component maps to do normalization
00150   // separately here for each channel?
00151 
00152   Mat frame_saliency = getSaliencyMap(frame);
00153   Mat depth_i = getIntensityMap(depth_map);
00154   Mat depth_o = getOrientationMap(depth_map);
00155 
00156   // Combine depth_saliency_i and _o
00157   float bar_max = 0;
00158   for (int r = 0; r < depth_i.rows; ++r)
00159   {
00160     for (int c = 0; c < depth_i.cols; ++c)
00161     {
00162       if (depth_i.at<float>(r,c) > bar_max)
00163         bar_max = depth_i.at<float>(r,c);
00164       if (depth_o.at<float>(r,c) > bar_max)
00165         bar_max = depth_o.at<float>(r,c);
00166       if (frame_saliency.at<float>(r,c) > bar_max)
00167         bar_max = frame_saliency.at<float>(r,c);
00168     }
00169   }
00170 
00171 
00172   Mat saliency_map(frame_saliency.rows, frame_saliency.cols, CV_32FC1);
00173   saliency_map = (normalize(frame_saliency, bar_max)*0.75 +
00174                   normalize(depth_i, bar_max)*0.125 +
00175                   normalize(depth_o, bar_max)*0.125);
00176   Mat scaled;
00177   scaled = scaleMap(saliency_map);
00178   Mat display_scaled = upSampleResponse(scaled, min_delta_, frame.size());
00179 #ifdef DISPLAY_SALIENCY_MAPS
00180   Mat display_fs = upSampleResponse(frame_saliency, min_delta_, frame.size());
00181   Mat display_i = upSampleResponse(depth_i, min_delta_, frame.size());
00182   Mat display_o = upSampleResponse(depth_o, min_delta_, frame.size());
00183   cv::imshow("Frame Saliency", display_fs);
00184   cv::imshow("Depth i", display_i);
00185   cv::imshow("Depth o", display_o);
00186   cv::imshow("Combined Saliency", display_scaled);
00187   // cv::waitKey();
00188 #endif // DISPLAY_SALIENCY_MAPS
00189 
00190   return display_scaled;
00191 }
00192 
00193 //
00194 // Core Saliency Functions
00195 //
00196 
00197 Mat CenterSurroundMapper::getSaliencyMap(Mat& frame)
00198 {
00199   Mat I(frame.rows, frame.cols, CV_32FC1);
00200   Mat R(frame.rows, frame.cols, CV_32FC1);
00201   Mat G(frame.rows, frame.cols, CV_32FC1);
00202   Mat B(frame.rows, frame.cols, CV_32FC1);
00203   Mat Y(frame.rows, frame.cols, CV_32FC1);
00204 
00205   vector<Mat> I_scales;
00206   vector<Mat> R_scales;
00207   vector<Mat> G_scales;
00208   vector<Mat> B_scales;
00209   vector<Mat> Y_scales;
00210   vector<vector<Mat> > O_n_theta; // First index scale, second index relative
00211                                   // orientation (default 0 - 3)
00212 
00213   // Get the component color channels
00214   vector<Mat> channels;
00215   Mat to_split(frame.rows, frame.cols, CV_32FC3);
00216   frame.convertTo(to_split, to_split.type(), 1.0/255.0);
00217   split(to_split, channels);
00218 
00219   // Get hue indepenednt intensity channel
00220   I = channels[R_INDEX]/3.0 + channels[G_INDEX]/3.0 + channels[B_INDEX]/3.0;
00221 
00222   float max_i = 0;
00223   for(int r = 0; r < I.rows; ++r)
00224   {
00225     for(int c = 0; c < I.cols; ++c)
00226     {
00227       if (I.at<float>(r,c) > max_i)
00228         max_i = I.at<float>(r,c);
00229     }
00230   }
00231 
00232   // Normalize r, g, and b at locations greater than 1/10 the image
00233   for(int r = 0; r < I.rows; ++r)
00234   {
00235     for(int c = 0; c < I.cols; ++c)
00236     {
00237       if (I.at<float>(r,c) > max_i*0.10)
00238       {
00239         float intensity = I.at<float>(r,c);
00240         channels[R_INDEX].at<float>(r,c) = (channels[R_INDEX].at<float>(r,c) /
00241                                             intensity);
00242         channels[G_INDEX].at<float>(r,c) = (channels[G_INDEX].at<float>(r,c) /
00243                                             intensity);
00244         channels[B_INDEX].at<float>(r,c) = (channels[B_INDEX].at<float>(r,c) /
00245                                             intensity);
00246       }
00247       else
00248       {
00249         channels[0].at<float>(r,c) = 0;
00250         channels[1].at<float>(r,c) = 0;
00251         channels[2].at<float>(r,c) = 0;
00252       }
00253     }
00254   }
00255 
00256   // Get intensity independent hue channels
00257   R = channels[R_INDEX] - (channels[G_INDEX]/2.0 + channels[B_INDEX]/2.0);
00258   G = channels[G_INDEX] - (channels[R_INDEX]/2.0 + channels[B_INDEX]/2.0);
00259   B = channels[B_INDEX] - (channels[R_INDEX]/2.0 + channels[G_INDEX]/2.0);
00260 
00261   Mat Y_abs(Y.rows, Y.cols, Y.type());
00262   for(int r = 0; r < Y.rows; ++r)
00263   {
00264     for(int c = 0; c < Y.cols; ++c)
00265     {
00266       Y_abs.at<float>(r,c) = std::abs(channels[R_INDEX].at<float>(r,c)/2.0 -
00267                                       channels[G_INDEX].at<float>(r,c)/2.0);
00268     }
00269   }
00270 
00271   Y = (channels[R_INDEX]/2.0 + channels[G_INDEX]/2.0) - Y_abs-channels[B_INDEX];
00272 
00273 #ifdef DISPLAY_COLOR_MAPS
00274   cv::imshow("I",I);
00275   cv::imshow("R",R);
00276   cv::imshow("G",G);
00277   cv::imshow("B",B);
00278   cv::imshow("Y",Y);
00279   cv::waitKey();
00280 #endif // DISPLAY_COLOR_MAPS
00281 
00282   // Get copies of the four feature maps at all scales
00283   cv::buildPyramid(I, I_scales, num_scales_);
00284   cv::buildPyramid(R, R_scales, num_scales_);
00285   cv::buildPyramid(G, G_scales, num_scales_);
00286   cv::buildPyramid(B, B_scales, num_scales_);
00287   cv::buildPyramid(Y, Y_scales, num_scales_);
00288 
00289   //
00290   // Build Gabor orientation Pyramid
00291   //
00292 
00293   // Get laplacians at all scales (TODO: do this while building pyr via DoG)
00294   for(unsigned int i = 0; i < I_scales.size(); ++i)
00295   {
00296     Mat lap;
00297     cv::Laplacian(I_scales[i], lap, I_scales[i].depth());
00298     vector<Mat> O_theta;
00299 
00300     // Get the N orientation maps for each scale
00301     for (int a = 0; a < N_; a++)
00302     {
00303       Mat m_a(lap.rows, lap.cols, lap.type());
00304       cv::filter2D(lap, m_a, -1, gabor_filters_[a]);
00305 
00306       // For each of the N orientation maps smooth and downsample
00307       O_theta.push_back(m_a);
00308       // std::stringstream m_name;
00309       // m_name << "O_" << i << "_" << a;
00310       // cv::imshow(m_name.str(), m_a);
00311       // cv::waitKey();
00312     }
00313     O_n_theta.push_back(O_theta);
00314   }
00315 
00316   //
00317   // Get multi-scale maps of the features
00318   //
00319   vector<Mat> I_cs;
00320   vector<Mat> RG_cs;
00321   vector<Mat> BY_cs;
00322   vector<Mat> C_cs;
00323   vector<vector<Mat> > O_theta_cs;
00324 
00325   // Initialize the vectors for the different orientations
00326   for (int a = 0; a < N_; ++a)
00327   {
00328     vector<Mat> O_cs;
00329     O_cs.clear();
00330     O_theta_cs.push_back(O_cs);
00331   }
00332 
00333   // Here we build the multi-scale maps
00334   for (int c = min_c_; c <= max_c_; c++)
00335   {
00336     for (int s = c + min_delta_; s <= c + max_delta_; s++)
00337     {
00338       // Build intensity maps
00339       I_cs.push_back(mapDifference(I_scales[c], I_scales[s], c, s));
00340 
00341       // Build red-green opponent color maps
00342       Mat RG_c = R_scales[c] - G_scales[c];
00343       Mat GR_s = G_scales[s] - R_scales[s];
00344       RG_cs.push_back(mapDifference(RG_c, GR_s, c, s));
00345 
00346       // Build blue-yellow opponent color maps
00347       Mat BY_c = B_scales[c] - Y_scales[c];
00348       Mat YB_s = Y_scales[s] - B_scales[s];
00349       BY_cs.push_back(mapDifference(BY_c, YB_s, c, s));
00350 
00351       // Build orientation maps for each of the N_ orientations
00352       for (int a = 0; a < N_; ++a)
00353       {
00354         O_theta_cs[a].push_back(mapDifference(O_n_theta[c][a],
00355                                               O_n_theta[s][a], c, s));
00356       }
00357     }
00358   }
00359 
00360   //
00361   // Normalize all maps, based on feature type
00362   //
00363 
00364   // Find max values to normalize maps by feature type
00365   float I_max = 0;
00366   float C_max = 0;
00367   vector<float> O_theta_max(N_, 0);
00368 
00369   for (unsigned int i = 0; i < I_cs.size(); ++i)
00370   {
00371     for(int r = 0; r < I_cs[i].rows; ++r)
00372     {
00373       for(int c = 0; c < I_cs[i].cols; ++c)
00374       {
00375         if (I_cs[i].at<float>(r,c) > I_max)
00376           I_max = I_cs[i].at<float>(r,c);
00377         if (RG_cs[i].at<float>(r,c) > C_max)
00378           C_max = RG_cs[i].at<float>(r,c);
00379         if (BY_cs[i].at<float>(r,c) > C_max)
00380           C_max = BY_cs[i].at<float>(r,c);
00381 
00382         for (int a = 0; a < N_; ++a)
00383         {
00384           if (O_theta_cs[a][i].at<float>(r,c) > O_theta_max[a])
00385           {
00386             O_theta_max[a] = O_theta_cs[a][i].at<float>(r,c);
00387           }
00388         }
00389       }
00390     }
00391   }
00392 
00393   // Perform the normalization
00394   for (unsigned int i = 0; i < I_cs.size(); ++i)
00395   {
00396     // Test for max value for normalization
00397     I_cs[i] = normalize(I_cs[i], I_max);
00398     cv::Mat rg_norm = normalize(RG_cs[i], C_max);
00399     cv::Mat by_norm = normalize(BY_cs[i], C_max);
00400     cv::Mat c_csi = rg_norm + by_norm;
00401     C_cs.push_back(c_csi);
00402     for (int a = 0; a < N_; ++a)
00403     {
00404       O_theta_cs[a][i] = normalize(O_theta_cs[a][i], O_theta_max[a]);
00405     }
00406   }
00407 
00408   // Combine conspicuity maps into feature maps
00409   Mat I_bar;
00410   Mat C_bar;
00411   vector<Mat> O_bars;
00412   I_bar = mapSum(I_cs);
00413   C_bar = mapSum(C_cs);
00414 
00415   // For the orientations we sum each orientation separately, then combined
00416   // normalized forms of those four
00417   for (int a = 0; a < N_; ++a)
00418   {
00419     Mat O_bar_a = mapSum(O_theta_cs[a]);
00420     O_bars.push_back(O_bar_a);
00421   }
00422 
00423   //
00424   // Normalize and sum the O_bars
00425   //
00426   Mat O_bar(I_bar.rows, I_bar.cols, I_bar.type(), 0.0);
00427   float O_max = 0;
00428 
00429   // Get the max of the o_bars
00430   for (int a = 0; a < N_; ++a)
00431   {
00432     for(int r = 0; r < O_bars[a].rows; ++r)
00433     {
00434       for(int c = 0; c < O_bars[a].cols; ++c)
00435       {
00436         if (O_bars[a].at<float>(r,c) > O_max)
00437           O_max = O_bars[a].at<float>(r,c);
00438       }
00439     }
00440   }
00441   for (int a = 0; a < N_; ++a)
00442   {
00443     cv::Mat norm_o_bar = normalize(O_bars[a], O_max);
00444     O_bar += norm_o_bar;
00445   }
00446 
00447   //
00448   // Normalize and sum the different feature channels
00449   //
00450 
00451   // Get the max of the different feature channel conspicuicy maps
00452   float bar_max = 0;
00453   for (int r = 0; r < I_bar.rows; ++r)
00454   {
00455     for (int c = 0; c < I_bar.cols; ++c)
00456     {
00457       if (I_bar.at<float>(r,c) > bar_max)
00458         bar_max = I_bar.at<float>(r,c);
00459       if (C_bar.at<float>(r,c) > bar_max)
00460         bar_max = C_bar.at<float>(r,c);
00461       if (O_bar.at<float>(r,c) > bar_max)
00462         bar_max = O_bar.at<float>(r,c);
00463     }
00464   }
00465 
00466   // Build the saliency map as the combination of the feature maps
00467   Mat saliency_map(I_bar.rows, I_bar.cols, CV_32FC1);
00468   Mat I_bar_norm = normalize(I_bar, bar_max);
00469   Mat C_bar_norm = normalize(C_bar, bar_max);
00470   Mat O_bar_norm = normalize(O_bar, bar_max);
00471   saliency_map = (I_bar_norm*(1/3.0)+C_bar_norm*(1/3.0)+O_bar_norm*(1/3.0));
00472 
00473 #ifdef DISPLAY_SALIENCY_MAPS
00474   Mat display_I_bar = upSampleResponse(I_bar, min_delta_, frame.size());
00475   Mat display_C_bar = upSampleResponse(C_bar, min_delta_, frame.size());
00476   Mat display_O_bar = upSampleResponse(O_bar, min_delta_, frame.size());
00477   cv::imshow("C bar small", C_bar);
00478   // cv::imshow("I bar", display_I_bar);
00479   cv::imshow("C bar", display_C_bar);
00480   cv::imshow("O bar", display_O_bar);
00481   double min_val = 0;
00482   double max_val = 0;
00483   cv::minMaxLoc(display_O_bar, &min_val, &max_val);
00484   std::cout << "O_bar Min_val: " << min_val << std::endl;
00485   std::cout << "O_bar Max_val: " << max_val << std::endl;
00486   cv::minMaxLoc(display_C_bar, &min_val, &max_val);
00487   std::cout << "C_bar Min_val: " << min_val << std::endl;
00488   std::cout << "C_bar Max_val: " << max_val << std::endl;
00489   cv::minMaxLoc(display_I_bar, &min_val, &max_val);
00490   std::cout << "I_bar Min_val: " << min_val << std::endl;
00491   std::cout << "I_bar Max_val: " << max_val << std::endl;
00492   cv::minMaxLoc(saliency_map, &min_val, &max_val);
00493   std::cout << "sal_map Min_val: " << min_val << std::endl;
00494   std::cout << "sal_map Max_val: " << max_val << std::endl;
00495 #endif // DISPLAY_SALIENCY_MAPS
00496   return saliency_map;
00497 }
00498 
00499 Mat CenterSurroundMapper::getIntensityMap(Mat& frame)
00500 {
00501   Mat I(frame.rows, frame.cols, CV_32FC1);
00502   vector<Mat> I_scales;
00503 
00504   // Get hue indepenednt intensity channel if frame is a color image
00505   if(frame.channels() == 3)
00506   {
00507     // Get the component color channels
00508     vector<Mat> channels;
00509     Mat to_split(frame.rows, frame.cols, CV_32FC3);
00510     frame.convertTo(to_split, to_split.type(), 1.0/255.0);
00511     split(to_split, channels);
00512     I = channels[R_INDEX]/3.0 + channels[G_INDEX]/3.0 + channels[B_INDEX]/3.0;
00513   }
00514   else
00515   {
00516     frame.convertTo(I, CV_32FC1);
00517   }
00518 
00519   float max_i = 0;
00520   for(int r = 0; r < I.rows; ++r)
00521   {
00522     for(int c = 0; c < I.cols; ++c)
00523     {
00524       if (I.at<float>(r,c) > max_i)
00525         max_i = I.at<float>(r,c);
00526     }
00527   }
00528 
00529   // Get copies of the feature maps at all scales
00530   cv::buildPyramid(I, I_scales, num_scales_);
00531 
00532   //
00533   // Get multi-scale maps of the features
00534   //
00535   vector<Mat> I_cs;
00536 
00537   // Here we build the multi-scale maps
00538   for (int c = min_c_; c <= max_c_; c++)
00539   {
00540     for (int s = c + min_delta_; s <= c + max_delta_; s++)
00541     {
00542       // Build intensity maps
00543       I_cs.push_back(mapDifference(I_scales[c], I_scales[s], c, s));
00544     }
00545   }
00546 
00547   // Find max values to normalize maps by feature type
00548   float I_max = 0;
00549   for (unsigned int i = 0; i < I_cs.size(); ++i)
00550   {
00551     for(int r = 0; r < I_cs[i].rows; ++r)
00552     {
00553       for(int c = 0; c < I_cs[i].cols; ++c)
00554       {
00555         if (I_cs[i].at<float>(r,c) > I_max)
00556           I_max = I_cs[i].at<float>(r,c);
00557       }
00558     }
00559   }
00560 
00561   // Perform the normalization
00562   for (unsigned int i = 0; i < I_cs.size(); ++i)
00563   {
00564     // Test for max value for normalization
00565     I_cs[i] = normalize(I_cs[i], I_max);
00566   }
00567 
00568   // Combine conspicuity maps into feature maps
00569   Mat I_bar;
00570   I_bar = mapSum(I_cs);
00571   return I_bar;
00572 }
00573 
00574 Mat CenterSurroundMapper::getOrientationMap(Mat& frame)
00575 {
00576   vector<vector<Mat> > O_n_theta; // First index scale, second index relative
00577                                   // orientation (default 0 - 3)
00578   Mat I(frame.rows, frame.cols, CV_32FC1);
00579   vector<Mat> I_scales;
00580 
00581   // Get hue indepenednt intensity channel if frame is a color image
00582   if(frame.channels() == 3)
00583   {
00584     // Get the component color channels
00585     vector<Mat> channels;
00586     Mat to_split(frame.rows, frame.cols, CV_32FC3);
00587     frame.convertTo(to_split, to_split.type(), 1.0/255.0);
00588     split(to_split, channels);
00589     I = channels[R_INDEX]/3.0 + channels[G_INDEX]/3.0 + channels[B_INDEX]/3.0;
00590   }
00591   else
00592   {
00593     frame.convertTo(I, CV_32FC1);
00594   }
00595 
00596   float max_i = 0;
00597   for(int r = 0; r < I.rows; ++r)
00598   {
00599     for(int c = 0; c < I.cols; ++c)
00600     {
00601       if (I.at<float>(r,c) > max_i)
00602         max_i = I.at<float>(r,c);
00603     }
00604   }
00605 
00606   // Get copies of the feature maps at all scales
00607   cv::buildPyramid(I, I_scales, num_scales_);
00608 
00609   // Get laplacians at all scales (TODO: do this while building pyr via DoG)
00610   for(unsigned int i = 0; i < I_scales.size(); ++i)
00611   {
00612     Mat lap;
00613     cv::Laplacian(I_scales[i], lap, I_scales[i].depth());
00614     vector<Mat> O_theta;
00615 
00616     // Get the N orientation maps for each scale
00617     for (int a = 0; a < N_; a++)
00618     {
00619       Mat m_a(lap.rows, lap.cols, lap.type());
00620       cv::filter2D(lap, m_a, -1, gabor_filters_[a]);
00621 
00622       // For each of the N orientation maps smooth and downsample
00623       O_theta.push_back(m_a);
00624 
00625 #ifdef DISPLAY_GABOR_FILTERS
00626       std::stringstream m_name;
00627       m_name << "O_" << i << "_" << a;
00628       cv::imshow(m_name.str(), m_a);
00629       cv::waitKey();
00630 #endif // DISPLAY_GABOR_FILTERS
00631 
00632     }
00633     O_n_theta.push_back(O_theta);
00634   }
00635 
00636   //
00637   // Get multi-scale maps of the features
00638   //
00639 
00640   vector<vector<Mat> > O_theta_cs;
00641 
00642   // Initialize the vectors for the different orientations
00643   for (int a = 0; a < N_; ++a)
00644   {
00645     vector<Mat> O_cs;
00646     O_cs.clear();
00647     O_theta_cs.push_back(O_cs);
00648   }
00649 
00650   // Here we build the multi-scale maps
00651   for (int c = min_c_; c <= max_c_; c++)
00652   {
00653     for (int s = c + min_delta_; s <= c + max_delta_; s++)
00654     {
00655       // Build orientation maps for each of the N_ orientations
00656       for (int a = 0; a < N_; ++a)
00657       {
00658         O_theta_cs[a].push_back(mapDifference(O_n_theta[c][a],
00659                                               O_n_theta[s][a], c, s));
00660       }
00661     }
00662   }
00663 
00664   // Find max values to normalize maps by feature type
00665   vector<float> O_theta_max(N_, 0);
00666   for (unsigned int i = 0; i < O_theta_cs[0].size(); ++i)
00667   {
00668     for(int r = 0; r < O_theta_cs[0][i].rows; ++r)
00669     {
00670       for(int c = 0; c < O_theta_cs[0][i].cols; ++c)
00671       {
00672         for (int a = 0; a < N_; ++a)
00673         {
00674           if (O_theta_cs[a][i].at<float>(r,c) > O_theta_max[a])
00675             O_theta_max[a] = O_theta_cs[a][i].at<float>(r,c);
00676         }
00677       }
00678     }
00679   }
00680 
00681   // Perform the normalization
00682   for (unsigned int i = 0; i < O_theta_cs[0].size(); ++i)
00683   {
00684     for (int a = 0; a < N_; ++a)
00685     {
00686       O_theta_cs[a][i] = normalize(O_theta_cs[a][i], O_theta_max[a]);
00687     }
00688   }
00689 
00690   // Combine conspicuity maps into feature maps
00691   vector<Mat> O_bars;
00692 
00693   // For the orientations we sum each orientation separately, then combined
00694   // normalized forms of those four
00695   for (int a = 0; a < N_; ++a)
00696   {
00697     Mat O_bar_a = mapSum(O_theta_cs[a]);
00698     O_bars.push_back(O_bar_a);
00699   }
00700 
00701   //
00702   // Normalize and sum the O_bars
00703   //
00704   int min_rows = 10000;
00705   int min_cols = 10000;
00706 
00707   // Find the smallest scale of the images
00708   for (unsigned int i = 0; i < O_bars.size(); ++i)
00709   {
00710     if (O_bars[i].rows < min_rows)
00711     {
00712       min_rows = O_bars[i].rows;
00713       min_cols = O_bars[i].cols;
00714     }
00715   }
00716 
00717   Mat O_bar = Mat::zeros(min_rows, min_cols, CV_32FC1);
00718   float O_max = 0;
00719 
00720   // Get the max of the o_bars
00721   for (int a = 0; a < N_; ++a)
00722   {
00723     for(int r = 0; r < O_bars[a].rows; ++r)
00724     {
00725       for(int c = 0; c < O_bars[a].cols; ++c)
00726       {
00727         if (O_bars[a].at<float>(r,c) > O_max)
00728           O_max = O_bars[a].at<float>(r,c);
00729       }
00730     }
00731   }
00732 
00733   for (int a = 0; a < N_; ++a)
00734   {
00735     O_bar += normalize(O_bars[a], O_max);
00736   }
00737   return O_bar;
00738 }
00739 
00740 Mat CenterSurroundMapper::getColorMap(Mat& frame)
00741 {
00742   Mat R(frame.rows, frame.cols, CV_32FC1);
00743   Mat G(frame.rows, frame.cols, CV_32FC1);
00744   Mat B(frame.rows, frame.cols, CV_32FC1);
00745   Mat Y(frame.rows, frame.cols, CV_32FC1);
00746   Mat I(frame.rows, frame.cols, CV_32FC1);
00747   vector<Mat> R_scales;
00748   vector<Mat> G_scales;
00749   vector<Mat> B_scales;
00750   vector<Mat> Y_scales;
00751 
00752   // Get the component color channels
00753   vector<Mat> channels;
00754   Mat to_split(frame.rows, frame.cols, CV_32FC3);
00755   frame.convertTo(to_split, to_split.type(), 1.0/255.0);
00756 
00757   split(to_split, channels);
00758 
00759   // Get hue indepenednt intensity channel
00760   I = channels[R_INDEX]/3.0 + channels[G_INDEX]/3.0 + channels[B_INDEX]/3.0;
00761 
00762   float max_i = 0;
00763   for(int r = 0; r < I.rows; ++r)
00764   {
00765     for(int c = 0; c < I.cols; ++c)
00766     {
00767       if (I.at<float>(r,c) > max_i)
00768         max_i = I.at<float>(r,c);
00769     }
00770   }
00771 
00772   // Normalize r, g, and b at locations greater than 1/10 the image
00773   for(int r = 0; r < I.rows; ++r)
00774   {
00775     for(int c = 0; c < I.cols; ++c)
00776     {
00777       if (I.at<float>(r,c) > max_i*0.10)
00778       {
00779         float intensity = I.at<float>(r,c);
00780         channels[R_INDEX].at<float>(r,c) = (channels[R_INDEX].at<float>(r,c) /
00781                                             intensity);
00782         channels[G_INDEX].at<float>(r,c) = (channels[G_INDEX].at<float>(r,c) /
00783                                             intensity);
00784         channels[B_INDEX].at<float>(r,c) = (channels[B_INDEX].at<float>(r,c) /
00785                                             intensity);
00786       }
00787       else
00788       {
00789         channels[0].at<float>(r,c) = 0;
00790         channels[1].at<float>(r,c) = 0;
00791         channels[2].at<float>(r,c) = 0;
00792       }
00793     }
00794   }
00795 
00796   // Get intensity independent hue channels
00797   R = channels[R_INDEX] - (channels[G_INDEX]/2.0 + channels[B_INDEX]/2.0);
00798   G = channels[G_INDEX] - (channels[R_INDEX]/2.0 + channels[B_INDEX]/2.0);
00799   B = channels[B_INDEX] - (channels[R_INDEX]/2.0 + channels[G_INDEX]/2.0);
00800 
00801   Mat Y_abs(Y.rows, Y.cols, Y.type());
00802   for(int r = 0; r < Y.rows; ++r)
00803   {
00804     for(int c = 0; c < Y.cols; ++c)
00805     {
00806       Y_abs.at<float>(r,c) = std::abs(channels[R_INDEX].at<float>(r,c)/2.0 -
00807                                       channels[G_INDEX].at<float>(r,c)/2.0);
00808     }
00809   }
00810 
00811   Y = (channels[R_INDEX]/2.0 + channels[G_INDEX]/2.0) - Y_abs-channels[B_INDEX];
00812 
00813 #ifdef DISPLAY_COLOR_MAPS
00814   cv::imshow("I",I);
00815   cv::imshow("R",R);
00816   cv::imshow("G",G);
00817   cv::imshow("B",B);
00818   cv::imshow("Y",Y);
00819   cv::waitKey();
00820 #endif // DISPLAY_COLOR_MAPS
00821 
00822   cv::buildPyramid(R, R_scales, num_scales_);
00823   cv::buildPyramid(G, G_scales, num_scales_);
00824   cv::buildPyramid(B, B_scales, num_scales_);
00825   cv::buildPyramid(Y, Y_scales, num_scales_);
00826 
00827   //
00828   // Get multi-scale maps of the features
00829   //
00830   vector<Mat> RG_cs;
00831   vector<Mat> BY_cs;
00832   vector<Mat> C_cs;
00833 
00834   // Here we build the multi-scale maps
00835   for (int c = min_c_; c <= max_c_; c++)
00836   {
00837     for (int s = c + min_delta_; s <= c + max_delta_; s++)
00838     {
00839       // Build red-green opponent color maps
00840       Mat RG_c = R_scales[c] - G_scales[c];
00841       Mat GR_s = G_scales[s] - R_scales[s];
00842       RG_cs.push_back(mapDifference(RG_c, GR_s, c, s));
00843 
00844       // Build blue-yellow opponent color maps
00845       Mat BY_c = B_scales[c] - Y_scales[c];
00846       Mat YB_s = Y_scales[s] - B_scales[s];
00847       BY_cs.push_back(mapDifference(BY_c, YB_s, c, s));
00848     }
00849   }
00850 
00851   //
00852   // Normalize all maps, based on feature type
00853   //
00854 
00855   // Find max values to normalize maps by feature type
00856   float C_max = 0;
00857   for (unsigned int i = 0; i < RG_cs.size(); ++i)
00858   {
00859     for(int r = 0; r < RG_cs[i].rows; ++r)
00860     {
00861       for(int c = 0; c < RG_cs[i].cols; ++c)
00862       {
00863         if (RG_cs[i].at<float>(r,c) > C_max)
00864           C_max = RG_cs[i].at<float>(r,c);
00865         if (BY_cs[i].at<float>(r,c) > C_max)
00866           C_max = BY_cs[i].at<float>(r,c);
00867       }
00868     }
00869   }
00870 
00871   // Perform the normalization
00872   for (unsigned int i = 0; i < RG_cs.size(); ++i)
00873   {
00874     // Test for max value for normalization
00875     C_cs.push_back(normalize(RG_cs[i], C_max) + normalize(BY_cs[i], C_max));
00876   }
00877 
00878   Mat C_bar;
00879   C_bar = mapSum(C_cs);
00880   return C_bar;
00881 }
00882 
00883 //
00884 // Helper methods
00885 //
00886 Mat CenterSurroundMapper::mapDifference(Mat& m_c, Mat& m_s, int c, int s)
00887 {
00888   // Upsample c to s resolution
00889   Mat m_s_prime(m_s.rows, m_s.cols, CV_32FC1);
00890   Mat temp;
00891   m_s.copyTo(temp);
00892 
00893   // Calculate the correct sizes to up sample to
00894   vector<cv::Size> sizes;
00895   cv::Size current_size(m_c.cols, m_c.rows);
00896   for (int i = c; i < s; i++)
00897   {
00898     sizes.push_back(current_size);
00899     current_size.width /= 2;
00900     current_size.height /= 2;
00901   }
00902 
00903   for (int i = c; i < s; i++)
00904   {
00905     cv::Size up_size;
00906     up_size = sizes.back();
00907     sizes.pop_back();
00908     cv::pyrUp(temp, m_s_prime, up_size);
00909     temp = m_s_prime;
00910   }
00911 
00912   // Take pixelwise difference
00913   Mat diff = abs(m_c - m_s_prime);
00914 
00915   return diff;
00916 }
00917 
00918 Mat CenterSurroundMapper::mapSum(vector<Mat>& maps)
00919 {
00920   int min_rows = 10000;
00921   int min_cols = 10000;
00922 
00923   // Find the smallest scale of the images
00924   for (unsigned int i = 0; i < maps.size(); ++i)
00925   {
00926     if (maps[i].rows < min_rows)
00927     {
00928       min_rows = maps[i].rows;
00929       min_cols = maps[i].cols;
00930     }
00931   }
00932 
00933   Mat sum = Mat::zeros(min_rows, min_cols, CV_32FC1);
00934 
00935   for (unsigned int i = 0; i < maps.size(); ++i)
00936   {
00937     Mat m_prime = maps[i];
00938     Mat temp = maps[i];
00939 
00940     while (temp.cols > min_cols)
00941     {
00942       cv::pyrDown(temp, m_prime);
00943       temp = m_prime;
00944     }
00945     sum += m_prime;
00946   }
00947 
00948   return sum;
00949 }
00950 
00951 Mat CenterSurroundMapper::normalize(Mat& map, float M)
00952 {
00953   // Normalize the values to the range 0 to max_val...
00954   float cur_max_val = 0;
00955   for(int r = 0; r < map.rows; ++r)
00956   {
00957     for (int c = 0; c < map.cols; ++c)
00958     {
00959       if (map.at<float>(r,c) > cur_max_val)
00960         cur_max_val = map.at<float>(r,c);
00961     }
00962   }
00963 
00964   if (cur_max_val == 0)
00965   {
00966     std::cout << "Max value is 0, don't want to get nans..." << std::endl;
00967     std::cout << "M value is: " << M << std::endl;
00968     cur_max_val = M;
00969     // cur_max_val = 1;
00970   }
00971 
00972   Mat normalized = map * (M/float(cur_max_val));
00973   float thresh = M*0.20;
00974   // float thresh = 0;
00975 
00976   // Find the local maxima
00977   vector<float> maxima;
00978   int nan_count = 0;
00979   for(int r = 0; r < map.rows; ++r)
00980   {
00981     for(int c = 0; c < map.cols; ++c)
00982     {
00983       float val = map.at<float>(r,c);
00984       if (isnan(val))
00985       {
00986         // std::cout << "have nan value" << std::endl;
00987         nan_count++;
00988         val = 0;
00989         map.at<float>(r,c) = 0;
00990       }
00991       if (val > thresh)
00992       {
00993         // Test if maximal over the 3 by 3 window
00994         // TODO: Make this test uneven (to remove ties or over-supression)
00995         if (c > 0) // Has left
00996         {
00997           if( val < map.at<float>(r,c-1))
00998             continue;
00999           if (r > 0 && val < map.at<float>(r - 1, c - 1))
01000             continue;
01001           if (r < map.rows - 1 && val < map.at<float>(r + 1, c - 1))
01002             continue;
01003         }
01004         if (c < map.cols - 1) // Has Right
01005         {
01006           if( val < map.at<float>(r,c + 1))
01007             continue;
01008           if (r > 0 && val < map.at<float>(r - 1, c + 1))
01009             continue;
01010           if (r < map.rows - 1 && val < map.at<float>(r + 1, c + 1))
01011             continue;
01012         }
01013         if (r > 0 && val < map.at<float>(r - 1, c)) // Has above
01014           continue;
01015         if (r < map.rows - 1 && val < map.at<float>(r + 1, c)) // Has below
01016           continue;
01017 
01018         // Store the local maxima value
01019         maxima.push_back(val);
01020       }
01021     }
01022   }
01023 
01024   if (nan_count > 0)
01025   {
01026     std::cout << "Image has: " << nan_count << " nans" << std::endl;
01027   }
01028   // Get mean of the local maxima
01029   float m_bar = 0;
01030 
01031   for (unsigned int i = 0; i < maxima.size(); ++i)
01032   {
01033     m_bar += maxima[i];
01034   }
01035 
01036   m_bar /= maxima.size();
01037 
01038   // Finally perform the normalization based on the difference between the
01039   // maximum value and the mean of the local maxima
01040   float fact = (M - m_bar)*(M - m_bar);
01041 
01042   // Scale it to be within the image range
01043   fact *= 1.0/(fact*M);
01044   normalized *= fact;
01045   return normalized;
01046 }
01047 
01048 Mat CenterSurroundMapper::scaleMap(Mat saliency_map)
01049 {
01050   Mat scaled;
01051   Mat saliency_int(saliency_map.rows, saliency_map.cols, CV_8UC1);
01052   saliency_map *= 255;
01053   saliency_map.convertTo(saliency_int, saliency_int.type());
01054   cv::equalizeHist(saliency_int, scaled);
01055   return scaled;
01056 }
01057 
01058 cv::Mat CenterSurroundMapper::upSampleResponse(cv::Mat& m_s, int s, cv::Size size0)
01059 {
01060   // Upsample from scale s to scale 0
01061   cv::Mat m_s_prime = m_s.clone();
01062   cv::Mat temp = m_s.clone();
01063 
01064   // Calculate the correct sizes to up sample to
01065   std::vector<cv::Size> sizes;
01066   cv::Size current_size = size0;
01067   while (true)
01068   {
01069     sizes.push_back(current_size);
01070     current_size.width /= 2;
01071     current_size.height /= 2;
01072     if (current_size.width == m_s.cols)
01073     {
01074       break;
01075     }
01076     else if (current_size.width-1 == m_s.cols)
01077     {
01078       std::cout << "Breaking on -1" << std::endl;
01079     }
01080     else if (current_size.width+1 == m_s.cols)
01081     {
01082       std::cout << "Breaking on +1" << std::endl;
01083     }
01084   }
01085   int num_ups = sizes.size();
01086   for (int i = 0; i < num_ups; i++)
01087   {
01088     cv::Size up_size;
01089     up_size = sizes.back();
01090     sizes.pop_back();
01091     cv::pyrUp(temp, m_s_prime, up_size);
01092     temp = m_s_prime;
01093   }
01094 
01095   return m_s_prime;
01096 }
01097 
01098 }


cpl_visual_features
Author(s): Tucker Hermans
autogenerated on Wed Nov 27 2013 11:52:35