Xmipp  v3.23.11-Nereus
Enumerations | Functions
Routines to work with fringes
Collaboration diagram for Routines to work with fringes:

Enumerations

enum  FP_TYPE {
  SIMPLY_OPEN_FRINGES, SIMPLY_CLOSED_FRINGES, COMPLEX_OPEN_FRINGES, COMPLEX_CLOSED_FRINGES,
  SIMPLY_CLOSED_FRINGES_MOD
}
 

Functions

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. More...
 
void SPTH (FourierTransformer &ftrans, MultidimArray< double > &im, MultidimArray< std::complex< double > > &imProc)
 
void orMinDer (const MultidimArray< double > &im, MultidimArray< double > &orMap, MultidimArray< double > &orModMap, int wSize, MultidimArray< bool > &ROI)
 
void normalize (FourierTransformer &ftrans, MultidimArray< double > &im, MultidimArray< double > &imN, MultidimArray< double > &imModMap, double R, double S, MultidimArray< bool > &ROI)
 
void normalizeWB (MultidimArray< double > &im, MultidimArray< double > &imN, MultidimArray< double > &imModMap, double rmax, double rmin, MultidimArray< bool > &ROI)
 
void normalizeWB2 (MultidimArray< double > &im, MultidimArray< double > &imN, MultidimArray< double > &imModMap, double rmax, double rmin, MultidimArray< bool > &ROI)
 
void direction (const MultidimArray< double > &orMap, MultidimArray< double > &qualityMap, double lambda, int size, MultidimArray< double > &dirMap, int x, int y)
 
void unwrapping (const MultidimArray< double > &wrappedPhase, MultidimArray< double > &qualityMap, double lambda, int size, MultidimArray< double > &unwrappedPhase)
 
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=0)
 
void demodulate2 (MultidimArray< double > &im, double lambda, int size, int rmin, int rmax, Matrix1D< double > &coeffs, int verbose=0)
 
void firsPSDZero (MultidimArray< double > &enhancedPSD, Matrix1D< double > &xPoints, Matrix1D< double > &yPoints, double rmin, double rmax, int numAngles, int verbose=0)
 
void fitEllipse (Matrix1D< double > &xPts, Matrix1D< double > &yPts, double &x0, double &y0, double &majorAxis, double &minorAxis, double &ellipseAngle)
 
void fitEllipse (MultidimArray< double > &normImag, double &x0, double &y0, double &majorAxis, double &minorAxis, double &ellipseAngle, size_t &area)
 
void calculateDefocus (double &defocusU, double &defocusV, double majorAxis, double minorAxis, double Q0, double lambda, double Cs, double imgSize, double Ts)
 

Detailed Description

Enumeration Type Documentation

◆ FP_TYPE

enum FP_TYPE

Types of fringes. This class groups all the fringe processing methods that will be used to process the CTF in the non parametric CTF estimation approach

Enumerator
SIMPLY_OPEN_FRINGES 
SIMPLY_CLOSED_FRINGES 
COMPLEX_OPEN_FRINGES 
COMPLEX_CLOSED_FRINGES 
SIMPLY_CLOSED_FRINGES_MOD 

Definition at line 44 of file fringe_processing.h.

Function Documentation

◆ calculateDefocus()

void calculateDefocus ( double &  defocusU,
double &  defocusV,
double  majorAxis,
double  minorAxis,
double  Q0,
double  lambda,
double  Cs,
double  imgSize,
double  Ts 
)

Definition at line 1235 of file fringe_processing.cpp.

1237 {
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);
1241 
1242  R = minorAxis / (imgSize * Ts);
1243  a = lambda * R * R;
1244  defocusV = -std::atan(-Q0) / (PI * a) - (1 / a) - 0.5 * (Cs * a * lambda);
1245 
1246 }
double * lambda
#define PI
Definition: tools.h:43
doublereal * a

◆ demodulate()

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 = 0 
)

Definition at line 671 of file fringe_processing.cpp.

673 {
674  //Initial Setup
675  MultidimArray< double > In, orMap, orModMap, dir, wphase;
678 
679  In.resizeNoCopy(im);
680  orMap.resizeNoCopy(im);
681  orModMap.resizeNoCopy(im);
682  dir.resizeNoCopy(im);
683  wphase.resizeNoCopy(im);
684  sph.resizeNoCopy(im);
685  ROI.resizeNoCopy(im);
686  mod.resizeNoCopy(im);
687  phase.resizeNoCopy(im);
688 
689  //First we define the Region of Interesting to perform all the operations inside these regions
690  ROI.setXmippOrigin();
691  dir.setXmippOrigin();
692 
694  {
695  double temp = std::sqrt(i*i+j*j);
696  if ( (temp > rmin) && (temp < rmax) )
697  {
698  A2D_ELEM(ROI,i,j)= true;
699  A2D_ELEM(dir,i,j)= -std::atan2(i,j);
700  }
701  else
702  A2D_ELEM(ROI,i,j)= false;
703  }
704 
705  STARTINGX(dir)=STARTINGY(dir)=0;
706 
707  //We obtain previous necessary maps
708  //Normalized version of im and modulation map
709  normalizeWB(im,In,mod, rmax, rmin, ROI);
710 
711  int imax, jmax;
712  mod.maxIndex(imax,jmax);
713  mod = mod/A2D_ELEM(mod,imax,jmax);
714 
715  Image<double> save;
716  if (verbose == 1)
717  {
718  save()=In;
719  save.write("PPP1.xmp");
720  }
721 
722  STARTINGX(mod)=STARTINGY(mod)=0;
723  //Orientation map of the fringes
724 
725  /*orMinDer(In, orMap, orModMap, size, ROI);
726  if (verbose == 2)
727 {
728  save()=orMap;
729  save.write("PPP21.xmp");
730  save()=orModMap;
731  save.write("PPP22.xmp");
732 }
733  */
734 
735  //We obtain the direction from the orientation map
736  //direction(orMap, orModMap, lambda, size, dir, x, y);
737  if (verbose == 2)
738  {
739  save()=dir;
740  save.write("PPP2.xmp");
741  }
742  //Spiral transform of the normalized image
743  FourierTransformer ftrans(FFTW_BACKWARD);
744  SPTH(ftrans,In,sph);
745  STARTINGX(sph)=STARTINGY(sph)=0;
746  STARTINGX(In)=STARTINGY(In)=0;
747  STARTINGX(ROI)=STARTINGY(ROI)=0;
748 
749  // We obtain the wrapped phase
750  auto ci = std::complex<double>(0,1.0);
751  std::complex<double> temp;
752 
753  //tempTheta is the Theta coordinate in polar. This is because we want to subtract in ROI the crux of the CTF
754  double tempTheta=0;
755  //val is a number that establish the quantity we want to subtract in the phase to avoid the crux
756  double val = 0.1;
758  {
759  if (A2D_ELEM(ROI,i,j))
760  {
761  temp = -ci*std::exp(ci*A2D_ELEM(dir,i,j))*A2D_ELEM(sph,i,j);
762  A2D_ELEM(wphase,i,j) = std::atan2(temp.real(),A2D_ELEM(In,i,j));
763  A2D_ELEM(orModMap,i,j) = std::sqrt( temp.real()*temp.real() + A2D_ELEM(In,i,j)*A2D_ELEM(In,i,j));
764 
765  //tempTheta = std::atan2(j-(int)((float) (XSIZE(ROI)) / 2.0), i-(int)((float) (XSIZE(ROI)) / 2.0));
766  tempTheta = A2D_ELEM(dir,i,j);
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))) )
769  {
770  A2D_ELEM(ROI,i,j) = false;
771  }
772 
773  }
774  else
775  {
776  A2D_ELEM(orModMap,i,j) = 0;
777  }
778  }
779 
780  if (verbose == 3)
781  {
782  save()=wphase;
783  save.write("PPP3.xmp");
784  Image<bool> saveROI;
785  saveROI() = ROI;
786  saveROI.write("PPP_roi.xmp");
787  }
788 
789  //unwrapping(wphase, orModMap, lambda, size, phase);
790  unwrapping(wphase, mod, lambda, size, phase);
791 
792  if (verbose == 4)
793  {
794  save()=phase;
795  save.write("PPP4.xmp");
796  }
797 
798  ROI.setXmippOrigin();
799  mod.setXmippOrigin();
800  Matrix1D<int> coefsInit(VEC_XSIZE(coeffs));
801 
802  for (size_t i=0; i<VEC_XSIZE(coeffs); i++)
803  VEC_ELEM(coefsInit,i) = (int)VEC_ELEM(coeffs,i);
804 
805  PolyZernikes polynom;
807  mod2.resizeNoCopy(im);
808  mod2.initConstant(1.0);
809  polynom.fit(coefsInit,phase,mod2,ROI,0);
810  polynom.fit(coefsInit,phase,mod2,ROI,verbose);
811 
812  int index = 0;
813  for (size_t i=0; i<VEC_XSIZE(coeffs); i++)
814  {
815  if ( VEC_ELEM(coeffs,i) != 0)
816  {
817  VEC_ELEM(coeffs,i) = VEC_ELEM(polynom.fittedCoeffs,index);
818  index++;
819  }
820 
821  else
822  VEC_ELEM(coeffs,i) = 0;
823  }
824 
825  STARTINGX(im)=STARTINGY(im)=0;
826  STARTINGX(phase)=STARTINGY(phase)=0;
827  STARTINGX(mod)=STARTINGY(mod)=0;
828 
829  Image<bool> save2;
830  if (verbose > 5)
831  {
832  save()=im;
833  save.write("PPP0.xmp");
834  save()=In;
835  save.write("PPP1.xmp");
836  save()=orMap;
837  save.write("PPP21.xmp");
838  save()=orModMap;
839  save.write("PPP22.xmp");
840  save()=dir;
841  save.write("PPP3.xmp");
842  save()=wphase;
843  save.write("PPP4.xmp");
844  save()=phase;
845  save.write("PPP5.xmp");
846  save()=mod;
847  save.write("PPP6.xmp");
848  save2()= ROI;
849  save2.write("PPP7.xmp");
850  }
851 }
#define A2D_ELEM(v, i, j)
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
#define FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(m)
#define VEC_XSIZE(m)
Definition: matrix1d.h:77
void resizeNoCopy(const MultidimArray< T1 > &v)
void sqrt(Image< double > &op)
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 maxIndex(int &jmax) const
void initConstant(T val)
#define STARTINGX(v)
#define i
#define FOR_ALL_ELEMENTS_IN_ARRAY2D(m)
#define STARTINGY(v)
Matrix1D< double > fittedCoeffs
double * lambda
void normalizeWB(MultidimArray< double > &im, MultidimArray< double > &imN, MultidimArray< double > &imModMap, double rmax, double rmin, MultidimArray< bool > &ROI)
viol index
void mod(const MultidimArray< T > &x, MultidimArray< T > &m, double y)
#define j
void unwrapping(const MultidimArray< double > &wrappedPhase, MultidimArray< double > &qualityMap, double lambda, int size, MultidimArray< double > &unwrappedPhase)
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)
#define PI
Definition: tools.h:43

◆ demodulate2()

void demodulate2 ( MultidimArray< double > &  im,
double  lambda,
int  size,
int  rmin,
int  rmax,
Matrix1D< double > &  coeffs,
int  verbose = 0 
)

Definition at line 853 of file fringe_processing.cpp.

855 {
856  //Initial Setup :
857  // In : Normalized image after using methods normalizeWB or normalize
858  // orModMap : modulation map of im using
859  //TODO : subtract orModMap
860  // dir : direction map. It can be obtained using direction method but in our case, we impose the direction map shape
861  // wphase : wrapped phase
862  // phase : absolute phase computed from wphase and method unwrapping
863  // mod : modulation map
864 
865  MultidimArray< double > In, mod, dir, wphase, phase;
868 
869  In.resizeNoCopy(im);
870  mod.resizeNoCopy(im);
871  dir.resizeNoCopy(im);
872  wphase.resizeNoCopy(im);
873  sph.resizeNoCopy(im);
874  ROI.resizeNoCopy(im);
875  phase.resizeNoCopy(im);
876 
877  //First we define the Region of Interesting to perform all the operations inside these regions
878  ROI.setXmippOrigin();
879  dir.setXmippOrigin();
880 
881  //We obtain the ROI from rmax and rmin parameters and the direction map
883  {
884  double temp = std::sqrt(i*i+j*j);
885  if ( (temp > rmin) && (temp < rmax) )
886  {
887  A2D_ELEM(ROI,i,j)= true;
888  A2D_ELEM(dir,i,j)= -std::atan2(i,j);
889  }
890  else
891  A2D_ELEM(ROI,i,j)= false;
892  }
893  STARTINGX(dir)=STARTINGY(dir)=0;
894 
895  //We obtain the normalized version of im and the modulation map mod
896  normalizeWB(im,In,mod, rmax, rmin, ROI);
897  STARTINGX(mod)=STARTINGY(mod)=0;
898 
899  //Here we make the max of mod to be 1
900  int imax, jmax;
901  mod.maxIndex(imax,jmax);
902  mod = mod/A2D_ELEM(mod,imax,jmax);
903 
904  Image<double> save;
905  if (verbose == 1)
906  {
907  save()=In;
908  save.write("PPP1.xmp");
909  }
910 
911  if (verbose == 2)
912  {
913  save()=dir;
914  save.write("PPP2.xmp");
915  }
916  //Spiral transform of the normalized image
917  FourierTransformer ftrans(FFTW_BACKWARD);
918  SPTH(ftrans,In,sph);
919  STARTINGX(sph)=STARTINGY(sph)=0;
920  STARTINGX(In)=STARTINGY(In)=0;
921  STARTINGX(ROI)=STARTINGY(ROI)=0;
922 
923  // We obtain the wrapped phase
924  auto ci = std::complex<double>(0,1.0);
925  std::complex<double> temp;
926 
927  //tempTheta is the Theta magnitude in polar coordinates. This is because we want to subtract in ROI the crux of the CTF
928  double tempTheta=0;
929  //val is a number that establish the quantity we want to subtract in the phase to avoid the crux
930  double val = 0.05;
932  {
933  if (A2D_ELEM(ROI,i,j))
934  {
935  temp = -ci*std::exp(ci*A2D_ELEM(dir,i,j))*A2D_ELEM(sph,i,j);
936  A2D_ELEM(wphase,i,j) = std::atan2(temp.real(),A2D_ELEM(In,i,j));
937  //A2D_ELEM(mod,i,j) = std::sqrt( temp.real()*temp.real() + A2D_ELEM(In,i,j)*A2D_ELEM(In,i,j));
938 
939  tempTheta = std::atan2(j-(int)((float) (XSIZE(ROI)) / 2.0), i-(int)((float) (XSIZE(ROI)) / 2.0));
940 
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))) )
943  {
944  A2D_ELEM(ROI,i,j) = false;
945  }
946  }
947  else
948  {
949  A2D_ELEM(mod,i,j) = 0;
950  }
951  }
952 
953  if (verbose == 3)
954  {
955  save()=wphase;
956  save.write("PPP3.xmp");
957  }
958 
959  unwrapping(wphase, mod, lambda, size, phase);
960 
961  if (verbose == 4)
962  {
963  save()=phase;
964  save.write("PPP4.xmp");
965  }
966 
967  ROI.setXmippOrigin();
968  mod.setXmippOrigin();
969  Matrix1D<int> coefsInit(VEC_XSIZE(coeffs));
970 
971  for (size_t i=0; i<VEC_XSIZE(coeffs); i++)
972  VEC_ELEM(coefsInit,i) = (int)VEC_ELEM(coeffs,i);
973 
974  PolyZernikes polynom;
975  MultidimArray< double > onesMatrix;
976  onesMatrix.resizeNoCopy(im);
977  onesMatrix.initConstant(1.0);
978  polynom.fit(coefsInit,phase,onesMatrix,ROI,0);
979  polynom.fit(coefsInit,phase,onesMatrix,ROI,verbose);
980 
981  int index = 0;
982  for (size_t i=0; i<VEC_XSIZE(coeffs); i++)
983  {
984  if ( VEC_ELEM(coeffs,i) != 0)
985  {
986  VEC_ELEM(coeffs,i) = VEC_ELEM(polynom.fittedCoeffs,index);
987  index++;
988  }
989 
990  else
991  VEC_ELEM(coeffs,i) = 0;
992  }
993 
994 
995  STARTINGX(im)=STARTINGY(im)=0;
996  STARTINGX(phase)=STARTINGY(phase)=0;
997  STARTINGX(mod)=STARTINGY(mod)=0;
998 
999  if (verbose > 5)
1000  {
1001  save()=im;
1002  save.write("PPP0.xmp");
1003  save()=In;
1004  save.write("PPP1.xmp");
1005  //save()=orModMap;
1006  //save.write("PPP2.xmp");
1007  save()=dir;
1008  save.write("PPP3.xmp");
1009  save()=wphase;
1010  save.write("PPP4.xmp");
1011  save()=phase;
1012  save.write("PPP5.xmp");
1013  save()=mod;
1014  save.write("PPP6.xmp");
1015 
1016  Image<bool> save2;
1017  save2()= ROI;
1018  save2.write("PPP7.xmp");
1019  }
1020 }
#define A2D_ELEM(v, i, j)
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
#define FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(m)
#define VEC_XSIZE(m)
Definition: matrix1d.h:77
void resizeNoCopy(const MultidimArray< T1 > &v)
void sqrt(Image< double > &op)
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 maxIndex(int &jmax) const
void initConstant(T val)
#define STARTINGX(v)
#define i
#define FOR_ALL_ELEMENTS_IN_ARRAY2D(m)
#define STARTINGY(v)
Matrix1D< double > fittedCoeffs
double * lambda
void normalizeWB(MultidimArray< double > &im, MultidimArray< double > &imN, MultidimArray< double > &imModMap, double rmax, double rmin, MultidimArray< bool > &ROI)
viol index
#define XSIZE(v)
void mod(const MultidimArray< T > &x, MultidimArray< T > &m, double y)
#define j
void unwrapping(const MultidimArray< double > &wrappedPhase, MultidimArray< double > &qualityMap, double lambda, int size, MultidimArray< double > &unwrappedPhase)
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)
#define PI
Definition: tools.h:43

◆ direction()

void direction ( const MultidimArray< double > &  orMap,
MultidimArray< double > &  qualityMap,
double  lambda,
int  size,
MultidimArray< double > &  dirMap,
int  x,
int  y 
)

Definition at line 436 of file fringe_processing.cpp.

437 {
438  //First we perform some setup stuff
439  double minQuality = 0;
440  int imax, jmax, max_levels=10;
441 
442 
443  //We look for the maximun value of qualityMap
444  qualityMap.selfABS();
445  qualityMap.maxIndex(imax,jmax);
446 
447  double maxQualityMapValue = (A2D_ELEM(qualityMap,imax,jmax));
448 
449  MultidimArray<double> qualityMapInt = ((qualityMap/maxQualityMapValue)*(max_levels-1));
450  qualityMapInt.selfROUND();
451 
452  MultidimArray<double> ds, dc, px, py;
453  sincos(orMap,ds,dc);
454  px.resizeNoCopy(orMap);
455  py.resizeNoCopy(orMap);
456 
457  //In hx and hy we store the pixels to process attending to their quality value
458  Matrix2D<double> G(2,2), invG(2,2);
459  Matrix1D<double> b(2),R(2);
460 
461  MultidimArray<bool> processed;
462  processed.initZeros(orMap);
463 
464  // Process the first point
465  size_t i=x;
466  size_t j=y;
467 
468  A2D_ELEM(dirMap,i,j) = std::atan2(A2D_ELEM(ds,i,j),A2D_ELEM(dc,i,j));
469  A2D_ELEM(processed,i,j) = true;
470  A2D_ELEM(px,i,j) = -A2D_ELEM(ds,i,j);
471  A2D_ELEM(py,i,j) = A2D_ELEM(dc,i,j);
472 
473  PointQuality p;
474  p.i=i;
475  p.j=j;
476  p.quality=A2D_ELEM(qualityMapInt,i,j);
477  std::priority_queue<PointQuality> queueToProcess;
478  queueToProcess.push(p);
479 
480  while (!queueToProcess.empty())
481  {
482  p=queueToProcess.top();
483  i=p.i;
484  j=p.j;
485  queueToProcess.pop();
486 
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};
489 
490  for(int k = 0; k< 8 ; k++)
491  {
492 
493  size_t ni=indi[k];
494  size_t nj=indj[k];
495 
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) )
501  {
502  double G11=0, G12=0, G22=0, b1=0, b2=0;
503 
504  for (int li=-size; li <= size; li++)
505  {
506  size_t nli=ni+li;
507  for (int lj=-size; lj <= size; lj++)
508  {
509  size_t nlj=nj+lj;
510  if (A2D_ELEM(processed,nli,nlj))
511  {
512  double tempx = A2D_ELEM(ds,nli,nlj);
513  double tempy = A2D_ELEM(dc,nli,nlj);
514 
515  G11 += tempx*tempx+lambda;
516  G12 += tempx*tempy;
517  G22 += tempy*tempy+lambda;
518 
519  b1 += lambda*(A2D_ELEM(px,nli,nlj));
520  b2 += lambda*(A2D_ELEM(py,nli,nlj));
521  }
522  }
523  }
524 
525  dMij(G,0,0) = G11;
526  dMij(G,0,1) = G12;
527  dMij(G,1,0) = G12;
528  dMij(G,1,1) = G22;
529 
530  VEC_ELEM(b,0) = b1;
531  VEC_ELEM(b,1) = b2;
532 
533  G.inv(invG);
534  R = invG*b;
535 
536  // Process this point
537  A2D_ELEM(dirMap,ni,nj) = std::atan2(-VEC_ELEM(R,1),VEC_ELEM(R,0));
538  A2D_ELEM(processed,ni,nj) = true;
539  A2D_ELEM(px,ni,nj) = VEC_ELEM(R,0);
540  A2D_ELEM(py,ni,nj) = VEC_ELEM(R,1);
541  p.i=ni;
542  p.j=nj;
543  p.quality=A2D_ELEM(qualityMapInt,ni,nj);
544  queueToProcess.push(p);
545  }
546 
547  }
548  }
549 }
#define YSIZE(v)
#define A2D_ELEM(v, i, j)
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
void resizeNoCopy(const MultidimArray< T1 > &v)
static double * y
void maxIndex(int &jmax) const
doublereal * x
#define i
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
doublereal * b
double * lambda
#define dMij(m, i, j)
Definition: matrix2d.h:139
#define XSIZE(v)
void sincos(const MultidimArray< double > &x, MultidimArray< double > &s, MultidimArray< double > &c)
int ni
#define j
void initZeros(const MultidimArray< T1 > &op)

◆ firsPSDZero()

void firsPSDZero ( MultidimArray< double > &  enhancedPSD,
Matrix1D< double > &  xPoints,
Matrix1D< double > &  yPoints,
double  rmin,
double  rmax,
int  numAngles,
int  verbose = 0 
)

Definition at line 1022 of file fringe_processing.cpp.

1024 {
1025  enhancedPSD.setXmippOrigin();
1026  Histogram1D hist;
1027  compute_hist(enhancedPSD, hist, 200);
1028 
1029  double eff0 = hist.percentil(0.1);
1030  double effF = hist.percentil(98);
1031  double thrs = (eff0 + effF)*0.5;
1032 
1033  int index = 0;
1034  double nAngles = 0;
1035 
1036  while (nAngles < 2*PI)
1037  {
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;
1043  double x = minX;
1044  double y = 0;
1045  double step = (maxX-minX)/numPts;
1046  double z;
1047 
1048  for (int i = 0; i < numPts; i++)
1049  {
1050  y = ( ((maxY-minY)/(maxX-minX))*(x-minX)+minY );
1051  z = A2D_ELEM(enhancedPSD,int(y),int(x));
1052 
1053  if (z < thrs)
1054  {
1055  VEC_ELEM(xPoints,index)=x;
1056  VEC_ELEM(yPoints,index)=y;
1057  continue;
1058  }
1059 
1060  x = (x+step);
1061  }
1062 
1063  index++;
1064  nAngles += 2*PI/numAngles;
1065  }
1066 
1067 
1068  double x0, y0, majorAxis, minorAxis, ellipseAngle;
1069  fitEllipse(xPoints, yPoints, x0, y0, majorAxis,minorAxis,ellipseAngle);
1070 
1071  //EXTERNAL FORCE FIELD CALCULATION
1072  //We perform the gradient of enhancedPSD to calculate the external force field (balloon force)
1073  //Calculation of the external force field normalized to maximum magnitude 1
1074 
1075  size_t NR, NC,NZ, Ndim;
1076  enhancedPSD.getDimensions(NR,NC,NZ,Ndim);
1077  MultidimArray<double > fx, fy;
1078  fx.resizeNoCopy(enhancedPSD);
1079  fy.resizeNoCopy(enhancedPSD);
1080 
1081  FOR_ALL_ELEMENTS_IN_ARRAY2D(enhancedPSD)
1082  {
1083  if ( (i==STARTINGX(enhancedPSD)) || (i== (FINISHINGX(enhancedPSD)-1)) || (j == STARTINGY(enhancedPSD)) || (j== (FINISHINGY(enhancedPSD)-1)) )
1084  {
1085  A2D_ELEM(fx,i,j) = 0;
1086  A2D_ELEM(fy,i,j) = 0;
1087  }
1088  else
1089  {
1090  A2D_ELEM(fx,i,j) = std::sqrt(2)*std::abs(A2D_ELEM(enhancedPSD,i+1,j)-A2D_ELEM(enhancedPSD,i-1,j));
1091  A2D_ELEM(fy,i,j) = std::sqrt(2)*std::abs(A2D_ELEM(enhancedPSD,i,j+1)-A2D_ELEM(enhancedPSD,i,j-1));
1092  }
1093  }
1094 
1095  //We calculate accumulative derivatives in the area 2*wSize+1
1096  int wSize = 5;
1097  MultidimArray<double> mask;
1098  mask.resizeNoCopy(enhancedPSD);
1099  for (int i = -std::abs(wSize); i <= std::abs(wSize); i++ )
1100  for (int j = -std::abs(wSize); j <= std::abs(wSize); j++ )
1101  A2D_ELEM(mask,i,j) = double(1)/ double(wSize*wSize+1);
1102 
1103  convolutionFFTStack(fx,mask,fx);
1104  convolutionFFTStack(fy,mask,fy);
1105 
1106  //From the derivatives we calculate kappa:
1107  MultidimArray<double> absFx;
1108  absFx.resizeNoCopy(fx);
1109  //absFx= fx.selfABS();
1110  /*absFx.maxIndex(imax,jmax);
1111  double maxFx = A2D_ELEM(fx,imax,jmax);
1112 
1113  MultidimArray<double> absFy = (fy.selfABS());
1114  absFy.maxIndex(imax,jmax);
1115  double maxFy = A2D_ELEM(fy,imax,jmax);
1116  */
1118  //
1119 
1120  STARTINGX(enhancedPSD)=STARTINGY(enhancedPSD)=0;
1121  //enhancedPSD
1122 
1123  if (verbose != 0)
1124  {
1125  Image<double> save;
1126 
1127  save() = enhancedPSD;
1128  save.write("PPP0.xmp");
1129 
1130  save() = fx;
1131  save.write("PPP1.xmp");
1132 
1133  save() = fy;
1134  save.write("PPP2.xmp");
1135  }
1136 }
void fitEllipse(Matrix1D< double > &xPts, Matrix1D< double > &yPts, double &x0, double &y0, double &majorAxis, double &minorAxis, double &ellipseAngle)
#define A2D_ELEM(v, i, j)
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
#define FINISHINGX(v)
void resizeNoCopy(const MultidimArray< T1 > &v)
void sqrt(Image< double > &op)
static double * y
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 abs(Image< double > &op)
double percentil(double percent_mass)
Definition: histogram.cpp:160
void convolutionFFTStack(const MultidimArray< double > &img, const MultidimArray< double > &kernel, MultidimArray< double > &result)
Definition: xmipp_fftw.cpp:429
#define STARTINGX(v)
doublereal * x
#define i
#define FOR_ALL_ELEMENTS_IN_ARRAY2D(m)
#define STARTINGY(v)
#define y0
void compute_hist(const MultidimArrayGeneric &array, Histogram1D &hist, int no_steps)
Definition: histogram.cpp:572
#define x0
viol index
double z
#define j
#define FINISHINGY(v)
#define PI
Definition: tools.h:43
void getDimensions(size_t &Xdim, size_t &Ydim, size_t &Zdim, size_t &Ndim) const

◆ fitEllipse() [1/2]

void fitEllipse ( Matrix1D< double > &  xPts,
Matrix1D< double > &  yPts,
double &  x0,
double &  y0,
double &  majorAxis,
double &  minorAxis,
double &  ellipseAngle 
)

Definition at line 1139 of file fringe_processing.cpp.

1141 {
1142  Matrix2D<double> B(VEC_XSIZE(xPts),5), A(2,2);
1143  Matrix1D<double> V(5), squareXPts(VEC_XSIZE(xPts));
1144 
1145  for (size_t nRow = 0; nRow < VEC_XSIZE(xPts); nRow++)
1146  {
1147  double tempX = VEC_ELEM(xPts,nRow);
1148  double tempY = VEC_ELEM(yPts,nRow);
1149 
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;
1154  dMij(B,nRow,4) = 1;
1155  VEC_ELEM(squareXPts,nRow) = -tempX*tempX;
1156 
1157  }
1158 
1159  Matrix1D<double> bv(2);
1160  Matrix2D<double> C =(B.inv());
1161  V = C*squareXPts;
1162 
1163  dMij(A,0,0) = 1-VEC_ELEM(V,1);
1164  dMij(A,0,1) = VEC_ELEM(V,0);
1165  dMij(A,1,0) = VEC_ELEM(V,0);
1166  dMij(A,1,1) = VEC_ELEM(V,1);
1167  VEC_ELEM(bv,0) = VEC_ELEM(V,2);
1168  VEC_ELEM(bv,1) = VEC_ELEM(V,3);
1169 
1170  double c = VEC_ELEM(V,4);
1171 
1172  Matrix2D<double> Ainv(2,2), u, v;
1175  M2x2_INV(Ainv,A);
1176  svdcmp(A,u,w,v);
1177 
1178  Matrix1D<double> t = -0.5*Ainv*bv;
1179 
1180  Matrix1D<double> temp1 = (t*(A * t));//c_h = t.transpose() * A * t + bv.transpose() * t + c;
1181  Matrix1D<double> temp2 = bv.transpose() * t;
1182  double c_h = temp1.sum()+temp2.sum()+c;
1183  x0 = VEC_ELEM(t,0);
1184  y0 = VEC_ELEM(t,1);
1185  majorAxis = std::sqrt(-c_h / VEC_ELEM(w,0));
1186  minorAxis = std::sqrt(-c_h / VEC_ELEM(w,1));
1187  ellipseAngle = std::atan2(dMij(u,0,0),-dMij(u,1,0));
1188 
1189  if ( (VEC_ELEM(w,0)==0) || (VEC_ELEM(w,1)==0) || (VEC_ELEM(w,0)*c_h>0) || (VEC_ELEM(w,1)*c_h>0) )
1190  {
1191  majorAxis=0;
1192  minorAxis=0;
1193  }
1194 
1195 
1196  double sEllipseAngle, cEllipseAngle;
1197  sincos(-ellipseAngle,&sEllipseAngle,&cEllipseAngle);
1198 
1199  double angle = 0, deltaAngle=(2*PI)/VEC_XSIZE(xPts);
1200  for (size_t nPoint = 0; nPoint < VEC_XSIZE(xPts); nPoint++)
1201  {
1202  double sAngle, cAngle;
1203  sincos(angle,&sAngle,&cAngle);
1204  //We impose that the origin always is zero and because this whe do not sum it
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;
1210  }
1211 }
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
#define SPEED_UP_temps0
Definition: xmipp_macros.h:394
#define VEC_XSIZE(m)
Definition: matrix1d.h:77
doublereal * c
void sqrt(Image< double > &op)
doublereal * w
double temp2
#define y0
#define x0
#define dMij(m, i, j)
Definition: matrix2d.h:139
void sincos(const MultidimArray< double > &x, MultidimArray< double > &s, MultidimArray< double > &c)
double sum(bool average=false) const
Definition: matrix1d.cpp:652
void svdcmp(const Matrix2D< T > &a, Matrix2D< double > &u, Matrix1D< double > &w, Matrix2D< double > &v)
Definition: matrix2d.cpp:125
doublereal * u
#define M2x2_INV(Ainv, A)
Definition: matrix2d.h:286
#define PI
Definition: tools.h:43
double temp1

◆ fitEllipse() [2/2]

void fitEllipse ( MultidimArray< double > &  normImag,
double &  x0,
double &  y0,
double &  majorAxis,
double &  minorAxis,
double &  ellipseAngle,
size_t &  area 
)

Definition at line 1213 of file fringe_processing.cpp.

1215 {
1216  area = (size_t)normImag.sum();
1217  Matrix1D<double> xPts(area);
1218  Matrix1D<double> yPts(area);
1219  int index = 0;
1220 
1221  FOR_ALL_ELEMENTS_IN_ARRAY2D(normImag)
1222  {
1223  if ( A2D_ELEM(normImag,i,j)==1 )
1224  {
1225  VEC_ELEM(xPts,index)=j;
1226  VEC_ELEM(yPts,index)=i;
1227  index++;
1228  }
1229  }
1230 
1231  fitEllipse(xPts,yPts,x0,y0,majorAxis,minorAxis,ellipseAngle);
1232 }
void fitEllipse(Matrix1D< double > &xPts, Matrix1D< double > &yPts, double &x0, double &y0, double &majorAxis, double &minorAxis, double &ellipseAngle)
#define A2D_ELEM(v, i, j)
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
#define i
#define FOR_ALL_ELEMENTS_IN_ARRAY2D(m)
#define y0
#define x0
viol index
#define j
double sum() const

◆ normalize()

void normalize ( FourierTransformer ftrans,
MultidimArray< double > &  im,
MultidimArray< double > &  imN,
MultidimArray< double > &  imModMap,
double  R,
double  S,
MultidimArray< bool > &  ROI 
)

Definition at line 233 of file fringe_processing.cpp.

235 {
236  // H is an Annular filter with radius=R and sigma=S and a Gaussian DC filter with sigma=1
238  H.resizeNoCopy(im);
239 
240  im.setXmippOrigin();
241  H.setXmippOrigin();
242 
243  MultidimArray<std::complex<double> > fftIm, imComplex;
244  typeCast(im, imComplex);
245 
246  //Fourier Transformer
247  ftrans.FourierTransform(imComplex, fftIm, false);
248 
249  double temp = 0;
250  double K=1.0/(2*S*S);
251  for (int i=STARTINGY(im); i<=FINISHINGY(im); ++i)
252  {
253  double di=i;
254  double i2=di*di;
255  for (int j=STARTINGX(im); j<=FINISHINGX(im); ++j)
256  {
257  double dj=j;
258  double j2=dj*dj;
259  double i2_j2=i2+j2;
260  // temp= std::exp(-std::pow((std::sqrt(std::pow((double)i,2)+std::pow((double)j,2))-R),2)/(2*std::pow(S,2)))*(1-(std::exp((-1)*(std::pow(double(i),2) + std::pow(double(j),2)) /(2*1))));
261  temp= std::exp(-std::pow((std::sqrt(i2_j2)-R),2)*K)*(1-(std::exp(-0.5*i2_j2)));
262  auto *ptr=(double*)&A2D_ELEM(H,i,j);
263  *ptr=*(ptr+1)=temp;
264  }
265  }
266 
267  CenterFFT(H,false);
268  fftIm *= H;
269  ftrans.inverseFourierTransform();
270 
271  //output of the program
272  imN.setXmippOrigin();
273  imModMap.setXmippOrigin();
274 
276  A2D_ELEM(imN,i,j) = A2D_ELEM(imComplex,i,j).real();
277 
278  SPTH(ftrans,imN,H);
279 
281  {
282  if (A2D_ELEM(ROI,i,j))
283  {
284  temp = std::abs(A2D_ELEM(H,i,j));
285  A2D_ELEM(imModMap,i,j) = std::sqrt(std::pow(temp,2)+std::pow(A2D_ELEM(imN,i,j),2));
286  A2D_ELEM(imN,i,j) = std::cos(std::atan2(temp, A2D_ELEM(imN,i,j)));
287  }
288  else
289  {
290  A2D_ELEM(imModMap,i,j) = 0;
291  A2D_ELEM(imN,i,j) = 0;
292  }
293  }
294 
295  STARTINGX(imN)=STARTINGY(imN)=0;
296 }
#define A2D_ELEM(v, i, j)
void inverseFourierTransform()
Definition: xmipp_fftw.cpp:329
#define FINISHINGX(v)
void resizeNoCopy(const MultidimArray< T1 > &v)
void sqrt(Image< double > &op)
void abs(Image< double > &op)
double * di
#define STARTINGX(v)
#define i
#define FOR_ALL_ELEMENTS_IN_ARRAY2D(m)
#define STARTINGY(v)
void CenterFFT(MultidimArray< T > &v, bool forward)
Definition: xmipp_fft.h:291
void FourierTransform(T &v, T1 &V, bool getCopy=true)
Definition: xmipp_fftw.h:166
#define j
#define FINISHINGY(v)
void typeCast(const Matrix1D< T1 > &v1, Matrix1D< T2 > &v2)
Definition: matrix1d.h:1227
void SPTH(FourierTransformer &ftrans, MultidimArray< double > &im, MultidimArray< std::complex< double > > &imProc)
constexpr int K

◆ normalizeWB()

void normalizeWB ( MultidimArray< double > &  im,
MultidimArray< double > &  imN,
MultidimArray< double > &  imModMap,
double  rmax,
double  rmin,
MultidimArray< bool > &  ROI 
)

Definition at line 298 of file fringe_processing.cpp.

299 {
300 
301  // H is an Annular bandpass filter.
303  H.resizeNoCopy(im);
304 
305  im.setXmippOrigin();
306  H.setXmippOrigin();
307 
308  MultidimArray<std::complex<double> > fftIm, imComplex;
309  typeCast(im, imComplex);
310 
311  //Fourier Transformer
312  FourierTransformer ftrans(FFTW_BACKWARD);
313  ftrans.FourierTransform(imComplex, fftIm, false);
314 
315  double rang = (rmax-rmin)/2;
316  //Inside rang we assume that there will a range of fringes per field from 2 to 10.
317  double freq2 = XSIZE(im)/(rang/1);
318  double freq1 = XSIZE(im)/(rang/15);
319 
321  {
322  double r2=i*i+j*j;
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);
326  }
327 
328  CenterFFT(H,false);
329  fftIm *= H;
330  ftrans.inverseFourierTransform();
331 
332  //output of the program
333  imN.setXmippOrigin();
334  imModMap.setXmippOrigin();
335 
337  {
338  A2D_ELEM(imN,i,j) = A2D_ELEM(imComplex,i,j).real();
339  }
340 
341  SPTH(ftrans,imN,H);
342 
344  {
345  if (A2D_ELEM(ROI,i,j))
346  {
347  double temp = std::abs(A2D_ELEM(H,i,j));
348  A2D_ELEM(imModMap,i,j) = std::sqrt(std::pow(temp,2)+std::pow(A2D_ELEM(imN,i,j),2));
349  A2D_ELEM(imN,i,j) = std::cos(std::atan2(temp, A2D_ELEM(imN,i,j)));
350  }
351  else
352  {
353  A2D_ELEM(imModMap,i,j) = 0;
354  A2D_ELEM(imN,i,j) = 0;
355  }
356  }
357 
358  STARTINGX(imN)=STARTINGY(imN)=0;
359 }
#define A2D_ELEM(v, i, j)
void resizeNoCopy(const MultidimArray< T1 > &v)
void sqrt(Image< double > &op)
void abs(Image< double > &op)
#define STARTINGX(v)
#define i
#define FOR_ALL_ELEMENTS_IN_ARRAY2D(m)
#define STARTINGY(v)
void CenterFFT(MultidimArray< T > &v, bool forward)
Definition: xmipp_fft.h:291
#define XSIZE(v)
#define j
void typeCast(const Matrix1D< T1 > &v1, Matrix1D< T2 > &v2)
Definition: matrix1d.h:1227
void SPTH(FourierTransformer &ftrans, MultidimArray< double > &im, MultidimArray< std::complex< double > > &imProc)
float r2

◆ normalizeWB2()

void normalizeWB2 ( MultidimArray< double > &  im,
MultidimArray< double > &  imN,
MultidimArray< double > &  imModMap,
double  rmax,
double  rmin,
MultidimArray< bool > &  ROI 
)

Definition at line 361 of file fringe_processing.cpp.

362 {
363 
364  // H is an Annular bandpass filter.
366  H.resizeNoCopy(im);
367 
368  im.setXmippOrigin();
369  H.setXmippOrigin();
370 
371  MultidimArray<std::complex<double> > fftIm, imComplex;
372  typeCast(im, imComplex);
373 
374  //Fourier Transformer
375  FourierTransformer ftrans(FFTW_BACKWARD);
376  ftrans.FourierTransform(imComplex, fftIm, false);
377 
378  double rang = (rmax-rmin)/2;
379  //Inside rang we assume that there will a range of fringes per field from 2 to 10.
380  double freq2 = XSIZE(im)/(rang/1);
381  double freq1 = XSIZE(im)/(rang/15);
382 
384  {
385  double r2=i*i+j*j;
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);
389  }
390 
391  CenterFFT(H,false);
392  fftIm *= H;
393  ftrans.inverseFourierTransform();
394 
395  //output of the program
396  imN.setXmippOrigin();
397  imModMap.setXmippOrigin();
398 
400  {
401  A2D_ELEM(imN,i,j) = A2D_ELEM(imComplex,i,j).real();
402  }
403 
404  SPTH(ftrans,imN,H);
405 
407  {
408  if (A2D_ELEM(ROI,i,j))
409  {
410  double temp = std::abs(A2D_ELEM(H,i,j));
411  A2D_ELEM(imModMap,i,j) = std::sqrt(std::pow(temp,2)+std::pow(A2D_ELEM(imN,i,j),2));
412  A2D_ELEM(imN,i,j) = std::cos(std::atan2(temp, A2D_ELEM(imN,i,j)));
413  }
414  else
415  {
416  A2D_ELEM(imModMap,i,j) = 0;
417  A2D_ELEM(imN,i,j) = 0;
418  }
419  }
420 
421  STARTINGX(imN)=STARTINGY(imN)=0;
422 }
#define A2D_ELEM(v, i, j)
void resizeNoCopy(const MultidimArray< T1 > &v)
void sqrt(Image< double > &op)
void abs(Image< double > &op)
#define STARTINGX(v)
#define i
#define FOR_ALL_ELEMENTS_IN_ARRAY2D(m)
#define STARTINGY(v)
void CenterFFT(MultidimArray< T > &v, bool forward)
Definition: xmipp_fft.h:291
#define XSIZE(v)
#define j
void typeCast(const Matrix1D< T1 > &v1, Matrix1D< T2 > &v2)
Definition: matrix1d.h:1227
void SPTH(FourierTransformer &ftrans, MultidimArray< double > &im, MultidimArray< std::complex< double > > &imProc)
float r2

◆ orMinDer()

void orMinDer ( const MultidimArray< double > &  im,
MultidimArray< double > &  orMap,
MultidimArray< double > &  orModMap,
int  wSize,
MultidimArray< bool > &  ROI 
)

Definition at line 149 of file fringe_processing.cpp.

150 {
151  size_t NR, NC,NZ, NDim;
152  im.getDimensions(NC,NR,NZ,NDim);
153 
154  if ( (NZ!=1) || (NDim!=1) )
155  REPORT_ERROR(ERR_MULTIDIM_DIM,(std::string)"ZDim and NDim has to be equals to one");
156 
157  MultidimArray<double > d0,d45,d90,d135,mask,orn,ornMod;
158  d0.resizeNoCopy(im);
159  d45.resizeNoCopy(im);
160  d90.resizeNoCopy(im);
161  d135.resizeNoCopy(im);
162  //mask is the region to perform average in d0, d90, d90, d135 maps
163  //ROI is the region of interested. If ROI is 1 then process in this px if ROI is 0 not process
164  mask.resizeNoCopy(im);
165  orn.resizeNoCopy(im);
166  ornMod.resizeNoCopy(im);
167 
169  {
170  A2D_ELEM(mask,i,j) = 0;
171 
172  if ( (i==0) || (i== (NR-1)) || (j == 0) || (j== (NC-1)) )
173  {
174  DIRECT_A2D_ELEM(d0,i,j) = 0;
175  DIRECT_A2D_ELEM(d45,i,j) = 0;
176  DIRECT_A2D_ELEM(d90,i,j) = 0;
177  DIRECT_A2D_ELEM(d135,i,j)= 0;
178  }
179  else
180  {
181  DIRECT_A2D_ELEM(d0,i,j) = std::sqrt(2)*std::abs(A2D_ELEM(im,i+1,j)-A2D_ELEM(im,i-1,j));
182  DIRECT_A2D_ELEM(d45,i,j) = std::abs(A2D_ELEM(im,i+1,j-1)-A2D_ELEM(im,i-1,j+1));
183  DIRECT_A2D_ELEM(d90,i,j) = std::sqrt(2)*std::abs(A2D_ELEM(im,i,j+1)-A2D_ELEM(im,i,j-1));
184  DIRECT_A2D_ELEM(d135,i,j)= std::abs(A2D_ELEM(im,i+1,j+1)-A2D_ELEM(im,i-1,j-1));
185  }
186  }
187 
188  mask.setXmippOrigin();
189 
190  for (int i = -std::abs(wSize); i <= std::abs(wSize); i++ )
191  for (int j = -std::abs(wSize); j <= std::abs(wSize); j++ )
192  A2D_ELEM(mask,i,j) = double(1)/ double(wSize*wSize+1);
193 
194 
195  convolutionFFT(d0,mask,d0);
196  convolutionFFT(d45,mask,d45);
197  convolutionFFT(d135,mask,d135);
198  convolutionFFT(d90,mask,d90);
199 
200  /*
201  * b=0.5*(D0-D90);
202  * c=-0.5*(D45-D135); %hay que tener en cuenta que la "y" esta downwards
203  * Or=0.5*atan2(c,b);
204  * orn = mod(Or+pi/2,pi);
205  * ornMod=mat2gray(abs(c+1i*b));
206  */
207 
208  STARTINGX(ROI)=STARTINGY(ROI)=0;
209  double tempx = 0;
210  double tempy = 0;
212  {
213  if (A2D_ELEM(ROI,i,j))
214  {
215  tempx = 0.5*(A2D_ELEM(d0,i,j)-A2D_ELEM(d90,i,j));
216  tempy = -0.5*(A2D_ELEM(d45,i,j)-A2D_ELEM(d135,i,j));
217  A2D_ELEM(orn,i,j) = std::fmod((0.5*std::atan2(tempy,tempx))+double(PI)/2,double(PI));
218  A2D_ELEM(ornMod,i,j) = std::sqrt(std::pow(tempx,2)+std::pow(tempy,2));
219  }
220  else
221  {
222  A2D_ELEM(orn,i,j) = 0;
223  A2D_ELEM(ornMod,i,j) = 0;
224  }
225  }
226 
227  //we set again the ROI in the xmipp origin
228  ROI.setXmippOrigin();
229  orMap = orn;
230  orModMap = ornMod;
231 }
#define A2D_ELEM(v, i, j)
#define REPORT_ERROR(nerr, ErrormMsg)
Definition: xmipp_error.h:211
#define FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(m)
void resizeNoCopy(const MultidimArray< T1 > &v)
void sqrt(Image< double > &op)
#define DIRECT_A2D_ELEM(v, i, j)
void abs(Image< double > &op)
#define STARTINGX(v)
#define i
#define FOR_ALL_ELEMENTS_IN_ARRAY2D(m)
#define STARTINGY(v)
void convolutionFFT(MultidimArray< double > &img, MultidimArray< double > &kernel, MultidimArray< double > &result)
Definition: xmipp_fftw.cpp:457
double * d0
#define j
Incorrect MultidimArray dimensions.
Definition: xmipp_error.h:173
#define PI
Definition: tools.h:43
void getDimensions(size_t &Xdim, size_t &Ydim, size_t &Zdim, size_t &Ndim) const

◆ simulPattern()

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.

Definition at line 42 of file fringe_processing.cpp.

43 {
44 
45  //FIRST WE TRAVEL ALL ELEMENTS fn im
46  int ndim = 1;
47  int zdim = 1;
48  im.resize(ndim,zdim,ydim,xdim,false);
49  im.setXmippOrigin();
50 
51  double iMaxDim2 = 2./std::max(xdim,ydim);
52  PolyZernikes polynom;
53 
55  {
57  ROI.resizeNoCopy(im);
58 
59  ROI.setXmippOrigin();
61  A2D_ELEM(ROI,i,j) = true;
62 
63  polynom.zernikePols(coefs,im,ROI);
64  im.setXmippOrigin();
65  }
66 
68  {
69  switch (type)
70  {
71 
72  case SIMPLY_OPEN_FRINGES :
73  {
74  A2D_ELEM(im,i,j) = std::cos(j*iMaxDim2*freq)+rnd_gaus(0,noiseLevel);
75  break;
76  }
78  {
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);
80  break;
81  }
82 
84  {
85  A2D_ELEM(im,i,j) = std::cos(j*iMaxDim2*freq+A2D_ELEM(im,i,j))+rnd_gaus(0,noiseLevel);
86  break;
87  }
89  {
90  A2D_ELEM(im,i,j) = std::cos(50*std::exp(-0.5*(std::pow(i*iMaxDim2,2)+std::pow(j*iMaxDim2,2)))+A2D_ELEM(im,i,j))+rnd_gaus(0,noiseLevel);
91  break;
92  }
93 
95  {
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));
97  break;
98  }
99  }
100  }
101 
102  STARTINGX(im)=STARTINGY(im)=0;
103 }
#define A2D_ELEM(v, i, j)
void resize(size_t Ndim, size_t Zdim, size_t Ydim, size_t Xdim, bool copy=true)
void resizeNoCopy(const MultidimArray< T1 > &v)
void zernikePols(const Matrix1D< int > coef, MultidimArray< double > &im, const MultidimArray< bool > &ROI, int verbose=0)
#define STARTINGX(v)
#define i
#define FOR_ALL_ELEMENTS_IN_ARRAY2D(m)
#define STARTINGY(v)
viol type
void max(Image< double > &op1, const Image< double > &op2)
#define j
double rnd_gaus()

◆ SPTH()

void SPTH ( FourierTransformer ftrans,
MultidimArray< double > &  im,
MultidimArray< std::complex< double > > &  imProc 
)

Function to simulate some test fringe patterns. sd=SPHT(c) computes the quadrture term of c still affected by the direction phase factor. Therefore for a real c=b*cos(phi) sd=SPHT(c)=i*exp(i*dir)*b*sin(phi) Ref: Kieran G. Larkin, Donald J. Bone, and Michael A. Oldfield, "Natural demodulation of two-dimensional fringe patterns. I. General background of the spiral phase quadrature transform," J. Opt. Soc. Am. A 18, 1862-1870 (2001)

Definition at line 111 of file fringe_processing.cpp.

112 {
113  im.setXmippOrigin();
114  MultidimArray<std::complex<double> > H, fftIm, imComplex;
115  typeCast(im, imComplex);
116 
117  // Fourier Transformer
118  ftrans.FourierTransform(imComplex, fftIm, false);
119 
120  H.resizeNoCopy(fftIm);
121  H.setXmippOrigin();
122 
123  //std::complex<double> compTemp(0, 0);
124  //i -> for exterior filas o Y
125  //j -> for interior columns o X
126  for (int i=STARTINGY(H); i<=FINISHINGY(H); i++)
127  {
128  double i2=((double)i)*i;
129  for (int j=STARTINGX(H); j<=FINISHINGX(H); j++)
130  {
131  double j2=((double)j)*j;
132  double iR;
133  if ( (i!=0) || (j!=0))
134  iR=1.0/std::sqrt(i2+j2);
135  else
136  iR=0.;
137  auto *ptr=(double*)&A2D_ELEM(H,i,j);
138  *ptr=j*iR;
139  *(ptr+1)=i*iR;
140  }
141  }
142 
143  fftIm *= H;
144  ftrans.inverseFourierTransform();
145  //Here in the Matlab code there is a complex conjugate s
146  imProc = imComplex;
147 }
#define A2D_ELEM(v, i, j)
void inverseFourierTransform()
Definition: xmipp_fftw.cpp:329
#define FINISHINGX(v)
void resizeNoCopy(const MultidimArray< T1 > &v)
void sqrt(Image< double > &op)
#define STARTINGX(v)
#define i
#define STARTINGY(v)
void FourierTransform(T &v, T1 &V, bool getCopy=true)
Definition: xmipp_fftw.h:166
#define j
#define FINISHINGY(v)
void typeCast(const Matrix1D< T1 > &v1, Matrix1D< T2 > &v2)
Definition: matrix1d.h:1227

◆ unwrapping()

void unwrapping ( const MultidimArray< double > &  wrappedPhase,
MultidimArray< double > &  qualityMap,
double  lambda,
int  size,
MultidimArray< double > &  unwrappedPhase 
)

Definition at line 552 of file fringe_processing.cpp.

553 {
554  //First we perform some setup stuff
555  double minQuality = 0.05;
556 
557  int imax, jmax, max_levels=10;
558  size_t i, j;
559 
560  //We look for the maximun value of qualityMap
561  qualityMap.selfABS();
562  qualityMap.maxIndex(imax,jmax);
563 
564  double maxQualityMapValue = A2D_ELEM(qualityMap,imax,jmax);
565  MultidimArray<double> qualityMapInt = ((qualityMap/maxQualityMapValue)*(max_levels-1));
566  qualityMapInt.selfROUND();
567 
568  Matrix2D<double> gaussian;
569  gaussian.initGaussian(2*size+1,1.5);
570 
571  MultidimArray<bool> processed;
572  processed.resizeNoCopy(wrappedPhase);
573  processed.initZeros();
574 
575  //Here appears the real processing
576  A2D_ELEM(unwrappedPhase,imax,jmax) = A2D_ELEM(wrappedPhase,imax,jmax);
577  A2D_ELEM(processed,imax,jmax) = true;
578 
579  PointQuality p;
580  p.i=imax;
581  p.j=jmax;
582  p.quality=A2D_ELEM(qualityMapInt,imax,jmax);
583  std::priority_queue<PointQuality> queueToProcess;
584  queueToProcess.push(p);
585 
586  //predictor and corrector
587  double pred = 0;
588  double cor = 0;
589  double uw = 0;
590  double wp = 0;
591  double g = 0;
592  int n = 0;
593  double t = 0;
594  double norm = 0;
595  double up = 0;
596 
597  while(!queueToProcess.empty())
598  {
599  p=queueToProcess.top();
600  if ((p.i <=0) || (p.j <= 0)) {
601  REPORT_ERROR(ERR_VALUE_INCORRECT,(std::string)"Access to these elements will cause size_t underflow");
602  }
603  i=(size_t)p.i;
604  j=(size_t)p.j;
605 
606  queueToProcess.pop();
607 
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};
610 
611  for(int k = 0; k< 8 ; k++)
612  {
613  size_t ni=indi[k];
614  size_t nj=indj[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))
617  {
618  wp = A2D_ELEM(wrappedPhase,ni,nj);
619 
620  for (int li=-size; li <= size; li++)
621  {
622  size_t nli=ni+li;
623  for (int lj=-size; lj <= size; lj++)
624  {
625  size_t nlj=nj+lj;
626  if (A2D_ELEM(processed,nli,nlj))
627  {
628  uw = A2D_ELEM(unwrappedPhase,nli,nlj);
629  //q = A2D_ELEM(qualityMap,nli,nlj);
630  g = dMij(gaussian,li+size,lj+size);
631 
632  t = (wp - uw);
633 
634  if ( t > 0 )
635  n = (int)(std::floor(t+3.14159265)/(2*3.14159265));
636  else
637  n = (int)(std::ceil(t-3.14159265)/(2*3.14159265));
638 
639  up = t - (2*3.14159265)*n;
640 
641  //pred += (uw*q*g);
642  //cor += (up*q*g);
643  //norm += (q*g);
644 
645  pred += (uw*g);
646  cor += (up*g);
647  norm += g;
648 
649  }
650  }
651 
652  }
653 
654  A2D_ELEM(unwrappedPhase,ni,nj) = pred/norm + (lambda*cor)/norm;
655  A2D_ELEM(processed,ni,nj) = true;
656 
657  norm = 0;
658  pred = 0;
659  cor = 0;
660 
661  p.i=ni;
662  p.j=nj;
663  p.quality=A2D_ELEM(qualityMapInt,ni,nj);
664  queueToProcess.push(p);
665 
666  }
667  }
668  }
669 }
#define YSIZE(v)
#define A2D_ELEM(v, i, j)
__host__ __device__ float2 floor(const float2 v)
#define REPORT_ERROR(nerr, ErrormMsg)
Definition: xmipp_error.h:211
doublereal * g
void resizeNoCopy(const MultidimArray< T1 > &v)
void maxIndex(int &jmax) const
T norm(const std::vector< T > &v)
Definition: vector_ops.h:399
#define i
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
double * lambda
#define dMij(m, i, j)
Definition: matrix2d.h:139
#define XSIZE(v)
int ni
#define j
void initZeros(const MultidimArray< T1 > &op)
void initGaussian(int Ydim, int Xdim, double op1=0., double op2=1.)
Definition: matrix2d.cpp:1118
Incorrect value received.
Definition: xmipp_error.h:195
int * n