00001 /*!@file SIFT/FeatureVector.C Feature vector for SIFT obj recognition */ 00002 00003 // //////////////////////////////////////////////////////////////////// // 00004 // The iLab Neuromorphic Vision C++ Toolkit - Copyright (C) 2001 by the // 00005 // University of Southern California (USC) and the iLab at USC. // 00006 // See http://iLab.usc.edu for information about this project. // 00007 // //////////////////////////////////////////////////////////////////// // 00008 // Major portions of the iLab Neuromorphic Vision Toolkit are protected // 00009 // under the U.S. patent ``Computation of Intrinsic Perceptual Saliency // 00010 // in Visual Environments, and Applications'' by Christof Koch and // 00011 // Laurent Itti, California Institute of Technology, 2001 (patent // 00012 // pending; application number 09/912,225 filed July 23, 2001; see // 00013 // http://pair.uspto.gov/cgi-bin/final/home.pl for current status). // 00014 // //////////////////////////////////////////////////////////////////// // 00015 // This file is part of the iLab Neuromorphic Vision C++ Toolkit. // 00016 // // 00017 // The iLab Neuromorphic Vision C++ Toolkit is free software; you can // 00018 // redistribute it and/or modify it under the terms of the GNU General // 00019 // Public License as published by the Free Software Foundation; either // 00020 // version 2 of the License, or (at your option) any later version. // 00021 // // 00022 // The iLab Neuromorphic Vision C++ Toolkit is distributed in the hope // 00023 // that it will be useful, but WITHOUT ANY WARRANTY; without even the // 00024 // implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // 00025 // PURPOSE. See the GNU General Public License for more details. // 00026 // // 00027 // You should have received a copy of the GNU General Public License // 00028 // along with the iLab Neuromorphic Vision C++ Toolkit; if not, write // 00029 // to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, // 00030 // Boston, MA 02111-1307 USA. // 00031 // //////////////////////////////////////////////////////////////////// // 00032 // 00033 // Primary maintainer for this file: James Bonaiuto <bonaiuto@usc.edu> 00034 // $HeadURL: svn://isvn.usc.edu/software/invt/trunk/saliency/src/SIFT/FeatureVector.C $ 00035 // $Id: FeatureVector.C 9412 2008-03-10 23:10:15Z farhan $ 00036 // 00037 00038 #include "SIFT/FeatureVector.H" 00039 #include "Util/Assert.H" 00040 #include "Util/Promotions.H" // for clamped_convert<T>() 00041 #include <cmath> 00042 #include "Image/DrawOps.H" 00043 00044 // ###################################################################### 00045 FeatureVector::FeatureVector(int xSize, int ySize, int zSize, bool wrapZ) : 00046 itsFV(xSize*ySize*zSize, 0.0F), 00047 itsXSize(xSize), itsYSize(ySize), itsZSize(zSize), itsWrapZ(wrapZ) 00048 { } 00049 00050 // ###################################################################### 00051 FeatureVector::~FeatureVector() 00052 { } 00053 00054 // ######################################################################## 00055 void FeatureVector::addValue(const float x, const float y, 00056 const float z, const float value) 00057 { 00058 int xi0, xi1, yi0, yi1, zi0, zi1; // bins 00059 float wx0, wy0, wz0, wx1, wy1, wz1; // corresponding weights 00060 00061 // if close to bounds then the values go fully into the end bins, 00062 // otherwise they split between two adjacent bins. Note: a value of 00063 // 2.0 should equally split between bins 1 and 2: 00064 if (x <= 0.5F) 00065 { xi0 = 0; xi1 = 0; wx0 = 0.5F; wx1 = 0.5F; } 00066 else if (x >= (float)itsXSize-0.5F) 00067 { xi0 = itsXSize-1; xi1 = itsXSize-1; wx0 = 0.5F; wx1 = 0.5F; } 00068 else 00069 { 00070 const float xx = x - 0.5F; 00071 xi0 = int(xx); xi1 = xi0 + 1; 00072 wx1 = xx - float(xi0); wx0 = 1.0F - wx1; 00073 } 00074 00075 if (y <= 0.5F) 00076 { yi0 = 0; yi1 = 0; wy0 = 0.5F; wy1 = 0.5F; } 00077 else if (y >= (float)itsYSize-0.5F) 00078 { yi0 = itsYSize-1; yi1 = itsYSize-1; wy0 = 0.5F; wy1 = 0.5F; } 00079 else 00080 { 00081 const float yy = y - 0.5F; 00082 yi0 = int(yy); yi1 = yi0 + 1; 00083 wy1 = yy - float(yi0); wy0 = 1.0F - wy1; 00084 } 00085 00086 00087 // the situation is different for orientation as we wrap around: 00088 // orientation are now labeld 'z' for more general purpose 00089 00090 if (itsWrapZ){ 00091 //Wrap the Z around the itsZSize 00092 if (z <= 0.5F) 00093 { 00094 zi0 = 0; zi1 = itsZSize-1; 00095 wz0 = 0.5F + z; wz1 = 1.0F - wz0; 00096 } 00097 else if (z >= itsZSize-0.5F) 00098 { 00099 zi0 = itsZSize-1; zi1 = 0; 00100 wz0 = ((float)itsZSize+0.5F) - z; wz1 = 1.0F - wz0; 00101 } 00102 else 00103 { 00104 const float zz = z - 0.5F; 00105 zi0 = int(zz); zi1 = zi0 + 1; 00106 wz1 = zz - float(zi0); wz0 = 1.0F - wz1; 00107 } 00108 } else { 00109 //Dont wrap z bin 00110 if (z <= 0.5F) 00111 { 00112 zi0 = 0; zi1 = 0; 00113 wz0 = 0.5F; wz1 =0.5F; 00114 } 00115 else if (z >= (float)itsZSize-0.5F) 00116 { 00117 zi0 = itsZSize-1; zi1 = itsZSize-1; 00118 wz0 = 0.5F; wz1 = 0.5F; 00119 } 00120 else 00121 { 00122 const float zz = z - 0.5F; 00123 zi0 = int(zz); zi1 = zi0 + 1; 00124 wz1 = zz - float(zi0); wz0 = 1.0F - wz1; 00125 } 00126 } 00127 00128 // convention: we add 1 for each unit of o (our fastest varying 00129 // index), then zSize for each unit of y, finally zSize*ySize for each unit of 00130 // x. Let's populate our 8 bins: 00131 00132 //No more opt calc unless we informace the size of the bins 00133 //Hopfully the compiler will optemize powers of 2 00134 //xi0 <<= 5; xi1 <<= 5; yi0 <<= 3; yi1 <<= 3; 00135 00136 xi0 *= itsZSize*itsYSize; xi1 *= itsZSize*itsYSize; 00137 yi0 *= itsZSize; yi1 *= itsZSize; 00138 00139 00140 itsFV[xi0 + yi0 + zi0] += value * wx0 * wy0 * wz0; 00141 itsFV[xi1 + yi0 + zi0] += value * wx1 * wy0 * wz0; 00142 itsFV[xi0 + yi1 + zi0] += value * wx0 * wy1 * wz0; 00143 itsFV[xi1 + yi1 + zi0] += value * wx1 * wy1 * wz0; 00144 itsFV[xi0 + yi0 + zi1] += value * wx0 * wy0 * wz1; 00145 itsFV[xi1 + yi0 + zi1] += value * wx1 * wy0 * wz1; 00146 itsFV[xi0 + yi1 + zi1] += value * wx0 * wy1 * wz1; 00147 itsFV[xi1 + yi1 + zi1] += value * wx1 * wy1 * wz1; 00148 00149 //LINFO("%f,%f,%f -> %d-%d(%f) %d-%d(%f) %d-%d(%f)", 00150 // x,y,o,xi0,xi1,wx0,yi0,yi1,wy0,oi0,oi1,wo0); 00151 } 00152 00153 // ######################################################################## 00154 void FeatureVector::normalize() 00155 { 00156 std::vector<float>::iterator ptr = itsFV.begin(), stop = itsFV.end(); 00157 00158 // compute sum of squares: 00159 float sq = 0.0f; 00160 while(ptr != stop) { sq += (*ptr) * (*ptr); ++ ptr; } 00161 00162 // if too small to normalize, forget it: 00163 if (sq < 1.0e-10) return; 00164 00165 // compute normalization factor: 00166 sq = 1.0F / sqrtf(sq); 00167 00168 // normalize it: 00169 ptr = itsFV.begin(); 00170 while(ptr != stop) *ptr++ *= sq; 00171 } 00172 00173 // ######################################################################## 00174 void FeatureVector::threshold(const float limit) 00175 { 00176 bool changed = false; 00177 00178 std::vector<float>::iterator ptr = itsFV.begin(), stop = itsFV.end(); 00179 00180 while(ptr != stop) 00181 { 00182 if (*ptr > limit) { *ptr = limit; changed = true; } 00183 ++ ptr; 00184 } 00185 00186 if (changed) normalize(); 00187 } 00188 00189 // ######################################################################## 00190 float FeatureVector::getValue(const uint index) const 00191 { 00192 ASSERT(index < itsFV.size()); 00193 return itsFV[index]; 00194 } 00195 00196 // ###################################################################### 00197 void FeatureVector::toByteKey(std::vector<byte>& dest, float thresh, bool norm) 00198 { 00199 dest.resize(itsFV.size()); 00200 00201 // do normalization and thresholding: 00202 00203 if (norm) normalize(); 00204 if (thresh > -1) threshold(thresh); //TODO: is -1 a good value to check for? 00205 if (norm) normalize(); 00206 00207 std::vector<float>::const_iterator src = itsFV.begin(), stop = itsFV.end(); 00208 std::vector<byte>::iterator dst = dest.begin(); 00209 00210 while (src != stop) *dst++ = clamped_convert<byte>((*src++) * 512.0F); 00211 } 00212 00213 // ###################################################################### 00214 Image<PixRGB<byte> > FeatureVector::getFeatureVectorImage(std::vector<byte> &fv) 00215 { 00216 00217 00218 int WIDTH = 256; 00219 int HEIGHT = 256; 00220 00221 Image<PixRGB<byte> > fvImg(WIDTH, HEIGHT, NO_INIT); 00222 fvImg.clear(PixRGB<byte>(255, 255, 255)); 00223 int xBins = int((float)WIDTH/4); 00224 int yBins = int((float)HEIGHT/4); 00225 00226 drawGrid(fvImg, xBins, yBins, 1, 1, PixRGB<byte>(0, 0, 0)); 00227 00228 for (int xx=0; xx<4; xx++){ 00229 for (int yy=0; yy<4; yy++){ 00230 for (int oo=0; oo<8; oo++){ 00231 Point2D<int> loc(xBins/2+(xBins*xx), yBins/2+(yBins*yy)); 00232 uint fv_i = xx*32+yy*8+oo; 00233 byte mag = 0; 00234 if (fv_i > fv.size()) 00235 LFATAL("Invalid feture vector"); 00236 else 00237 mag = fv[xx*32+yy*8+oo]/4; 00238 drawDisk(fvImg, loc, 2, PixRGB<byte>(255, 0, 0)); 00239 drawLine(fvImg, loc, 00240 Point2D<int>(int(loc.i + mag*cosf(oo*M_PI/4)), 00241 int(loc.j + mag*sinf(oo*M_PI/4))), 00242 PixRGB<byte>(255, 0, 0)); 00243 } 00244 } 00245 } 00246 00247 return fvImg; 00248 } 00249 00250 00251 double FeatureVector::getMag() 00252 { 00253 00254 double mag = 0; 00255 std::vector<float>::iterator ptr = itsFV.begin(), stop = itsFV.end(); 00256 00257 while(ptr != stop) 00258 { 00259 mag += (*ptr * *ptr); 00260 ++ ptr; 00261 } 00262 00263 return sqrt(mag); 00264 00265 } 00266 00267 00268 00269 00270 // ###################################################################### 00271 /* So things look consistent in everyone's emacs... */ 00272 /* Local Variables: */ 00273 /* indent-tabs-mode: nil */ 00274 /* End: */