00001
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
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
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
00060
00061 vector<Vec3f> BoxingTools::colors = vector<Vec3f>();
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
00137
00138
00139
00140 unsigned int local_columns = mColumns;
00141 if ( mRows < mColumns )
00142 {
00143
00144
00145
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
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
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
00246 gsl_matrix * ref_coords = gsl_matrix_calloc( mClasses, mColumns );
00247
00248
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
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
00288
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
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
00318 gsl_matrix * ref_coords = gsl_matrix_calloc( mClasses, mColumns );
00319
00320
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
00334
00335
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
00352
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
00369
00370
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
00394
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
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
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
00440 mappings_tally[*it] = 0;
00441 }
00442
00443
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
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);
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
00522 int xx = x+j;
00523 int yy = y+k;
00524
00525
00526 if ( xx >= image->get_xsize() || xx < 0 ) continue;
00527 if ( yy >= image->get_ysize() || yy < 0 ) continue;
00528
00529
00530 if ( xx == x && yy == y) continue;
00531
00532
00533 int square_length = k*k + j*j;
00534 if (square_length > radius_squared ) continue;
00535
00536
00537
00538 int idx = -1;
00539
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
00546
00547 idx -= 1;
00548
00549 if ( mode == SWARM_DIFFERENCE ) {
00550
00551 float val = peakval - image->get_value_at(xx,yy);
00552
00553
00554 if ( profile[idx] > val || profile[idx] == 0 ) profile[idx] = val;
00555 }
00556 else if (mode == SWARM_RATIO) {
00557
00558 float val = (peakval - image->get_value_at(xx,yy) ) / peakval;
00559
00560
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
00590 int xx = x+j;
00591 int yy = y+k;
00592
00593
00594 if ( xx >= image->get_xsize() || xx < 0 ) continue;
00595 if ( yy >= image->get_ysize() || yy < 0 ) continue;
00596
00597
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 {
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
00697 int xx = x+j;
00698 int yy = y+k;
00699
00700
00701 if ( xx >= image->get_xsize() || xx < 0 ) continue;
00702 if ( yy >= image->get_ysize() || yy < 0 ) continue;
00703
00704
00705 if ( xx == x && yy == y) continue;
00706
00707
00708 int square_length = k*k + j*j;
00709 if (square_length > radius_squared ) continue;
00710
00711
00712 if ( image->get_value_at(xx,yy) > peakval) return false;
00713
00714
00715
00716 int idx = -1;
00717
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
00724
00725 idx -= 1;
00726
00727 if (mode == SWARM_DIFFERENCE) {
00728
00729 float val = peakval - image->get_value_at(xx,yy);
00730
00731
00732 if ( profile[idx] > val || profile[idx] == 0 ) profile[idx] = val;
00733 }
00734 else if (mode == SWARM_RATIO) {
00735
00736 float val = (peakval - image->get_value_at(xx,yy) ) / peakval;
00737
00738
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
00770 int xx = x+j;
00771 int yy = y+k;
00772
00773 if ((k*k+j*j)>radius_squared) continue;
00774
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
00793 int xx = x+j;
00794 int yy = y+k;
00795
00796
00797 if ( xx >= map->get_xsize() || xx < 0 ) continue;
00798 if ( yy >= map->get_ysize() || yy < 0 ) continue;
00799
00800
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
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