1 /* im_maxpos_avg.c
2 *
3 * Copyright: 2006, The Nottingham Trent University
4 * Copyright: 2006, Tom Vajzovic
5 *
6 * Author: Tom Vajzovic
7 *
8 * Written on: 2006-09-25
9 * 15/10/07 JC
10 * - changed spelling of occurrences
11 * - check for !occurrences before using val
12 * - renamed avg as sum, a bit clearer
13 * 2/9/09
14 * - gtkdoc comment
15 * 8/9/08
16 * - rewrite from im_maxpos()
17 * - now handles many bands, complex, faster
18 * 27/7/14
19 * - fix a race ... did not merge states if max was equal
20 * 26/3/15
21 * - avoid NaN, thanks Paul
22 */
23
24 /*
25
26 This file is part of VIPS.
27
28 VIPS is free software; you can redistribute it and/or modify
29 it under the terms of the GNU Lesser General Public License as published by
30 the Free Software Foundation; either version 2 of the License, or
31 (at your option) any later version.
32
33 This program is distributed in the hope that it will be useful,
34 but WITHOUT ANY WARRANTY; without even the implied warranty of
35 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
36 GNU Lesser General Public License for more details.
37
38 You should have received a copy of the GNU Lesser General Public License
39 along with this program; if not, write to the Free Software
40 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
41 02110-1301 USA
42
43 */
44
45 /*
46
47 These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk
48
49 */
50
51 /*
52 #define DEBUG
53 */
54
55 #ifdef HAVE_CONFIG_H
56 #include <config.h>
57 #endif /*HAVE_CONFIG_H*/
58 #include <vips/intl.h>
59
60 #include <stdio.h>
61 #include <stdlib.h>
62 #include <math.h>
63
64 #include <vips/vips.h>
65 #include <vips/vips7compat.h>
66 #include <vips/internal.h>
67
68 /* A position and maximum.
69 */
70 typedef struct _Maxposavg {
71 int xpos;
72 int ypos;
73 double max;
74
75 /* occurences == 0 means we found no points, or we are uninitialised.
76 */
77 int occurences;
78 } Maxposavg;
79
80 /* New sequence value.
81 */
82 static void *
maxposavg_start(IMAGE * in,void * a,void * b)83 maxposavg_start( IMAGE *in, void *a, void *b )
84 {
85 Maxposavg *global_maxposavg = (Maxposavg *) b;
86 Maxposavg *maxposavg;
87
88 if( !(maxposavg = IM_NEW( NULL, Maxposavg )) )
89 return( NULL );
90 *maxposavg = *global_maxposavg;
91
92 return( (void *) maxposavg );
93 }
94
95 /* Merge the sequence value back into the per-call state.
96 */
97 static int
maxposavg_stop(void * seq,void * a,void * b)98 maxposavg_stop( void *seq, void *a, void *b )
99 {
100 Maxposavg *global_maxposavg = (Maxposavg *) b;
101 Maxposavg *maxposavg = (Maxposavg *) seq;
102
103 /* Merge.
104 */
105 if( maxposavg->occurences == 0 ) {
106 }
107 else if( maxposavg->max > global_maxposavg->max )
108 *global_maxposavg = *maxposavg;
109 else if( maxposavg->max == global_maxposavg->max ) {
110 global_maxposavg->xpos += maxposavg->xpos;
111 global_maxposavg->ypos += maxposavg->ypos;
112 global_maxposavg->occurences += maxposavg->occurences;
113 }
114
115 im_free( seq );
116
117 return( 0 );
118 }
119
120 /* int loop.
121 */
122 #define ILOOP( TYPE ) { \
123 TYPE *p = (TYPE *) in; \
124 TYPE m; \
125 \
126 m = max; \
127 \
128 for( x = 0; x < sz; x++ ) { \
129 TYPE v = p[x]; \
130 \
131 if( occurences == 0 || v > m ) { \
132 m = v; \
133 xpos = r->left + x / reg->im->Bands; \
134 ypos = r->top + y; \
135 occurences = 1; \
136 } \
137 else if( v == m ) { \
138 xpos += r->left + x / reg->im->Bands; \
139 ypos += r->top + y; \
140 occurences += 1; \
141 } \
142 } \
143 \
144 max = m; \
145 }
146
147 /* float/double loop ... avoid NaN.
148 */
149 #define FLOOP( TYPE ) { \
150 TYPE *p = (TYPE *) in; \
151 TYPE m; \
152 \
153 m = max; \
154 \
155 for( x = 0; x < sz; x++ ) { \
156 TYPE v = p[x]; \
157 \
158 if( isnan( v ) ) { \
159 } \
160 else if( occurences == 0 || v > m ) { \
161 m = v; \
162 xpos = r->left + x / reg->im->Bands; \
163 ypos = r->top + y; \
164 occurences = 1; \
165 } \
166 else if( v == m ) { \
167 xpos += r->left + x / reg->im->Bands; \
168 ypos += r->top + y; \
169 occurences += 1; \
170 } \
171 } \
172 \
173 max = m; \
174 }
175
176 /* complex/dpcomplex loop ... avoid NaN.
177 */
178 #define CLOOP( TYPE ) { \
179 TYPE *p = (TYPE *) in; \
180 \
181 for( x = 0; x < sz; x++ ) { \
182 double mod, re, im; \
183 \
184 re = p[0]; \
185 im = p[1]; \
186 p += 2; \
187 mod = re * re + im * im; \
188 \
189 if( isnan( mod ) ) { \
190 } \
191 else if( occurences == 0 || mod > max ) { \
192 max = mod; \
193 xpos = r->left + x / reg->im->Bands; \
194 ypos = r->top + y; \
195 occurences = 1; \
196 } \
197 else if( mod == max ) { \
198 xpos += r->left + x / reg->im->Bands; \
199 ypos += r->top + y; \
200 occurences += 1; \
201 } \
202 } \
203 }
204
205 /* Loop over region, adding to seq.
206 */
207 static int
maxposavg_scan(REGION * reg,void * seq,void * a,void * b,gboolean * stop)208 maxposavg_scan( REGION *reg, void *seq, void *a, void *b, gboolean *stop )
209 {
210 const Rect *r = ®->valid;
211 const int sz = IM_REGION_N_ELEMENTS( reg );
212 Maxposavg *maxposavg = (Maxposavg *) seq;
213
214 int x, y;
215 double max;
216 int xpos, ypos, occurences;
217
218 xpos = maxposavg->xpos;
219 ypos = maxposavg->ypos;
220 max = maxposavg->max;
221 occurences = maxposavg->occurences;
222
223 for( y = 0; y < r->height; y++ ) {
224 VipsPel *in = VIPS_REGION_ADDR( reg, r->left, r->top + y );
225
226 switch( reg->im->BandFmt ) {
227 case IM_BANDFMT_UCHAR: ILOOP( unsigned char ); break;
228 case IM_BANDFMT_CHAR: ILOOP( signed char ); break;
229 case IM_BANDFMT_USHORT: ILOOP( unsigned short ); break;
230 case IM_BANDFMT_SHORT: ILOOP( signed short ); break;
231 case IM_BANDFMT_UINT: ILOOP( unsigned int ); break;
232 case IM_BANDFMT_INT: ILOOP( signed int ); break;
233 case IM_BANDFMT_FLOAT: FLOOP( float ); break;
234 case IM_BANDFMT_DOUBLE: FLOOP( double ); break;
235 case IM_BANDFMT_COMPLEX: CLOOP( float ); break;
236 case IM_BANDFMT_DPCOMPLEX: CLOOP( double ); break;
237
238 default:
239 g_assert( 0 );
240 }
241 }
242
243 maxposavg->xpos = xpos;
244 maxposavg->ypos = ypos;
245 maxposavg->max = max;
246 maxposavg->occurences = occurences;
247
248 return( 0 );
249 }
250
251 /**
252 * im_maxpos_avg:
253 * @im: image to scan
254 * @xpos: returned X position
255 * @ypos: returned Y position
256 * @out: returned value
257 *
258 * Function to find the maximum of an image. Returns coords and value at
259 * double precision. In the event of a draw, returns average of all
260 * drawing coords.
261 *
262 * See also: im_maxpos(), im_min(), im_stats().
263 *
264 * Returns: 0 on success, -1 on error
265 */
266 int
im_maxpos_avg(IMAGE * in,double * xpos,double * ypos,double * out)267 im_maxpos_avg( IMAGE *in, double *xpos, double *ypos, double *out )
268 {
269 Maxposavg *global_maxposavg;
270
271 if( im_pincheck( in ) ||
272 im_check_uncoded( "im_maxpos_avg", in ) )
273 return( -1 );
274
275 if( !(global_maxposavg = IM_NEW( in, Maxposavg )) )
276 return( -1 );
277 global_maxposavg->occurences = 0;
278
279 if( vips_sink( in, maxposavg_start, maxposavg_scan, maxposavg_stop,
280 in, global_maxposavg ) )
281 return( -1 );
282
283 if( global_maxposavg->occurences == 0 ) {
284 *xpos = nan("");
285 *ypos = nan("");
286 *out = nan("");
287 }
288 else {
289 /* Back to modulus.
290 */
291 if( vips_band_format_iscomplex( in->BandFmt ) )
292 global_maxposavg->max = sqrt( global_maxposavg->max );
293
294 if( xpos )
295 *xpos = (double) global_maxposavg->xpos /
296 global_maxposavg->occurences;
297 if( ypos )
298 *ypos = (double) global_maxposavg->ypos /
299 global_maxposavg->occurences;
300 if( out )
301 *out = global_maxposavg->max;
302 }
303
304 return( 0 );
305 }
306
307