Line data Source code
1 : !> \file mo_spatialsimilarity.f90
2 : !> \brief \copybrief mo_spatialsimilarity
3 : !> \details \copydetails mo_spatialsimilarity
4 :
5 : !> \brief Routines for bias insensitive comparison of spatial patterns.
6 : !> \details These routines are based on the idea that spatial similarity can be assessed by comparing
7 : !! the magnitude of neighboring pixels (e.g. is the neighboring pixel larger or smaller).
8 : !> \author Matthias Zink
9 : !> \date Mar 2013
10 : !> \copyright Copyright 2005-\today, the CHS Developers, Sabine Attinger: All rights reserved.
11 : !! FORCES is released under the LGPLv3+ license \license_note
12 : MODULE mo_spatialsimilarity
13 :
14 : USE mo_kind, ONLY : i4, sp, dp
15 :
16 : IMPLICIT NONE
17 :
18 : PUBLIC :: NNDV ! number of neighboring dominating values
19 : PUBLIC :: PD ! patter dissimilarity measure
20 :
21 : ! ------------------------------------------------------------------
22 :
23 : ! NAME
24 : ! NNDV - number of neighboring dominating values
25 :
26 : ! PURPOSE
27 : !> \brief Calculates the number of neighboring dominating values, a measure for spatial dissimilarity.
28 : !> \details
29 : !> NNDV = 1 - sum(abs(dominating_neighbors(mat1) - dominating_neighbors(mat2))) / count(mask)
30 : !> dominating_neighbors(mat1) = comparison if pixel is larger than its neighbouring values
31 : !>
32 : !> An array element value is compared with its 8 neighbouring cells to check if these cells are larger
33 : !> than the array element value. The result is a 3x3 matrix in which larger cells are indicated by a
34 : !> true value. For comparison this is done with both input arrays. The resulting array is the sum of
35 : !> substraction of the 3x3 matrices for each of the both arrays. The resulting matrix is afterwards
36 : !> normalized to its available neighbors.
37 : !> Furthermore an average over the entire field is calculated. The valid interval of
38 : !> the values for NNDV is [0..1]. In which 1 indicates full agreement and 0 full dismatching.
39 : !> <pre>
40 : !> EXAMPLE:\n
41 : !> mat1 = | 12 17 1 | , mat2 = | 7 9 12 |
42 : !> | 4 10 11 | | 12 11 11 |
43 : !> | 15 2 20 | | 5 13 7 |
44 : !> booleans determined for every grid cell following fortran array scrolling
45 : !> i.e. (/col1_row1, col1_row2, col1_row3, col2_row1, .. ,col3_row3/),(/3,3/)
46 : !>
47 : !> comp1 = | FFF FFF FTF, FFF FFF FFF, FTT FFT FFF |
48 : !> | FFF TFT TTF, TFT TFF FTT, TFF FFT FFF |
49 : !> | FFF FFF FFF, TTF TFF TTF, FFF FFF FFF |
50 : !>
51 : !> comp2 = | FFF FFT FTT, FFT FFT FTT, FFF FFF FFF |
52 : !> | FFF FFF FFT, FTF FFT TFF, FFT TFF FFF |
53 : !> | FFF TFF TTF, FFF FFF FFF, TTF TFF FFF |
54 : !>
55 : !> NNDVMatrix =
56 : !> abs( count(comp1) - count(comp2) ) = | 1-3, 0-4, 3-0 | = | 2, 4, 3 |
57 : !> | 4-1, 5-3, 2-2 | | 3, 2, 0 |
58 : !> | 0-3, 5-0, 0-3 | | 3, 5, 3 |
59 : !>
60 : !> DISSIMILAR / VALID NEIGH CELLS
61 : !> NNDVMatrix / VALID NEIGH CELLS = | 2, 4, 3 | / | 3, 5, 3 |
62 : !> | 3, 2, 0 | | 5, 8, 5 |
63 : !> | 3, 5, 3 | | 3, 5, 3 |
64 : !>
65 : !> = | 0.66, 0.80, 1.00 |
66 : !> | 0.60, 0.25, 0.00 |
67 : !> | 1.0, 1.00, 1.00 |
68 : !>
69 : !> NNDV = 1 - sum(NNDVMatrix) / count(mask) = 1 - (6.31666666 / 9) = 0.2981
70 : !> </pre>
71 : !>
72 : !> If an optinal mask is given, the calculations are over those locations that correspond to true values in the mask.
73 : !> mat1 and mat2 can be single or double precision. The result will have the same numerical precision.
74 :
75 : ! CALLING SEQUENCE
76 : ! out = NNDV(mat1, mat2, mask=mask, valid=valid)
77 :
78 : ! INDENT(IN)
79 : !> \param[in] "real(sp/dp), dimension(:,:) :: mat1" 2D-array with input numbers
80 : !> \param[in] "real(sp/dp), dimension(:,:) :: mat2" 2D-array with input numbers
81 :
82 : ! INDENT(INOUT)
83 : ! None
84 :
85 : ! INDENT(OUT)
86 : ! None
87 :
88 : ! INDENT(IN), OPTIONAL
89 : !> \param[in] "logical,dimension(:,:),optinal :: mask" 2D-array of logical values with size(mat1/mat2).
90 : !> If present, only those locations in mat1/mat2 having true values in mask are evaluated.
91 :
92 : ! INDENT(INOUT), OPTIONAL
93 : ! None
94 :
95 : ! INDENT(OUT), OPTIONAL
96 : !> \param[out] "logical :: valid" indicates if the function could determine a valid value
97 : !> result can be unvalid if entire mask is .false. for ex.
98 : !> in this case PatternDissim is set to 0 (worst case)
99 :
100 : ! RETURN
101 : !> \return real(sp/dp) :: NNDV — Number of neighboring dominating values
102 :
103 : ! RESTRICTIONS
104 : ! Input values must be floating points.
105 :
106 : ! EXAMPLE
107 : ! mat1 = reshape(/ 1., 2, 3., -999., 5., 6. /, (/3,3/))
108 : ! mat2 = reshape(/ 1., 2, 3., -999., 5., 6. /, (/3,3/))
109 : ! out = NNDV(mat1, mat2, mask=(mat1 >= 0. .and. mat2 >= 0.))
110 : ! -> see also example in test directory
111 :
112 : ! LITERATURE
113 : !> \note routine based on algorithm by Luis Samaniego 2009
114 :
115 : ! HISTORY
116 : !> \author Matthias Zink
117 : !> \date Nov 2012
118 : ! update May 2015 created documentation
119 : INTERFACE NNDV
120 : MODULE PROCEDURE NNDV_sp, NNDV_dp
121 : END INTERFACE NNDV
122 :
123 : ! ------------------------------------------------------------------
124 :
125 : ! NAME
126 : ! PD
127 :
128 : ! PURPOSE
129 : !> \brief Calculates pattern dissimilarity (PD) measure
130 : !> \details
131 : !> PD = 1 - sum(dissimilarity(mat1, mat2)) / count(mask)
132 : !> dissimilarity(mat1, mat2) = comparison if pixel is larger than its neighbouring values
133 : !>
134 : !> An array element value is compared with its 8 neighbouring cells to check if these cells are larger
135 : !> than the array element value. The result is a 3x3 matrix in which larger cells are indicated by a
136 : !> true value. For comparison this is done with both input arrays. The resulting array is the sum of
137 : !> xor values of the 3x3 matrices for each of the both arrays. This means only neighbourhood comparisons
138 : !> which are different in the 2 matrices are counted. This resulting matrix is afterwards normalized to its
139 : !> available neighbors. Furthermore an average over the entire field is calculated. The valid interval of
140 : !> the values for PD is [0..1]. In which 1 indicates full agreement and 0 full dismatching.
141 : !>
142 : !> <pre>
143 : !> EXAMPLE:\n
144 : !> mat1 = | 12 17 1 | , mat2 = | 7 9 12 |
145 : !> | 4 10 11 | | 12 11 11 |
146 : !> | 15 2 20 | | 5 13 7 |
147 : !>
148 : !> booleans determined for every grid cell following fortran array scrolling
149 : !> i.e. (/col1_row1, col1_row2, col1_row3, col2_row1, .. ,col3_row3/),(/3,3/)
150 : !>
151 : !> comp1 = | FFF FFF FTF, FFF FFF FFF, FTT FFT FFF |
152 : !> | FFF TFT TTF, TFT TFF FTT, TFF FFT FFF |
153 : !> | FFF FFF FFF, TTF TFF TTF, FFF FFF FFF |
154 : !>
155 : !> comp2 = | FFF FFT FTT, FFT FFT FTT, FFF FFF FFF |
156 : !> | FFF FFF FFT, FTF FFT TFF, FFT TFF FFF |
157 : !> | FFF TFF TTF, FFF FFF FFF, TTF TFF FFF |
158 : !>
159 : !>
160 : !> xor=neq = | FFF FFT FFT, FFT FFT FTT, FTT FFT FFF |
161 : !> | FFF TFT TTT, TTT TFT TTT, TFT TFT FFF |
162 : !> | FFF TFF TTF, TTF TFF TTF, TTF TFF FFF |
163 : !>
164 : !> DISSIMILAR / VALID NEIGH CELLS
165 : !> PDMatrix = | 2, 4, 3 | / | 3, 5, 3 | = | 0.66, 0.80, 1.00 |
166 : !> | 5, 8, 4 | | 5, 8, 5 | | 1.00, 1.00, 0.80 |
167 : !> | 3, 5, 3 | | 3, 5, 3 | | 1.00, 1.00, 1.00 |
168 : !>
169 : !> PD = 1 - sum(PDMatrix) / count(mask) = 1 - (8.2666666 / 9) = 0.08148
170 : !> </pre>
171 : !>
172 : !> If an optinal mask is given, the calculations are over those locations that correspond to true values in the mask.
173 : !> mat1 and mat2 can be single or double precision. The result will have the same numerical precision.
174 :
175 : ! CALLING SEQUENCE
176 : ! out = PD(mat1, mat2, mask=mask, valid=valid)
177 :
178 : ! INDENT(IN)
179 : !> \param[in] "real(sp/dp), dimension(:,:) :: mat1" 2D-array with input numbers
180 : !> \param[in] "real(sp/dp), dimension(:,:) :: mat2" 2D-array with input numbers
181 :
182 : ! INDENT(INOUT)
183 : ! None
184 :
185 : ! INDENT(OUT)
186 : ! None
187 :
188 : ! INDENT(IN), OPTIONAL
189 : !> \param[in] "logical,dimension(:,:),optinal :: mask" 2D-array of logical values with size(mat1/mat2)
190 : !> If present, only those locations in mat1/mat2 having true values in mask are evaluated.
191 :
192 : ! INDENT(INOUT), OPTIONAL
193 : ! None
194 :
195 : ! INDENT(OUT), OPTIONAL
196 : !> \param[out] "logical,optinal :: valid" indicates if the function could determine a valid value
197 : !> result can be unvalid if entire mask is .false. for ex.
198 : !> in this case PD is set to 0 (worst case)
199 :
200 : ! RETURN
201 : !> \return real(sp/dp) :: PD — pattern dissimilarity measure
202 :
203 : ! RESTRICTIONS
204 : ! Input values must be floating points.
205 :
206 : ! EXAMPLE
207 : ! mat1 = reshape(/ 1., 2, 3., -999., 5., 6. /, (/3,3/))
208 : ! mat2 = reshape(/ 1., 2, 3., -999., 5., 6. /, (/3,3/))
209 : ! out = PD(mat1, mat2, mask=(mat1 >= 0. .and. mat2 >= 0.))
210 : ! -> see also example in test directory
211 :
212 : ! LITERATURE
213 : ! None
214 :
215 : ! HISTORY
216 : !> \author Matthias Zink and Juliane Mai
217 : !> \date Jan 2013
218 : INTERFACE PD
219 : MODULE PROCEDURE PD_sp, PD_dp
220 : END INTERFACE PD
221 :
222 : ! ------------------------------------------------------------------
223 :
224 : CONTAINS
225 :
226 4 : FUNCTION NNDV_sp(mat1, mat2, mask, valid)
227 :
228 : IMPLICIT NONE
229 :
230 : REAL(sp), DIMENSION(:, :), INTENT(IN) :: mat1, mat2
231 : LOGICAL, DIMENSION(:, :), OPTIONAL, INTENT(IN) :: mask
232 : LOGICAL, OPTIONAL, INTENT(OUT) :: valid
233 : REAL(sp) :: NNDV_sp
234 :
235 : INTEGER(i4) :: iCo, iRo
236 : INTEGER(i4) :: noValidPixels
237 : INTEGER(i4), DIMENSION(size(shape(mat1))) :: shapemask
238 0 : INTEGER(i4), DIMENSION(size(mat1, dim = 1), size(mat1, dim = 2)) :: validcount
239 126 : REAL(sp), DIMENSION(size(mat1, dim = 1) + 2_i4, size(mat1, dim = 2) + 2_i4) :: bufferedMat1, bufferedMat2
240 28 : REAL(sp), DIMENSION(size(mat1, dim = 1), size(mat1, dim = 2)) :: NNDVMatrix
241 2 : LOGICAL, DIMENSION(size(mat1, dim = 1) + 2_i4, size(mat1, dim = 2) + 2_i4) :: maske
242 2 : LOGICAL, DIMENSION(size(mat1, dim = 1), size(mat1, dim = 2)) :: validmaske
243 :
244 : ! check if input has all the same dimensions
245 2 : if (present(mask)) then
246 3 : shapemask = shape(mask)
247 : else
248 3 : shapemask = shape(mat1)
249 : end if
250 : !
251 6 : if (any(shape(mat1) .NE. shape(mat2))) &
252 0 : stop 'NNDV_sp: shapes of input matrix 1 and input matrix 2 are not matching'
253 6 : if (any(shape(mat1) .NE. shapemask)) &
254 0 : stop 'NNDV_sp: shapes of input matrices and mask are not matching'
255 : !
256 : ! additional 2 rows and 2 cols added for checking the border cells without crashing the search agorithm
257 : ! so the search windows can always cover 9 cells even if it is checking a border cell
258 : ! buffer rows are initialized as false values within mask, i.e. they are not considered for the calculation
259 : ! of the criteria
260 : !
261 : ! initialize mask with default=.false.
262 62 : maske = .false.
263 2 : if (present(mask)) then
264 14 : maske(2 : (size(maske, dim = 1) - 1_i4), 2 : (size(maske, dim = 2) - 1_i4)) = mask
265 : else
266 14 : maske(2 : (size(maske, dim = 1) - 1_i4), 2 : (size(maske, dim = 2) - 1_i4)) = .true.
267 : end if
268 : !
269 : ! initialize bufferedMat1 & bufferedMat2
270 62 : bufferedMat1 = 0.0_sp
271 62 : bufferedMat2 = 0.0_sp
272 28 : bufferedMat1(2 : (size(maske, dim = 1) - 1), 2 : (size(maske, dim = 2) - 1)) = mat1
273 28 : bufferedMat2(2 : (size(maske, dim = 1) - 1), 2 : (size(maske, dim = 2) - 1)) = mat2
274 : !
275 : ! initialize NNDV
276 26 : NNDVMatrix = 0.0_sp
277 : !
278 2 : NNDV_sp = 0.0_sp
279 8 : do iCo = 2_i4, size(bufferedMat1, dim = 2) - 1
280 26 : do iRo = 2_i4, size(bufferedMat1, dim = 1) - 1
281 18 : if (.NOT. maske(iRo, iCo)) cycle
282 30 : NNDVMatrix(iRo - 1_i4, iCo - 1_i4) = &
283 0 : real(abs(count((bufferedMat1(iRo - 1 : iRo + 1, iCo - 1 : iCo + 1) &
284 : - bufferedMat1(iRo, iCo) > epsilon(0.0_sp)) .AND. &
285 : (maske(iRo - 1 : iRo + 1, iCo - 1 : iCo + 1))) - &
286 : count((bufferedMat2(iRo - 1 : iRo + 1, iCo - 1 : iCo + 1) &
287 : - bufferedMat2(iRo, iCo) > epsilon(0.0_sp)) .AND. &
288 : (maske(iRo - 1 : iRo + 1, iCo - 1 : iCo + 1))) &
289 : ), &
290 405 : sp)
291 : ! count - 1 to exclude referendce cell (iRo, iCo)
292 204 : validcount(iRo - 1_i4, iCo - 1_i4) = count(maske(iRo - 1_i4 : iRo + 1_i4, iCo - 1_i4 : iCo + 1_i4)) - 1_i4
293 : end do
294 : end do
295 : !
296 : ! normalize every pixel to number of valid neighbouring cells (defined by maske) [0..8] --> [0..1]
297 28 : validmaske = (maske(2 : size(maske, dim = 1) - 1_i4, 2 : size(maske, dim = 2) - 1_i4) .and. (validcount > 0_i4))
298 26 : noValidPixels = count(validmaske)
299 2 : if (noValidPixels .GT. 0_i4) then
300 26 : NNDVMatrix = merge(NNDVMatrix / real(validcount, sp), NNDVMatrix, validmaske)
301 : ! average over all pixels
302 26 : NNDV_sp = 1.0_sp - sum(NNDVMatrix, mask = validmaske) / noValidPixels
303 2 : if (present(valid)) valid = .TRUE.
304 : ! case if maske is full of .false. or no valid neighbouring cells available for every pixel (validcount(:,:)=0)
305 : else
306 0 : NNDV_sp = 0.0_sp
307 0 : if (present(valid)) valid = .FALSE.
308 : end if
309 : !
310 2 : END FUNCTION NNDV_sp
311 :
312 6 : FUNCTION NNDV_dp(mat1, mat2, mask, valid)
313 :
314 : IMPLICIT NONE
315 :
316 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: mat1, mat2
317 : LOGICAL, DIMENSION(:, :), OPTIONAL, INTENT(IN) :: mask
318 : LOGICAL, OPTIONAL, INTENT(OUT) :: valid
319 : REAL(dp) :: NNDV_dp
320 :
321 : INTEGER(i4) :: iCo, iRo
322 : INTEGER(i4) :: noValidPixels
323 : INTEGER(i4), DIMENSION(size(shape(mat1))) :: shapemask
324 0 : INTEGER(i4), DIMENSION(size(mat1, dim = 1), size(mat1, dim = 2)) :: validcount
325 189 : REAL(dp), DIMENSION(size(mat1, dim = 1) + 2_i4, size(mat1, dim = 2) + 2_i4) :: bufferedMat1, bufferedMat2
326 42 : REAL(dp), DIMENSION(size(mat1, dim = 1), size(mat1, dim = 2)) :: NNDVMatrix
327 3 : LOGICAL, DIMENSION(size(mat1, dim = 1) + 2_i4, size(mat1, dim = 2) + 2_i4) :: maske
328 3 : LOGICAL, DIMENSION(size(mat1, dim = 1), size(mat1, dim = 2)) :: validmaske
329 :
330 : ! check if input has all the same dimensions
331 3 : if (present(mask)) then
332 6 : shapemask = shape(mask)
333 : else
334 3 : shapemask = shape(mat1)
335 : end if
336 : !
337 9 : if (any(shape(mat1) .NE. shape(mat2))) &
338 0 : stop 'NNDV_dp: shapes of input matrix 1 and input matrix 2 are not matching'
339 9 : if (any(shape(mat1) .NE. shapemask)) &
340 0 : stop 'NNDV_dp: shapes of input matrices and mask are not matching'
341 : !
342 : ! additional 2 rows and 2 cols added for checking the border cells without crashing the search agorithm
343 : ! so the search windows can always cover 9 cells even if it is checking a border cell
344 : ! buffer rows are initialized as false values within mask, i.e. they are not considered for the calculation
345 : ! of the criteria
346 : !
347 : ! initialize mask with default=.false.
348 93 : maske = .false.
349 3 : if (present(mask)) then
350 28 : maske(2 : (size(maske, dim = 1) - 1_i4), 2 : (size(maske, dim = 2) - 1_i4)) = mask
351 : else
352 14 : maske(2 : (size(maske, dim = 1) - 1_i4), 2 : (size(maske, dim = 2) - 1_i4)) = .true.
353 : end if
354 : !
355 : ! initialize bufferedMat1 & bufferedMat2
356 93 : bufferedMat1 = 0.0_dp
357 93 : bufferedMat2 = 0.0_dp
358 42 : bufferedMat1(2 : (size(maske, dim = 1) - 1), 2 : (size(maske, dim = 2) - 1)) = mat1
359 42 : bufferedMat2(2 : (size(maske, dim = 1) - 1), 2 : (size(maske, dim = 2) - 1)) = mat2
360 : !
361 : ! initialize NNDV
362 39 : NNDVMatrix = 0.0_dp
363 : !
364 3 : NNDV_dp = 0.0_dp
365 12 : do iCo = 2_i4, size(bufferedMat1, dim = 2) - 1
366 39 : do iRo = 2_i4, size(bufferedMat1, dim = 1) - 1
367 27 : if (.NOT. maske(iRo, iCo)) cycle
368 30 : NNDVMatrix(iRo - 1_i4, iCo - 1_i4) = &
369 0 : real(abs(count((bufferedMat1(iRo - 1 : iRo + 1, iCo - 1 : iCo + 1) - &
370 : bufferedMat1(iRo, iCo) > epsilon(0.0_dp)) .AND. &
371 : (maske(iRo - 1 : iRo + 1, iCo - 1 : iCo + 1))) - &
372 : count((bufferedMat2(iRo - 1 : iRo + 1, iCo - 1 : iCo + 1) - &
373 : bufferedMat2(iRo, iCo) > epsilon(0.0_dp)) .AND. &
374 : (maske(iRo - 1 : iRo + 1, iCo - 1 : iCo + 1))) &
375 : ), &
376 405 : dp)
377 : ! count - 1 to exclude referendce cell (iRo, iCo)
378 216 : validcount(iRo - 1_i4, iCo - 1_i4) = count(maske(iRo - 1_i4 : iRo + 1_i4, iCo - 1_i4 : iCo + 1_i4)) - 1_i4
379 : end do
380 : end do
381 : !
382 : ! normalize every pixel to number of valid neighbouring cells (defined by maske) [0..8] --> [0..1]
383 42 : validmaske = (maske(2 : size(maske, dim = 1) - 1_i4, 2 : size(maske, dim = 2) - 1_i4) .and. (validcount > 0_i4))
384 39 : noValidPixels = count(validmaske)
385 3 : if (noValidPixels .GT. 0_i4) then
386 26 : NNDVMatrix = merge(NNDVMatrix / real(validcount, dp), NNDVMatrix, validmaske)
387 : ! average over all pixels
388 26 : NNDV_dp = 1.0_dp - sum(NNDVMatrix, mask = validmaske) / noValidPixels
389 2 : if (present(valid)) valid = .TRUE.
390 : ! case if maske is full of .false. or no valid neighbouring cells available for every pixel (validcount(:,:)=0)
391 : else
392 1 : NNDV_dp = 0.0_dp
393 1 : if (present(valid)) valid = .FALSE.
394 : end if
395 : !
396 2 : END FUNCTION NNDV_dp
397 :
398 : ! ----------------------------------------------------------------------------------------------------------------
399 :
400 4 : FUNCTION PD_sp(mat1, mat2, mask, valid)
401 :
402 : IMPLICIT NONE
403 :
404 : REAL(sp), DIMENSION(:, :), INTENT(IN) :: mat1, mat2
405 : LOGICAL, DIMENSION(:, :), OPTIONAL, INTENT(IN) :: mask
406 : LOGICAL, OPTIONAL, INTENT(OUT) :: valid
407 : REAL(sp) :: PD_sp
408 :
409 : INTEGER(i4) :: iCo, iRo
410 : INTEGER(i4) :: noValidPixels
411 : INTEGER(i4), DIMENSION(size(shape(mat1))) :: shapemask
412 0 : INTEGER(i4), DIMENSION(size(mat1, dim = 1), size(mat1, dim = 2)) :: validcount
413 126 : REAL(sp), DIMENSION(size(mat1, dim = 1) + 2_i4, size(mat1, dim = 2) + 2_i4) :: bufferedMat1, bufferedMat2
414 28 : REAL(sp), DIMENSION(size(mat1, dim = 1), size(mat1, dim = 2)) :: PDMatrix
415 2 : LOGICAL, DIMENSION(size(mat1, dim = 1) + 2_i4, size(mat1, dim = 2) + 2_i4) :: maske
416 2 : LOGICAL, DIMENSION(size(mat1, dim = 1), size(mat1, dim = 2)) :: validmaske
417 : ! check if input has all the same dimensions
418 2 : if (present(mask)) then
419 6 : shapemask = shape(mask)
420 : else
421 0 : shapemask = shape(mat1)
422 : end if
423 : !
424 6 : if (any(shape(mat1) .NE. shape(mat2))) &
425 0 : stop 'PD_sp: shapes of input matrix 1 and input matrix 2 are not matching'
426 6 : if (any(shape(mat1) .NE. shapemask)) &
427 0 : stop 'PD_sp: shapes of input matrices and mask are not matching'
428 : !
429 : ! additional 2 rows and 2 cols added for checking the border cells without crashing the search agorithm
430 : ! so the search windows can always cover 9 cells even if it is checking a border cell
431 : ! buffer rows are initialized as false values within mask, i.e. they are not considered for the calculation
432 : ! of the criteria
433 : !
434 : ! initialize mask with default=.false.
435 62 : maske = .false.
436 2 : if (present(mask)) then
437 28 : maske(2 : (size(maske, dim = 1) - 1_i4), 2 : (size(maske, dim = 2) - 1_i4)) = mask
438 : else
439 0 : maske(2 : (size(maske, dim = 1) - 1_i4), 2 : (size(maske, dim = 2) - 1_i4)) = .true.
440 : end if
441 : !
442 : ! initialize bufferedMat1 & bufferedMat2
443 62 : bufferedMat1 = 0.0_sp
444 62 : bufferedMat2 = 0.0_sp
445 28 : bufferedMat1(2 : (size(maske, dim = 1) - 1), 2 : (size(maske, dim = 2) - 1)) = mat1
446 28 : bufferedMat2(2 : (size(maske, dim = 1) - 1), 2 : (size(maske, dim = 2) - 1)) = mat2
447 : !
448 : ! initialize PD
449 26 : PDMatrix = 0.0_sp
450 8 : do iCo = 2_i4, size(bufferedMat1, dim = 2) - 1_i4
451 26 : do iRo = 2_i4, size(bufferedMat1, dim = 1) - 1_i4
452 : ! no calculation at unmasked values
453 18 : if (.NOT. maske(iRo, iCo)) cycle
454 : ! .NEQV. is the fortran standard for .XOR.
455 : ! result is written to -1 column and row because of buffer
456 30 : PDMatrix(iRo - 1_i4, iCo - 1_i4) = &
457 : real(count(( &
458 : ! determine higher neighbouring values in mat1
459 0 : (bufferedMat1(iRo - 1_i4 : iRo + 1_i4, iCo - 1_i4 : iCo + 1_i4) - &
460 : bufferedMat1(iRo, iCo) > epsilon(0.0_sp)) .NEQV. &
461 : ! determine higher neighbouring values in mat2
462 : (bufferedMat2(iRo - 1_i4 : iRo + 1_i4, iCo - 1_i4 : iCo + 1_i4) - &
463 : bufferedMat2(iRo, iCo) > epsilon(0.0_sp)) &
464 : ) &
465 : ! exclude unmasked values
466 : .and. (maske(iRo - 1_i4 : iRo + 1_i4, iCo - 1_i4 : iCo + 1_i4)) &
467 : ), &
468 225 : sp)
469 : ! count - 1 to exclude reference cell / center pixel (iRo, iCo)
470 204 : validcount(iRo - 1_i4, iCo - 1_i4) = count(maske(iRo - 1_i4 : iRo + 1_i4, iCo - 1_i4 : iCo + 1_i4)) - 1_i4
471 : !
472 : end do
473 : end do
474 : !
475 : ! normalize every pixel to number of valid neighbouring cells (defined by maske) [0..8] --> [0..1]
476 28 : validmaske = (maske(2 : size(maske, dim = 1) - 1_i4, 2 : size(maske, dim = 2) - 1_i4) .and. (validcount > 0_i4))
477 26 : noValidPixels = count(validmaske)
478 2 : if (noValidPixels .GT. 0_i4) then
479 26 : PDMatrix = merge(PDMatrix / real(validcount, sp), PDMatrix, validmaske)
480 : ! average over all pixels
481 26 : PD_sp = 1.0_sp - sum(PDMatrix, mask = validmaske) / noValidPixels
482 2 : if (present(valid)) valid = .TRUE.
483 : ! case if maske is full of .false. or no valid neighbouring cells available for every pixel (validcount(:,:)=0)
484 : else
485 0 : PD_sp = 0.0_sp
486 0 : if (present(valid)) valid = .FALSE.
487 : end if
488 : !
489 3 : END FUNCTION PD_sp
490 :
491 3 : FUNCTION PD_dp(mat1, mat2, mask, valid)
492 :
493 : IMPLICIT NONE
494 :
495 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: mat1, mat2
496 : LOGICAL, DIMENSION(:, :), OPTIONAL, INTENT(IN) :: mask
497 : LOGICAL, OPTIONAL, INTENT(OUT) :: valid
498 : REAL(dp) :: PD_dp
499 :
500 : INTEGER(i4) :: iCo, iRo
501 : INTEGER(i4) :: noValidPixels
502 : INTEGER(i4), DIMENSION(size(shape(mat1))) :: shapemask
503 0 : INTEGER(i4), DIMENSION(size(mat1, dim = 1), size(mat1, dim = 2)) :: validcount
504 189 : REAL(dp), DIMENSION(size(mat1, dim = 1) + 2_i4, size(mat1, dim = 2) + 2_i4) :: bufferedMat1, bufferedMat2
505 42 : REAL(dp), DIMENSION(size(mat1, dim = 1), size(mat1, dim = 2)) :: PDMatrix
506 3 : LOGICAL, DIMENSION(size(mat1, dim = 1) + 2_i4, size(mat1, dim = 2) + 2_i4) :: maske
507 3 : LOGICAL, DIMENSION(size(mat1, dim = 1), size(mat1, dim = 2)) :: validmaske
508 : ! check if input has all the same dimensions
509 3 : if (present(mask)) then
510 9 : shapemask = shape(mask)
511 : else
512 0 : shapemask = shape(mat1)
513 : end if
514 : !
515 9 : if (any(shape(mat1) .NE. shape(mat2))) &
516 0 : stop 'PD_dp: shapes of input matrix 1 and input matrix 2 are not matching'
517 9 : if (any(shape(mat1) .NE. shapemask)) &
518 0 : stop 'PD_dp: shapes of input matrices and mask are not matching'
519 : !
520 : ! additional 2 rows and 2 cols added for checking the border cells without crashing the search agorithm
521 : ! so the search windows can always cover 9 cells even if it is checking a border cell
522 : ! buffer rows are initialized as false values within mask, i.e. they are not considered for the calculation
523 : ! of the criteria
524 : !
525 : ! initialize mask with default=.false.
526 93 : maske = .false.
527 3 : if (present(mask)) then
528 42 : maske(2 : (size(maske, dim = 1) - 1_i4), 2 : (size(maske, dim = 2) - 1_i4)) = mask
529 : else
530 0 : maske(2 : (size(maske, dim = 1) - 1_i4), 2 : (size(maske, dim = 2) - 1_i4)) = .true.
531 : end if
532 : !
533 : ! initialize bufferedMat1 & bufferedMat2
534 93 : bufferedMat1 = 0.0_dp
535 93 : bufferedMat2 = 0.0_dp
536 42 : bufferedMat1(2 : (size(maske, dim = 1) - 1), 2 : (size(maske, dim = 2) - 1)) = mat1
537 42 : bufferedMat2(2 : (size(maske, dim = 1) - 1), 2 : (size(maske, dim = 2) - 1)) = mat2
538 : !
539 : ! initialize PD
540 39 : PDMatrix = 0.0_dp
541 12 : do iCo = 2_i4, size(bufferedMat1, dim = 2) - 1_i4
542 39 : do iRo = 2_i4, size(bufferedMat1, dim = 1) - 1_i4
543 : ! no calculation at unmasked values
544 27 : if (.NOT. maske(iRo, iCo)) cycle
545 : ! .NEQV. is the fortran standard for .XOR.
546 : ! result is written to -1 column and row because of buffer
547 30 : PDMatrix(iRo - 1_i4, iCo - 1_i4) = &
548 : real(count(( &
549 : ! determine higher neighbouring values in mat1
550 0 : (bufferedMat1(iRo - 1_i4 : iRo + 1_i4, iCo - 1_i4 : iCo + 1_i4) - &
551 : bufferedMat1(iRo, iCo) > epsilon(0.0_dp)) .NEQV. &
552 : ! determine higher neighbouring values in mat2
553 : (bufferedMat2(iRo - 1_i4 : iRo + 1_i4, iCo - 1_i4 : iCo + 1_i4) - &
554 : bufferedMat2(iRo, iCo) > epsilon(0.0_dp)) &
555 : ) &
556 : ! exclude unmasked values
557 : .and. (maske(iRo - 1_i4 : iRo + 1_i4, iCo - 1_i4 : iCo + 1_i4)) &
558 : ), &
559 225 : dp)
560 : ! count - 1 to exclude referendce cell (iRo, iCo)
561 216 : validcount(iRo - 1_i4, iCo - 1_i4) = count(maske(iRo - 1_i4 : iRo + 1_i4, iCo - 1_i4 : iCo + 1_i4)) - 1_i4
562 : !
563 : end do
564 : end do
565 : !
566 : ! normalize every pixel to number of valid neighbouring cells (defined by maske) [0..8] --> [0..1]
567 42 : validmaske = (maske(2 : size(maske, dim = 1) - 1_i4, 2 : size(maske, dim = 2) - 1_i4) .and. (validcount > 0_i4))
568 39 : noValidPixels = count(validmaske)
569 3 : if (noValidPixels .GT. 0_i4) then
570 26 : PDMatrix = merge(PDMatrix / real(validcount, dp), PDMatrix, validmaske)
571 : ! average over all pixels
572 26 : PD_dp = 1.0_dp - sum(PDMatrix, mask = validmaske) / noValidPixels
573 2 : if (present(valid)) valid = .TRUE.
574 : ! case if maske is full of .false. or no valid neighbouring cells available for every pixel (validcount(:,:)=0)
575 : else
576 1 : PD_dp = 0.0_dp
577 1 : if (present(valid)) valid = .FALSE.
578 : end if
579 : !
580 2 : END FUNCTION PD_dp
581 : !
582 : END MODULE mo_spatialsimilarity
|
