PLearn 0.1
ChemicalICP.cc
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // ChemicalICP.cc
00004 //
00005 // Copyright (C) 2006 Pascal Lamblin
00006 //
00007 // Redistribution and use in source and binary forms, with or without
00008 // modification, are permitted provided that the following conditions are met:
00009 //
00010 //  1. Redistributions of source code must retain the above copyright
00011 //     notice, this list of conditions and the following disclaimer.
00012 //
00013 //  2. Redistributions in binary form must reproduce the above copyright
00014 //     notice, this list of conditions and the following disclaimer in the
00015 //     documentation and/or other materials provided with the distribution.
00016 //
00017 //  3. The name of the authors may not be used to endorse or promote
00018 //     products derived from this software without specific prior written
00019 //     permission.
00020 //
00021 // THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR
00022 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
00023 // OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN
00024 // NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
00025 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
00026 // TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00027 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00028 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00029 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00030 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00031 //
00032 // This file is part of the PLearn library. For more information on the PLearn
00033 // library, go to the PLearn Web site at www.plearn.org
00034 
00035 // Authors: Pascal Lamblin
00036 
00039 #include "ChemicalICP.h"
00040 #include <plearn/base/stringutils.h>
00041 #include <plearn/io/openFile.h>
00042 #include <plearn/var/UnaryVariable.h>
00043 #include <plearn/var/VarColumnsVariable.h>
00044 #include <plearn/io/pl_log.h>
00045 #include "geometry.h"
00046 
00047 namespace PLearn {
00048 using namespace std;
00049 
00050 PLEARN_IMPLEMENT_OBJECT(
00051     ChemicalICP,
00052     "Iterative Closest Points (ICP) aligner.",
00053     ""
00054 );
00055 
00056 ChemicalICP::ChemicalICP():
00057       mol_feat_indices( new UnaryVariable() ),
00058       template_feat_indices( new UnaryVariable() ),
00059       matching_neighbors( new UnaryVariable() ),
00060       weighting_method( "features_sigmoid" ),
00061       weighting_params( Vec(2,1) ),
00062       matching_method( "exhaustive" ),
00063       memory( "none" ),
00064       initial_angles_step( 0 ),
00065       max_iter( 50 ),
00066       error_t( 0 ),
00067       angle_t( 0.5 ),
00068       trans_t( 0 ),
00069       rotation( 3, 3 ),
00070       translation( 3 ),
00071       all_mol_features( 0 ),
00072       all_template_features( 0 ),
00073       all_template_feat_dev( 0 ),
00074       mol_coordinates( 0 ),
00075       used_mol_features( new VarColumnsVariable( all_mol_features,
00076                                                  mol_feat_indices ) ),
00077       template_coordinates( 0 ),
00078       template_geom_dev( 0 ),
00079       used_template_features( new VarColumnsVariable(all_template_features,
00080                                                      template_feat_indices) ),
00081       used_template_feat_dev( new VarColumnsVariable(all_template_feat_dev,
00082                                                      template_feat_indices) )
00083 {}
00084 
00085 ChemicalICP::ChemicalICP( const MolTemplate& the_template,
00086                           const Mol& the_molecule,
00087                           const TVec<string>& the_feature_names,
00088                           string the_weighting_method,
00089                           const Var& the_weighting_params,
00090                           string the_matching_method ):
00091       mol_feat_indices( new UnaryVariable() ),
00092       template_feat_indices( new UnaryVariable() ),
00093       matching_neighbors( new UnaryVariable() ),
00094       mol_template(the_template),
00095       molecule(the_molecule),
00096       feature_names(the_feature_names),
00097       weighting_method(the_weighting_method),
00098       weighting_params(the_weighting_params),
00099       matching_method( the_matching_method ),
00100       memory( "none" ),
00101       initial_angles_step( 0 ),
00102       max_iter( 50 ),
00103       error_t( 0 ),
00104       angle_t( 0.5 ),
00105       trans_t( 0 ),
00106       rotation( 3, 3 ),
00107       translation( 3 ),
00108       all_mol_features( 0 ), // calls SourceVariable()
00109       all_template_features( 0 ),
00110       all_template_feat_dev( 0 ),
00111       mol_coordinates( 0 ),
00112       used_mol_features( new VarColumnsVariable( all_mol_features,
00113                                                  mol_feat_indices ) ),
00114       template_coordinates( 0 ),
00115       template_geom_dev( 0 ),
00116       used_template_features( new VarColumnsVariable(all_template_features,
00117                                                      template_feat_indices) ),
00118       used_template_feat_dev( new VarColumnsVariable(all_template_feat_dev,
00119                                                      template_feat_indices) )
00120 {
00121 //pout << "Big constructor called" << endl;
00122     build();
00123 }
00124 
00125 void ChemicalICP::setMolecule( const Mol& the_molecule )
00126 {
00127     molecule = the_molecule;
00128     computeUsedFeatures();
00129     computeVariables();
00130 
00131     if( matching_method == "exhaustive" )
00132         cacheFeatureDistances();
00133 }
00134 
00135 void ChemicalICP::computeUsedFeatures()
00136 {
00137     // compute intersection between mol_template->feature_names and
00138     // feature_names
00139     TVec<string> template_names = mol_template->feature_names;
00140     TVec<string> common_names;
00141     TVec<int> common_indices;
00142 
00143     if( feature_names.length() == 0 )
00144     {
00145         common_names = template_names.copy();
00146         common_indices = TVec<int>( 0, template_names.length() - 1, 1 );
00147     }
00148     else if( feature_names[0] == "none" )
00149     {
00150         used_feat_names.resize( 0 );
00151         mol_feat_indices->resize( 0, 0 );
00152         template_feat_indices->resize( 0, 0 );
00153         return;
00154     }
00155     else
00156     {
00157         common_names.resize(0, feature_names.length()); // to have some space
00158         common_indices.resize(0, feature_names.length());
00159         TVec<int> indices = feature_names.find( template_names );
00160 
00161         for( int i=0 ; i<indices.length() ; i++ )
00162         {
00163             if( indices[i] >= 0 )
00164             {
00165                 // template_names[i] is present in feature_names
00166                 common_names.append( template_names[i] );
00167                 common_indices.append( i );
00168             }
00169         }
00170     }
00171 
00172     // then compute intersection with molecule->feature_names
00173     TVec<string> mol_names = molecule->feature_names;
00174     used_feat_names.resize(0, common_names.length());
00175 
00176     // storage for indices
00177     Vec template_feat_indices_;
00178     Vec mol_feat_indices_;
00179 
00180     TVec<int> indices = mol_names.find( common_names );
00181     for( int i=0 ; i<indices.length() ; i++ )
00182     {
00183         if( indices[i] >= 0 )
00184         {
00185             // common_names[i] is present in mol_names at position indices[i]
00186             used_feat_names.append( common_names[i] );
00187             template_feat_indices_.append( common_indices[i] );
00188             mol_feat_indices_.append( indices[i] );
00189         }
00190     }
00191 
00192     // update Var's value
00193     template_feat_indices->defineValueLocation(
00194         template_feat_indices_.toMat( template_feat_indices_.size(), 1 ) );
00195 
00196     mol_feat_indices->defineValueLocation(
00197         mol_feat_indices_.toMat( mol_feat_indices_.size(), 1 ) );
00198 }
00199 
00200 void ChemicalICP::computeVariables()
00201 {
00202     // make mol_coordinates point to the new molecule's coordinates
00203     mol_coordinates->defineValueLocation( molecule->coordinates );
00204 
00205     // make mol_features point to the new molecule's features
00206     all_mol_features->defineValueLocation( molecule->features );
00207 
00208     // update used_... Var's so they have the right size
00209     used_mol_features->sizefprop();
00210     used_template_features->sizefprop();
00211     used_template_feat_dev->sizefprop();
00212 }
00213 
00214 void ChemicalICP::cacheFeatureDistances()
00215 {
00216 //pout << "cacheFeatureDistances()" << endl;
00217     Mat t_features = used_template_features->matValue;
00218     Mat t_feat_dev = used_template_feat_dev->matValue;
00219     Mat m_features = used_mol_features->matValue;
00220 
00221     int n_template_points = t_features.length();
00222     int n_mol_points = m_features.length();
00223     int n_features = m_features.width();
00224 
00225     feat_distances2.resize( n_template_points, n_mol_points );
00226     feat_distances2.fill(0);
00227 
00228     for( int i=0 ; i<n_template_points ; i++ )
00229     {
00230         for( int j=0 ; j<n_mol_points ; j++ )
00231         {
00232             for( int k=0 ; k<n_features ; k++ )
00233             {
00234                 real diff = (t_features(i, k) - m_features(j, k))
00235                                 / t_feat_dev(i, k);
00236                 feat_distances2(i, j) += diff * diff;
00237             }
00238         }
00239     }
00240 /*pout << "t_features(0) = " << t_features(0) << endl;
00241 pout << "t_feat_dev(0) = " << t_feat_dev(0) << endl;
00242 pout << "feat_distances2(0) = " << feat_distances2(0) << endl;*/
00243 }
00244 
00246 // forgetMemorizedAlignments //
00248 void ChemicalICP::forgetMemorizedAlignments()
00249 {
00250     mem_aligned_molecules.resize(0);
00251     mem_error.resize(0);
00252     mem_rotation.resize(0);
00253     mem_translation.resize(0);
00254     mem_matching.resize(0);
00255     mem_count.resize(0);
00256 }
00257 
00259 // run //
00261 void ChemicalICP::run()
00262 {
00263     if( initial_angles_list.length() > 0 )
00264     {
00265         if( initial_angles_step > 0 )
00266         {
00267             PLWARNING( "ChemicalICP::run - both 'initial_angles_step' and"
00268                        " 'initial_angles_list'\n"
00269                        "are provided. Setting 'initial_angles_step' to 0.\n" );
00270             initial_angles_step = 0;
00271         }
00272     }
00273     else
00274     {
00275         if( fast_is_equal( initial_angles_step, 0. ) )
00276             initial_angles_step = 360;
00277 
00278         initial_angles_list.resize( 0, 3 );
00279         for( real rx=0. ; rx<360. ; rx += initial_angles_step )
00280             for( real ry=0. ; ry<360. ; ry += initial_angles_step )
00281                 for( real rz=0. ; rz<180. ; rz += initial_angles_step )
00282                 {
00283                     Vec angles( 3 );
00284                     angles[0] = rx;
00285                     angles[1] = ry;
00286                     angles[2] = rz;
00287                     initial_angles_list.appendRow( angles );
00288                 }
00289     }
00290     EXTREME_NAMED_LOG("ICP+NNet") << "begin run()" << endl;
00291 
00292     // Only perform the alignment if we do not just obtain the result from a
00293     // previously memorized alignment.
00294     bool realign = true;
00295     int find_mol = -1;
00296     if (memory != "none") {
00297         find_mol = mem_aligned_molecules.find(molecule);
00298         if (find_mol == -1) {
00299             // This molecule had not been previously aligned.
00300             find_mol = mem_aligned_molecules.length();
00301             mem_aligned_molecules.append(molecule);
00302             // Dummy additions just to resize the vectors correctly.
00303             mem_error.append(0);
00304             mem_rotation.append(Mat());
00305             mem_translation.append(Vec());
00306             mem_matching.append(TVec<int>());
00307             mem_count.append(0);
00308         }
00309     }
00310 
00311     if (every_x > 0) {
00312         PLASSERT( find_mol >= 0 );
00313         if (mem_count[find_mol] % every_x != 0)
00314             realign = false;
00315     }
00316 
00317     if (realign) {
00318 
00319     int n_points = mol_template->n_points();
00320     if( n_points < 3 )
00321         PLERROR( "ChemicalICP::run() - not enough points in template (%d).\n",
00322                  n_points );
00323 
00324     Mat best_rotation(3, 3);
00325     Vec best_translation(3);
00326     TVec<int> best_matching( n_points );
00327     Vec best_weights( n_points );
00328     real best_error = REAL_MAX;
00329 
00330     // transformed template coordinates
00331     Mat tr_template_coords( n_points, 3 );
00332 
00333     // coordinates of molecule points matched to template ones
00334     Mat matched_mol_coords( n_points, 3 );
00335 
00336     int n_initial_angles = initial_angles_list.length();
00337     EXTREME_NAMED_LOG("ICP+NNet") << n_initial_angles << " initial angles to try" << endl;
00338     for( int i=0 ; i<n_initial_angles ; i++ )
00339     {
00340                 EXTREME_NAMED_LOG("ICP+NNet") << "global iteration number " << i << endl;
00341         // initialization
00342         rotation = rotationMatrixFromAngles( initial_angles_list(i) );
00343         translation.fill(0);
00344         matching = TVec<int>( n_points, -1 );
00345         weights = Vec( n_points );
00346 
00347         int n_iter = 0;
00348         real delta_rot_length = REAL_MAX;
00349         real delta_trans_length = REAL_MAX;
00350 
00351         applyGeomTransformation( rotation, translation,
00352                                  template_coordinates->matValue,
00353                                  tr_template_coords );
00354                 EXTREME_NAMED_LOG("ICPDEBUG") << "beginning ICP:" << endl
00355                         << "rotation = " << endl << rotation << endl
00356                 << "translation = " << translation << endl
00357                 << "template_coordinates->matValue = " << endl
00358                 << template_coordinates->matValue << endl
00359                 << "tr_template_coords = " << endl << tr_template_coords << endl
00360                 << endl;
00361         // main loop
00362         do
00363         {
00364                         EXTREME_NAMED_LOG("ICPDEBUG") << "    beginning of main loop" << endl;
00365                         EXTREME_NAMED_LOG("ICPDEBUG") << "    matchNearestNeighbors()" << endl;
00366             matchNearestNeighbors( tr_template_coords, matched_mol_coords );
00367                         EXTREME_NAMED_LOG("ICPDEBUG") << "    minimizeWeightedDistance()" << endl;
00368             minimizeWeightedDistance( tr_template_coords, matched_mol_coords,
00369                                       delta_rot_length, delta_trans_length );
00370                         EXTREME_NAMED_LOG("ICPDEBUG") << "    applyGeomTransformation()" << endl;
00371             applyGeomTransformation( rotation, translation,
00372                                      template_coordinates->matValue,
00373                                      tr_template_coords );
00374                         EXTREME_NAMED_LOG("ICPDEBUG") << "tr_template_coords = " << endl << tr_template_coords << endl;
00375             error = computeWeightedDistance( tr_template_coords,
00376                                              matched_mol_coords );
00377             n_iter++;
00378                         EXTREME_NAMED_LOG("ICPDEBUG") << "end of main loop" << endl;
00379                         EXTREME_NAMED_LOG("ICPDEBUG") << "    iteration = " << n_iter << " / " << max_iter << endl
00380                                 << "    error = " << error << " / " << error_t << endl
00381                                 << "    delta_rot_length = " << delta_rot_length << " / "
00382                                 << angle_t << endl
00383                                 << "    delta_trans_length = " << delta_trans_length << " / "
00384                                 << trans_t << endl
00385                                 << endl;
00386         }
00387         while( n_iter < max_iter &&
00388                error > error_t &&
00389                delta_rot_length > angle_t &&
00390                delta_trans_length > trans_t  );
00391 
00392         // keep the best one
00393         if( error < best_error )
00394         {
00395             best_error = error;
00396             best_rotation << rotation;
00397             best_translation << translation;
00398             best_matching << matching;
00399         }
00400                 EXTREME_NAMED_LOG("ICP+NNet") << "end global iteration number " << i << endl;
00401     }
00402 
00403     // get best parameters
00404     error = best_error;
00405     rotation = best_rotation;
00406     translation = best_translation;
00407     matching = best_matching;
00408 
00409     if (memory != "none") {
00410         // Remember this alignment.
00411         PLASSERT( find_mol >= 0 );
00412         mem_error[find_mol] = error;
00413         mem_rotation[find_mol].resize(rotation.length(), rotation.width());
00414         mem_rotation[find_mol] << rotation;
00415         mem_translation[find_mol].resize(translation.length());
00416         mem_translation[find_mol] << translation;
00417         mem_matching[find_mol].resize(matching.length());
00418         mem_matching[find_mol] << matching;
00419     }
00420 
00421     } else {
00422         PLASSERT( find_mol >= 0 );
00423         error = mem_error[find_mol];
00424         // Copies make the code simpler, though it would be slightly more
00425         // efficient to use resize and the << operator.
00426         rotation = mem_rotation[find_mol].copy();
00427         translation = mem_translation[find_mol].copy();
00428         matching = mem_matching[find_mol].copy();
00429     }
00430 
00431     if (find_mol >= 0)
00432         // Increment counter of alignments for this molecule.
00433         mem_count[find_mol]++;
00434 
00435     // Update the 'matching_neighbors' variable.
00436     matching_neighbors->resize(matching.length(), 1);
00437     for (int i = 0; i < matching.length(); i++)
00438         matching_neighbors->value[i] = matching[i];
00439 
00440     if( !fast_is_equal( initial_angles_step, 0. ) )
00441         initial_angles_list.resize( 0, 3 );
00442 }
00443 
00445 // matchNearestNeighbors //
00447 void ChemicalICP::matchNearestNeighbors( const Mat& tr_template_coords,
00448                                          const Mat& matched_mol_coords )
00449 {
00450 //pout << "matchNearestNeighbors()" << endl;
00451     Mat mol_coords = mol_coordinates->matValue;
00452     int n_template_points = tr_template_coords.length();
00453     int n_mol_points = mol_coords.length();
00454 
00455     if( matching_method == "exhaustive" )
00456     {
00457         // bruteforce searche
00458         for( int i=0 ; i<n_template_points ; i++ )
00459         {
00460             Vec t_point = tr_template_coords( i );
00461 //pout << "tr_template_coords(" << i  << ") = " << t_point << endl;
00462             Vec dists( n_template_points );
00463             real closest_dist2 = REAL_MAX;
00464 //pout << "closest_dist2 = " << closest_dist2 << endl;
00465 
00466             for( int j=0 ; j<n_mol_points ; j++ )
00467             {
00468                 // compute distance
00469                 Vec m_point = mol_coords( j );
00470                 real dist2 = powdistance( t_point, m_point, 2 )
00471                                 + feat_distances2(i, j);
00472 /*pout << "m_point = " << m_point << endl;
00473 pout << "powdistance( t_point, m_point, 2 ) = "
00474     << powdistance( t_point, m_point, 2 ) << endl;
00475 pout << "feat_distances2(" << i << "," << j << ") = " <<
00476     feat_distances2(i, j) << endl;
00477 pout << "dist2 = " << dist2 << endl;*/
00478 
00479                 // keep the smallest
00480                 if( dist2 < closest_dist2 )
00481                 {
00482 //pout << "beep " << endl;
00483 //pout << " matching[" << i << "] = " << j << endl;
00484                     closest_dist2 = dist2;
00485                     matching[i] = j;
00486                 }
00487             }
00488             matched_mol_coords( i ) << mol_coords( matching[i] );
00489         }
00490     }
00491 }
00492 
00493 void ChemicalICP::minimizeWeightedDistance( const Mat& tr_template_coords,
00494                                             const Mat& matched_mol_coords,
00495                                             real& delta_rot_length,
00496                                             real& delta_trans_length )
00497 {
00498     Mat delta_rot( 3, 3 );
00499     Vec delta_trans( 3 );
00500     real err = REAL_MAX;
00501 
00502 //pout << "before computeWeights()" << endl
00503 //    << tr_template_coords(0) << endl;
00504     computeWeights( tr_template_coords, matched_mol_coords );
00505 //pout << "after computeWeights() / before transformationFromWeightedMatchedPoints" << endl
00506 //    << tr_template_coords(0) << endl;
00507     transformationFromWeightedMatchedPoints( tr_template_coords,
00508                                              matched_mol_coords,
00509                                              weights,
00510                                              delta_rot, delta_trans,
00511                                              err );
00512 /*pout << "after transformationFromWeightedMatchedPoints" << endl
00513     << tr_template_coords(0) << endl;
00514 pout << "weights = " << weights << endl
00515      << "delta_rot = " << endl << delta_rot << endl
00516      << "delta_trans = " << delta_trans << endl;*/
00517 
00518     delta_rot_length = norm( anglesFromRotationMatrix( delta_rot ), 2 );
00519     delta_trans_length = norm( delta_trans );
00520 
00521     // accumulate transformation ensuring normalization
00522     Vec angles = anglesFromRotationMatrix( product( rotation, delta_rot ) );
00523     rotation << rotationMatrixFromAngles( angles );
00524 
00525     // translation = delta_trans + delta_rot * translation
00526     productAcc( delta_trans, delta_rot, translation );
00527     translation << delta_trans;
00528 //pout << "rotation = " << endl << rotation << endl
00529 //     << "translation = " << translation << endl;
00530 }
00531 
00532 real ChemicalICP::computeWeightedDistance( const Mat& tr_template_coords,
00533                                            const Mat& matched_mol_coords )
00534 {
00535 //pout << "computeWeightedDistance()" << endl;
00536     real err = 0;
00537     int n_points = tr_template_coords.length();
00538     if( matching_method == "exhaustive" )
00539     {
00540         for( int i=0 ; i<n_points ; i++ )
00541         {
00542             real diff2 = feat_distances2( i, matching[i] ) +
00543                 powdistance( tr_template_coords(i), matched_mol_coords(i), 2 );
00544             err += weights[i] * sqrt( diff2 );
00545 /*pout << "i = " << i << endl
00546     << "tr_template_coords = " << tr_template_coords(i) << endl
00547     << "matched_mol_coords = " << matched_mol_coords(i) << endl
00548     << "weights = " <<  weights[i] << endl
00549     << "delta_err = " << weights[i] * sqrt( diff2 ) << endl
00550     << endl;*/
00551         }
00552     }
00553     return err;
00554 }
00555 
00556 void ChemicalICP::computeWeights( const Mat& tr_template_coords,
00557                                   const Mat& matched_mol_coords )
00558 {
00559 //pout << "computeWeights()" << endl;
00560     int n_points = tr_template_coords.length();
00561 //pout << "n_points = " << n_points << endl;
00562 
00563     if( weighting_method == "none" )
00564     {
00565 //pout << "pouet!" << endl;
00566         weights.fill( 1. / real( n_points ) );
00567 //pout << "1. / real( n_points ) = " << 1. / real( n_points ) << endl;
00568 //pout << "weights[0] = " << weights[0] << endl;
00569     }
00570     else if( weighting_method == "features_sigmoid" )
00571     {
00572 //pout << "plouf!" << endl;
00573         real total_weight = 0;
00574         real mid = weighting_params->value[0];
00575         real slope = weighting_params->value[1];
00576 
00577         if( matching_method == "exhaustive" )
00578         {
00579             for( int i=0 ; i<n_points ; i++ )
00580             {
00581                 real diff2 = feat_distances2( i, matching[i] );
00582                 weights[i] = sigmoid( slope * (mid - sqrt(diff2)) );
00583                 total_weight += weights[i];
00584 /*pout << "i = " << i << endl
00585     << "diff2 = " << diff2 << endl
00586     << "slope = " << slope << endl
00587     << "mid = " << mid << endl
00588     << "weights = " <<  weights[i] << endl
00589     << endl;*/
00590             }
00591             weights /= total_weight;
00592 //pout << "total_weight = " << total_weight << endl;
00593         }
00594     }
00595 }
00596 
00597 void ChemicalICP::saveMatch( const PPath& filename )
00598 {
00599     PStream file = openFile( filename, PStream::plearn_ascii, "w" );
00600     file<< "rotation =" << endl
00601         << rotation << endl
00602         << "translation =" << endl
00603         << translation << endl
00604         << "matching =" << endl
00605         << matching << endl
00606         << "error = " << error << endl;
00607 }
00608 
00609 // ### Nothing to add here, simply calls build_
00610 void ChemicalICP::build()
00611 {
00612     inherited::build();
00613     build_();
00614 }
00615 
00617 // makeDeepCopyFromShallowCopy //
00619 void ChemicalICP::makeDeepCopyFromShallowCopy(CopiesMap& copies)
00620 {
00621     inherited::makeDeepCopyFromShallowCopy(copies);
00622 
00623     varDeepCopyField(mol_feat_indices, copies);
00624     varDeepCopyField(template_feat_indices, copies);
00625     varDeepCopyField(matching_neighbors, copies);
00626     deepCopyField(mol_template, copies);
00627     deepCopyField(molecule, copies);
00628     deepCopyField(feature_names, copies);
00629     varDeepCopyField(weighting_params, copies);
00630     deepCopyField(initial_angles_list, copies);
00631     deepCopyField(rotation, copies);
00632     deepCopyField(translation, copies);
00633     deepCopyField(matching, copies);
00634     deepCopyField(weights, copies);
00635     deepCopyField(used_properties, copies);
00636     deepCopyField(other_base_properties, copies);
00637 
00638     deepCopyField(mem_aligned_molecules,    copies);
00639     deepCopyField(mem_error,                copies);
00640     deepCopyField(mem_rotation,             copies);
00641     deepCopyField(mem_translation,          copies);
00642     deepCopyField(mem_matching,             copies);
00643     deepCopyField(mem_count,                copies);
00644     
00645     deepCopyField(used_feat_names, copies);
00646     deepCopyField(feat_distances2, copies);
00647 
00648     varDeepCopyField(all_mol_features, copies);
00649     varDeepCopyField(all_template_features, copies);
00650     varDeepCopyField(all_template_feat_dev, copies);
00651     varDeepCopyField(mol_coordinates, copies);
00652     varDeepCopyField(used_mol_features, copies);
00653     varDeepCopyField(template_coordinates, copies);
00654     varDeepCopyField(template_geom_dev, copies);
00655     varDeepCopyField(used_template_features, copies);
00656     varDeepCopyField(used_template_feat_dev, copies);
00657 }
00658 
00660 // declareOptions //
00662 void ChemicalICP::declareOptions(OptionList& ol)
00663 {
00664     // declareOption(ol, "myoption", &ChemicalICP::myoption, OptionBase::buildoption,
00665     //               "Help text describing this option");
00666     // ...
00667 
00668     declareOption(ol, "mol_template", &ChemicalICP::mol_template,
00669                   OptionBase::buildoption,
00670                   "The template we try to align on the molecule");
00671 
00672     declareOption(ol, "molecule", &ChemicalICP::molecule,
00673                   OptionBase::buildoption,
00674                   "The molecule");
00675 
00676     declareOption(ol, "feature_names", &ChemicalICP::feature_names,
00677                   OptionBase::buildoption,
00678                   "Names of features to use during alignment.\n"
00679                   "Empty TVec means 'use all available features.\n"
00680                   "Use '[ \"none\" ]' if you don't want to use any feature.\n"
00681                   );
00682 
00683     declareOption(ol, "weighting_method", &ChemicalICP::weighting_method,
00684                   OptionBase::buildoption,
00685                   "Method used to compute the weight of a pair of point."
00686                   " One of:\n"
00687                   "    - \"features_sigmoid\": sigmoid of feature distance,\n"
00688                   "    - \"none\": same weight for each pair.\n"
00689                   );
00690 
00691     declareOption(ol, "weighting_params", &ChemicalICP::weighting_params,
00692                   OptionBase::buildoption,
00693                   "Var containing parameters used during weighting.\n"
00694                   "Size and meaning depends on value of 'weighting_method'.\n"
00695                   );
00696 
00697     declareOption(ol, "matching_method", &ChemicalICP::matching_method,
00698                   OptionBase::buildoption,
00699                   "Method used to find the nearest neighbors. For the moment,"
00700                   " only one:\n"
00701                   "    - \"exhaustive\": exhaustive search (caching feature"
00702                   " distances).\n"
00703                   );
00704 
00705     declareOption(ol, "memory", &ChemicalICP::memory, OptionBase::buildoption,
00706         "The kind of memory allowing ICP to re-use previous alignments:\n"
00707         "- 'none'   : no memory\n"
00708         "- 'every-X': the alignment for a given molecule is re-computed only\n"
00709         "             once every X times it is asked\n"
00710     );
00711 
00712     declareOption(ol, "initial_angles_step", &ChemicalICP::initial_angles_step,
00713                   OptionBase::buildoption,
00714                   "Tries initial rotations every \"initial_angles_step\""
00715                   " degrees");
00716 
00717     declareOption(ol, "initial_angles_list", &ChemicalICP::initial_angles_list,
00718                   OptionBase::buildoption,
00719                   "Explicit list of initial rotations angles");
00720 
00721     declareOption(ol, "max_iter", &ChemicalICP::max_iter,
00722                   OptionBase::buildoption,
00723                   "Maximum number of iterations to perform during alignment");
00724 
00725     declareOption(ol, "error_t", &ChemicalICP::error_t,
00726                   OptionBase::buildoption,
00727                   "Stop alignment if error falls below this threshold");
00728 
00729     declareOption(ol, "angle_t", &ChemicalICP::angle_t,
00730                   OptionBase::buildoption,
00731                   "Stop alignment if angles falls below this threshold");
00732 
00733     declareOption(ol, "trans_t", &ChemicalICP::trans_t,
00734                   OptionBase::buildoption,
00735                   "Stop alignment if translation falls below this threshold");
00736 
00737     declareOption(ol, "rotation", &ChemicalICP::rotation,
00738                   OptionBase::learntoption,
00739                   "Learned rotation matrix");
00740 
00741     declareOption(ol, "translation", &ChemicalICP::translation,
00742                   OptionBase::learntoption,
00743                   "Learned translation vector");
00744 
00745     declareOption(ol, "matching", &ChemicalICP::matching,
00746                   OptionBase::learntoption,
00747                   "matching[i] is the index of the molecule point being\n"
00748                   "the nearest neighbor of template point i.\n");
00749 
00750     declareOption(ol, "weights", &ChemicalICP::weights,
00751                   OptionBase::learntoption,
00752                   "Weight of the pair of points (i, matching[i])");
00753 
00754     declareOption(ol, "error", &ChemicalICP::error,
00755                   OptionBase::learntoption,
00756                   "Weigted error of the alignment");
00757 
00758     // Now call the parent class' declareOptions
00759     inherited::declareOptions(ol);
00760 }
00761 
00763 // build_ //
00765 void ChemicalICP::build_()
00766 {
00767 
00768 #ifdef BOUNDCHECK
00769     // Variable names for debugging.
00770     used_mol_features->setName("used_mol_features");
00771     used_template_features->setName("used_template_features");
00772     used_template_feat_dev->setName("used_template_feat_dev");
00773 #endif
00774 
00775     if (memory == "none") {
00776         // A negative value for 'every_x' means it is not used.
00777         every_x = -1;
00778     } else if (memory.find("every-") == 0) {
00779         string x = memory.substr(6);
00780         every_x = atoi(x.c_str());
00781         if (every_x <= 0)
00782             PLERROR("In ChemicalICP::build_ - The X value in 'every-X' (option"
00783                     " 'memory') must be positive");
00784     } else
00785         PLERROR("In ChemicalICP::build_ - Invalid value for option 'memory'");
00786 
00787     if( feature_names.size() > 0 &&
00788         lowerstring( feature_names[0] ) == "none" )
00789     {
00790         // no feature will be used during the alignment nor score computation
00791         feature_names[0] = "none";
00792         if( feature_names.size() > 1 )
00793         {
00794             PLWARNING("First element of 'feature_names' is 'none', but"
00795                       " other features are present.\n"
00796                       "Resizing 'feature_names' to 1.\n");
00797             feature_names.resize( 1 );
00798         }
00799     }
00800 
00801     if( lowerstring( weighting_method ) == "none" || weighting_method == "" )
00802     {
00803         weighting_method = "none";
00804         weighting_params->resize(0, 0);
00805     }
00806     else if( lowerstring( weighting_method ) == "features_sigmoid" )
00807     {
00808         weighting_method = "features_sigmoid";
00809         weighting_params->resize(2, 1);
00810     }
00811     else
00812         PLERROR( "ChemicalICP::build_ - weighting_method '%s' is unknown.\n",
00813                  weighting_method.c_str() );
00814 
00815 
00816     if( lowerstring( matching_method ) == "exhaustive" )
00817         matching_method = "exhaustive";
00818     else
00819         PLERROR( "ChemicalICP::build_ - matching_method '%s' is unknown.\n",
00820                  matching_method.c_str() );
00821 
00822     if( mol_template )
00823     {
00824         // make the Var's relative to the template have the right storage
00825         template_coordinates->defineValueLocation( mol_template->coordinates );
00826         template_geom_dev->defineValueLocation(
00827             mol_template->geom_dev.toMat( mol_template->n_points(), 1 ) );
00828         all_template_features->defineValueLocation( mol_template->features );
00829         all_template_feat_dev->defineValueLocation( mol_template->feat_dev );
00830 
00831         if( molecule )
00832         {
00833             // make as if 'setMolecule' were called
00834             computeUsedFeatures();
00835             computeVariables();
00836 
00837             if( matching_method == "exhaustive" )
00838                 cacheFeatureDistances();
00839         }
00840 
00841         // Resize 'matching_neighbors' so that it has the correct size (i.e.
00842         // the number of points in the template).
00843         matching_neighbors->resize(mol_template->n_points(), 1);
00844     }
00845 
00846     // build VarArray
00847     used_properties = mol_coordinates
00848         & template_coordinates & template_geom_dev
00849         & used_mol_features & used_template_features & used_template_feat_dev;
00850 
00851     other_base_properties = all_mol_features
00852         & all_template_features & all_template_feat_dev;
00853 //pout << "end build_()" << endl;
00854 
00855 }
00856 
00857 
00858 } // end of namespace PLearn
00859 
00860 
00861 /*
00862   Local Variables:
00863   mode:c++
00864   c-basic-offset:4
00865   c-file-style:"stroustrup"
00866   c-file-offsets:((innamespace . 0)(inline-open . 0))
00867   indent-tabs-mode:nil
00868   fill-column:79
00869   End:
00870 */
00871 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=79 :
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines