Xmipp  v3.23.11-Nereus
Classes | Functions
Spectrum modelling by ARMA filters
Collaboration diagram for Spectrum modelling by ARMA filters:

Classes

class  ARMA_parameters
 

Functions

void CausalARMA (MultidimArray< double > &Img, ARMA_parameters &prm)
 
void ARMAFilter (MultidimArray< double > &Img, MultidimArray< double > &Filter, ARMA_parameters &prm)
 

Detailed Description

Function Documentation

◆ ARMAFilter()

void ARMAFilter ( MultidimArray< double > &  Img,
MultidimArray< double > &  Filter,
ARMA_parameters prm 
)

ARMAFilter. This function returns the ARMA Filter associated to an ARMA model. This filter should be applied to a random matrix of independent identically distributed variables. In Xmipp the program to do that is Fourierfilter

PARAMETERS: Img - The matrix - Here it's supposed that it comes from an input image Filter - The matrix that will contain the filter. ARParameters - The matrix with the AR model coeficients, as is returned by CausalARMA or NonCausalARMA MAParameters - The matrix with the MA model coeficients, as is returned by CausalARMA or NonCausalARMA dSigma - The Sigma Coeficient for the ARMA model

OUTPUT: The function stores the output in Filter.

DATE: 26-3-2001

Definition at line 194 of file ctf_estimate_psd_with_arma.cpp.

196 {
197  bool apply_final_median_filter = false;
198  Matrix1D<double> dDigitalFreq(2);
199 
200  // Resize the Filter to the image dimensions
201  Filter.initZeros(YSIZE(Img), XSIZE(Img));
202 
203  // Compute the filter (only half the values are computed)
204  // The other half is computed based in symmetry.
205  int sizeX = XSIZE(Img);
206  int sizeY = YSIZE(Img);
207  long NumberOfMAParameters = YSIZE(prm.MAParameters);
208  long NumberOfARParameters = YSIZE(prm.ARParameters);
209  MultidimArray<int> iMAParameters(NumberOfMAParameters, 2),
210  iARParameters(NumberOfARParameters, 2);
211  for (long n = 0 ; n < NumberOfMAParameters; n++)
212  {
213  DIRECT_A2D_ELEM(iMAParameters, n, 0) = (int)DIRECT_A2D_ELEM(prm.MAParameters, n, 0);
214  DIRECT_A2D_ELEM(iMAParameters, n, 1) = (int)DIRECT_A2D_ELEM(prm.MAParameters, n, 1);
215  }
216  for (long n = 0 ; n < NumberOfARParameters; n++)
217  {
218  DIRECT_A2D_ELEM(iARParameters, n, 0) = (int)DIRECT_A2D_ELEM(prm.ARParameters, n, 0);
219  DIRECT_A2D_ELEM(iARParameters, n, 1) = (int)DIRECT_A2D_ELEM(prm.ARParameters, n, 1);
220  }
221  for (int i = 0;i < sizeY;i++)
222  for (int j = 0;j < (sizeX / 2);j++)
223  {
224  // Compute dDigitalFreq
225  XX(dDigitalFreq) = j / (double)sizeX;
226  YY(dDigitalFreq) = i / (double)sizeY;
227 
228  // Compute B
229  double B = 0;
230  for (long n = 0 ; n < NumberOfMAParameters; n++)
231  {
232  int p = DIRECT_A2D_ELEM(iMAParameters, n, 0);
233  int q = DIRECT_A2D_ELEM(iMAParameters, n, 1);
234  B += DIRECT_A2D_ELEM(prm.MAParameters, n, 2) * 2 *
235  cos((-2) * PI * (p * YY(dDigitalFreq) + q * XX(dDigitalFreq)));
236  }
237  B = B + 1.0;
238 
239  // Compute A
240  double A = 0;
241  for (long n = 0 ; n < NumberOfARParameters; n++)
242  {
243  int p = DIRECT_A2D_ELEM(iARParameters, n, 0);
244  int q = DIRECT_A2D_ELEM(iARParameters, n, 1);
245  A += DIRECT_A2D_ELEM(prm.ARParameters, n, 2) * 2 *
246  cos((-2) * PI * (p * YY(dDigitalFreq) + q * XX(dDigitalFreq)));
247  }
248  A = 1.0 - A;
249 
250  // Check for posible problems calculating the value of A that could lead
251  // to ARMA filters with erroneous values due to a value of A near 0.
252  if (A < 0)
253  apply_final_median_filter = true;
254 
255  // This is the filter proposed by original paper
256  // As the filter values are symmetrical respect the center (frequency zero),
257  // take advantage of this fact.
258  double val2 = prm.dSigma * B / A;
259  //if (val2>=0) val=sqrt(val2);
260  //else {val=abs(sqrt(std::complex<double>(val2,0)));}
261  A2D_ELEM(Filter,sizeY - 1 - i, sizeX - 1 - j) = A2D_ELEM(Filter,i, j) = fabs(val2);
262  }
263 
264  // Apply final median filter
265  if (apply_final_median_filter)
266  {
268  medianFilter3x3(Filter, aux);
269  // Copy all but the borders
270  for (size_t i = 1; i < YSIZE(Filter) - 1; i++)
271  for (size_t j = 1; j < XSIZE(Filter) - 1; j++)
272  DIRECT_A2D_ELEM(Filter, i, j) = DIRECT_A2D_ELEM(aux, i, j);
273  }
274 }
#define YSIZE(v)
MultidimArray< double > MAParameters
#define A2D_ELEM(v, i, j)
#define DIRECT_A2D_ELEM(v, i, j)
#define i
MultidimArray< double > ARParameters
#define XX(v)
Definition: matrix1d.h:85
#define XSIZE(v)
#define j
#define YY(v)
Definition: matrix1d.h:93
void medianFilter3x3(MultidimArray< T > &m, MultidimArray< T > &out)
Definition: filters.h:1088
void initZeros(const MultidimArray< T1 > &op)
#define PI
Definition: tools.h:43
int * n

◆ CausalARMA()

void CausalARMA ( MultidimArray< double > &  Img,
ARMA_parameters prm 
)

CausalARMA.

This function determines the coeficients of an 2D - ARMA model

Img(y,x)=sum( AR(p,q)*Img(y+p,x+q)) + sqrt(sigma) * e(y,x)

Where: (p,q) is a vector belonging to a support region N1 AR(p,q) is the matrix with the AR part coeficients of the model sigma is the variance of the random correlated input e(x,y) e(x,y) is random correlated input with zero mean and cross-spectral density:

E[e(x,y)*Img(x+a,y+b)]= sqrt(sigma) if (p,q)=(0,0) sqrt(sigma)*MA(p,q) if (p,q) belongs to N2 0 otherwise

N1 - The support region for the AR part is the first quadrant region defined by N_AR and M_AR N2 - The support region the the MA part is the second quadrant region defined by N_MA and M_MA

This model is determined following the method depicted in:

R. L. Kashyap, "Characterization and Estimation of Two-Dimensional ARMA Models, " IEEE Transactions on Information Theory, Vol. IT-30, No. 5, pp. 736-745, September 1984.

PARAMETERS: Img - The matrix - Here it's supposed that it comes from an image N_AR, M_AR - The order in Rows and Columns directions of the AR part of the model. N_AR, M_AR - The order in Rows and Columns directions of the MA part of the model. ARParameters - The matrix to return the resulting parameteres for the AR part of the model. MAParameters - The matrix to return the resulting parameteres for the MA part of the model.

OUTPUT: The function stores the ARMA parameters into ARParameters and MAParameters Every row of this output matrices has 3 values: 1nd and 2nd- indicate the support point (p,q) for the coeficient 3th column - indicates the value of the coeficient

DATE: 26-3-2001

Definition at line 92 of file ctf_estimate_psd_with_arma.cpp.

93 {
94  // Calculate the autocorrelation matrix
97  R.setXmippOrigin();
98 
99  /**********************************************************************/
100  // Set equation system for AR part of the model
101  /**********************************************************************/
102  Matrix2D<double> Coeficients;
103  Matrix1D<double> Indep_terms, ARcoeficients;
105 
106  // Assign the support region for the AR part of the model (N1)
108  // Assign the support region for the MA part of the model (N2)
110  // Assign the support region for the AR equations (N3)
111  // Here is the same of N1, but it has not to be
112  First_Quadrant_Neighbors(prm.N_AR, prm.M_AR, N3);
113 
114  long NumberOfARParameters = YSIZE(prm.ARParameters);
115  long NumberOfMAParameters = YSIZE(prm.MAParameters);
116 
117  Coeficients.resizeNoCopy(NumberOfARParameters, NumberOfARParameters);
118  Indep_terms.resizeNoCopy(NumberOfARParameters);
119  ARcoeficients.resizeNoCopy(NumberOfARParameters);
120 
121  // Generate matrix (eq stands for equation number and co for coeficents)
122  for (long eq = 0 ; eq < NumberOfARParameters; eq++)
123  {
124  // take the independet term from the correlation matrix (or calculate it
125  // if it was not calculated before).
126  auto l = (int)A2D_ELEM(N3,eq, 0);
127  auto m = (int)A2D_ELEM(N3,eq, 1);
128  VEC_ELEM(Indep_terms,eq) = A2D_ELEM(R, l, m);
129 
130  // take the coeficients
131  for (long co = 0 ; co < NumberOfARParameters; co++)
132  {
133  // Take the pertinent coeficient form the correlation matrix (or calculate it)
134  auto alpha1 = (int)(A2D_ELEM(N3,eq, 0) - A2D_ELEM(prm.ARParameters,co, 0));
135  auto alpha2 = (int)(A2D_ELEM(N3,eq, 1) - A2D_ELEM(prm.ARParameters,co, 1));
136  auto beta1 = (int)(A2D_ELEM(N3,eq, 0) + A2D_ELEM(prm.ARParameters,co, 0));
137  auto beta2 = (int)(A2D_ELEM(N3,eq, 1) + A2D_ELEM(prm.ARParameters,co, 1));
138  MAT_ELEM(Coeficients,eq, co) = A2D_ELEM(R,alpha1, alpha2) + A2D_ELEM(R,beta1, beta2);
139  }
140  }
141 
142  /**********************************************************************/
143  // Solve the equation system to determine the AR model coeficients and sigma.
144  /**********************************************************************/
145  solveViaCholesky(Coeficients, Indep_terms, ARcoeficients);
146 
147  // Assign the results to the matrix given as parameter
148  for (long n = 0 ; n < NumberOfARParameters; n++)
149  A2D_ELEM(prm.ARParameters, n, 2) = VEC_ELEM(ARcoeficients,n);
150 
151  /**********************************************************************/
152  // Determine the sigma coeficient from the equation for (p,q)=(0,0)
153  /**********************************************************************/
154  double dSum = 0;
155  for (long n = 0 ; n < NumberOfARParameters; n++)
156  {
157  auto p = (int)A2D_ELEM(prm.ARParameters,n, 0);
158  auto q = (int)A2D_ELEM(prm.ARParameters,n, 1);
159  dSum += A2D_ELEM(prm.ARParameters, n, 2) * A2D_ELEM(R, p, q);
160  }
161 
162  // And calculate sigma
163  prm.dSigma = (A2D_ELEM(R, 0, 0) - 2 * dSum);
164  double idSigma = 1.0 / prm.dSigma;
165 
166  /**********************************************************************/
167  // Determine the MA parameters of the model using the AR parameters and
168  // sigma
169  /**********************************************************************/
170  for (long n = 0 ; n < NumberOfMAParameters; n++)
171  {
172  dSum = 0;
173  double MAn0 = A2D_ELEM(prm.MAParameters, n, 0);
174  double MAn1 = A2D_ELEM(prm.MAParameters, n, 1);
175  for (long m = 0 ; m < NumberOfARParameters; m++)
176  {
177  double ARm0 = A2D_ELEM(prm.ARParameters, m, 0);
178  double ARm1 = A2D_ELEM(prm.ARParameters, m, 1);
179  auto alpha1 = (int)(MAn0 - ARm0);
180  auto alpha2 = (int)(MAn1 - ARm1);
181  auto beta1 = (int)(MAn0 + ARm0);
182  auto beta2 = (int)(MAn1 + ARm1);
183  dSum += A2D_ELEM(prm.ARParameters, m, 2) * (
184  A2D_ELEM(R, alpha1, alpha2) + A2D_ELEM(R, beta1, beta2));
185  }
186 
187  auto p = (int)MAn0;
188  auto q = (int)MAn1;
189  A2D_ELEM(prm.MAParameters, n, 2) = (A2D_ELEM(R, p, q) - dSum) * idSigma;
190  }
191 }
#define YSIZE(v)
MultidimArray< double > MAParameters
#define A2D_ELEM(v, i, j)
#define VEC_ELEM(v, i)
Definition: matrix1d.h:245
void Second_Quadrant_Neighbors(int N, int M, MultidimArray< double > &Neighbors)
void First_Quadrant_Neighbors(int N, int M, MultidimArray< double > &Neighbors)
void auto_correlation_matrix(const MultidimArray< double > &Img, MultidimArray< double > &R, CorrelationAux &aux)
Definition: xmipp_fftw.cpp:957
void resizeNoCopy(int Ydim, int Xdim)
Definition: matrix2d.h:534
MultidimArray< double > ARParameters
#define MAT_ELEM(m, i, j)
Definition: matrix2d.h:116
void eq(Image< double > &op1, const Image< double > &op2)
Be careful with integer images for relational operations...due to double comparisons.
int m
void resizeNoCopy(int Xdim)
Definition: matrix1d.h:458
void solveViaCholesky(const Matrix2D< double > &A, const Matrix1D< double > &b, Matrix1D< double > &result)
int * n