w1 = (double)(i1)*mf1; w2 = (double)(i2-m2)*mf2; theta = atan2(w1,w2); c = (float)cos(theta-thetan); c = TWO_THIRDS*c*c;
if (eui<evi) eui = evi; if (theta!=null) theta[i3][i2][i1] = acos(u1i); if (phi!=null) phi[i3][i2][i1] = atan2(u3i,u2i); if (u1!=null) u1[i3][i2][i1] = u1i; if (u2!=null) u2[i3][i2][i1] = u2i;
/** * Finds extrema of output values for a 2nd-order steerable filter. * The filter is comprised of three steering filters for angles 0, * PI/3 and 2*PI/3. * <p> * Given three values output from the three steering filters, this method * computes two angles and corresponding steered output values such that * the derivative of steered output with respect to angle is zero. * Returned angles are in the range [0,PI]. * @param f0 the value for steering filter with angle 0*PI/3. * @param f1 the value for steering filter with angle 1*PI/3. * @param f2 the value for steering filter with angle 2*PI/3. * @return array {{theta1,value1},{theta2,value2}}. The angle theta1 * corresponds to value1 with maximum absolute value; the angle theta2 * corresponds to value2 with minimum absolute value. Both theta1 and * theta2 are in the range [0,PI] and the absolute difference in these * angles is PI/2. (This method by Dave Hale 2007.12.19.) */ private static double[][] findExtrema(double f0, double f1, double f2) { double theta1 = 0.5*(PI+atan2(SQRT3*(f1-f2),(2.0*f0-f1-f2))); double value1 = eval0(f0,f1,f2,theta1); double theta2 = modPi(theta1+0.5*PI); double value2 = eval0(f0,f1,f2,theta2); if (abs(value1)>=abs(value2)) { return new double[][]{{theta1,value1},{theta2,value2}}; } else { return new double[][]{{theta2,value2},{theta1,value1}}; } }
/** * Computes eigenvalues of a symmetric 3x3 matrix using Cardano's * analytical method. */ private static void getEigenvaluesSymmetric33(double[][] a, double[] d) { double a00 = a[0][0], a01 = a[0][1], a11 = a[1][1], a02 = a[0][2], a12 = a[1][2], a22 = a[2][2]; double de = a01*a12; double dd = a01*a01; double ee = a12*a12; double ff = a02*a02; double c2 = a00+a11+a22; double c1 = (a00*a11+a00*a22+a11*a22)-(dd+ee+ff); double c0 = a22*dd+a00*ee+a11*ff-a00*a11*a22-2.0*a02*de; double p = c2*c2-3.0*c1; double q = c2*(p-1.5*c1)-13.5*c0; // 13.5 = 27/2 double t = 27.0*(0.25*c1*c1*(p-c1)+c0*(q+6.75*c0)); // 6.75 = 27/4 double phi = ONE_THIRD*atan2(sqrt(abs(t)),q); double sqrtp = sqrt(abs(p)); double c = sqrtp*cos(phi); double s = ONE_OVER_SQRT3*sqrtp*sin(phi); double dt = ONE_THIRD*(c2-c); d[0] = dt+c; d[1] = dt+s; d[2] = dt-s; }
private static Cdouble[] adjustPoles(double sigma, Cdouble[] poles) { // Simple search for scale factor q that yields the desired sigma. double q = sigma; double s = computeSigma(q,poles); for (int iter=0; abs(sigma-s)>sigma*1.0e-8; ++iter) { //System.out.println("sigma="+sigma+" s="+s+" q="+q); Check.state(iter<100,"number of iterations less than 100"); s = computeSigma(q,poles); q *= sigma/s; } // Adjust poles. int npole = poles.length; Cdouble[] apoles = new Cdouble[npole]; for (int ipole=0; ipole<npole; ++ipole) { Cdouble pi = poles[ipole]; double a = pow(pi.abs(),2.0/q); double t = atan2(pi.i,pi.r)*2.0/q; apoles[ipole] = Cdouble.polar(a,t).inv(); } return apoles; }
private static double computeSigma(double sigma, Cdouble[] poles) { int npole = poles.length; double q = sigma/2.0; Cdouble c1 = new Cdouble(1.0); Cdouble cs = new Cdouble(); for (int ipole=0; ipole<npole; ++ipole) { Cdouble pi = poles[ipole]; double a = pow(pi.abs(),-1.0/q); double t = atan2(pi.i,pi.r)/q; Cdouble b = Cdouble.polar(a,t); Cdouble c = c1.minus(b); Cdouble d = c.times(c); cs.plusEquals(b.times(2.0).over(d)); } return sqrt(cs.r); } }