MPIopenvision3.C

Go to the documentation of this file.
00001 /*!@file INVT/MPIopenvision3.C  MPI version of openvision3.C */
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: Laurent Itti <itti@usc.edu>
00034 // $HeadURL: svn://isvn.usc.edu/software/invt/trunk/saliency/src/INVT/MPIopenvision3.C $
00035 // $Id: MPIopenvision3.C 10845 2009-02-13 08:49:12Z itti $
00036 //
00037 
00038 #include "Channels/ChannelOpts.H"
00039 #include "Channels/ConvolveChannel.H"
00040 #include "Component/ModelManager.H"
00041 #include "Image/MathOps.H"
00042 #include "Image/MatrixOps.H"
00043 #include "Image/Pixels.H"
00044 #include "Media/FrameSeries.H"
00045 #include "Channels/RawVisualCortex.H"
00046 #include "Raster/Raster.H"
00047 #include "Util/sformat.H"
00048 
00049 #include <fstream>
00050 #ifdef HAVE_MPI_H
00051 #include <mpi.h>
00052 #endif
00053 
00054 #define ITMAX 200
00055 #define FTOL 0.000001f
00056 #define MASTER 0
00057 #define MSG_STOP 10
00058 #define MSG_PIC 1
00059 #define MSG_DATA 2
00060 #define MSG_RESULT 200
00061 #define PATH "/home/hpc-26/itti/openvision"
00062 #define INPUTFILE "/home/hpc-26/itti/openvision/start.txt"
00063 #define NB_PICS 10
00064 #define NB_FILTERS 3
00065 #define NCOEFFS 8
00066 #define MCOEFFS NCOEFFS * NCOEFFS * NB_FILTERS
00067 #define NB_EVALS 9
00068 #define NB_WAVES 5
00069 #define NB_PROCS NB_PICS * NB_EVALS * 2 + 1
00070 #define RAD 64
00071 
00072 int main(int argc, char **argv)
00073 {
00074 #ifndef HAVE_MPI_H
00075 
00076   LFATAL("<mpi.h> must be installed to use this program");
00077 
00078 #else
00079 
00080   // Set verbosity:
00081   int loglev = LOG_ERR;
00082   MYLOGVERB = loglev;
00083 
00084   // Number of coefficients:
00085   int n = NCOEFFS;
00086   int m = MCOEFFS;
00087 
00088   // Initialize the MPI system and get the number of processes:
00089   int myrank, size;
00090   MPI_Init(&argc, &argv);
00091   MPI_Comm_size(MPI_COMM_WORLD, &size);
00092 
00093   // Check for proper number of processes:
00094   if(size < NB_PROCS)
00095     {
00096       LERROR("*** Error: %i Processes needed. ***", NB_PROCS);
00097       MPI_Finalize();
00098       return 1;
00099     }
00100 
00101   // Get the rank of the process:
00102   MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
00103 
00104   if (myrank == MASTER)
00105     {
00106       // ################ MASTER PROCESS ################
00107 
00108       // **************** INITIALIZATION ****************
00109 
00110       MPI_Status status;
00111 
00112       // Load start coefficients:
00113       float p[m];
00114 //       std::ifstream inputfile (INPUTFILE);
00115 //       if (inputfile.is_open())
00116 //         {
00117 //           for (int j = 0; j < m; j++)
00118 //             inputfile >> p[j];
00119 //           inputfile.close();
00120 //         }
00121 //       else
00122 //         {
00123 //           LERROR("*** Cannot open input file !");
00124 //           return 1;
00125 //         }
00126       for (int j = 0; j < m; j++)
00127         p[j] = ((rand()%200) - 50) / 200.0f;
00128 
00129       // Generate the starting point, directions and function value:
00130       float xi[m][m];
00131       for(int i = 0; i < m; i++)
00132         {
00133           for(int j = 0; j < m; j++)
00134             xi[j][i] = 0.0f;
00135           xi[i][i] = 1.0f;
00136         }
00137       float fp = 10.0f;
00138 
00139       // **************** OPTIMIZATION ****************
00140 
00141       for (int iter = 0; ; iter++)
00142         {
00143           LERROR("### ITERATION %i ###", iter);
00144 
00145           int stop = 0;
00146           float del = 0.0f;
00147           int ibig = 0;
00148           float fstart = fp;
00149           float start[m];
00150           for (int i = 0; i < m; i++)
00151             start[m] = p[m];
00152           float fpt = fp;
00153           int mubig = 0;
00154 
00155           // Loop over all directions:
00156           for (int i = 0; i < m; i++)
00157             {
00158               mubig = 0;
00159               for (int wave = 0; wave < NB_WAVES; wave++)
00160                 {
00161                   int muplus = NB_EVALS * wave;
00162                   // Send data to slaves:
00163                   for (int part = 0; part < 2; part++)
00164                     for (int mu = 1; mu <= NB_EVALS; mu++)
00165                       {
00166                         int mu2 = (1 - 2 * part) * (mu + muplus);
00167                         float ptry[m];
00168                         for (int j = 0; j < m; j++)
00169                           ptry[j] = p[j] + (mu2 * 0.02f * xi[j][i]);
00170                         for (int pic = 0; pic < NB_PICS; pic++)
00171                           {
00172                             int proc = (NB_EVALS * NB_PICS * part);
00173                             proc += (NB_PICS * (mu - 1)) + pic + 1;
00174                             MPI_Send(&stop, 1, MPI_INT, proc,
00175                                      MSG_STOP, MPI_COMM_WORLD);
00176                             MPI_Send(&pic, 1, MPI_INT, proc,
00177                                      MSG_PIC, MPI_COMM_WORLD);
00178                             MPI_Send(ptry, m, MPI_FLOAT, proc,
00179                                      MSG_DATA, MPI_COMM_WORLD);
00180 
00181                           }
00182                       }
00183 
00184                   // Receive result from slaves and keep best decrease:
00185                   for (int part = 0; part < 2; part++)
00186                     for (int mu = 1; mu <= NB_EVALS; mu++)
00187                       {
00188                         float fptry = 0.0f;
00189                         for (int pic = 0; pic < NB_PICS; pic++)
00190                           {
00191                             float value;
00192                             int proc = (NB_EVALS * NB_PICS * part);
00193                             proc += (NB_PICS * (mu - 1)) + pic + 1;
00194                             MPI_Recv(&value, 1, MPI_FLOAT, proc,
00195                                      MSG_RESULT, MPI_COMM_WORLD,
00196                                      &status);
00197                             fptry += value;
00198                           }
00199                         fptry /= NB_PICS;
00200                         if (fptry < fpt)
00201                           {
00202                             fpt = fptry;
00203                             mubig = (1 - 2 * part) * (mu + muplus);
00204                           }
00205                       }
00206                 }
00207 
00208               // Keep the greatest decrease:
00209               if (fp - fpt > del)
00210                 {
00211                   del = fp - fpt;
00212                   ibig = i;
00213                 }
00214               fp = fpt;
00215               for (int j = 0; j < m; j++)
00216                 p[j] += mubig * 0.02f * xi[j][i];
00217 
00218               LERROR("*** Direction %i : mu = %i", i, mubig);
00219               LERROR("*** New value : %f", fp);
00220 
00221 //               float maxcoeff1 = 0;
00222 //               float maxcoeff2 = 0;
00223 //               float maxcoeff3 = 0;
00224 //               for (int j = 0; j < m / 3; j++)
00225 //                 {
00226 //                   if (fabs(p[j]) > maxcoeff1)
00227 //                     maxcoeff1 = fabs(p[j]);
00228 //                   if (fabs(p[j + m / 3]) > maxcoeff2)
00229 //                     maxcoeff2 = fabs(p[j + m / 3]);
00230 //                   if (fabs(p[j + 2 * m / 3]) > maxcoeff3)
00231 //                     maxcoeff3 = fabs(p[j + 2 * m / 3]);
00232 //                 }
00233 //               for (int j = 0; j < m / 3; j++)
00234 //                 {
00235 //                   p[j] /= maxcoeff1;
00236 //                   p[j + m / 3] /= maxcoeff2;
00237 //                   p[j + 2 * m / 3] /= maxcoeff3;
00238 //                 }
00239              }
00240 
00241           LERROR("### Total decrease for iteration : %f", fstart - fp);
00242 
00243           // Save current coefficients:
00244           const std::string filename =
00245             sformat("%s/results/step%03i.txt", PATH, iter);
00246           std::ofstream resultfile (filename.c_str());
00247           if (resultfile.is_open())
00248             {
00249               for (int j = 0; j < m; j++)
00250                 resultfile << p[j] << "\n";
00251               resultfile.close();
00252             }
00253 
00254           // Exit the loop if ITMAX is reach or if the decrease is too low:
00255           if (2.0 * fabs(fstart - fp) <= FTOL * (fabs(fstart) + fabs(fp)))
00256             {
00257               stop = 1;
00258               for (int i = 1; i < NB_PROCS; i++)
00259                 MPI_Send(&stop, 1, MPI_INT, i,
00260                          MSG_STOP, MPI_COMM_WORLD);
00261               LERROR("### Low decrease ! ###");
00262               sleep(10);
00263               // MPI_Finalize();
00264               return 0;
00265             }
00266           if (iter == ITMAX)
00267             {
00268               stop = 1;
00269               for (int i = 1; i < NB_PROCS; i++)
00270                 MPI_Send(&stop, 1, MPI_INT, i,
00271                          MSG_STOP, MPI_COMM_WORLD);
00272               LERROR("### Maximum iterations exceeded ! ###");
00273               sleep(10);
00274               // MPI_Finalize();
00275               return 1;
00276             }
00277 
00278           // Update data:
00279           float xit[m];
00280           float norm = 0;
00281           for (int j = 0; j < m; j++)
00282             {
00283               xit[j] = p[j] - start[j];
00284               norm += xit[j] * xit[j];
00285             }
00286           norm = sqrt(norm);
00287           for (int j = 0; j < m; j++)
00288             {
00289               xi[j][ibig] = xit[j] / norm;
00290             }
00291         }
00292     }
00293   else
00294     {
00295       // ################ SLAVES PROCESS ################
00296 
00297       // **************** INITIALIZATION ****************
00298 
00299       MPI_Status status;
00300 
00301       // Generate the haar transform matrix:
00302       Image<float> hmat(n, n, ZEROS);
00303       for(int i = 0; i < n; i++)
00304         {
00305           hmat.setVal(i, 0, 1.0f);
00306         }
00307       for(int i = 0; i < n / 2; i++)
00308         {
00309           hmat.setVal(i, 1, 1.0f);
00310           hmat.setVal(i + n / 2, 1, -1.0f);
00311           if (i - 2 < 0)
00312             {
00313               hmat.setVal(i, 2, 1.0f);
00314               hmat.setVal(i + 2, 2, -1.0f);
00315             }
00316           else
00317             {
00318               hmat.setVal(i + 2, 3, 1.0f);
00319               hmat.setVal(i + 4, 3, -1.0f);
00320             }
00321           hmat.setVal(2 * i, i + n / 2, 1.0f);
00322           hmat.setVal(2 * i + 1, i + n / 2, -1.0f);
00323         }
00324 
00325       // Load the input images:
00326       ImageSet< PixRGB<byte> > input(NB_PICS);
00327       Point2D<int> ch[20][NB_PICS];
00328       int nch[NB_PICS];
00329       Point2D<int> kp[20][NB_PICS];
00330       int nkp[NB_PICS];
00331       Point2D<int> pw[20][NB_PICS];
00332       int npw[NB_PICS];
00333       Point2D<int> wkm[20][NB_PICS];
00334       int nwkm[NB_PICS];
00335       for (int i = 0; i < NB_PICS; i++)
00336         {
00337           input[i] =
00338             Raster::ReadRGB(sformat("%s/pictures/sat_%03i.ppm", PATH, i));
00339           int count = 0;
00340           std::ifstream chfile(sformat("%s/pictures/sat_%03i.ch.dat", PATH, i).c_str());
00341           bool eofile = false;
00342           while (!chfile.eof())
00343             {
00344               int px, py;
00345               chfile >> px >> py;
00346               if (!chfile.eof())
00347                 {
00348                   ch[count][i].i = px;
00349                   ch[count][i].j = py;
00350                   count++;
00351                 }
00352               else
00353                 eofile = true;
00354             }
00355           chfile.close();
00356           nch[i] = count;
00357           count = 0;
00358           std::ifstream kpfile(sformat("%s/pictures/sat_%03i.kp.dat", PATH, i).c_str());
00359           eofile = false;
00360           while (!eofile)
00361             {
00362               int px, py;
00363               kpfile >> px >>py;
00364               if (!kpfile.eof())
00365                 {
00366                   kp[count][i].i = px;
00367                   kp[count][i].j = py;
00368                   count++;
00369                 }
00370               else
00371                 eofile = true;
00372             }
00373           kpfile.close();
00374           nkp[i] = count;
00375           count = 0;
00376           std::ifstream pwfile(sformat("%s/pictures/sat_%03i.pw.dat", PATH, i).c_str());
00377           eofile = false;
00378           while (!eofile)
00379             {
00380               int px, py;
00381               pwfile >> px >> py;
00382               if (!pwfile.eof())
00383                 {
00384                   pw[count][i].i = px;
00385                   pw[count][i].j = py;
00386                   count++;
00387                 }
00388               else
00389                 eofile = true;
00390             }
00391           pwfile.close();
00392           npw[i] = count;
00393           count = 0;
00394           std::ifstream wkmfile(sformat("%s/pictures/sat_%03i.wkm.dat", PATH, i).c_str());
00395           eofile = false;
00396           while (!eofile)
00397             {
00398               int px, py;
00399               wkmfile >> px >> py;
00400               if (!wkmfile.eof())
00401                 {
00402                   wkm[count][i].i = px;
00403                   wkm[count][i].j = py;
00404                   count++;
00405                 }
00406               else
00407                 eofile = true;
00408             }
00409           wkmfile.close();
00410           nwkm[i] = count;
00411         }
00412 
00413       // **************** OPTIMIZATION ****************
00414 
00415       for (int iter = 0; ; ++iter)
00416         {
00417           // Receive the stop message:
00418           int stop;
00419           MPI_Recv(&stop, 1, MPI_INT, MASTER,
00420                    MSG_STOP, MPI_COMM_WORLD, &status);
00421           if (stop == 1)
00422             {
00423               // MPI_Finalize();
00424               return 0;
00425             }
00426 
00427           // Receive data from master:
00428           int pic;
00429           MPI_Recv(&pic, 1, MPI_INT, MASTER,
00430                    MSG_PIC, MPI_COMM_WORLD, &status);
00431           float data[m];
00432           MPI_Recv(data, m, MPI_FLOAT, MASTER,
00433                    MSG_DATA, MPI_COMM_WORLD, &status);
00434 
00435           // Convert data into filters:
00436           ImageSet<float> trans(NB_FILTERS);
00437           for (int i = 0; i < NB_FILTERS; i++)
00438             trans[i] = Image<float>(data + (n * n * i), n, n);
00439           ImageSet<float> filter(NB_FILTERS);
00440           for (int i = 0; i < NB_FILTERS; i++)
00441             filter[i] = matrixMult(transpose(hmat),
00442                                    matrixMult(trans[i], hmat));
00443 
00444           // Compute the dot product of the filters:
00445           int nb_prod = (NB_FILTERS * (NB_FILTERS - 1)) / 2;
00446           ImageSet<float> prod(nb_prod);
00447           int k = 0;
00448           for (int j = 0; j < NB_FILTERS - 1; j++)
00449             for (int i = j + 1; i < NB_FILTERS; i++)
00450               {
00451                 prod[k] = filter[j] * filter[i];
00452                 k++;
00453               }
00454           float dotprod = 0.0f;
00455           for (int i = 0; i < nb_prod; i++)
00456             dotprod += sum(prod[i]);
00457 
00458           // Instantiate a ModelManager:
00459           ModelManager manager("Open Attention Model");
00460 
00461           // Instantiate our various ModelComponents:
00462           nub::soft_ref<RawVisualCortex> vcx(new RawVisualCortex(manager));
00463           manager.addSubComponent(vcx);
00464           manager.setOptionValString(&OPT_MaxNormType, "Maxnorm");
00465 
00466           for (int i = 0; i < NB_FILTERS; i++)
00467             {
00468               // Create a channel attached to each filter:
00469               nub::soft_ref<ConvolveChannel>
00470                 channel(new ConvolveChannel(manager));
00471               channel->setDescriptiveName(sformat("Convolve%d", i));
00472               channel->setTagName(sformat("conv%d", i));
00473 
00474               // Assign the filter to the channel:
00475               channel->setFilter(filter[i], CONV_BOUNDARY_ZERO);
00476 
00477               // Attach the channel to our visual cortex:
00478               vcx->addSubChan(channel);
00479             }
00480 
00481           // Let's get all our ModelComponent instances started:
00482           manager.exportOptions(MC_RECURSE);
00483           manager.setModelParamString("UsingFPE", "false");
00484           manager.setModelParamString("TestMode", "true");
00485 
00486           // reduce amount of log messages:
00487           MYLOGVERB = loglev;
00488 
00489           manager.start();
00490 
00491           // Process the image through the visual cortex:
00492           vcx->input(InputFrame::fromRgb(&input[pic]));
00493 
00494           // Get the resulting saliency map:
00495           Image<float> sm = vcx->getOutput();
00496 
00497           // Normalize the saliency map:
00498           inplaceNormalize(sm, 0.0f, 255.0f);
00499 
00500           // Get the average saliency:
00501           double avgsal = mean(sm);
00502 
00503           // Get the map level to scale things down:
00504           const LevelSpec lspec =
00505             vcx->getModelParamVal<LevelSpec>("LevelSpec");
00506           int sml = lspec.mapLevel();
00507 
00508           // Scale the radius
00509           int radius = RAD >> sml;
00510 
00511           // Get the saliency on each end of saccade:
00512           float chsal = 0;
00513           int sc = 1 << sml;
00514           for (int i = 0; i < nch[pic]; i++)
00515             chsal += getLocalMax(sm, ch[i][pic] / sc, radius);
00516           chsal /= nch[pic];
00517           float kpsal = 0;
00518           for (int i = 0; i < nkp[pic]; i++)
00519             kpsal += getLocalMax(sm, kp[i][pic] / sc, radius);
00520           kpsal /= nkp[pic];
00521           float pwsal = 0;
00522           for (int i = 0; i < npw[pic]; i++)
00523             pwsal += getLocalMax(sm, pw[i][pic] / sc, radius);
00524           pwsal /= npw[pic];
00525           float wkmsal = 0;
00526           for (int i = 0; i < nwkm[pic]; i++)
00527             wkmsal += getLocalMax(sm, wkm[i][pic] / sc, radius);
00528           wkmsal /= nwkm[pic];
00529 
00530           float goodsal = (chsal + kpsal + pwsal + wkmsal) / 4;
00531 
00532           // Compute the error:
00533           float error = (1 + avgsal) / (1 + goodsal);
00534           float result = error + 0.00001f * fabs(dotprod);
00535 
00536           // Stop all our ModelComponents
00537           manager.stop();
00538 
00539           // Send the result to the master :
00540           MPI_Send(&result, 1, MPI_FLOAT, MASTER,
00541                    MSG_RESULT, MPI_COMM_WORLD);
00542         }
00543     }
00544 
00545 #endif // HAVE_MPI_H
00546 
00547 }
00548 
00549 // ######################################################################
00550 /* So things look consistent in everyone's emacs... */
00551 /* Local Variables: */
00552 /* indent-tabs-mode: nil */
00553 /* End: */
00554 
Generated on Sun May 8 08:40:57 2011 for iLab Neuromorphic Vision Toolkit by  doxygen 1.6.3