39 const FileName &model_name,
const std::vector<MDRowSql> &vMicList)
138 size_t(
YSIZE(inputMicrograph)),
139 size_t(
XSIZE(inputMicrograph)));
150 filter.
w2=(filter.
w1)+0.025;
155 Iaux=inputMicrograph;
163 for (
size_t i=0;
i<MD.size();
i++)
175 FileName micFile, posFile, preMicFile;
179 for (
size_t i=0;
i<MD.size();
i++)
182 if (strcmp(micFile.c_str(),preMicFile.c_str()))
189 if (
i==(MD.size()-1))
244 for (
int i=0;
i<num_part;
i++)
256 for (
size_t j=0;
j<
NSIZE(IpolarCorr);
j++)
289 std::vector<Particle2> negativeSamples;
290 size_t numNegInv = 0;
291 size_t numNegPart = 0;
299 std::random_device rd;
300 std::mt19937 gen(rd());
301 std::uniform_real_distribution<> dist(0.0, 1.0);
302 randomValues.
resize(1,1,1,negativeSamples.size());
313 numNegatives=num_part*2;
326 for (
int i=0;
i<numNegatives;
i++)
334 for (
size_t j=0;
j<
NSIZE(IpolarCorr);
j++)
444 for (
int j=0;
j<12;
j++)
488 std::cout <<
"Read the data for PCA, Rot PCA, Average ..." << std::endl;
498 std::cout <<
"Read rotational PCA model ..." << std::endl;
527 std::ofstream fhTrain;
549 std::vector<Particle2> positionArray;
562 positionArray.clear();
574 auto num=(int)(positionArray.size()*(proc_prec/100.0));
578 for (
int k=0;
k<num;
k++)
587 int j=positionArray[
k].x;
588 int i=positionArray[
k].y;
657 std::vector<Particle2> positionArray;
662 auto num=(int)(positionArray.size()*(proc_prec/100.0));
664 for (
int k=0;
k<num;
k++)
673 int j=positionArray[
k].x;
674 int i=positionArray[
k].y;
738 auto num=(int)(positionArray.size()*(proc_prec/100.0));
740 for (
int k=0;
k<num;
k++)
744 int j=positionArray[
k].x;
745 int i=positionArray[
k].y;
750 buildVector(IpolarCorr,staticVec,featVec,pieceImage);
759 this->picker = picker;
760 waitingForResults =
false;
771 this->fnmicrograph = fnMic;
791 positionArray.clear();
794 if (waitingForResults)
804 this->picker->generateFeatVec(fnmicrograph,
proc_prec, positionArray);
811 setMicrograph(fnMic, proc_prec);
821 waitingForResults =
true;
830 size_t lastObjId =
micList.size();
836 if (!strcmp(currentMic.c_str(),fnmicrograph.c_str()))
839 fnmicrograph=currentMic;
858 int idx=0,enabled,
x,
y;
864 for (
size_t i=0;
i<removedParticlesMD.size();
i++)
872 removedParticlesMD[
i].getValue(
MDL_COST,cost);
884 train(addedParticlesMD,
true,0,0,0,0);
899 for (
size_t objId : removedParticlesMD.
ids())
912 int limit = cntNeg + yDataSet;
913 for (
int n=yDataSet;
n<limit;
n++)
917 for (
size_t objId : removedParticlesMD.
ids())
939 const MultidimArray< std::complex< double > > &fourierPolarStack,
957 double dx=(p1.
x-p2.
x);
958 double dy= (p1.
y-p2.
y);
959 return sqrt(dx*dx+dy*dy);
971 for (
int i=1;
i<num_part;
i++)
1008 int nF =
NSIZE(fourierPolarStack);
1011 for (
int n=0;
n<nF;++
n)
1025 features.
resize(1,1,1,12);
1030 inputVec.
sort(sortedVec);
1032 for (
int i=2;
i<12;
i++)
1064 for (
int j=0;
j<12;
j++)
1107 size_t iBlockMax =
XSIZE(pcaBasis) / 4;
1108 for (
size_t i = 0;
i < iBlockMax;
i++)
1110 dotProduct += (*ptr2++) * (*ptr3++);
1111 dotProduct += (*ptr2++) * (*ptr3++);
1112 dotProduct += (*ptr2++) * (*ptr3++);
1113 dotProduct += (*ptr2++) * (*ptr3++);
1115 for (
size_t i = iBlockMax * 4;
i <
XSIZE(pcaBasis); ++
i)
1116 dotProduct += (*ptr2++) * (*ptr3++);
1126 FileName fnPositiveInvariatn=fnPositiveFeat+
"_invariant_Positive.stk";
1127 positiveInvariant.
read(fnPositiveInvariatn,
HEADER);
1136 positiveInvariant().getImage(pcaVec);
1167 if (numClassifier==1)
1195 if (!fn_Invariant.
exists())
1218 positiveInvariant.
readMapped(fn_Invariant,
k*num_correlation+
i+1);
1219 positiveInvariant().getImage(vec);
1230 positiveInvariant().getImage(vec);
1233 for (
int j=0;
j<12;
j++)
1251 int limit=0, size=0, newSize=0;
1309 FileName fnPositiveInvariatn=fnInvariantFeat+
"_Positive.stk";
1310 FileName fnPositiveParticles=fnParticles+
"_Positive.stk";
1313 for (
int i=0;
i<num_part;
i++)
1353 std::vector<Particle2> negativeSamples;
1354 std::vector<Micrograph::Point> positionVec;
1360 FileName fnNegativeInvariatn=fnInvariantFeat+
"_Negative.stk";
1361 FileName fnNegativeParticles=fnParticles+
"_Negative.stk";
1366 if (negativeSamples.size()==0)
1369 std::random_device rd;
1370 std::mt19937 gen(rd());
1371 std::uniform_real_distribution<> dist(0.0, 1.0);
1372 randomValues.
resize(1,1,1,negativeSamples.size());
1383 numNegatives=int(num_part/2);
1384 for (
int i=0;
i<numNegatives;
i++)
1411 int startX, startY, endX, endY;
1417 filter.
window(particleImage, startY, startX, endY, endX);
1424 int endX,endY,startX,startY;
1444 for (
int i=startY;
i<endY;
i=
i+gridStep)
1445 for (
int j= startX;
j<endX;
j=
j+gridStep)
1452 negativePosition.push_back(negSample);
1471 std::ifstream fhTrain;
1472 fhTrain.open(fn.c_str());
1477 for (
int i=0;
i<
y;
i++)
1480 for (
int j= 0;
j<
x;
j++)
1491 std::ofstream fhTrain;
1513 fhTrain.open(fn.c_str());
1533 FileName fnPcaModel=fn_root+
"_pca_model.stk";
1665 positionArray.push_back(p);
1667 for (
size_t i=0;
i<positionArray.size();++
i)
1668 for (
size_t j=0;
j<positionArray.size()-
i-1;
j++)
1669 if (positionArray[
j].cost<positionArray[
j + 1].cost)
1671 p=positionArray[
j+1];
1672 positionArray[
j+1]=positionArray[
j];
1705 for (
int deg=3;deg<360;deg+=3)
1710 avgRotatedLarge+=avgRotated;
1712 avgRotatedLarge/=120;
1730 for (
int deg=3;deg<360;deg+=3)
1735 avgRotatedLarge=avgRotated;
1756 mode = getParam(
"--mode");
1757 if (
mode ==
"buildinv")
1759 fn_train = getParam(
"--mode", 1);
1761 fn_root = getParam(
"--outputRoot");
1762 autoPicking = std::make_unique<AutoParticlePicking2>();
1763 autoPicking->readParams(
this);
1780 program->
addParamsLine(
" -i <micrograph> : Micrograph image");
1781 program->
addParamsLine(
" --outputRoot <rootname> : Output rootname");
1784 program->
addParamsLine(
" try : Try to autoselect within the training phase.");
1785 program->
addParamsLine(
" train : Train the classifier using the invariants features.");
1786 program->
addParamsLine(
" : <rootname>_auto_feature_vectors.txt contains the particle structure created by this program when used in automatic selection mode");
1787 program->
addParamsLine(
" : <rootname>_false_positives.xmd contains the list of false positives among the automatically picked particles");
1789 program->
addParamsLine(
" buildinv <posfile=\"\"> : posfile contains the coordinates of manually picked particles");
1790 program->
addParamsLine(
" --model <model_rootname> : Bayesian model of the particles to pick");
1791 program->
addParamsLine(
" --particleSize <size> : Particle size in pixels");
1792 program->
addParamsLine(
" [--thr <p=1>] : Number of threads for automatic picking");
1793 program->
addParamsLine(
" [--fast] : Perform a fast preprocessing of the micrograph (Fourier filter instead of Wavelet filter)");
1794 program->
addParamsLine(
" [--in_core] : Read the micrograph in memory");
1795 program->
addParamsLine(
" [--filter_num <n=6>] : The number of filters in filter bank");
1796 program->
addParamsLine(
" [--NPCA <n=4>] : The number of PCA components");
1797 program->
addParamsLine(
" [--NCORR <n=2>] : The number of PCA components");
1798 program->
addParamsLine(
" [--autoPercent <n=90>] : The number of PCA components");
1799 program->
addExampleLine(
"Automatically select particles during training:",
false);
1800 program->
addExampleLine(
"xmipp_micrograph_automatic_picking -i micrograph.tif --particleSize 100 --model model --thr 4 --outputRoot micrograph --mode try ");
1802 program->
addExampleLine(
"xmipp_micrograph_automatic_picking -i micrograph.tif --particleSize 100 --model model --thr 4 --outputRoot micrograph --mode train manual.pos");
1803 program->
addExampleLine(
"Automatically select particles after training:",
false);
1804 program->
addExampleLine(
"xmipp_micrograph_automatic_picking -i micrograph.tif --particleSize 100 --model model --thr 4 --outputRoot micrograph --mode autoselect");
1808 addUsageLine(
"Automatic particle picking for micrographs");
1809 addUsageLine(
"+The algorithm is designed to learn the particles from the user, as well as from its own errors.");
1810 addUsageLine(
"+The algorithm is fully described in [[http://www.ncbi.nlm.nih.gov/pubmed/19555764][this paper]].");
1823 autoPicking = std::make_unique<AutoParticlePicking2>(autoPicking->particle_size,autoPicking->filter_num,autoPicking->corr_num,autoPicking->NPCA,
fn_model, std::vector<MDRowSql>());
1824 autoPicking->automaticWithouThread(
fn_micrograph,proc_prec,fnAutoParticles);
double euclidean_distance(const Particle2 &p1, const Particle2 &p2)
void buildSearchSpace(std::vector< Particle2 > &positionArray, bool fast)
void extractPositiveInvariant()
void selfWindow(int n0, int z0, int y0, int x0, int nF, int zF, int yF, int xF, T init_value=0)
FeaturesThread(AutoParticlePicking2 *picker)
void min(Image< double > &op1, const Image< double > &op2)
MultidimArray< double > autoFeatVec
void resize(size_t Ndim, size_t Zdim, size_t Ydim, size_t Xdim, bool copy=true)
void normalize_NewXmipp(MultidimArray< double > &I, const MultidimArray< int > &bg_mask)
void addVector(const MultidimArray< float > &_v)
Add vector.
void LoadModel(const FileName &fnModel)
__host__ __device__ float2 floor(const float2 v)
void image_convertCartesianToPolar_ZoomAtCenter(const MultidimArray< double > &in, MultidimArray< double > &out, Matrix1D< double > &R, double zoomFactor, double Rmin, double Rmax, int NRSteps, float angMin, double angMax, int NAngSteps)
double alignImages(const MultidimArray< double > &Iref, const AlignmentTransforms &IrefTransforms, MultidimArray< double > &I, Matrix2D< double > &M, bool wrap, AlignmentAux &aux, CorrelationAux &aux2, RotationalCorrelationAux &aux3)
void sort(MultidimArray< T > &result) const
#define FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(m)
int readNextMic(FileName &fnmicrograph)
void sqrt(Image< double > &op)
void aliasImageInStack(const MultidimArray< T > &m, size_t select_image)
void extractNegativeInvariant()
void BinaryCircularMask(MultidimArray< int > &mask, double radius, int mode, double x0, double y0, double z0)
MultidimArray< double > dataSet
Image< double > microImage
void SaveModel(const FileName &fnModel)
#define DIRECT_A2D_ELEM(v, i, j)
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)
void subtractAvg()
Subtract average.
void correlation_matrix(const MultidimArray< double > &m1, const MultidimArray< double > &m2, MultidimArray< double > &R, CorrelationAux &aux, bool center)
void filterBankGenerator()
void saveAutoParticles(MetaData &md)
static const int NangSteps
FileName beforeLastOf(const String &str) const
SVMClassifier classifier2
FileName removeAllExtensions() const
Image< double > micrographStack
void correlationBetweenPolarChannels(int n1, int n2, int nF, const MultidimArray< std::complex< double > > &fourierPolarStack, MultidimArray< double > &mIpolarCorr, CorrelationAux &aux)
#define FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY1D(v)
int getParticlesThreshold()
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 convert2PolarFourier(MultidimArray< double > &particleImage, MultidimArray< std::complex< double > > &polar)
Convert an image to its polar form.
void open_micrograph(const FileName &fn_micrograph)
void batchBuildInvariant(const std::vector< MDRowSql > &MD)
void CenterFFT(MultidimArray< T > &v, bool forward)
Image< double > micrographStackPre
AutoParticlePicking2()=default
#define DIRECT_A1D_ELEM(v, i)
void readParams()
Read parameters.
void window(MultidimArray< T1 > &result, int n0, int z0, int y0, int x0, int nF, int zF, int yF, int xF, T1 init_value=0) const
void savePCAModel(const FileName &fn)
Save the PCA basis and average for each channel.
MultidimArray< double > classLabel1
void loadTrainingSet(const FileName &fn_root)
Load training set into the related array.
const char * getParam(const char *param, int arg=0)
void setParameters(double c, double gamma)
MultidimArray< double > dataSet1
void getDimensions(size_t &Xdim, size_t &Ydim, size_t &Zdim, size_t &Ndim) const
bool isLocalMaxima(MultidimArray< double > &inputArray, int x, int y)
MultidimArray< double > positiveParticleStack
int automaticWithouThread(FileName fnmicrograph, int proc_prec, const FileName &fn)
void setValue(const MDObject &object) override
std::vector< Particle2 > rejected_particles
void extractStatics(MultidimArray< double > &inputVec, MultidimArray< double > &features)
void applyConvolution(bool fast)
Convolve the micrograph with the different templates.
MultidimArray< double > pcaRotModel
void setMicrograph(const FileName &fnMic, int proc_prec)
void max(Image< double > &op1, const Image< double > &op2)
std::vector< Particle2 > accepted_particles
void buildInvariant(const std::vector< MDRowSql > &MD)
void addExampleLine(const char *example, bool verbatim=true)
void workOnMicrograph(const FileName &fnMic, int proc_prec)
void selfScaleToSizeFourier(int Zdim, int Ydim, int Xdim, MultidimArray< double > &mda, int nThreads)
std::vector< MultidimArray< double > > PCAbasis
void trainSVM(const FileName &fnModel, int numClassifier)
MultidimArray< std::complex< double > > FFT1
MultidimArray< double > positiveInvariatnStack
void normalizeDataset(int a, int b)
void polarCorrelation(const MultidimArray< std::complex< double > > &fourierPolarStack, MultidimArray< double > &IpolarCorr)
Calculate the correlation of different polar channels.
int automaticallySelectParticles(FileName fnmicrograph, int proc_prec, std::vector< MDRowSql > &md)
std::vector< Particle2 > positionArray
MultidimArray< double > dataSetNormal
double PCAProject(MultidimArray< double > &pcaBasis, MultidimArray< double > &vec)
Project a vector on one pca basis.
Image< double > microImagePrev
void readMic(const FileName &fn_micrograph, int keepPrev)
Read micrograph from the file.
MultidimArray< double > classLabel
void trainRotPCA(const FileName &fnAvgModel, const FileName &fnPCARotModel)
double predict(MultidimArray< double > &featVec, double &score)
int add_coord(int x, int y, int label, double cost)
MultidimArray< double > filterBankStack
static void defineParams(XmippProgram *program)
Define the parameters of the main program.
T dotProduct(const Matrix1D< T > &v1, const Matrix1D< T > &v2)
std::vector< Particle2 > auto_candidates
ProgImageRotationalPCA rotPcaAnalyzer
#define FIRST_XMIPP_INDEX(size)
bool checkDist(Particle2 &p)
MultidimArray< double > pcaModel
void applyMaskFourierSpace(const MultidimArray< double > &v, MultidimArray< std::complex< double > > &V)
String formatString(const char *format,...)
Particle_coords & coord(int n)
void defineParams()
Define Parameters.
MultidimArray< double > negativeParticleStack
double cost
Cost, scaled between 0 and 1.
void learnPCABasis(size_t NPCA, size_t Niter)
Learn basis.
bool checkParam(const char *param)
int read(const FileName &name, DataMode datamode=DATA, size_t select_img=ALL_IMAGES, bool mapData=false, int mode=WRITE_READONLY)
MultidimArray< double > avg
PCAMahalanobisAnalyzer pcaAnalyzer
void normalize_OldXmipp(MultidimArray< double > &I)
MultidimArray< double > particleAvg
double computeStddev() const
void initZeros(const MultidimArray< T1 > &op)
double computeAvg() const
~AutoParticlePicking2()
Destructor.
#define LAST_XMIPP_INDEX(size)
int getIntParam(const char *param, int arg=0)
int read(const FileName &name, DataMode datamode=DATA, size_t select_img=ALL_IMAGES, bool mapData=false)
void SVMTrain(MultidimArray< double > &trainSet, MultidimArray< double > &lable)
void generateMask(MultidimArray< double > &v)
void buildVector(MultidimArray< double > &inputVec, MultidimArray< double > &staticVec, MultidimArray< double > &featureVec, MultidimArray< double > &pieceImage)
void extractParticle(const int x, const int y, MultidimArray< double > &filter, MultidimArray< double > &particleImage, bool normal)
int readMapped(const FileName &name, size_t select_img=ALL_IMAGES, int mode=WRITE_READONLY)
MultidimArray< double > negativeInvariatnStack
MultidimArray< double > vec
void generateFeatVec(const FileName &fnmicrograph, int proc_prec, std::vector< Particle2 > &positionArray)
void addParamsLine(const String &line)
void indexSort(MultidimArray< int > &indx) const
void train(const std::vector< MDRowSql > &MD, bool corrFlag, int x, int y, int width, int height)
void readParams(XmippProgram *program)
Read the parmaeters of the main program.
void applyMaskSpace(MultidimArray< double > &v)
void getImage(size_t n, MultidimArray< T > &M, size_t n2=0) const
std::vector< MDRowSql > micList
void correction(const std::vector< MDRowSql > &addedParticlesMD, const std::vector< MDRowSql > &removedParticlesMD)
MultidimArray< double > convolveRes
void extractNonParticle(std::vector< Particle2 > &negativePosition)