Note
Function that returns p_acv which contains parameters on anisotropy.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=R8), | intent(in), | dimension(1:long, 1:larg) | :: | tabin |
surface acf array |
|
| integer(kind=I4), | intent(in) | :: | long |
surface length |
||
| integer(kind=I4), | intent(in) | :: | larg |
surface width |
||
| real(kind=R8), | intent(out), | dimension(1:8) | :: | p_acv |
vector containing anisotropy outputs |
|
| real(kind=R8), | intent(in), | optional | :: | cut |
cut height |
|
| real(kind=R8), | intent(in), | dimension(1:2) | :: | scale_xy |
lag along x and y in micrometers |
|
| logical(kind=I4), | intent(in) | :: | omp |
multithreaded ? |
subroutine ellipse_acf(tabin, long, larg, p_acv, cut, scale_xy, omp) !================================================================================================ !< @note Function that returns p_acv which contains parameters on anisotropy. !< !< - p_acv(1) = axe_a, ellipsis big axis !< - p_acv(2) = axe_b, ellipsis small axis !< - p_acv(3) = axe_a/axe_b another anisotropy factor !< - p_acv(4) = nint(angle/inc_a), main texture orientation !< - p_acv(5) = ray_pente, radius of greatest slope !< - p_acv(6) = max_pente, greatest slope !< - p_acv(7) = max_pente/min_pente slope anisotropy factor !< - p_acv(8) = highest curvature/smallest curvature, curvature anisotropy factor !< !< @endnote !------------------------------------------------------------------------------------------------ implicit none integer(kind=I4), intent(in ) :: long !! *surface length* integer(kind=I4), intent(in ) :: larg !! *surface width* logical(kind=I4), intent(in ) :: omp !! *multithreaded ?* real (kind=R8), intent(in ), optional :: cut !! *cut height* real (kind=R8), intent(in ), dimension(1:long, 1:larg) :: tabin !! *surface acf array* real (kind=R8), intent(in ), dimension(1:2) :: scale_xy !! *lag along x and y in micrometers* real (kind=R8), intent(out), dimension(1:8) :: p_acv !! *vector containing anisotropy outputs* integer(kind=I4) :: i, k, ll, p, q, qp, qm, nb_p, nx, ny, lo2, la2, funit, nb_th logical(kind=I4) :: verif real (kind=R8) :: x, y, xb, yb, xm, ym, xp, yp, inc_a real (kind=R8) :: r, r1, rc, theta, h0, hh, h1, h2, h3, h4, coupe, pente_locale, pente_locale_precede real (kind=R8) :: angle, axe_a, axe_b, ech_x, ech_y, ech_z, ech_r real (kind=R8) :: min_pente, max_pente, ray_pente real (kind=R8) :: c, s integer(kind=I4), dimension(0:179) :: p_min, e_angle integer(kind=I4), dimension(1:2) :: loc_max!, imax_acv real (kind=R8), dimension(0:179) :: pente_max real (kind=R8), dimension(0:359) :: ellipse real (kind=R8), allocatable, dimension(:,:) :: tabou, tab_tmp real (kind=R8), allocatable, dimension(:) :: courbure ech_x = scale_xy(1) !SCALE_IMG%dx * unit2IUf(SCALE_IMG%dx_unit) / 1.e-6 ! x lag in micron ech_y = scale_xy(2) !SCALE_IMG%dy * unit2IUf(SCALE_IMG%dy_unit) / 1.e-6 ! y lag in micron ech_z = 1. ! acf(0,0) = 1 loc_max(1:2) = maxloc( tabin(1:long, 1:larg) ) lo2 = loc_max(1) la2 = loc_max(2) ! on prend une surface carrée inscrite ll = min(lo2, la2) - 1 allocate( tabou(0:ll, 0:359) ) allocate( courbure( 0:359) ) tabou = EPS_R8 courbure = EPS_R8 !~ ! avec fftw, il faut mettre le max de l'acv en (lo2,la2) !~ imax_acv = maxloc( tabin(1:long, 1:larg) ) allocate( tab_tmp(1:long, 1:larg) ) tab_tmp(1:long, 1:larg) = tabin(1:long, 1:larg) verif = .false. if ( verif ) then ! sortie xyz call get_unit(funit) open(funit,file='out/test_pol.dat') endif ! angle increment for polar representation of the rosette inc_a = 2*PI_R8/360 ! determination of heights on a rosette diameter, obtained by linear interpolation : ! a point “falls” into a rectangular element [h1,h2,h3,h4], its height is determined by lin. interp. nb_p = 0 nb_th = 1 if (omp) nb_th = omp_get_num_procs() !$OMP PARALLEL DEFAULT(SHARED) NUM_THREADS(nb_th) IF (omp) !$OMP DO SCHEDULE (STATIC,(ll+1)/nb_th) PRIVATE(p,r,q,theta,x,y,xm,ym,h1,h2,h3,h4,hh) do p = 0, ll ! identifying a point on the diameter r = p ! corresponding algebraic radius do q = 0, 359 ! angular increment identification theta = q*inc_a ! projection on x and y of the point marked by its radius and angle, taking the lower integer, ! gives the number of the bottom line and left-hand column of the rectangle ! the remainder (x-nx) represents the abscissa of the point in the rectangle with sides=1 ! the 0.9999 coefficient is used to avoid falling right on an existing point x = lo2 + r * cos(theta) * 0.9999_R8 ; nx = floor(x) ; xb = x -nx y = la2 + r * sin(theta) * 0.9999_R8 ; ny = floor(y) ; yb = y -ny xm = UN -xb ; xp = xb ym = UN -yb ; yp = yb if ( nx+1 <= long .and. ny+1 <= larg .and. & nx >= 1 .and. ny >= 1) then ! attention r may be greater than lo2 or la2 h1 = tab_tmp(nx , ny ) h2 = tab_tmp(nx +1, ny ) h3 = tab_tmp(nx +1, ny +1) h4 = tab_tmp(nx , ny +1) hh = h1*xm*ym + h2*xp*ym + h3*xp*yp + h4*xm*yp tabou(p, q) = hh nb_p = nb_p +1 else hh = 0. endif if ( verif ) write(funit,*) real(x,kind=R4), real(y,kind=R4), real(hh,kind=R4) enddo enddo !$OMP END DO !$OMP END PARALLEL deallocate( tab_tmp ) if ( verif ) close(funit) if ( present(cut) ) then coupe = cut else coupe = 0.5 endif do q = 0, 359 ! angle en deg theta = q*inc_a ! angle en rad ech_r = sqrt( (ech_x*cos(theta))**2 + & ! (ech_y*sin(theta))**2 ) ! unit according to angle q ellipse(q) = (ll -1)*ech_r ! max value of ellipse radius do p = 0, ll -1 ! identifying a point on the diameter r1 = p ! algebraic radius h1 = tabou( p , q ) h2 = tabou( p +1, q ) if (abs(h2) < 10*EPS_R8) exit ! useful for ll = floor(sqrt(UN*(lo2**2 +la2**2))) -1 if (h1 > coupe .and. h2 < coupe) then rc = r1 +(h1 -coupe)/(h1 -h2) ellipse(q) = rc*ech_r exit ! if you don't pass here: no intersection with the cutting plane -> max value taken by default endif enddo ! curvature averaged over 3 points rc = 0 do p = 1, 3 h0 = tabou( p - 1, q ) h1 = tabou( p , q ) h2 = tabou( p + 1, q ) if (abs(h2) < 10*EPS_R8) exit ! useful for ll = floor(sqrt(UN*(lo2**2 +la2**2))) -1 rc = rc + ( h0 - 2 * h1 + h2 ) enddo courbure(q) = - ( rc / ( 2 * ech_r**2 ) ) / 3. enddo do q = 0, 179 e_angle(q) = 2 * q ! angle doubled for the continuity of sin and cos enddo call sort_array2( tab_inout = ellipse(0:179), & ! tab1 = e_angle(0:179), & ! n = 180 ) ! axe_b = sum(ellipse( 0: 2))/3. axe_a = sum(ellipse(177:179))/3. c = sum( cos(e_angle(175:179) * PI_R8 / 180) ) / 5. s = sum( sin(e_angle(175:179) * PI_R8 / 180) ) / 5. angle = atan2(s, c) * 180 / PI_R8 / 2 ! angle halfed (see above) p_acv(1) = axe_a p_acv(2) = axe_b p_acv(3) = axe_a/axe_b p_acv(4) = nint(angle) !----------------------------------------------- p_min(0:179) = ll -1 do q = 0, 179 ! incrément angulaire qp = q +1 ; if (qp>179) qp = 0 qm = q -1 ; if (qm<0 ) qm = 179 theta = q*inc_a ech_r = sqrt( (ech_x*cos(theta))**2 + & ! (ech_y*sin(theta))**2 ) ! unit according to angle q pente_locale_precede = 0._R8 do p = 1, ll-1 pente_locale = ( (tabou(p+1, q ) -tabou(p-1, q )) / 2 )**2 + & ! ( (tabou(p , qp) -tabou(p , qm)) / (2 * p * inc_a) )**2 ! pente_locale = sqrt( pente_locale ) * ech_z / ech_r if ( abs( pente_locale ) < abs( pente_locale_precede ) ) then p_min(q) = p ! the steepest slope distance is recorded exit else pente_locale_precede = pente_locale endif enddo enddo call sort_array2( tab_inout = p_min(0:179), & ! n = 180 ) ! k = int( sum( p_min(0:4) ) / 5., kind = I4 ) ! the minimum distance do q = 0, 179 ! angular increment qp = q +1 ; if (qp > 179) qp = 0 qm = q -1 ; if (qm < 0 ) qm = 179 theta = q * inc_a ech_r = sqrt( ( ech_x * cos(theta) )**2 + & ! (ech_y * sin(theta) )**2 ) ! unit according to angle q pente_locale = ( (tabou(k+1, q ) - tabou(k-1, q )) / 2 )**2 + & ! ( (tabou(k , qp) - tabou(k , qm)) / (2 * k * inc_a) )**2 ! pente_locale = sqrt( pente_locale ) * ech_z / ech_r pente_max(q) = abs( pente_locale ) ! for angle q, greater slope enddo do q = 0, 179 e_angle(q) = q enddo call sort_array2( tab_inout = pente_max(0:179), & ! tab1 = e_angle(0:179), & ! n = 180 ) ! angle = sum( e_angle(175:179) ) / 5. theta = angle*inc_a ech_r = sqrt( ( ech_x * cos(theta) )**2 + & ! ( ech_y * sin(theta) )**2 ) ! unit according to angle q min_pente = sum( pente_max( 0: 2) ) / 3. max_pente = sum( pente_max(177:179) ) / 3. ray_pente = k * ech_r p_acv(5) = ray_pente p_acv(6) = max_pente p_acv(7) = max_pente / min_pente ! anisotropy indicator p_acv(8) = maxval( courbure(0:359) ) / minval( courbure(0:359) ) deallocate( tabou, courbure ) return endsubroutine ellipse_acf