00001
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 <ctime>
00038 #include <memory>
00039 #include "emdata.h"
00040 #include "analyzer.h"
00041 #include "sparx/analyzer_sparx.h"
00042 #include "util.h"
00043 #include "cmp.h"
00044 #include "sparx/lapackblas.h"
00045 #include "sparx/varimax.h"
00046
00047 using namespace EMAN;
00048
00049 namespace EMAN {
00050
00051 const string PCAsmall::NAME = "pca";
00052 const string PCAlarge::NAME = "pca_large";
00053 const string varimax::NAME = "varimax";
00054 const string KMeansAnalyzer::NAME = "kmeans";
00055 const string SVDAnalyzer::NAME = "svd_gsl";
00056
00057 template <> Factory < Analyzer >::Factory()
00058 {
00059 force_add<PCAsmall>();
00060 force_add<PCAlarge>();
00061 force_add<varimax>();
00062 force_add<KMeansAnalyzer>();
00063 force_add<SVDAnalyzer>();
00064 }
00065
00066 }
00067
00068 int Analyzer::insert_images_list(vector<EMData *> image_list)
00069 {
00070 vector<EMData *>::const_iterator iter;
00071 for(iter=image_list.begin(); iter!=image_list.end(); ++iter) {
00072 insert_image(*iter);
00073 }
00074 return 0;
00075 }
00076
00077 void KMeansAnalyzer::set_params(const Dict & new_params)
00078 {
00079 params = new_params;
00080 if (params.has_key("ncls")) ncls = params["ncls"];
00081 if (params.has_key("maxiter"))maxiter = params["maxiter"];
00082 if (params.has_key("minchange"))minchange = params["minchange"];
00083 if (params.has_key("mininclass"))mininclass = params["mininclass"];
00084 if (params.has_key("slowseed"))slowseed = params["slowseed"];
00085 if (params.has_key("verbose"))verbose = params["verbose"];
00086 if (params.has_key("calcsigmamean")) calcsigmamean=params["calcsigmamean"];
00087
00088 }
00089
00090 vector<EMData *> KMeansAnalyzer::analyze()
00091 {
00092 if (ncls<=1) return vector<EMData *>();
00093
00094
00095
00096 int nptcl=images.size();
00097 int nclstot=ncls;
00098 if (calcsigmamean) centers.resize(nclstot*2);
00099 else centers.resize(nclstot);
00100 if (mininclass<1) mininclass=1;
00101
00102 for (int i=0; i<nptcl; i++) images[i]->set_attr("is_ok_center",(int)5);
00103
00104 if (slowseed) {
00105 if (ncls>25) slowseed=ncls/25+1;
00106
00107
00108 }
00109
00110 for (int i=0; i<ncls; i++) {
00111
00112 centers[i]=images[Util::get_irand(0,nptcl-1)]->copy();
00113
00114 }
00115
00116 if (calcsigmamean) {
00117 for (int i=nclstot; i<nclstot*2; i++) centers[i]=new EMData(images[0]->get_xsize(),images[0]->get_ysize(),images[0]->get_zsize());
00118 }
00119
00120
00121 for (int i=0; i<maxiter; i++) {
00122 nchanged=0;
00123 reclassify();
00124 if (verbose) printf("iter %d> %d (%d)\n",i,nchanged,ncls);
00125 if (nchanged<minchange && ncls==nclstot) break;
00126 update_centers();
00127
00128 if (slowseed && i%3==2 && ncls<nclstot) {
00129 for (int j=0; j<slowseed && ncls<nclstot; j++) {
00130 centers[ncls]=0;
00131 ncls++;
00132 }
00133 reseed();
00134 }
00135 }
00136 update_centers(calcsigmamean);
00137
00138 return centers;
00139 }
00140
00141 void KMeansAnalyzer::update_centers(int sigmas) {
00142 int nptcl=images.size();
00143
00144 int * repr = new int[ncls];
00145
00146 for (int i=0; i<ncls; i++) {
00147 centers[i]->to_zero();
00148 if (sigmas) centers[i+ncls]->to_zero();
00149 repr[i]=0;
00150 }
00151
00152
00153 for (int i=0; i<nptcl; i++) {
00154 int cid=images[i]->get_attr("class_id");
00155 if ((int)images[i]->get_attr("is_ok_center")>0) {
00156 centers[cid]->add(*images[i]);
00157 if (sigmas) centers[cid+ncls]->addsquare(*images[i]);
00158 repr[cid]++;
00159 }
00160 }
00161
00162 for (int i=0; i<ncls; i++) {
00163
00164 if (repr[i]<mininclass) {
00165
00166
00167 for (int j=0; j<nptcl; j++) {
00168 if ((int)images[j]->get_attr("class_id")==i) images[i]->set_attr("is_ok_center",(int)images[i]->get_attr("is_ok_center")-1);
00169 }
00170
00171 delete centers[i];
00172 centers[i]=0;
00173 repr[i]=0;
00174 }
00175
00176 else {
00177 centers[i]->mult((float)1.0/(float)(repr[i]));
00178 centers[i]->set_attr("ptcl_repr",repr[i]);
00179 if (sigmas) {
00180 centers[i+ncls]->mult((float)1.0/(float)(repr[i]));
00181 centers[i+ncls]->subsquare(*centers[i]);
00182 centers[i+ncls]->process("math.sqrt");
00183 centers[i+ncls]->mult((float)1.0/(float)sqrt((float)repr[i]));
00184 }
00185
00186 }
00187 if (verbose>1) printf("%d(%d)\t",i,(int)repr[i]);
00188 }
00189
00190 if (verbose>1) printf("\n");
00191
00192 reseed();
00193
00194 delete [] repr;
00195 }
00196
00197
00198 void KMeansAnalyzer::reseed() {
00199 int nptcl=images.size();
00200 int i,j;
00201
00202
00203 for (i=0; i<ncls; i++) {
00204 if (!centers[i]) break;
00205 }
00206 if (i==ncls) return;
00207
00208
00209 vector<int> goodcen;
00210 for (int i=0; i<nptcl; i++) if ((int)images[i]->get_attr("is_ok_center")>0) goodcen.push_back(i);
00211
00212 if (goodcen.size()==0) throw UnexpectedBehaviorException("Kmeans ran out of valid center particles with the provided parameters");
00213
00214
00215 for (i=0; i<ncls; i++) {
00216 if (centers[i]) continue;
00217 j=Util::get_irand(0,goodcen.size()-1);
00218 centers[i]=images[j]->copy();
00219 centers[i]->set_attr("ptcl_repr",1);
00220 printf("reseed %d -> %d\n",i,j);
00221 }
00222
00223
00224 }
00225
00226
00227 void KMeansAnalyzer::reclassify() {
00228 int nptcl=images.size();
00229
00230 Cmp *c = Factory < Cmp >::get("sqeuclidean");
00231 for (int i=0; i<nptcl; i++) {
00232 float best=1.0e38f;
00233 int bestn=0;
00234 for (int j=0; j<ncls; j++) {
00235 float d=c->cmp(images[i],centers[j]);
00236
00237 if (d<best) { best=d; bestn=j; }
00238 }
00239 int oldn=images[i]->get_attr_default("class_id",0);
00240 if (oldn!=bestn) nchanged++;
00241 images[i]->set_attr("class_id",bestn);
00242 }
00243 delete c;
00244 }
00245
00246 #define covmat(i,j) covmat[ ((j)-1)*nx + (i)-1 ]
00247 #define imgdata(i) imgdata[ (i)-1 ]
00248 int PCAsmall::insert_image(EMData * image)
00249 {
00250 if(mask==0)
00251 throw NullPointerException("Null mask image pointer, set_params() first");
00252
00253 EMData *maskedimage = Util::compress_image_mask(image,mask);
00254
00255 int nx = maskedimage->get_xsize();
00256 float *imgdata = maskedimage->get_data();
00257 if (nx != ncov) {
00258 fprintf(stderr,"insert_image: something is wrong...\n");
00259 exit(1);
00260 }
00261
00262
00263 nimages++;
00264 for (int j = 1; j <= nx; j++)
00265 for (int i = 1; i<=nx; i++) {
00266 covmat(i,j) += imgdata(i)*imgdata(j);
00267 }
00268
00269 EMDeletePtr(maskedimage);
00270 return 0;
00271 }
00272 #undef covmat
00273
00274 #define eigvec(i,j) eigvec[(j)*ncov + (i)]
00275 vector<EMData*> PCAsmall::analyze()
00276 {
00277 float *eigvec;
00278 int status = 0;
00279
00280 eigval = (float*)calloc(ncov,sizeof(float));
00281 eigvec = (float*)calloc(ncov*ncov,sizeof(float));
00282 status = Util::coveig(ncov, covmat, eigval, eigvec);
00283
00284
00285
00286
00287 EMData *eigenimage = new EMData();
00288 eigenimage->set_size(ncov,1,1);
00289 float *rdata = eigenimage->get_data();
00290 for (int j = 1; j<= nvec; j++) {
00291 for (int i = 0; i < ncov; i++) rdata[i] = eigvec(i,ncov-j);
00292
00293 EMData* recons_eigvec = Util::reconstitute_image_mask(eigenimage,mask);
00294 recons_eigvec->set_attr( "eigval", eigval[j-1] );
00295 images.push_back(recons_eigvec);
00296 }
00297
00298 free(eigvec);
00299 EMDeletePtr(eigenimage);
00300
00301 return images;
00302 }
00303 #undef eigvec
00304
00305 void PCAsmall::set_params(const Dict & new_params)
00306 {
00307 params = new_params;
00308 mask = params["mask"];
00309 nvec = params["nvec"];
00310
00311
00312
00313 EMData *dummy = new EMData();
00314
00315 int nx = mask->get_xsize();
00316 int ny = mask->get_ysize();
00317 int nz = mask->get_zsize();
00318
00319 dummy->set_size(nx,ny,nz);
00320
00321 EMData *dummy1d = Util::compress_image_mask(dummy,mask);
00322 ncov = dummy1d->get_xsize();
00323 EMDeletePtr(dummy);
00324 EMDeletePtr(dummy1d);
00325
00326
00327 nimages = 0;
00328 covmat = (float*)calloc(ncov*ncov,sizeof(float));
00329 }
00330
00331
00332
00333
00334 int PCAlarge::insert_image(EMData * image)
00335 {
00336 if(mask==0)
00337 throw NullPointerException("Null mask image pointer, set_params() first");
00338
00339 EMData *maskedimage = Util::compress_image_mask(image,mask);
00340
00341 FILE *fp;
00342 string scratchfile = string("maskedimages.scratch");
00343
00344 fp = fopen(scratchfile.c_str(),"ab");
00345
00346 int nx = maskedimage->get_xsize();
00347 float *imgdata = maskedimage->get_data();
00348 fwrite(imgdata, sizeof(float), nx, fp);
00349 nimages++;
00350
00351 fclose(fp);
00352
00353 EMDeletePtr(maskedimage);
00354
00355 return 0;
00356 }
00357
00358 void PCAlarge::set_params(const Dict & new_params)
00359 {
00360 params = new_params;
00361 mask = params["mask"];
00362 nvec = params["nvec"];
00363
00364
00365
00366 EMData *dummy = new EMData();
00367
00368 int nx = mask->get_xsize();
00369 int ny = mask->get_ysize();
00370 int nz = mask->get_zsize();
00371
00372 dummy->set_size(nx,ny,nz);
00373
00374 EMData *dummy1d = Util::compress_image_mask(dummy,mask);
00375
00376 ncov = dummy1d->get_xsize();
00377
00378 EMDeletePtr(dummy);
00379 EMDeletePtr(dummy1d);
00380
00381 nimages = 0;
00382 }
00383
00384 #define qmat(i,j) qmat[((j)-1)*kstep + (i) -1]
00385 #define diag(i) diag[(i)-1]
00386 #define rdata(i) rdata[(i)-1]
00387 #define eigvec(i,j) eigvec[((j)-1)*ncov + (i)-1]
00388 #define eigval(i) eigval[(i)-1]
00389
00390 vector<EMData*> PCAlarge::analyze()
00391 {
00392 int status = 0;
00393 int ione = 1;
00394 float one = 1.0, zero = 0.0;
00395 char trans;
00396 float *eigvec;
00397 string scratchfile = string("maskedimages.scratch");
00398 char command[100];
00399
00400 printf("start analyzing..., ncov = %d\n", ncov);
00401
00402 float resnrm = 0.0;
00403
00404 if ( nvec > nimages || nvec ==0 ) nvec = nimages;
00405 int nx = ncov;
00406
00407
00408 int kstep = nvec + 20;
00409 if (kstep > nimages) kstep = nimages;
00410
00411 float *diag = new float[kstep];
00412 float *subdiag = new float[kstep-1];
00413 float *vmat = new float[nx*kstep];
00414
00415
00416 status = Lanczos(scratchfile, &kstep, diag, subdiag,
00417 vmat, &resnrm);
00418
00419
00420 #ifdef _WIN32
00421 if (_unlink(scratchfile.c_str()) == -1) {
00422 fprintf(stderr,"PCAlarge: cannot remove scratchfile\n");
00423 }
00424 #else
00425 sprintf(command,"rm -f %s\n", scratchfile.c_str());
00426 status = system(command);
00427 if (status != 0) {
00428 fprintf(stderr,"PCAlarge: cannot remove scratchfile\n");
00429 }
00430 #endif //_WIN32
00431
00432 char jobz[2] = "V";
00433 float *qmat = new float[kstep*kstep];
00434
00435 int lwork = 100 + 4*kstep + kstep*kstep;
00436 int liwork = 3+5*kstep;
00437
00438 float *work = new float[lwork];
00439 int *iwork = new int[liwork];
00440 int info = 0;
00441
00442
00443 sstevd_(jobz, &kstep, diag, subdiag, qmat, &kstep, work, &lwork,
00444 iwork, &liwork, &info);
00445
00446
00447 eigval = (float*)calloc(ncov,sizeof(float));
00448 eigvec = (float*)calloc(ncov*nvec,sizeof(float));
00449
00450 for (int j = 0; j < nvec; j++) {
00451 eigval[j] = diag(kstep-j);
00452 }
00453
00454
00455
00456
00457
00458 for (int j=1; j<=nvec; j++) {
00459 trans = 'N';
00460 sgemv_(&trans, &nx, &kstep, &one, vmat, &nx, &qmat(1,kstep-j+1),
00461 &ione, &zero, &eigvec(1,j), &ione);
00462 }
00463
00464
00465 EMData *eigenimage = new EMData();
00466 eigenimage->set_size(ncov,1,1);
00467 float *rdata = eigenimage->get_data();
00468 for (int j = 1; j<= nvec; j++) {
00469 for (int i = 1; i <= ncov; i++)
00470 rdata(i) = eigvec(i,j);
00471
00472 EMData* recons_eigvec = Util::reconstitute_image_mask(eigenimage,mask);
00473
00474 recons_eigvec->set_attr( "eigval", eigval[j-1] );
00475
00476 images.push_back( recons_eigvec );
00477 }
00478
00479 free(eigvec);
00480 EMDeletePtr(eigenimage);
00481
00482 return images;
00483 }
00484 #undef qmat
00485 #undef diag
00486 #undef rdata
00487 #undef eigvec
00488 #undef eigval
00489
00490 #define TOL 1e-7
00491 #define V(i,j) V[((j)-1)*imgsize + (i) - 1]
00492 #define v0(i) v0[(i)-1]
00493 #define Av(i) Av[(i)-1]
00494 #define subdiag(i) subdiag[(i)-1]
00495 #define diag(i) diag[(i)-1]
00496 #define hvec(i) hvec[(i)-1]
00497
00498 int PCAlarge::Lanczos(const string &maskedimages, int *kstep,
00499 float *diag, float *subdiag, float *V,
00500 float *beta)
00501 {
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531 int i, iter;
00532
00533 float alpha;
00534 int ione = 1;
00535 float zero = 0.0, one = 1.0, mone = -1.0;
00536 int status = 0;
00537
00538 char trans;
00539 int imgsize = 0;
00540 float *v0, *Av, *hvec, *htmp, *imgdata;
00541 FILE *fp=NULL;
00542
00543 if (nimages <= 0) {
00544 status = 2;
00545 goto EXIT;
00546 }
00547
00548 imgsize = ncov;
00549 if (nimages <= 0) {
00550 status = 3;
00551 goto EXIT;
00552 }
00553
00554 v0 = new float[imgsize];
00555 Av = new float[imgsize];
00556 hvec = new float[*kstep];
00557 htmp = new float[*kstep];
00558 imgdata = new float[imgsize];
00559
00560 if (v0 == NULL || Av == NULL || hvec == NULL ||
00561 htmp == NULL || imgdata == NULL) {
00562 fprintf(stderr, "Lanczos: failed to allocate v0,Av,hvec,htmp\n");
00563 status = -1;
00564 goto EXIT;
00565 }
00566
00567
00568 for ( i = 1; i <= imgsize; i++)
00569 {
00570 v0(i) = 1.0;
00571 Av(i) = 0.0;
00572 }
00573
00574
00575 *beta = snrm2_(&imgsize, v0, &ione);
00576 for (i = 1; i<=imgsize; i++)
00577 V(i,1) = v0(i) / (*beta);
00578
00579
00580 fp = fopen(maskedimages.c_str(),"rb");
00581 if (fp==NULL) {
00582 fprintf(stderr,"Lanczos: cannot open %s\n", maskedimages.c_str());
00583 }
00584 for (i = 0; i < nimages; i++) {
00585 fread(imgdata, sizeof(float), imgsize, fp);
00586 alpha = sdot_(&imgsize, imgdata, &ione, V, &ione);
00587 saxpy_(&imgsize, &alpha, imgdata, &ione, Av, &ione);
00588 }
00589 fclose(fp);
00590
00591
00592
00593 diag(1) = sdot_(&imgsize, V, &ione, Av, &ione);
00594 alpha = -diag(1);
00595 saxpy_(&imgsize, &alpha, V, &ione, Av, &ione);
00596
00597
00598 for ( iter = 2 ; iter <= *kstep ; iter++ ) {
00599 *beta = snrm2_(&imgsize, Av, &ione);
00600
00601 if (*beta < TOL) {
00602
00603 *kstep = iter;
00604 break;
00605 }
00606
00607 subdiag(iter-1) = *beta;
00608 for ( i = 1 ; i <= imgsize ; i++ ) {
00609 V(i,iter) = Av(i) / (*beta);
00610 }
00611
00612
00613 for (i = 0; i < imgsize; i++) Av[i] = 0;
00614 fp = fopen(maskedimages.c_str(),"rb");
00615 for (i = 0; i < nimages; i++) {
00616 fread(imgdata, sizeof(float), imgsize, fp);
00617 alpha = sdot_(&imgsize, imgdata, &ione, &V(1,iter), &ione);
00618 saxpy_(&imgsize, &alpha, imgdata, &ione, Av, &ione);
00619 }
00620 fclose(fp);
00621
00622
00623 trans = 'T';
00624 status = sgemv_(&trans, &imgsize, &iter, &one, V, &imgsize, Av, &ione,
00625 &zero , hvec , &ione);
00626 trans = 'N';
00627 status = sgemv_(&trans, &imgsize, &iter, &mone, V, &imgsize, hvec,
00628 &ione , &one , Av, &ione);
00629
00630
00631 trans = 'T';
00632 status = sgemv_(&trans, &imgsize, &iter, &one, V, &imgsize, Av, &ione,
00633 &zero , htmp , &ione);
00634 saxpy_(&iter, &one, htmp, &ione, hvec, &ione);
00635 trans = 'N';
00636 status = sgemv_(&trans, &imgsize, &iter, &mone, V, &imgsize, htmp,
00637 &ione , &one , Av, &ione);
00638 diag(iter) = hvec(iter);
00639 }
00640
00641 EMDeleteArray(v0);
00642 EMDeleteArray(Av);
00643 EMDeleteArray(hvec);
00644 EMDeleteArray(htmp);
00645 EMDeleteArray(imgdata);
00646
00647 EXIT:
00648 return status;
00649
00650 }
00651 #undef v0
00652 #undef Av
00653 #undef V
00654 #undef hvec
00655 #undef diag
00656 #undef subdiag
00657 #undef TOL
00658
00659 void varimax::set_params(const Dict & new_params)
00660 {
00661 params = new_params;
00662 m_mask = params["mask"];
00663
00664
00665
00666 EMData *dummy = new EMData();
00667
00668 int nx = m_mask->get_xsize();
00669 int ny = m_mask->get_ysize();
00670 int nz = m_mask->get_zsize();
00671
00672 dummy->set_size(nx,ny,nz);
00673
00674 EMData *dummy1d = Util::compress_image_mask(dummy,m_mask);
00675
00676 m_nlen = dummy1d->get_xsize();
00677 m_nfac = 0;
00678
00679 EMDeletePtr(dummy);
00680 EMDeletePtr(dummy1d);
00681 }
00682
00683 int varimax::insert_image(EMData* image)
00684 {
00685 if(m_mask==0)
00686 throw NullPointerException("Null mask image pointer, set_params() first");
00687
00688 EMData* img1d = Util::compress_image_mask(image,m_mask);
00689
00690 m_data.insert( m_data.end(), img1d->get_data(), img1d->get_data() + m_nlen );
00691
00692 m_nfac++;
00693
00694 Assert( (int)m_data.size() == m_nfac*m_nlen);
00695
00696 return 0;
00697 }
00698
00699 vector<EMData*> varimax::analyze()
00700 {
00701 int itmax = 10000;
00702 float eps = 1e-4f;
00703 int verbose = 1;
00704 float params[4];
00705 params[0] = 1.0;
00706 varmx( &m_data[0], m_nlen, m_nfac, IVARIMAX, params, NULL, itmax, eps, verbose);
00707
00708 vector<EMData*> images;
00709
00710 EMData* img1d = new EMData();
00711 img1d->set_size(m_nlen, 1, 1);
00712 for( int i=0; i < m_nfac; ++i )
00713 {
00714 float* imgdata = img1d->get_data();
00715
00716 int offset = i * m_nlen;
00717 for( int i=0; i < m_nlen; ++i )
00718 {
00719 imgdata[i] = m_data[offset+i];
00720 }
00721
00722 EMData* img = Util::reconstitute_image_mask(img1d,m_mask);
00723 images.push_back(img);
00724 }
00725
00726 EMDeletePtr(img1d);
00727
00728 return images;
00729 }
00730
00731 int SVDAnalyzer::insert_image(EMData * image)
00732 {
00733 if (mask==0)
00734 throw NullPointerException("Null mask image pointer, set_params() first");
00735
00736
00737 size_t totpix=mask->get_xsize()*mask->get_ysize()*mask->get_zsize();
00738 float *d=image->get_data();
00739 float *md=mask ->get_data();
00740 for (size_t i=0,j=0; i<totpix; ++i) {
00741 if (md[i]) {
00742 gsl_matrix_set(A,j,nsofar,d[i]);
00743 j++;
00744 }
00745 }
00746 nsofar++;
00747
00748 return 0;
00749 }
00750 #undef covmat
00751
00752 #define eigvec(i,j) eigvec[(j)*ncov + (i)]
00753 vector<EMData*> SVDAnalyzer::analyze()
00754 {
00755
00756 gsl_vector *work=gsl_vector_alloc(nimg);
00757 gsl_vector *S=gsl_vector_alloc(nimg);
00758 gsl_matrix *V=gsl_matrix_alloc(nimg,nimg);
00759 gsl_matrix *X=gsl_matrix_alloc(nimg,nimg);
00760
00761
00762
00763 gsl_linalg_SV_decomp_mod (A,X, V, S, work);
00764
00765
00766 vector<EMData*> ret;
00767
00768 float *md=mask->get_data();
00769 size_t totpix=mask->get_xsize()*mask->get_ysize()*mask->get_zsize();
00770 for (int k=0; k<nvec; k++) {
00771 EMData *img = new EMData;
00772 img->set_size(mask->get_xsize(),mask->get_ysize(),mask->get_zsize());
00773
00774 float *d=img->get_data();
00775 for (size_t i=0,j=0; i<totpix; ++i) {
00776 if (md[i]) {
00777 d[i]=(float)gsl_matrix_get(A,j,k);
00778 j++;
00779 }
00780 }
00781 img->set_attr( "eigval", gsl_vector_get(S,k));
00782 ret.push_back(img);
00783 }
00784
00785 gsl_vector_free(work);
00786 gsl_vector_free(S);
00787 gsl_matrix_free(V);
00788 gsl_matrix_free(X);
00789
00790 gsl_matrix_free(A);
00791 A=NULL;
00792 mask=NULL;
00793
00794 return ret;
00795 }
00796
00797 void SVDAnalyzer::set_params(const Dict & new_params)
00798 {
00799 params = new_params;
00800 mask = params["mask"];
00801 nvec = params["nvec"];
00802 nimg = params["nimg"];
00803
00804
00805 pixels=0;
00806 size_t totpix=mask->get_xsize()*mask->get_ysize()*mask->get_zsize();
00807 float *d=mask->get_data();
00808 for (size_t i=0; i<totpix; ++i) if (d[i]) ++pixels;
00809
00810 printf("%d,%d\n",pixels,nimg);
00811 A=gsl_matrix_alloc(pixels,nimg);
00812 nsofar=0;
00813 }
00814
00815
00816 void EMAN::dump_analyzers()
00817 {
00818 dump_factory < Analyzer > ();
00819 }
00820
00821 map<string, vector<string> > EMAN::dump_analyzers_list()
00822 {
00823 return dump_factory_list < Analyzer > ();
00824 }
00825
00826
00827
00828
00829
00830
00831