40 this->
read(argc, argv);
73 if (limitRclass < -1. || limitRclass > 1.)
75 "limitRclass should be a percentage: provide values between -100 and 100.");
89 if (limitRper < -100. || limitRper > 100.)
91 "limitRper should be a percentage: provide values between -100 and 100.");
126 addUsageLine(
"Make class average images and corresponding selfiles from angular_projection_matching docfiles.");
127 addSeeAlsoLine(
"angular_project_library, angular_projection_matching");
130 addParamsLine(
" -i <doc_file> : Docfile with assigned angles for all experimental particles");
131 addParamsLine(
" --lib <doc_file> : Docfile with angles used to generate the projection matching library");
132 addParamsLine(
" -o <root_name> : Output rootname for class averages and selfiles");
133 addParamsLine(
" [--split ] : Also output averages of random halves of the data");
134 addParamsLine(
" [--wien <img=\"\"> ] : Apply this Wiener filter to the averages");
135 addParamsLine(
" [--pad <factor=1.> ] : Padding factor for Wiener correction");
136 addParamsLine(
" [--save_images_assigned_to_classes] : Save images assigned te each class in output metadatas");
138 addParamsLine(
"==+ IMAGE SELECTION BASED ON INPUT DOCFILE (select one between: limit 0, F and R ==");
139 addParamsLine(
" [--select <col=\"maxCC\">] : Column to use for image selection (limit0, limitF or limitR)");
140 addParamsLine(
" [--limit0 <l0>] : Discard images below <l0>");
141 addParamsLine(
" [--limitF <lF>] : Discard images above <lF>");
142 addParamsLine(
" [--limitRclass <lRc>] : if (lRc>0 && lRc< 100): discard lowest <lRc> % in each class");
143 addParamsLine(
" : if (lRc<0 && lR>-100): discard highest <lRc> % in each class");
144 addParamsLine(
" [--limitRper <lRp>] : if (lRp>0 && lRp< 100): discard lowest <lRa> %");
145 addParamsLine(
" : if (lRp<0 && lRp>-100): discard highest <lRa> %");
147 addParamsLine(
" [--pcaSorting ] : Perform PCA sorting to obtain the average classes");
150 addParamsLine(
" [--iter <nr_iter=0>] : Number of iterations for re-alignment");
151 addParamsLine(
" [--Ri <ri=1>] : Inner radius to limit rotational search");
152 addParamsLine(
" [--Ro <r0=-1>] : Outer radius to limit rotational search");
154 addParamsLine(
" [--mpi_job_size <size=10>] : Number of images sent to a cpu in a single job ");
157 addExampleLine(
"Sample at default values and calculating output averages of random halves of the data",
false);
158 addExampleLine(
"xmipp_angular_class_average -i proj_match.doc --lib ref_angles.doc -o out_dir --split");
177 size_t jobId = 0, size;
178 size_t finishedNodes = 1;
179 bool whileLoop =
true;
182 int ctfGroup, ref3d, ref2d;
187 MPI_Probe(MPI_ANY_SOURCE, MPI_ANY_TAG, MPI_COMM_WORLD, &
status);
214 MPI_Recv(
nullptr, 0, MPI_INT, MPI_ANY_SOURCE,
TAG_I_AM_FREE, MPI_COMM_WORLD,
223 std::cerr <<
"Sending job to worker " <<
status.MPI_SOURCE <<std::endl;
231 MPI_Send(
nullptr, 0, MPI_INT,
status.MPI_SOURCE,
TAG_STOP, MPI_COMM_WORLD);
233 if (finishedNodes >=
node->size)
247 std::cerr <<
"Blocking. lockIndex: " << lockIndex <<
" | status.MPI_SOURCE: " <<
status.MPI_SOURCE
248 <<
" | lockArray[" << order_index<<
"," <<ref3d_index<<
"]: " <<
dAij(
lockArray,order_index,ref3d_index) << std::endl;
253 MPI_Send(
nullptr, 0, MPI_INT,
status.MPI_SOURCE,
272 std::cerr <<
"[" <<
node->rank <<
"] TAG_YES_YOU_MAY_WRITE.lockIndex: " << lockIndex <<std::endl;
273 std::cerr <<
"lockWeightIndexes[index_weight]" << lockWeightIndexes[
index_weight]<<std::endl;
290 std::cerr <<
"Unblocking. lockIndex: " << lockIndex <<
" | status.MPI_SOURCE: " <<
status.MPI_SOURCE << std::endl;
305 std::cerr <<
"WRONG TAG RECEIVED" << std::endl;
312 bool whileLoop =
true;
317 std::cerr <<
"[" <<
node->rank <<
"] Asking for a job " <<std::endl;
320 MPI_Send(
nullptr, 0, MPI_INT, 0,
TAG_I_AM_FREE, MPI_COMM_WORLD);
322 MPI_Probe(0, MPI_ANY_TAG, MPI_COMM_WORLD, &
status);
327 MPI_Recv(
nullptr, 0, MPI_INT, 0,
TAG_STOP,
359 int size =
ROUND(Def_3Dref_2Dref_JobNo[0]);
361 for(
int i=0;
i<size;
i++)
372 std::cerr<<
"["<<
node->rank<<
"]" 379 <<
" rot: " << Def_3Dref_2Dref_JobNo[
index_Rot]
380 <<
" tilt: " << Def_3Dref_2Dref_JobNo[
index_Tilt]
381 <<
" Sat node: " <<
node->rank
390 double psi, xshift, yshift,
w, w1, w2, scale;
392 int ref_number, this_image, ref3d, defGroup;
393 static int defGroup_last = 0;
394 int isplit, lockIndex;
404 order_number =
ROUND(Def_3Dref_2Dref_JobNo[index_Order]);
405 ref_number =
ROUND(Def_3Dref_2Dref_JobNo[index_2DRef]);
406 defGroup =
ROUND(Def_3Dref_2Dref_JobNo[index_DefGroup]);
407 ref3d =
ROUND(Def_3Dref_2Dref_JobNo[index_3DRef]);
409 lockIndex =
ROUND(Def_3Dref_2Dref_JobNo[index_jobId]);
411 if (
fn_wien !=
"" && defGroup_last != defGroup)
418 auxImg.
read(fn_wfilter);
420 defGroup_last = defGroup;
438 "Program should never execute this line, something went wrong");
446 if (_DF_temp.
size() > 1)
448 "Program should never execute this line, something went wrong");
449 if (_DF_temp.
size() == 0)
454 for (
size_t objId : _DF_temp.
ids())
457 img_ref.
read(fn_img);
462 std::vector<int> exp_number, exp_split;
463 std::vector<Image<double> > exp_imgs;
474 pcaAnalyzerSplit1.
clear();
475 pcaAnalyzerSplit2.
clear();
478 for (
size_t objId : _DF.
ids())
490 img().setXmippOrigin();
505 exp_imgs.push_back(img);
506 exp_number.push_back(this_image);
514 exp_split.push_back(isplit);
577 int static static_i=0;
588 std::cerr <<
"static_i:" << static_i << std::endl;
612 if (_DF.
size() <= 5 )
615 for (
size_t objId : _DF.
ids())
620 avg1() += img()*(max_cc);
628 for (
size_t objId : _DF.
ids())
641 img().setXmippOrigin();
654 double score = pcaAnalyzer.
getZscore(index);
658 avg1() += img()*(5-score)*max_cc;
659 w1 +=(5-score)*max_cc;
677 for (
size_t objId : _DF.
ids())
679 isplit = exp_split[
index];
692 img().setXmippOrigin();
707 if (pcaAnalyzerSplit1.
v.size() > 6 )
709 double score = pcaAnalyzerSplit1.
getZscore(index1);
712 avg1() += img()*(5-score)*max_cc;
713 w1 +=(5-score)*max_cc;
727 avg1() += img()*(w1);
734 if (pcaAnalyzerSplit2.
v.size() > 6 )
736 double score = pcaAnalyzerSplit2.
getZscore(index2);
739 avg2() += img()*(5-score)*max_cc;
740 w2 +=(5-score)*max_cc;
754 avg2() += img()*(w2);
768 avg() = avg1() + avg2();
776 reAlignClass(avg1, avg2, SFclass1, SFclass2, exp_imgs, exp_split,
777 exp_number, order_number, my_output);
796 avg() = avg1() + avg2();
805 img_ref().setXmippOrigin();
853 SFclassDiscarded,_DF, w1, w2, w, lockIndex);
910 double weight_old, weights1_old, weights2_old;
913 bool whileLoop =
true;
917 std::cerr <<
"[" <<
node->rank <<
"] May I write. lockIndex: " << lockIndex <<std::endl;
921 MPI_Probe(0, MPI_ANY_TAG, MPI_COMM_WORLD, &
status);
925 MPI_Recv(
nullptr, 0, MPI_INT, 0, MPI_ANY_TAG,
938 std::cerr <<
"[" <<
node->rank <<
"] TAG_YES_YOU_MAY_WRITE.lockIndex: " << lockIndex <<std::endl;
939 std::cerr <<
"lockWeightIndexes[index_weight]" << lockWeightIndexes[
index_weight]<<std::endl;
947 mpi_write(dirno, ref3dIndex, avg, avg1, avg2, w1, w2, weight_old,
948 weights1_old, weights2_old);
952 std::cerr <<
"process WRONG TAG RECEIVED at " <<
node->rank << std::endl;
957 std::cerr <<
"[" <<
node->rank <<
"] Finish writer " <<std::endl;
990 fn_tmp.
compose(dirno, fileNameStk);
1022 node->barrierWait();
1027 Iempty().setXmippOrigin();
1028 node->barrierWait();
1067 MPI_Bcast(&
master_seed,1,MPI_UNSIGNED ,0,MPI_COMM_WORLD);
1092 MPI_Bcast(&
ctfNum,1,MPI_INT,0,MPI_COMM_WORLD);
1102 node->barrierWait();
1157 int size = auxDF.
size();
1159 auxF1.
sort(auxDF,codifyLabel,asc,limit,0);
1177 const MDLabel myGroupByLabels[] =
1181 std::vector<MDLabel> groupbyLabels(myGroupByLabels,myGroupByLabels+6);
1199 int ref3d, defgroup;
1200 size_t order, jobCount, jobCount2;
1202 for (
size_t objId : auxMdJobList.
ids())
1216 dAkij(multiCounter, order, defgroup, ref3d) = jobCount2;
1229 for (
size_t objId : auxF1.
ids())
1235 if (
dAkij(multiCounter, order, defgroup, ref3d) > 0)
1238 dAkij(multiCounter, order, defgroup, ref3d)--;
1247 auxF1.
write(
"/tmp/inputfileAfterREmove.xmd");
1256 std::stringstream comment;
1260 comment <<
"Discarded images";
1262 comment <<
". Min value = " <<
limit0;
1264 comment <<
". Max value = " <<
limitF;
1270 comment <<
". Drop " <<
limitRclass*100 <<
"% of images (per class) with lower " << col_select
1271 <<
". If the ROUND(num_images_per_class * limitRclass / 100) == 0 then images are randomly dropped" 1272 <<
" so the percentage is satisfied";
1274 comment <<
". Drop " <<
limitRclass*100 <<
"% of images (per class) with higher " << col_select
1275 <<
". If the ROUND(num_images_per_class * limitRclass / 100) == 0 then images are randomly dropped" 1276 <<
" so the percentage is satisfied";
1283 auxDFsort.
write(fileNameXmd);
1314 "Incompatible padding factor for this Wiener filter");
1337 unlink((
fn_out + fn_tmp +
".xmd").c_str());
1338 unlink((
fn_out + fn_tmp +
".stk").c_str());
1343 unlink((
fn_out1 + fn_tmp +
".xmd").c_str());
1344 unlink((
fn_out1 + fn_tmp +
".stk").c_str());
1347 unlink((
fn_out2 + fn_tmp +
".xmd").c_str());
1348 unlink((
fn_out2 + fn_tmp +
".stk").c_str());
1352 unlink((
fn_out + fn_tmp +
"_discarded.xmd").c_str());
1362 FileName imageName, fileNameXmd, blockNameXmd;
1363 FileName imageNames1, fileNameXmds1, imageNames2, fileNameXmds2;
1365 size_t order_number;
1367 double weights1, weights2, weights;
1374 (
String)
"This file contains a list of class averages with direction projections and weights.";
1384 "Ref3D_%03lu@%s_Ref3D_%03lu.xmd",
i,
fn_out.c_str(),
i);
1390 const MDLabel myGroupByLabels[] =
1394 std::vector<MDLabel> groupbyLabels(myGroupByLabels,myGroupByLabels+5);
1401 for (
size_t objId : auxMd2.
ids())
1417 auxMd2.
write(fileNameXmd);
1421 for (
size_t j= 1;
j <=
Ndim;
j++)
1432 if(auxMd3.
size() != 0)
1448 "Ref3D_%03lu@%s_Ref3D_%03lu.xmd",
i,
fn_out1.c_str(),
i);
1451 "Ref3D_%03lu@%s_Ref3D_%03lu.xmd",
i,
fn_out2.c_str(),
i);
1453 auto idIter1 = auxMds1.
ids().begin();
1454 auto idIter2 = auxMds2.
ids().begin();
1456 const auto totalSize = auxMds1.
ids().end();
1457 for (; idIter1 != totalSize; ++idIter1, ++idIter2)
1494 auxMds1.
write(fileNameXmds1);
1496 auxMds2.
write(fileNameXmds2);
1503 const MDLabel myGroupByLabels[] =
1507 std::vector<MDLabel> groupbyLabels(myGroupByLabels,myGroupByLabels+6);
1521 Polar<std::complex<double> > &fP,
bool conjugated,
float xoff,
1536 std::vector<int> numbers,
size_t dirno,
double * my_output)
1539 std::vector<double> ccfs(splits.size());
1542 double maxcorr, new_xoff=0., new_yoff=0.;
1543 double w1=0., w2=0., opt_flip = 0., opt_psi = 0.;
1548 Mref = avg1() + avg2();
1554 auxImg.
write(
"ref.xmp");
1568 std::cerr<<
" entering iter "<<
iter<<std::endl;
1571 for (
size_t imgno = 0; imgno < imgs.size(); imgno++)
1576 getPolar(imgs[imgno](), fPimg,
false, (
float) -imgs[imgno].Xoff(),
1577 (
float) -imgs[imgno].Yoff());
1582 if (
corr(
k) > maxcorr)
1594 if (
corr(
k) > maxcorr)
1597 opt_psi =
realWRAP(360. - ang(
k), -180., 180.);
1613 xmipp_transformation::DONT_WRAP);
1623 imgs[imgno].setPsi(opt_psi);
1624 imgs[imgno].setFlip(opt_flip);
1625 imgs[imgno].setShifts(new_xoff, new_yoff);
1627 imgs[imgno].setShifts(-new_xoff, new_yoff);
1630 if (splits[imgno] == 0)
1635 else if (splits[imgno] == 1)
1647 Mref = avg1() + avg2();
1656 for (
size_t imgno = 0; imgno < imgs.size(); imgno++)
1658 if (splits[imgno] < 0)
1671 if (splits[imgno] == 0)
1676 else if (splits[imgno] == 1)
void mpi_produceSideInfo()
void saveDiscardedImages()
void selfWindow(int n0, int z0, int y0, int x0, int nF, int zF, int yF, int xF, T init_value=0)
void mpi_process_loop(double *Def_3Dref_2Dref_JobNo)
void addVector(const MultidimArray< float > &_v)
Add vector.
void setScale(double scale, const size_t n=0)
double getDoubleParam(const char *param, int arg=0)
static MDLabel str2Label(const String &labelName)
const int & getValue2(int) const
void addAndQuery(MDQuery &query)
#define REPORT_ERROR(nerr, ErrormMsg)
constexpr int index_2DRef
constexpr int index_lockIndex
#define FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(m)
void mpi_write(size_t dirno, int ref3dIndex, Image< double > avg, Image< double > avg1, Image< double > avg2, double w1, double w2, double old_w, double old_w1, double old_w2)
constexpr int index_DefGroup
void resizeNoCopy(const MultidimArray< T1 > &v)
MpiProgAngularClassAverage()
Just for debugging, situation that can't happens.
void mpi_writeFile(Image< double > avg, size_t dirno, FileName fileNameStk, double w_old)
int getSampleNoOuterRing() const
constexpr int AVG_OUPUT_SIZE
constexpr int TAG_YES_YOU_MAY_WRITE
constexpr int index_weights1
void rotationalCorrelation(const Polar< std::complex< double > > &M1, const Polar< std::complex< double > > &M2, MultidimArray< double > &angles, RotationalCorrelationAux &aux)
void reAlignClass(Image< double > &avg1, Image< double > &avg2, MetaData &SFclass1, MetaData &SFclass2, std::vector< Image< double > > imgs, std::vector< int > splits, std::vector< int > numbers, size_t dirno, double *my_output)
void write(const FileName &name="", size_t select_img=ALL_IMAGES, bool isStack=false, int mode=WRITE_OVERWRITE, CastWriteMode castMode=CW_CAST, int _swapWrite=0)
std::vector< MultidimArray< float > > v
double correlationIndex(const MultidimArray< T > &x, const MultidimArray< T > &y, const MultidimArray< int > *mask=NULL, MultidimArray< double > *Contributions=NULL)
RotationalCorrelationAux rotAux
void compose(const String &str, const size_t no, const String &ext="")
Polar_fftw_plans global_plans
FileName removeAllExtensions() const
void InverseFourierTransform(const MultidimArray< std::complex< double > > &in, MultidimArray< double > &out)
double rot(const size_t n=0) const
MultidimArray< bool > lockArray
MultidimArray< double > Mwien
constexpr int index_Count
double weight(const size_t n=0) const
constexpr int TAG_DO_NOT_DARE_TO_WRITE
void init_random_generator(int seed)
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
#define rotate(a, i, j, k, l)
void addSeeAlsoLine(const char *seeAlso)
double tilt(const size_t n=0) const
constexpr int index_Order
void addKeywords(const char *keywords)
String getExtension() const
int argc
Original command line arguments.
void setShifts(double xoff, double yoff, double zoff=0., const size_t n=0)
MultidimArray< double > weightArray
MultidimArray< double > weightArrays1
const char * getParam(const char *param, int arg=0)
void getPolarFromCartesianBSpline(const MultidimArray< T > &M1, int first_ring, int last_ring, int BsplineOrder=3, double xoff=0., double yoff=0., double oversample1=1., int mode1=FULL_CIRCLES)
void evaluateZScore(int NPCA, int Niter, bool trained=false, const char *filename="temp.txt", int numdesc=0)
constexpr int index_weight
void getPolar(MultidimArray< double > &img, Polar< std::complex< double > > &fP, bool conjugated=false, float xoff=0., float yoff=0.)
constexpr int lockWeightIndexesSize
void setFlip(bool flip, const size_t n=0)
void mpi_writeController(size_t dirno, Image< double > avg, Image< double > avg1, Image< double > avg2, const MetaDataDb &SFclass, const MetaDataDb &SFclass1, const MetaDataDb &SFclass2, const MetaDataDb &SFclassDiscarded, const MetaDataDb &_DF, double w1, double w2, double w, int lockIndex)
std::shared_ptr< MpiNode > node
#define FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(v)
#define dAkij(V, k, i, j)
constexpr int index_3DRef
#define DIRECT_MULTIDIM_ELEM(v, n)
void addExampleLine(const char *example, bool verbatim=true)
void createEmptyFile(const FileName &filename, int xdim, int ydim, int Zdim, size_t select_img, bool isStack, int mode, int _swapWrite, const MDRowVec *md)
constexpr int TAG_MAY_I_WRITE
MultidimArray< double > corr
void FourierTransform(const MultidimArray< double > &in, MultidimArray< std::complex< double > > &out)
constexpr int TAG_I_FINISH_WRITTING
MultidimArray< double > weightArrays2
const char * name() const
void fourierTransformRings(Polar< double > &in, Polar< std::complex< double > > &out, Polar_fftw_plans &plans, bool conjugated)
constexpr int index_ref3d
void mpi_process(double *Def_3Dref_2Dref_JobNo)
constexpr int TAG_I_AM_FREE
void setEulerAngles(double rot, double tilt, double psi, const size_t n=0)
constexpr int index_weights2
void applyWienerFilter(MultidimArray< double > &img)
void filterInputMetadata()
FileName removeBlockName() const
FileName withoutExtension() const
FourierTransformer local_transformer
#define FIRST_XMIPP_INDEX(size)
double psi(const double x)
String formatString(const char *format,...)
#define realWRAP(x, x0, xF)
bool checkParam(const char *param)
void getTransformationMatrix(Matrix2D< double > &A, bool only_apply_shifts=false, const size_t n=0)
int read(const FileName &name, DataMode datamode=DATA, size_t select_img=ALL_IMAGES, bool mapData=false, int mode=WRITE_READONLY)
bool do_save_images_assigned_to_classes
static String label2Str(const MDLabel &label)
void addUsageLine(const char *line, bool verbatim=false)
void initZeros(const MultidimArray< T1 > &op)
unsigned int randomize_random_generator()
void setWeight(double weight, const size_t n=0)
#define LAST_XMIPP_INDEX(size)
int getIntParam(const char *param, int arg=0)
void calculateFftwPlans(Polar_fftw_plans &out)
Incorrect value received.
void formatStringFast(String &str, const char *format,...)
void read(int argc, char **argv)
void addParamsLine(const String &line)
constexpr int index_jobId