GAMPIvision.C

Go to the documentation of this file.
00001 /*!@file GA/GAMPIvision.C Learning filters with genetic algorithm */
00002 
00003 // //////////////////////////////////////////////////////////////////// //
00004 // The iLab Neuromorphic Vision C++ Toolkit - Copyright (C) 2000-2002   //
00005 // by the 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/GA/GAMPIvision.C $
00035 // $Id: GAMPIvision.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 "GA/GAPopulation.H"
00042 #include "Image/MathOps.H"
00043 #include "Image/MatrixOps.H"
00044 #include "Image/Pixels.H"
00045 #include "Channels/RawVisualCortex.H"
00046 #include "Raster/Raster.H"
00047 
00048 #include <fstream>
00049 #ifdef HAVE_MPI_H
00050 #include <mpi.h>
00051 #endif
00052 
00053 #define ITMAX 400
00054 #define MASTER 0
00055 #define MSG_DATA 0
00056 #define MSG_RESULT 1
00057 #define MSG_STOP 2
00058 #define PATH "/home/hpc-26/itti/openvision"
00059 #define NB_PICS 10
00060 #define NB_FILTERS 3
00061 #define NCOEFFS 3
00062 #define POP_SIZE 196
00063 #define RAD 64
00064 
00065 int main(const int argc, const char **argv)
00066 {
00067 #ifndef HAVE_MPI_H
00068 
00069   LFATAL("<mpi.h> must be installed to use this program");
00070 
00071 #else
00072 
00073   // Set verbosity :
00074   int loglev = LOG_ERR;
00075   MYLOGVERB = loglev;
00076 
00077   // Number of coefficients :
00078   int n = 1 << NCOEFFS;
00079   int csize = n * n * NB_FILTERS;
00080   int psize = POP_SIZE;
00081 
00082   // Initialize the MPI system and get the number of processes :
00083   int myrank, nodes;
00084   MPI_Init(const_cast<int *>(&argc), const_cast<char ***>(&argv));
00085   MPI_Comm_size(MPI_COMM_WORLD, &nodes);
00086 
00087   // Check for proper number of processes :
00088   if(nodes < 2)
00089     {
00090       LERROR("*** At least two nodes needed ***");
00091       MPI_Finalize();
00092       return 1;
00093     }
00094 
00095   // Get the rank of the process :
00096   MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
00097 
00098   if (myrank == MASTER)
00099     {
00100       // ################ MASTER PROCESS ################
00101 
00102       MPI_Status status;
00103       int stop = 0;
00104 
00105       // Create a new population :
00106       LERROR("### START ! ###");
00107       GAPopulation pop(psize, csize);
00108 
00109       for (int iter = 1; ; iter++)
00110         {
00111           // Compute the fitness of every chromosome :
00112           int current = 0;
00113           int done = 0;
00114           float low_fit = 1000.0f;
00115           float high_fit = 0.0f;
00116           while (psize - done > 0)
00117             {
00118               int max;
00119               if (psize - done > nodes - 1)
00120                 max = nodes;
00121               else
00122                 max = psize - done + 1;
00123               for (int i = MASTER + 1; i < max; i++)
00124                 {
00125                   GAChromosome c = pop.get_chromosome(current);
00126                   int genes[csize];
00127                   for (int j = 0; j < csize; j++)
00128                     genes[j] = c.get_gene(j);
00129                   MPI_Send(&stop, 1, MPI_INT, i, MSG_STOP,
00130                            MPI_COMM_WORLD);
00131                   MPI_Send(genes, csize, MPI_INT, i, MSG_DATA,
00132                            MPI_COMM_WORLD);
00133                   current++;
00134                 }
00135               for (int i = MASTER + 1; i < max; i++)
00136                 {
00137                   float fit;
00138                   GAChromosome c = pop.get_chromosome(done);
00139                   MPI_Recv(&fit, 1, MPI_FLOAT, i, MSG_RESULT,
00140                            MPI_COMM_WORLD, &status);
00141                   c.set_fitness(fit);
00142                   pop.set_chromosome(done, c);
00143                   done++;
00144                   if (fit > high_fit)
00145                     high_fit = fit;
00146                   if (fit < low_fit)
00147                     low_fit = fit;
00148                 }
00149             }
00150 
00151           // Compute mean fitness and standard deviation :
00152           pop.compute_pop_fitness();
00153           pop.compute_sigma();
00154 
00155           // Display various stuff :
00156           LERROR("*** Highest fitness : %f", high_fit);
00157           LERROR("*** Lowest fitness : %f", low_fit);
00158           LERROR("*** Mean fitness : %f", pop.get_mean_fitness());
00159           LERROR("*** Standard deviation : %f", pop.get_sigma());
00160 
00161           // Select population :
00162           pop.linear_scaling();
00163           pop.selection();
00164 
00165           // Save current population :
00166           char filename[1024];
00167           sprintf(filename, "%s/results/pop%03i.txt", PATH, iter);
00168           std::ofstream resultfile (filename);
00169           if (resultfile.is_open())
00170             {
00171               resultfile << pop;
00172               resultfile.close();
00173             }
00174 
00175           // Population evolution :
00176           LERROR("### EVOLUTION %i ###", iter);
00177           pop.crossover();
00178           pop.mutate();
00179           pop.update();
00180 
00181           // Exit the loop if ITMAX is reach :
00182           if (iter == ITMAX)
00183             {
00184               stop = 1;
00185               for (int i = 1; i < nodes; i++)
00186                 MPI_Send(&stop, 1, MPI_INT, i,
00187                          MSG_STOP, MPI_COMM_WORLD);
00188               LERROR("### Maximum evolutions exceeded ! ###");
00189               sleep(10);
00190               // MPI_Finalize();
00191               return 1;
00192             }
00193 
00194           // TODO : add a test to exit when the evolution is stalled.
00195 
00196         }
00197     }
00198   else
00199     {
00200       // ################ SLAVES PROCESS ################
00201 
00202       // **************** INITIALIZATION ****************
00203 
00204       MPI_Status status;
00205 
00206       // Generate the haar transform matrix:
00207       Image<float> hmat(n, n, ZEROS);
00208       for(int i = 0; i < n; i++)
00209         {
00210           hmat.setVal(i, 0, 1.0f);
00211         }
00212       for(int i = 0; i < n / 2; i++)
00213         {
00214           hmat.setVal(i, 1, 1.0f);
00215           hmat.setVal(i + n / 2, 1, -1.0f);
00216           if (i - 2 < 0)
00217             {
00218               hmat.setVal(i, 2, 1.0f);
00219               hmat.setVal(i + 2, 2, -1.0f);
00220             }
00221           else
00222             {
00223               hmat.setVal(i + 2, 3, 1.0f);
00224               hmat.setVal(i + 4, 3, -1.0f);
00225             }
00226           hmat.setVal(2 * i, i + n / 2, 1.0f);
00227           hmat.setVal(2 * i + 1, i + n / 2, -1.0f);
00228         }
00229 
00230       // Load the input images and the targets locations :
00231 
00232       // Set of pictures :
00233       ImageSet< PixRGB<byte> > input(NB_PICS);
00234 
00235       // Targets locations for each subjects and pictures :
00236       Point2D<int> ch[20][NB_PICS];
00237       Point2D<int> kp[20][NB_PICS];
00238       Point2D<int> pw[20][NB_PICS];
00239       Point2D<int> wkm[20][NB_PICS];
00240 
00241       // Numbers of targets for each subject and pictures :
00242       int nch[NB_PICS];
00243       int nkp[NB_PICS];
00244       int npw[NB_PICS];
00245       int nwkm[NB_PICS];
00246 
00247       // Data loading :
00248       for (int i = 0; i < NB_PICS; i++)
00249         {
00250           char framename[1024];
00251           char chname[1024];
00252           char kpname[1024];
00253           char pwname[1024];
00254           char wkmname[1024];
00255           sprintf(framename,
00256                   "%s/pictures/sat_%03i.ppm", PATH, i);
00257           sprintf(chname,
00258                   "%s/pictures/sat_%03i.ch.dat", PATH, i);
00259           sprintf(kpname,
00260                   "%s/pictures/sat_%03i.kp.dat", PATH, i);
00261           sprintf(pwname,
00262                   "%s/pictures/sat_%03i.pw.dat", PATH, i);
00263           sprintf(wkmname,
00264                   "%s/pictures/sat_%03i.wkm.dat", PATH, i);
00265           input[i] = Raster::ReadRGB(framename);
00266           int count = 0;
00267           std::ifstream chfile(chname);
00268           bool eofile = false;
00269           while (!chfile.eof())
00270             {
00271               int px, py;
00272               chfile >> px >> py;
00273               if (!chfile.eof())
00274                 {
00275                   ch[count][i].i = px;
00276                   ch[count][i].j = py;
00277                   count++;
00278                 }
00279               else
00280                 eofile = true;
00281             }
00282           chfile.close();
00283           nch[i] = count;
00284           count = 0;
00285           std::ifstream kpfile(kpname);
00286           eofile = false;
00287           while (!eofile)
00288             {
00289               int px, py;
00290               kpfile >> px >> py;
00291               if (!kpfile.eof())
00292                 {
00293                   kp[count][i].i = px;
00294                   kp[count][i].j = py;
00295                   count++;
00296                 }
00297               else
00298                 eofile = true;
00299             }
00300           kpfile.close();
00301           nkp[i] = count;
00302           count = 0;
00303           std::ifstream pwfile(chname);
00304           eofile = false;
00305           while (!eofile)
00306             {
00307               int px, py;
00308               pwfile >> px >> py;
00309               if (!pwfile.eof())
00310                 {
00311                   pw[count][i].i = px;
00312                   pw[count][i].j = py;
00313                   count++;
00314                 }
00315               else
00316                 eofile = true;
00317             }
00318           pwfile.close();
00319           npw[i] = count;
00320           count = 0;
00321           std::ifstream wkmfile(chname);
00322           eofile = false;
00323           while (!eofile)
00324             {
00325               int px, py;
00326               wkmfile >> px >> py;
00327               if (!wkmfile.eof())
00328                 {
00329                   wkm[count][i].i = px;
00330                   wkm[count][i].j = py;
00331                   count++;
00332                 }
00333               else
00334                 eofile = true;
00335             }
00336           wkmfile.close();
00337           nwkm[i] = count;
00338         }
00339 
00340       // **************** OPTIMIZATION ****************
00341 
00342       for (int iter = 0; ; ++iter)
00343         {
00344           // Receive the stop message :
00345           int stop;
00346           MPI_Recv(&stop, 1, MPI_INT, MASTER,
00347                    MSG_STOP, MPI_COMM_WORLD, &status);
00348           if (stop)
00349             {
00350               // MPI_Finalize();
00351               return 0;
00352             }
00353 
00354           // Receive data from master :
00355           int data[csize];
00356           MPI_Recv(data, csize, MPI_INT, MASTER, MSG_DATA,
00357                    MPI_COMM_WORLD, &status);
00358 
00359           // Reorder data in a picture :
00360           ImageSet<float> trans(NB_FILTERS);
00361           for (int i = 0; i < NB_FILTERS; i++)
00362             {
00363               trans[i] = Image<float>(n, n, NO_INIT);
00364               trans[i].setVal(0, data[0]);
00365               int range = 1;
00366               int current = 1;
00367               for (int j = 0; j < NCOEFFS; j++)
00368                 {
00369                   for (int k = 0; k < 3; k++)
00370                     {
00371                       int si = range;
00372                       int sj = range;
00373                       if (k == 0)
00374                         sj = 0;
00375                       if (k == 1)
00376                         si = 0;
00377                       for (int ii = si; ii < range + si; ii++)
00378                         for (int jj = sj; jj < range + sj; jj++)
00379                           {
00380                             trans[i].setVal(ii, jj, data[current]);
00381                             current++;
00382                           }
00383                     }
00384                   range <<= 1;
00385                 }
00386             }
00387 
00388           // Convert data into filters :
00389           ImageSet<float> filter(NB_FILTERS);
00390           for (int i = 0; i < NB_FILTERS; i++)
00391             filter[i] = matrixMult(transpose(hmat),
00392                                    matrixMult(trans[i], hmat));
00393 
00394           // Compute the dot product of the filters :
00395           int nb_prod = (NB_FILTERS * (NB_FILTERS - 1)) / 2;
00396           ImageSet<float> prod(nb_prod);
00397           int k = 0;
00398           for (int j = 0; j < NB_FILTERS - 1; j++)
00399             for (int i = j + 1; i < NB_FILTERS; i++)
00400               {
00401                 prod[k] = filter[j] * filter[i];
00402                 k++;
00403               }
00404           float dotprod = 0.0f;
00405           for (int i = 0; i < nb_prod; i++)
00406             dotprod += sum(prod[i]);
00407 
00408           // Instantiate a ModelManager :
00409           ModelManager manager("Open Attention Model");
00410 
00411           // Instantiate our various ModelComponents :
00412           nub::soft_ref<RawVisualCortex> vcx(new RawVisualCortex(manager));
00413           manager.addSubComponent(vcx);
00414           manager.setOptionValString(&OPT_MaxNormType, "Maxnorm");
00415 
00416           for (int i = 0; i < NB_FILTERS; i++)
00417             {
00418               // Create a channel attached to each filter :
00419               nub::soft_ref<ConvolveChannel>
00420                 channel(new ConvolveChannel(manager));
00421               char txt[100];
00422               sprintf(txt, "Convolve%d", i);
00423               channel->setDescriptiveName(txt);
00424               sprintf(txt, "conv%d", i);
00425               channel->setTagName(txt);
00426 
00427               // Assign the filter to the channel:
00428               channel->setFilter(filter[i], CONV_BOUNDARY_ZERO);
00429 
00430               // Attach the channel to our visual cortex:
00431               vcx->addSubChan(channel);
00432             }
00433 
00434           // Let's get all our ModelComponent instances started :
00435           manager.exportOptions(MC_RECURSE);
00436           manager.setModelParamString("UsingFPE", "false");
00437           manager.setModelParamString("TestMode", "true");
00438 
00439           // Reduce amount of log messages :
00440           MYLOGVERB = loglev;
00441 
00442           // Start the manager :
00443           manager.start();
00444 
00445           float error = 0.0f;
00446 
00447           // Compute error for each picture :
00448           for (int pic = 0; pic < NB_PICS; pic++)
00449             {
00450               // Process the image through the visual cortex :
00451               vcx->input(InputFrame::fromRgb(&input[pic]));
00452 
00453               // Get the resulting saliency map :
00454               Image<float> sm = vcx->getOutput();
00455 
00456               // Reset the visual cortex :
00457               vcx->reset(MC_RECURSE);
00458 
00459               // Normalize the saliency map :
00460               inplaceNormalize(sm, 0.0f, 255.0f);
00461 
00462               // Get the average saliency :
00463               double avgsal = mean(sm);
00464 
00465               // Get the map level to scale things down :
00466               const LevelSpec lspec =
00467                 vcx->getModelParamVal<LevelSpec>("LevelSpec");
00468               int sml = lspec.mapLevel();
00469 
00470               // Scale the radius :
00471               int radius = RAD >> sml;
00472 
00473               // Get the saliency on each end of saccade :
00474               float chsal = 0;
00475               int sc = 1 << sml;
00476               for (int i = 0; i < nch[pic]; i++)
00477                 chsal += getLocalMax(sm, ch[i][pic] / sc, radius);
00478               chsal /= nch[pic];
00479               float kpsal = 0;
00480               for (int i = 0; i < nkp[pic]; i++)
00481                 kpsal += getLocalMax(sm, kp[i][pic] / sc, radius);
00482               kpsal /= nkp[pic];
00483               float pwsal = 0;
00484               for (int i = 0; i < npw[pic]; i++)
00485                 pwsal += getLocalMax(sm, pw[i][pic] / sc, radius);
00486               pwsal /= npw[pic];
00487               float wkmsal = 0;
00488               for (int i = 0; i < nwkm[pic]; i++)
00489                 wkmsal += getLocalMax(sm, wkm[i][pic] / sc, radius);
00490               wkmsal /= nwkm[pic];
00491 
00492               float goodsal = (chsal + kpsal + pwsal + wkmsal) / 4;
00493 
00494               // Compute the error :
00495               error += (1 + avgsal) / (1 + goodsal);
00496             }
00497 
00498           // Stop all our ModelComponents :
00499           manager.stop();
00500 
00501           // Compute total result :
00502           float result = (error / NB_PICS) + 0.00001f * fabs(dotprod);
00503           result = 1.0f / result;
00504 
00505           // Send the result to the master :
00506           MPI_Send(&result, 1, MPI_FLOAT, MASTER, MSG_RESULT,
00507                    MPI_COMM_WORLD);
00508         }
00509     }
00510 
00511 #endif // HAVE_MPI_H
00512 
00513 }
00514 
00515 // ######################################################################
00516 /* So things look consistent in everyone's emacs... */
00517 /* Local Variables: */
00518 /* indent-tabs-mode: nil */
00519 /* End: */
Generated on Sun May 8 08:40:39 2011 for iLab Neuromorphic Vision Toolkit by  doxygen 1.6.3