EMAN2
boxingtools.cpp
Go to the documentation of this file.
00001 
00005 /*
00006  * Author: David Woolford, 04/15/2008 (woolford@bcm.edu)
00007  * Copyright (c) 2000-2006 Baylor College of Medicine
00008  *
00009  * This software is issued under a joint BSD/GNU license. You may use the
00010  * source code in this file under either license. However, note that the
00011  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00012  * so you are responsible for compliance with the licenses of these packages
00013  * if you opt to use BSD licensing. The warranty disclaimer below holds
00014  * in either instance.
00015  *
00016  * This complete copyright notice must be included in any revised version of the
00017  * source code. Additional authorship citations may be added, but existing
00018  * author citations must be preserved.
00019  *
00020  * This program is free software; you can redistribute it and/or modify
00021  * it under the terms of the GNU General Public License as published by
00022  * the Free Software Foundation; either version 2 of the License, or
00023  * (at your option) any later version.
00024  *
00025  * This program is distributed in the hope that it will be useful,
00026  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00027  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00028  * GNU General Public License for more details.
00029  *
00030  * You should have received a copy of the GNU General Public License
00031  * along with this program; if not, write to the Free Software
00032  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
00033  *
00034  * */
00035 
00036 
00037 #include "boxingtools.h"
00038 #include "exception.h"
00039 using namespace EMAN;
00040 
00041 
00042 #include <gsl/gsl_math.h>
00043 #include <gsl/gsl_vector.h>
00044 #include <gsl/gsl_linalg.h>
00045 
00046 // For random
00047 #include <cstdlib>
00048 #include <ctime>
00049 
00050 #include <iostream>
00051 using std::cerr;
00052 using std::cout;
00053 using std::endl;
00054 
00055 #include <string>
00056 using std::string;
00057 
00058 #include <algorithm>
00059 // find, min_element
00060 
00061 vector<Vec3f> BoxingTools::colors = vector<Vec3f>(); // static init
00062 BoxingTools::CmpMode BoxingTools::mode = SWARM_AVERAGE_RATIO;
00063 
00064 #define SVD_CLASSIFIER_DEBUG 0
00065 
00066 
00067 #if SVD_CLASSIFIER_DEBUG
00068  void printMatrix( const gsl_matrix * const A, const unsigned int rows, const unsigned int cols, const string& title = "" )
00069 {
00070         cout << "Printing matrix " << title << endl;
00071         cout << "It has " << rows << " rows and " << cols << " columns " << endl;
00072         for( unsigned int i = 0; i < rows; ++i )
00073         {
00074                 for( unsigned int j = 0; j < cols; ++j )
00075                 {
00076                         cout << gsl_matrix_get(A,i,j) << " ";
00077                 }
00078                 cout << endl;
00079         }
00080 
00081 }
00082 
00083 void printVector( const gsl_vector * const A, const unsigned int length, const string& title = "" )
00084 {
00085         cout << "Printing vector " << title << endl;
00086         for( unsigned int i = 0; i < length; ++i )
00087         {
00088                 cout << gsl_vector_get(A,i) << " ";
00089         }
00090         cout << endl;
00091 }
00092 
00093 void print_map( const map<unsigned int, unsigned int>& mapping )
00094 {
00095         for(map<unsigned int, unsigned int>::const_iterator it = mapping.begin(); it != mapping.end(); ++it)
00096         {
00097                 cout << it->first << " " << it->second << endl;
00098         }
00099 }
00100 #endif
00101 
00102 BoxSVDClassifier::BoxSVDClassifier(const vector<vector<float> >& data, const unsigned int& classes) :
00103                 mData(data), mClasses(classes)
00104 {
00105         setDims( mData );
00106 }
00107 
00108 
00109 BoxSVDClassifier::~BoxSVDClassifier()
00110 {
00111 
00112 }
00113 
00114 bool BoxSVDClassifier::setDims( const vector<vector<float> >& data )
00115 {
00116         mColumns = mData.size();
00117         vector<vector<float> >::const_iterator it = data.begin();
00118         mRows = it->size();
00119         it++;
00120         for( ; it != data.end(); ++it )
00121         {
00122                 if ( it->size() != mRows )
00123                 {
00124                         cerr << "ERROR: can not initial the BoxSVDClassifier with vectors of un-equal lengths " << endl;
00125                         cerr << "The vector lengths that did not agree were " <<  mRows << " and " << it->size() << endl;
00126                         return false;
00127                 }
00128         }
00129 
00130         return true;
00131 }
00132 
00133 
00134 map< unsigned int, unsigned int> BoxSVDClassifier::go()
00135 {
00136         //      This is done in the constructor
00137         //      setDims(mData);
00138 
00139 
00140         unsigned int local_columns = mColumns;
00141         if ( mRows < mColumns )
00142         {
00143 //              cerr << "Warning: gsl SVD works only when m > n, you have m = " << mRows << " and n = " << mColumns << endl;
00144                 // This local adaptation means things will proceed the same way even if there are more columns in A then rows
00145                 // Every input data is still classified, just the SVD eigenvectors are found using a subset of all the data
00146                 local_columns = mRows;
00147         }
00148 
00149         gsl_matrix * U = gsl_matrix_calloc( mRows, local_columns );
00150         gsl_matrix * A = gsl_matrix_calloc( mRows, mColumns );
00151         for ( unsigned int i = 0; i < mRows; ++i )
00152         {
00153                 for ( unsigned int j = 0; j < mColumns; ++j )
00154                 {
00155                         gsl_matrix_set( A, i, j, mData[j][i] );
00156                         if ( j < local_columns )
00157                                 gsl_matrix_set( U, i, j, mData[j][i] );
00158                 }
00159         }
00160 #if SVD_CLASSIFIER_DEBUG
00161         printMatrix( A, mRows, mColumns, "A" );
00162 #endif
00163 
00164         gsl_matrix * V = gsl_matrix_calloc( local_columns, local_columns );
00165         gsl_vector * S = gsl_vector_calloc( local_columns );
00166         gsl_vector * work = gsl_vector_calloc( local_columns );
00167 
00168         if ( gsl_linalg_SV_decomp (U, V, S, work) )
00169         {
00170                 cerr << "ERROR: gsl returned a non zero value on application of the SVD" << endl;
00171         }
00172 
00173 #if SVD_CLASSIFIER_DEBUG
00174         printMatrix( U, mRows, local_columns, "U" );
00175         printVector( S, local_columns, "S" );
00176         printMatrix( V, local_columns, local_columns, "V");
00177 #endif
00178 
00179         // normalize the columns of matrix A
00180         for ( unsigned int j = 0; j < mColumns; ++j )
00181         {
00182                 float norm = 0;
00183                 for ( unsigned int i = 0; i < mRows; ++i )
00184                 {
00185                         norm += (float)(gsl_matrix_get( A, i, j)*gsl_matrix_get( A, i, j));
00186                 }
00187                 norm = sqrtf(norm);
00188                 for ( unsigned int i = 0; i < mRows; ++i )
00189                 {
00190                         gsl_matrix_set( A, i, j, gsl_matrix_get(A,i,j)/norm);
00191                 }
00192         }
00193 
00194 #if SVD_CLASSIFIER_DEBUG
00195         for ( unsigned int j = 0; j < mColumns; ++j )
00196         {
00197                 double norm = 0;
00198                 for ( unsigned int i = 0; i < mRows; ++i )
00199                 {
00200                         norm += gsl_matrix_get( A, i, j)*gsl_matrix_get( A, i, j);
00201                 }
00202                 cout << "For column " << j << " the squared norm is " << norm << endl;
00203         }
00204 #endif
00205 
00206 
00207         gsl_matrix * svd_coords = gsl_matrix_calloc( mColumns, mColumns );
00208         // Correlate the columns of A with the columns of U and store the information in a martrix called svd_coords
00209         for ( unsigned int i = 0; i < mColumns; ++i )
00210         {
00211                 for ( unsigned int j = 0; j < local_columns; ++j )
00212                 {
00213                         double result = 0.0;
00214                         for ( unsigned int k = 0; k < mRows; ++k )
00215                         {
00216                                 result += gsl_matrix_get(A,k,i)*gsl_matrix_get(U,k,j);
00217                         }
00218                         gsl_matrix_set( svd_coords, i, j, result);
00219                 }
00220         }
00221 
00222 #if SVD_CLASSIFIER_DEBUG
00223         printMatrix( svd_coords, mColumns, mColumns, "svd_coords" );
00224 #endif
00225 
00226         map< unsigned int, unsigned int> grouping = randomSeedCluster(svd_coords, mColumns);
00227 
00228         for ( unsigned int i = 0; i < 20; ++ i )
00229         {
00230                 grouping = getIterativeCluster(svd_coords, grouping);
00231         }
00232 
00233         gsl_matrix_free(A);
00234         gsl_matrix_free(U);
00235         gsl_matrix_free(V);
00236         gsl_vector_free(S);
00237         gsl_vector_free(work);
00238         gsl_matrix_free(svd_coords);
00239 
00240         return grouping;
00241 }
00242 
00243 map< unsigned int, unsigned int> BoxSVDClassifier::getIterativeCluster(const gsl_matrix* const svd_coords, const map< unsigned int, unsigned int>& current_grouping)
00244 {
00245         // Space to store the reference vectors
00246         gsl_matrix * ref_coords = gsl_matrix_calloc( mClasses, mColumns );
00247 
00248         // Assumes there are a total of mClasses in the current_groupings mapping
00249         for(unsigned int i = 0; i < mClasses; ++i)
00250         {
00251                 unsigned int tally = 0;
00252                 for (map< unsigned int, unsigned int>::const_iterator it = current_grouping.begin(); it != current_grouping.end(); ++it )
00253                 {
00254                         if ( it->second == i )
00255                         {
00256                                 for( unsigned int j = 0; j < mColumns; ++j )
00257                                 {
00258                                         gsl_matrix_set(ref_coords, i, j, gsl_matrix_get( svd_coords, it->first, j ) + gsl_matrix_get( ref_coords, i, j));
00259                                 }
00260                                 ++tally;
00261                         }
00262 
00263                 }
00264                 // then normalize the the addition
00265                 if (tally != 0)
00266                         for( unsigned int j = 0; j < mColumns; ++j )
00267                 {
00268                         gsl_matrix_set(ref_coords, i, j, gsl_matrix_get( ref_coords, i, j )/((float) tally));
00269                 }
00270         }
00271 
00272         vector<vector<float> > distances = getDistances(svd_coords, ref_coords);
00273 
00274 #if SVD_CLASSIFIER_DEBUG
00275         cout << "The distance matrix is " << endl;
00276         for( unsigned int i = 0; i < distances.size(); ++i )
00277         {
00278                 for( unsigned int j = 0; j < distances[i].size(); ++j )
00279                 {
00280                         cout << distances[i][j] << " ";
00281                 }
00282                 cout << endl;
00283         }
00284 #endif
00285 
00286 
00287         // Finally decide which of the randomly chosen vectors is closest to each of the input vectors
00288         // and use that as the basis of the grouping
00289         map< unsigned int, unsigned int> return_map = getMapping(distances);
00290 
00291 #if SVD_CLASSIFIER_DEBUG
00292         cout << "Printing classification map" << endl;
00293         print_map(return_map);
00294 #endif
00295 
00296         gsl_matrix_free(ref_coords);
00297 
00298         return return_map;
00299 }
00300 
00301 
00302 map< unsigned int, unsigned int> BoxSVDClassifier::randomSeedCluster(const gsl_matrix* const svd_coords, unsigned int matrix_dims)
00303 {
00304         // Seed the random number generator
00305         srand(static_cast<unsigned int>(time(0)));
00306 
00307         vector<unsigned int> random_seed_indices;
00308         while ( random_seed_indices.size() < mClasses )
00309         {
00310                 unsigned int random_idx = static_cast<int>(((float)rand()/RAND_MAX)*matrix_dims);
00311                 if ( find( random_seed_indices.begin(), random_seed_indices.end(), random_idx ) == random_seed_indices.end() )
00312                 {
00313                         random_seed_indices.push_back( random_idx );
00314                 }
00315         }
00316 
00317         // Space to store the reference vectors
00318         gsl_matrix * ref_coords = gsl_matrix_calloc( mClasses, mColumns );
00319 
00320         // Put the reference vectors into a matrix to make the approach transparent to the reader
00321         for(unsigned int i = 0; i < random_seed_indices.size(); ++i)
00322         {
00323                 for( unsigned int j = 0; j < matrix_dims; ++j )
00324                 {
00325                         gsl_matrix_set(ref_coords, i, j, gsl_matrix_get( svd_coords, random_seed_indices[i], j ));
00326                 }
00327         }
00328 
00329 #if SVD_CLASSIFIER_DEBUG
00330         printMatrix( ref_coords, mClasses, matrix_dims, "Reference matrix in first grouping");
00331 #endif
00332 
00333         // accrue the distance data - this could be done more concisely, but there shouldn't be much cost
00334         // because the data should be fairl small. By more concisely I mean, the distance data would not need
00335         // to be stored, it could be determined without storing it in distances.
00336         vector<vector<float> > distances = getDistances(svd_coords, ref_coords);
00337 
00338 #if SVD_CLASSIFIER_DEBUG
00339         cout << "The distance matrix is " << endl;
00340         for( unsigned int i = 0; i < distances.size(); ++i )
00341         {
00342                 for( unsigned int j = 0; j < distances[i].size(); ++j )
00343                 {
00344                         cout << distances[i][j] << " ";
00345                 }
00346                 cout << endl;
00347         }
00348 #endif
00349 
00350 
00351         // Finally decide which of the randomly chosen vectors is closest to each of the input vectors
00352         // and use that as the basis of the grouping
00353         map< unsigned int, unsigned int> return_map = getMapping(distances);
00354 
00355 #if SVD_CLASSIFIER_DEBUG
00356         cout << "Printing classification map, randomly seeded" << endl;
00357         print_map(return_map);
00358 #endif
00359 
00360         gsl_matrix_free(ref_coords);
00361 
00362         return return_map;
00363 }
00364 
00365 
00366 vector<vector<float> > BoxSVDClassifier::getDistances( const gsl_matrix* const svd_coords, const gsl_matrix* const ref_coords)
00367 {
00368         // accrue the distance data - this could be done more concisely, but there shouldn't be much cost
00369         // because the data should be fairl small. By more concisely I mean, the distance data would not need
00370         // to be stored, it could be determined without storing it in distances.
00371         vector<vector<float> > distances;
00372         for (unsigned int i = 0; i < mColumns; ++i )
00373         {
00374                 vector<float> ith_distances;
00375                 for( unsigned int random_seed_idx = 0; random_seed_idx < mClasses; ++random_seed_idx )
00376                 {
00377                         float distance = 0;
00378                         for (unsigned int j = 0; j < mColumns; ++j )
00379                         {
00380                                 float value = (float)( (gsl_matrix_get( ref_coords, random_seed_idx, j) - gsl_matrix_get( svd_coords, i, j)) );
00381                                 distance += value * value;
00382                         }
00383                         ith_distances.push_back(distance);
00384                 }
00385                 distances.push_back(ith_distances);
00386         }
00387 
00388         return distances;
00389 }
00390 
00391 map< unsigned int, unsigned int> BoxSVDClassifier::getMapping(const vector<vector<float> >& distances)
00392 {
00393         // Finally decide which of the randomly chosen vectors is closest to each of the input vectors
00394         // and use that as the basis of the grouping
00395         map< unsigned int, unsigned int> return_map;
00396         unsigned int vector_idx = 0;
00397         for( vector<vector<float> >::const_iterator it = distances.begin(); it != distances.end(); ++it, ++vector_idx )
00398         {
00399                 vector<float>::const_iterator mIt = it->begin();
00400                 float min = *mIt;
00401                 unsigned int min_idx = 0;
00402                 for ( unsigned int current_idx = 0; mIt != it->end(); ++mIt, ++current_idx )
00403                 {
00404                         if ( *mIt < min )
00405                         {
00406                                 min = *mIt;
00407                                 min_idx = current_idx;
00408                         }
00409                 }
00410                 return_map[vector_idx] = min_idx;
00411         }
00412 
00413         return return_map;
00414 }
00415 
00416 map< unsigned int, unsigned int> BoxSVDClassifier::colorMappingByClassSize( const map< unsigned int, unsigned int>& grouping )
00417 {
00418 
00419         vector<unsigned int> current_mappings;
00420         // Get the extent of the current mappings
00421         for (map< unsigned int, unsigned int>::const_iterator it = grouping.begin(); it != grouping.end(); ++it )
00422         {
00423                 if ( find( current_mappings.begin(), current_mappings.end(), it->second ) == current_mappings.end() )
00424                 {
00425                         current_mappings.push_back( it->second );
00426                 }
00427         }
00428 
00429         if ( current_mappings.size() < 2 )
00430         {
00431                 cerr << "Error, cannot call colMappingByClassSize when less than 2 classes have been specified, I think you created " << current_mappings.size() << " classes " << endl;
00432                 throw;
00433         }
00434 
00435         // Record how many data points are in each class.
00436         map<unsigned int, unsigned int> mappings_tally;
00437         for( vector<unsigned int>::const_iterator it = current_mappings.begin(); it != current_mappings.end(); ++it )
00438         {
00439                 // First initialize each total to zero
00440                 mappings_tally[*it] = 0;
00441         }
00442 
00443         // Now do the actual counting
00444         for (map< unsigned int, unsigned int>::const_iterator it = grouping.begin(); it != grouping.end(); ++it )
00445         {
00446                 mappings_tally[it->second] += 1;
00447         }
00448 
00449         // find the largest tally
00450         unsigned int current_mapping_idx = 0;
00451         map< unsigned int, unsigned int> return_map;
00452         while ( mappings_tally.size() > 0 )
00453         {
00454 #if SVD_CLASSIFIER_DEBUG
00455                 cout << "Printing mappings_tally" << endl;
00456                 print_map(mappings_tally);
00457 #endif
00458 
00459                 map< unsigned int, unsigned int>::iterator it = mappings_tally.begin();
00460                 map< unsigned int, unsigned int>::iterator track_it = mappings_tally.begin();
00461                 unsigned int current_max = it->second;
00462                 unsigned int current_idx = it->first;
00463                 ++it;
00464                 for (; it != mappings_tally.end(); ++it )
00465                 {
00466                         if ( it->second > current_max )
00467                         {
00468                                 current_max = it->second;
00469                                 current_idx = it->first;
00470                                 track_it = it;
00471                         }
00472                 }
00473 
00474 #if SVD_CLASSIFIER_DEBUG
00475                 cout << "The mapping is " << current_idx << " to " << current_mapping_idx << endl;
00476 #endif
00477                 for (map< unsigned int, unsigned int>::const_iterator group_it = grouping.begin(); group_it != grouping.end(); ++group_it )
00478                 {
00479                         if ( group_it->second == current_idx )
00480                         {
00481                                 return_map[group_it->first] = current_mapping_idx;
00482                         }
00483                 }
00484 
00485                 mappings_tally.erase( current_idx );
00486 
00487                 current_mapping_idx++;
00488         }
00489 
00490 
00491 #if SVD_CLASSIFIER_DEBUG
00492         cout << "Printing adjusted classification map" << endl;
00493         print_map(return_map);
00494 #endif
00495 
00496 
00497         return return_map;
00498 }
00499 
00500 
00501 
00502 vector<float> BoxingTools::get_min_delta_profile(const EMData* const image, int x, int y, int radius)
00503 {
00504         float peakval = image->get_value_at(x,y);
00505 
00506         vector<float> profile(radius,0); // this sets the vectors size to radius, and the values to 0
00507         int radius_squared = radius*radius;
00508 
00509         static vector<float> squared_numbers;
00510         if ( (unsigned int)(radius+1) > squared_numbers.size() ) {
00511                 for(int i = squared_numbers.size(); i <= radius; ++i) {
00512                         squared_numbers.push_back((float)(i*i));
00513                 }
00514         }
00515 
00516         vector<unsigned int> tally;
00517         if (mode == SWARM_AVERAGE_RATIO) tally.resize(profile.size(),0);
00518 
00519         for(int k = -radius; k <= radius; ++k) {
00520                 for(int j = -radius; j <= radius; ++j) {
00521                         // Translate coordinates
00522                         int xx = x+j;
00523                         int yy = y+k;
00524 
00525                         // Protect against accessing pixels out of bounds
00526                         if ( xx >= image->get_xsize() || xx < 0 ) continue;
00527                         if ( yy >= image->get_ysize() || yy < 0 ) continue;
00528 
00529                         // We don't need to pay attention to the origin
00530                         if ( xx == x && yy == y) continue;
00531 
00532                         // Protect against vector accesses beyond the boundary
00533                         int square_length = k*k + j*j;
00534                         if (square_length > radius_squared ) continue;
00535 
00536                         // The idx is the radius, rounded down. This creates a certain type of pattern that
00537                         // can only really be explained visually...
00538                         int idx = -1;
00539                         // This little loop avoids the use of sqrtf
00540                         for(int i = 1; i < radius; ++i) {
00541                                 if ( square_length >= squared_numbers[i] && square_length <= squared_numbers[i+1] ) {
00542                                         idx = i;
00543                                 }
00544                         }
00545 //                      int idx = (int) sqrtf(k*k + j*j);
00546                         // decrement the idx, because the origin information is redundant
00547                         idx -= 1;
00548 
00549                         if ( mode == SWARM_DIFFERENCE ) {
00550                                 // Finally, get the drop relative to the origin
00551                                 float val = peakval - image->get_value_at(xx,yy);
00552 
00553                                 // Store it if the drop is smaller than the current value (or if there is no value)
00554                                 if ( profile[idx] > val || profile[idx] == 0 ) profile[idx] = val;
00555                         }
00556                         else if (mode == SWARM_RATIO) {
00557                                 // Finally, get the drop relative to the origin
00558                                 float val =  (peakval - image->get_value_at(xx,yy) ) / peakval;
00559 
00560                                 // Store it if the drop is smaller than the current value (or if there is no value)
00561                                 if ( profile[idx] > val || profile[idx] == 0 ) profile[idx] = val;
00562                         }
00563                         else if (mode == SWARM_AVERAGE_RATIO) {
00564                                 profile[idx] += image->get_value_at(xx,yy);
00565                                 tally[idx]++;
00566                         }
00567 
00568                 }
00569         }
00570 
00571         if (mode == SWARM_AVERAGE_RATIO) {
00572                 for(unsigned int i = 0; i < profile.size(); ++i) {
00573                         if (tally[i] != 0) {
00574                                 profile[i] /= static_cast<float>(tally[i]);
00575                                 profile[i] = (peakval - profile[i] ) / peakval;
00576                         }
00577                 }
00578         }
00579 
00580         return profile;
00581 }
00582 
00583 bool BoxingTools::is_local_maximum(const EMData* const image, int x, int y, int radius,EMData* exclusion_map)
00584 {
00585         float peakval = image->get_value_at(x,y);
00586         int radius_squared = radius*radius;
00587         for(int k = -radius; k <= radius; ++k) {
00588                 for(int j = -radius; j <= radius; ++j) {
00589                         // Translate coordinates
00590                         int xx = x+j;
00591                         int yy = y+k;
00592 
00593 //                      // Protect against accessing pixel out of bounds
00594                         if ( xx >= image->get_xsize() || xx < 0 ) continue;
00595                         if ( yy >= image->get_ysize() || yy < 0 ) continue;
00596 
00597                         // We don't need to pay attention to the origin
00598                         if ( xx == x && yy == y) continue;
00599 
00600                         if ((k*k+j*j)>radius_squared) continue;
00601 
00602                         if ( image->get_value_at(xx,yy) > peakval)  return false;
00603                 }
00604         }
00605 
00606         set_radial_non_zero(exclusion_map,x,y,radius);
00607 
00608         return true;
00609 
00610 }
00611 
00612 vector<IntPoint> BoxingTools::auto_correlation_pick(const EMData* const image, float threshold, int radius, const vector<float>& profile, EMData* const exclusion, const int cradius, int mode)
00613 {
00614         if (mode < 0 || mode > 2 ) {
00615                 throw InvalidValueException(mode,"Error, the mode can only be 0,1, or 2.");
00616         }
00617 
00618         if ( radius < 0) {
00619                 throw InvalidValueException(radius,"Radius must be greater than 1");
00620         }
00621 
00622         if ( cradius < 0) {
00623                 throw InvalidValueException(cradius,"CRadius must be greater than 1");
00624         }
00625 
00626 
00627         int nx = image->get_xsize();
00628         int ny = image->get_ysize();
00629 
00630         vector<IntPoint> solution;
00631 
00632         int r = radius+1;
00633 
00634         for(int j = r; j < ny-r;++j) {
00635                 for(int k = r; k < nx-r;++k) {
00636 
00637                         if (exclusion->get_value_at(k,j) != 0 ) continue;
00638 
00639                         if (image->get_value_at(k,j) < threshold) continue;
00640                         if ( mode == 0 ) {
00641                                 solution.push_back(IntPoint(k,j));
00642                                 set_radial_non_zero(exclusion,k,j,radius);
00643                                 continue;
00644                         }
00645 
00646                         vector<float> p(r,0);
00647 
00648                         if (hi_brid(image,k,j,r,exclusion,p)) {
00649                                 if ( mode == 1 ) {
00650                                         if (p[cradius] >= profile[cradius]) {
00651                                                 solution.push_back(IntPoint(k,j));
00652                                                 set_radial_non_zero(exclusion,k,j,radius);
00653                                                 continue;
00654                                         }
00655                                 }
00656                                 else /* mode == 2 */{
00657                                         bool bad = false;
00658                                         for (int ii = 0; ii <= cradius; ++ii) {
00659                                                 if (p[ii] < profile[ii]) {
00660                                                         bad = true;
00661                                                         break;
00662                                                 }
00663                                         }
00664                                         if (bad) continue;
00665                                         solution.push_back(IntPoint(k,j));
00666                                         set_radial_non_zero(exclusion,k,j,radius);
00667                                 }
00668 
00669 
00670                         }
00671                 }
00672         }
00673         return solution;
00674 }
00675 
00676 
00677 bool BoxingTools::hi_brid(const EMData* const image, int x, int y, int radius,EMData* const exclusion_map, vector<float>& profile)
00678 {
00679 
00680         float peakval = image->get_value_at(x,y);
00681 
00682         int radius_squared = radius*radius;
00683 
00684         static vector<float> squared_numbers;
00685         if ( (unsigned int)(radius+1) > squared_numbers.size() ) {
00686                 for(int i = squared_numbers.size(); i <= radius; ++i) {
00687                         squared_numbers.push_back((float)(i*i));
00688                 }
00689         }
00690 
00691         vector<unsigned int> tally;
00692         if (mode == SWARM_AVERAGE_RATIO) tally.resize(profile.size(),0);
00693 
00694         for(int k = -radius; k <= radius; ++k) {
00695                 for(int j = -radius; j <= radius; ++j) {
00696                         // Translate coordinates
00697                         int xx = x+j;
00698                         int yy = y+k;
00699 
00700                         // Protect against accessing pixels out of bounds
00701                         if ( xx >= image->get_xsize() || xx < 0 ) continue;
00702                         if ( yy >= image->get_ysize() || yy < 0 ) continue;
00703 
00704                         // We don't need to pay attention to the origin
00705                         if ( xx == x && yy == y) continue;
00706 
00707                         // Protect against vector accesses beyond the boundary
00708                         int square_length = k*k + j*j;
00709                         if (square_length > radius_squared ) continue;
00710 
00711                         // It's not a local maximum!
00712                         if ( image->get_value_at(xx,yy) > peakval)  return false;
00713 
00714                         // The idx is the radius, rounded down. This creates a certain type of pattern that
00715                         // can only really be explained visually...
00716                         int idx = -1;
00717                         // This little loop avoids the use of sqrtf
00718                         for(int i = 1; i < radius; ++i) {
00719                                 if ( square_length >= squared_numbers[i] && square_length <= squared_numbers[i+1] ) {
00720                                         idx = i;
00721                                 }
00722                         }
00723 //                      int idx = (int) sqrtf(k*k + j*j);
00724                         // decrement the idx, because the origin information is redundant
00725                         idx -= 1;
00726 
00727                         if (mode == SWARM_DIFFERENCE) {
00728                                 // Finally, get the drop relative to the origin
00729                                 float val = peakval - image->get_value_at(xx,yy);
00730 
00731                                 // Store it if the drop is smaller than the current value (or if there is no value)
00732                                 if ( profile[idx] > val || profile[idx] == 0 ) profile[idx] = val;
00733                         }
00734                         else if (mode == SWARM_RATIO) {
00735                                 // Finally, get the drop relative to the origin
00736                                 float val =  (peakval - image->get_value_at(xx,yy) ) / peakval;
00737 
00738                                 // Store it if the drop is smaller than the current value (or if there is no value)
00739                                 if ( profile[idx] > val || profile[idx] == 0 ) profile[idx] = val;
00740                         }
00741                         else if (mode == SWARM_AVERAGE_RATIO) {
00742                                 profile[idx] += image->get_value_at(xx,yy);
00743                                 tally[idx]++;
00744                         }
00745 
00746                 }
00747         }
00748 
00749         if (mode == SWARM_AVERAGE_RATIO) {
00750                 for(unsigned int i = 0; i < profile.size(); ++i) {
00751                         if (tally[i] != 0) {
00752                                 profile[i] /= static_cast<float>(tally[i]);
00753                                 profile[i] = (peakval - profile[i] ) / peakval;
00754                         }
00755                 }
00756         }
00757 
00758         set_radial_non_zero(exclusion_map,x,y,radius);
00759 
00760         return true;
00761 }
00762 
00763 
00764 void BoxingTools::set_radial_non_zero(EMData* const exclusion, int x, int y, int radius)
00765 {
00766         int radius_squared = radius*radius;
00767         for(int k = -radius; k <= radius; ++k) {
00768                 for(int j = -radius; j <= radius; ++j) {
00769                         // Translate coordinates
00770                         int xx = x+j;
00771                         int yy = y+k;
00772 
00773                         if ((k*k+j*j)>radius_squared) continue;
00774                         // Protect against accessing pixel out of bounds
00775                         if ( xx >= exclusion->get_xsize() || xx < 0 ) continue;
00776                         if ( yy >= exclusion->get_ysize() || yy < 0 ) continue;
00777 
00778                         exclusion->set_value_at(xx,yy,1);
00779                 }
00780         }
00781 }
00782 
00783 IntPoint BoxingTools::find_radial_max(const EMData* const map, int x, int y, int radius)
00784 {
00785         float currentmax = map->get_value_at(x,y);
00786 
00787         IntPoint soln(x,y);
00788 
00789         int radius_squared = radius*radius;
00790         for(int k = -radius; k <= radius; ++k) {
00791                 for(int j = -radius; j <= radius; ++j) {
00792                         // Translate coordinates
00793                         int xx = x+j;
00794                         int yy = y+k;
00795 
00796                         // Protect against accessing pixels out of bounds
00797                         if ( xx >= map->get_xsize() || xx < 0 ) continue;
00798                         if ( yy >= map->get_ysize() || yy < 0 ) continue;
00799 
00800                         // Protect against vector accesses beyond the boundary
00801                         int square_length = k*k + j*j;
00802                         if (square_length > radius_squared ) continue;
00803 
00804                         float val = map->get_value_at(xx,yy);
00805 
00806                         if (val > currentmax) {
00807                                 currentmax = val;
00808                                 soln[0] = xx;
00809                                 soln[1] = yy;
00810                         }
00811                 }
00812         }
00813 
00814         return soln;
00815 }
00816 
00817 
00818 void BoxingTools::set_region( EMData* const image, const EMData* const mask, const int x, const int y, const float& val) {
00819 
00820         // Works only in 2D
00821         int inx = image->get_xsize();
00822         int iny = image->get_ysize();
00823         int mnx = mask->get_xsize();
00824         int mny = mask->get_ysize();
00825 
00826         int startx = x-mnx/2;
00827         int endx =startx + mnx;
00828         int xoffset = 0;
00829         if (startx < 0) {
00830                 xoffset = abs(startx);
00831                 startx = 0;
00832         }
00833         if (endx > inx) endx = inx;
00834 
00835         int starty = y-mny/2;
00836         int endy =starty + mny;
00837         int yoffset = 0;
00838         if (starty < 0) {
00839                 yoffset = abs(starty);
00840                 starty = 0;
00841         }
00842         if (endy > iny) endy = iny;
00843 
00844 
00845         for (int j = starty; j < endy; ++j ) {
00846                 for (int i = startx; i < endx; ++i) {
00847                         if (mask->get_value_at(xoffset+i-startx,yoffset+j-starty) != 0 ) {
00848                                 image->set_value_at(i,j,val);
00849                         }
00850                 }
00851         }
00852 }
00853 
00854 map<unsigned int, unsigned int> BoxingTools::classify(const vector<vector<float> >& data, const unsigned int& classes)
00855 {
00856         BoxSVDClassifier classifier(data, classes);
00857         map< unsigned int, unsigned int> mapping = classifier.go();
00858 
00859         mapping = BoxSVDClassifier::colorMappingByClassSize( mapping );
00860 
00861         return mapping;
00862 
00863 }
00864 
00865 
00866 Vec3f BoxingTools::get_color( const unsigned int index )
00867 {
00868         if ( colors.size() == 0 ) {
00869                 colors.push_back(Vec3f(0,0,0));
00870                 colors.push_back(Vec3f(0,0,1));
00871                 colors.push_back(Vec3f(0,1,0));
00872                 colors.push_back(Vec3f(1,0,0));
00873                 colors.push_back(Vec3f(1,1,0));
00874                 colors.push_back(Vec3f(1,0,1));
00875                 colors.push_back(Vec3f(0,1,1));
00876                 colors.push_back(Vec3f(1,1,1));
00877         }
00878         if ( index >= colors.size() )
00879         {
00880                 while ( colors.size() <= index )
00881                 {
00882                         bool found = false;
00883                         while ( !found )
00884                         {
00885                                 unsigned int random_idx = rand() % colors.size();
00886                                 unsigned int random_idx2 = rand() % colors.size();
00887                                 while ( random_idx2 == random_idx )
00888                                 {
00889                                         random_idx2 = rand() % colors.size();
00890                                 }
00891 
00892                                 Vec3f result = (colors[random_idx] + colors[random_idx2])/2.0;
00893                                 if ( find( colors.begin(), colors.end(), result ) == colors.end() )
00894                                 {
00895                                         colors.push_back( result );
00896                                         found = true;
00897                                 }
00898                         }
00899                 }
00900         }
00901         return colors[index];
00902 }
00903