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