48 im.
resize(ndim,zdim,ydim,xdim,
false);
51 double iMaxDim2 = 2./
std::max(xdim,ydim);
79 A2D_ELEM(im,
i,
j) = std::cos(50*std::exp(-0.5*(std::pow(
i*iMaxDim2*freq,2)+std::pow(
j*iMaxDim2*freq,2))))+
rnd_gaus(0,noiseLevel);
96 A2D_ELEM(im,
i,
j) = (std::exp(-((
i*
i)+(
j*
j))/(2e3*freq*freq)))*(std::cos(50*std::exp(-0.5*(std::pow(i*iMaxDim2*freq,2)+std::pow(j*iMaxDim2*freq,2))))+
rnd_gaus(0,noiseLevel));
128 double i2=((double)
i)*
i;
131 double j2=((double)
j)*
j;
133 if ( (
i!=0) || (
j!=0))
151 size_t NR, NC,NZ, NDim;
154 if ( (NZ!=1) || (NDim!=1) )
172 if ( (
i==0) || (
i== (NR-1)) || (
j == 0) || (
j== (NC-1)) )
192 A2D_ELEM(mask,
i,
j) = double(1)/ double(wSize*wSize+1);
217 A2D_ELEM(orn,
i,
j) = std::fmod((0.5*std::atan2(tempy,tempx))+
double(
PI)/2,
double(
PI));
250 double K=1.0/(2*S*S);
261 temp= std::exp(-std::pow((
std::sqrt(i2_j2)-R),2)*K)*(1-(std::exp(-0.5*i2_j2)));
315 double rang = (rmax-rmin)/2;
317 double freq2 =
XSIZE(im)/(rang/1);
318 double freq1 =
XSIZE(im)/(rang/15);
323 double temp= (1/(1+std::exp((
std::sqrt(r2)-freq1)/10)))*
324 (1-(std::exp(-r2 /(2*freq2*freq2))));
325 A2D_ELEM(H,
i,j) = std::complex<double>(temp,temp);
378 double rang = (rmax-rmin)/2;
380 double freq2 =
XSIZE(im)/(rang/1);
381 double freq1 =
XSIZE(im)/(rang/15);
386 double temp= (1/(1+std::exp((
std::sqrt(r2)-freq1)/10)))*
387 (1-(std::exp(-r2 /(2*freq2*freq2))));
388 A2D_ELEM(H,
i,j) = std::complex<double>(temp,temp);
439 double minQuality = 0;
440 int imax, jmax, max_levels=10;
447 double maxQualityMapValue = (
A2D_ELEM(qualityMap,imax,jmax));
477 std::priority_queue<PointQuality> queueToProcess;
478 queueToProcess.push(p);
480 while (!queueToProcess.empty())
482 p=queueToProcess.top();
485 queueToProcess.pop();
487 size_t indi[8] = {i-1,i-1,i-1,
i,
i,i+1,i+1,i+1};
488 size_t indj[8] = {j-1,
j,j+1,j-1,j+1,j-1,
j,j+1};
490 for(
int k = 0;
k< 8 ;
k++)
496 if ( (!
A2D_ELEM(processed,ni,nj)) &&
497 (
A2D_ELEM(qualityMap,ni,nj) > minQuality) &&
498 ((ni-size)>0) && ((nj-size)>0) &&
499 ((ni+size)<
YSIZE(processed)-1) &&
500 ((nj+size)<
XSIZE(processed)-1) )
502 double G11=0, G12=0, G22=0, b1=0, b2=0;
504 for (
int li=-size; li <= size; li++)
507 for (
int lj=-size; lj <= size; lj++)
512 double tempx =
A2D_ELEM(ds,nli,nlj);
513 double tempy =
A2D_ELEM(dc,nli,nlj);
515 G11 += tempx*tempx+
lambda;
517 G22 += tempy*tempy+
lambda;
519 b1 += lambda*(
A2D_ELEM(px,nli,nlj));
520 b2 += lambda*(
A2D_ELEM(py,nli,nlj));
544 queueToProcess.push(p);
555 double minQuality = 0.05;
557 int imax, jmax, max_levels=10;
564 double maxQualityMapValue =
A2D_ELEM(qualityMap,imax,jmax);
577 A2D_ELEM(processed,imax,jmax) =
true;
583 std::priority_queue<PointQuality> queueToProcess;
584 queueToProcess.push(p);
597 while(!queueToProcess.empty())
599 p=queueToProcess.top();
600 if ((p.
i <=0) || (p.
j <= 0)) {
606 queueToProcess.pop();
608 size_t indi[8] = {i-1,i-1,i-1,
i,
i,i+1,i+1,i+1};
609 size_t indj[8] = {j-1,
j,j+1,j-1,j+1,j-1,
j,j+1};
611 for(
int k = 0;
k< 8 ;
k++)
615 if ( (!
A2D_ELEM(processed,ni,nj)) && (
A2D_ELEM(qualityMap,ni,nj) > minQuality) && ((ni-size)>0) && ((nj-size)>0) && ((ni+size)<
YSIZE(processed)-1)
616 && ((nj+size)<
XSIZE(processed)-1))
620 for (
int li=-size; li <= size; li++)
623 for (
int lj=-size; lj <= size; lj++)
628 uw =
A2D_ELEM(unwrappedPhase,nli,nlj);
630 g =
dMij(gaussian,li+size,lj+size);
635 n = (int)(
std::floor(t+3.14159265)/(2*3.14159265));
637 n = (int)(std::ceil(t-3.14159265)/(2*3.14159265));
639 up = t - (2*3.14159265)*
n;
654 A2D_ELEM(unwrappedPhase,ni,nj) = pred/norm + (lambda*cor)/norm;
664 queueToProcess.push(p);
696 if ( (temp > rmin) && (temp < rmax) )
719 save.
write(
"PPP1.xmp");
740 save.
write(
"PPP2.xmp");
750 auto ci = std::complex<double>(0,1.0);
751 std::complex<double> temp;
767 if ( !(((tempTheta > val) | (tempTheta < -val)) && !((tempTheta >
PI-val) || (tempTheta < -
PI+val)) &&
768 !( !(tempTheta >
PI/2+val) && (tempTheta >
PI/2-val)) && !( !(tempTheta > -
PI/2+val) && (tempTheta > -
PI/2-val))) )
783 save.
write(
"PPP3.xmp");
786 saveROI.
write(
"PPP_roi.xmp");
795 save.
write(
"PPP4.xmp");
809 polynom.
fit(coefsInit,phase,mod2,ROI,0);
810 polynom.
fit(coefsInit,phase,mod2,ROI,verbose);
833 save.
write(
"PPP0.xmp");
835 save.
write(
"PPP1.xmp");
837 save.
write(
"PPP21.xmp");
839 save.
write(
"PPP22.xmp");
841 save.
write(
"PPP3.xmp");
843 save.
write(
"PPP4.xmp");
845 save.
write(
"PPP5.xmp");
847 save.
write(
"PPP6.xmp");
849 save2.
write(
"PPP7.xmp");
885 if ( (temp > rmin) && (temp < rmax) )
908 save.
write(
"PPP1.xmp");
914 save.
write(
"PPP2.xmp");
924 auto ci = std::complex<double>(0,1.0);
925 std::complex<double> temp;
939 tempTheta = std::atan2(
j-(
int)((
float) (
XSIZE(ROI)) / 2.0),
i-(
int)((
float) (
XSIZE(ROI)) / 2.0));
941 if ( !(((tempTheta > val) | (tempTheta < -val)) && !((tempTheta >
PI-val) || (tempTheta < -
PI+val)) &&
942 !( !(tempTheta >
PI/2+val) && (tempTheta >
PI/2-val)) && !( !(tempTheta > -
PI/2+val) && (tempTheta > -
PI/2-val))) )
956 save.
write(
"PPP3.xmp");
964 save.
write(
"PPP4.xmp");
978 polynom.
fit(coefsInit,phase,onesMatrix,ROI,0);
979 polynom.
fit(coefsInit,phase,onesMatrix,ROI,verbose);
1002 save.
write(
"PPP0.xmp");
1004 save.
write(
"PPP1.xmp");
1008 save.
write(
"PPP3.xmp");
1010 save.
write(
"PPP4.xmp");
1012 save.
write(
"PPP5.xmp");
1014 save.
write(
"PPP6.xmp");
1018 save2.
write(
"PPP7.xmp");
1023 double rmin,
double rmax,
int numAngles,
int verbose)
1031 double thrs = (eff0 + effF)*0.5;
1036 while (nAngles < 2*
PI)
1038 double minX = (rmin/2)*std::cos(nAngles);
1039 double minY = (rmin/2)*std::sin(nAngles);
1040 double maxX = (rmax/2)*std::cos(nAngles);
1041 double maxY = (rmax/2)*std::sin(nAngles);
1042 int numPts = numAngles;
1045 double step = (maxX-minX)/numPts;
1048 for (
int i = 0;
i < numPts;
i++)
1050 y = ( ((maxY-minY)/(maxX-minX))*(x-minX)+minY );
1051 z =
A2D_ELEM(enhancedPSD,
int(y),
int(x));
1064 nAngles += 2*
PI/numAngles;
1068 double x0,
y0, majorAxis, minorAxis, ellipseAngle;
1069 fitEllipse(xPoints, yPoints, x0, y0, majorAxis,minorAxis,ellipseAngle);
1075 size_t NR, NC,NZ, Ndim;
1101 A2D_ELEM(mask,
i,
j) = double(1)/ double(wSize*wSize+1);
1127 save() = enhancedPSD;
1128 save.
write(
"PPP0.xmp");
1131 save.
write(
"PPP1.xmp");
1134 save.
write(
"PPP2.xmp");
1140 double & minorAxis,
double & ellipseAngle)
1145 for (
size_t nRow = 0; nRow <
VEC_XSIZE(xPts); nRow++)
1147 double tempX =
VEC_ELEM(xPts,nRow);
1148 double tempY =
VEC_ELEM(yPts,nRow);
1150 dMij(B,nRow,0) = 2*tempX*tempY;
1151 dMij(B,nRow,1) = tempY*tempY-tempX*tempX;
1152 dMij(B,nRow,2) = tempX;
1153 dMij(B,nRow,3) = tempY;
1155 VEC_ELEM(squareXPts,nRow) = -tempX*tempX;
1182 double c_h = temp1.
sum()+temp2.
sum()+
c;
1187 ellipseAngle = std::atan2(
dMij(u,0,0),-
dMij(u,1,0));
1196 double sEllipseAngle, cEllipseAngle;
1197 sincos(-ellipseAngle,&sEllipseAngle,&cEllipseAngle);
1199 double angle = 0, deltaAngle=(2*
PI)/
VEC_XSIZE(xPts);
1200 for (
size_t nPoint = 0; nPoint <
VEC_XSIZE(xPts); nPoint++)
1202 double sAngle, cAngle;
1203 sincos(angle,&sAngle,&cAngle);
1205 double K1=majorAxis*cAngle;
1206 double K2=minorAxis*sAngle;
1207 VEC_ELEM(xPts,nPoint) = K1*cEllipseAngle - K2*sEllipseAngle+
x0;
1208 VEC_ELEM(yPts,nPoint) = K1*sEllipseAngle + K2*cEllipseAngle+
y0;
1209 angle += deltaAngle;
1214 double & minorAxis,
double & ellipseAngle,
size_t & area)
1216 area = (size_t)normImag.
sum();
1231 fitEllipse(xPts,yPts,x0,y0,majorAxis,minorAxis,ellipseAngle);
1236 double Cs,
double imgSize,
double Ts)
1238 double R = majorAxis / (imgSize * Ts);
1239 double a = lambda * R * R;
1240 defocusU = -std::atan(-Q0) / (
PI *
a) - (1 / a) - 0.5 * (Cs * a *
lambda);
1242 R = minorAxis / (imgSize * Ts);
1244 defocusV = -std::atan(-Q0) / (
PI *
a) - (1 / a) - 0.5 * (Cs * a *
lambda);
void fitEllipse(Matrix1D< double > &xPts, Matrix1D< double > &yPts, double &x0, double &y0, double &majorAxis, double &minorAxis, double &ellipseAngle)
#define A2D_ELEM(v, i, j)
void resize(size_t Ndim, size_t Zdim, size_t Ydim, size_t Xdim, bool copy=true)
__host__ __device__ float2 floor(const float2 v)
#define REPORT_ERROR(nerr, ErrormMsg)
#define FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(m)
void resizeNoCopy(const MultidimArray< T1 > &v)
void sqrt(Image< double > &op)
void orMinDer(const MultidimArray< double > &im, MultidimArray< double > &orMap, MultidimArray< double > &orModMap, int wSize, MultidimArray< bool > &ROI)
#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 normalize(FourierTransformer &ftrans, MultidimArray< double > &im, MultidimArray< double > &imN, MultidimArray< double > &imModMap, double R, double S, MultidimArray< bool > &ROI)
void maxIndex(int &jmax) const
void abs(Image< double > &op)
double percentil(double percent_mass)
void zernikePols(const Matrix1D< int > coef, MultidimArray< double > &im, const MultidimArray< bool > &ROI, int verbose=0)
void demodulate(MultidimArray< double > &im, double lambda, int size, int x, int y, int rmin, int rmax, MultidimArray< double > &phase, MultidimArray< double > &mod, Matrix1D< double > &coeffs, int verbose)
void convolutionFFTStack(const MultidimArray< double > &img, const MultidimArray< double > &kernel, MultidimArray< double > &result)
T norm(const std::vector< T > &v)
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 FOR_ALL_ELEMENTS_IN_ARRAY2D(m)
void CenterFFT(MultidimArray< T > &v, bool forward)
void demodulate2(MultidimArray< double > &im, double lambda, int size, int rmin, int rmax, Matrix1D< double > &coeffs, int verbose)
Matrix1D< double > fittedCoeffs
void compute_hist(const MultidimArrayGeneric &array, Histogram1D &hist, int no_steps)
void normalizeWB(MultidimArray< double > &im, MultidimArray< double > &imN, MultidimArray< double > &imModMap, double rmax, double rmin, MultidimArray< bool > &ROI)
void calculateDefocus(double &defocusU, double &defocusV, double majorAxis, double minorAxis, double Q0, double lambda, double Cs, double imgSize, double Ts)
void max(Image< double > &op1, const Image< double > &op2)
void sincos(const MultidimArray< double > &x, MultidimArray< double > &s, MultidimArray< double > &c)
void convolutionFFT(MultidimArray< double > &img, MultidimArray< double > &kernel, MultidimArray< double > &result)
double sum(bool average=false) const
void direction(const MultidimArray< double > &orMap, MultidimArray< double > &qualityMap, double lambda, int size, MultidimArray< double > &dirMap, int x, int y)
void mod(const MultidimArray< T > &x, MultidimArray< T > &m, double y)
void simulPattern(MultidimArray< double > &im, enum FP_TYPE type, int xdim, int ydim, double noiseLevel, const double freq, Matrix1D< int > coefs)
Function to simulate some test fringe patterns.
void svdcmp(const Matrix2D< T > &a, Matrix2D< double > &u, Matrix1D< double > &w, Matrix2D< double > &v)
const bool operator<(const PointQuality &p) const
void unwrapping(const MultidimArray< double > &wrappedPhase, MultidimArray< double > &qualityMap, double lambda, int size, MultidimArray< double > &unwrappedPhase)
void typeCast(const Matrix1D< T1 > &v1, Matrix1D< T2 > &v2)
void normalizeWB2(MultidimArray< double > &im, MultidimArray< double > &imN, MultidimArray< double > &imModMap, double rmax, double rmin, MultidimArray< bool > &ROI)
void SPTH(FourierTransformer &ftrans, MultidimArray< double > &im, MultidimArray< std::complex< double > > &imProc)
void fit(const Matrix1D< int > &coef, MultidimArray< double > &im, const MultidimArray< double > &weight, const MultidimArray< bool > &ROI, int verbose=0)
Incorrect MultidimArray dimensions.
void initZeros(const MultidimArray< T1 > &op)
#define M2x2_INV(Ainv, A)
void initGaussian(int Ydim, int Xdim, double op1=0., double op2=1.)
Incorrect value received.
void firsPSDZero(MultidimArray< double > &enhancedPSD, Matrix1D< double > &xPoints, Matrix1D< double > &yPoints, double rmin, double rmax, int numAngles, int verbose)
void getDimensions(size_t &Xdim, size_t &Ydim, size_t &Zdim, size_t &Ndim) const