1 /* VipsDeviate
2 *
3 * Copyright: 1990, J. Cupitt
4 *
5 * Author: J. Cupitt
6 * Written on: 02/08/1990
7 * Modified on:
8 * 5/5/93 JC
9 * - now does partial images
10 * - less likely to overflow
11 * - adapted from im_deviate
12 * 1/7/93 JC
13 * - adapted for partial v2
14 * - ANSIfied
15 * 21/2/95 JC
16 * - modernised again
17 * 11/5/95 JC
18 * - oops! return( NULL ) in im_deviate(), instead of return( -1 )
19 * 20/6/95 JC
20 * - now returns double, not float
21 * 13/1/05
22 * - use 64 bit arithmetic
23 * 8/12/06
24 * - add liboil support
25 * 2/9/09
26 * - gtk-doc comment
27 * - minor reformatting
28 * 4/9/09
29 * - use im__wrapscan()
30 * 31/7/10
31 * - remove liboil
32 * 6/11/11
33 * - rewrite as a class
34 */
35
36 /*
37
38 This file is part of VIPS.
39
40 VIPS is free software; you can redistribute it and/or modify
41 it under the terms of the GNU Lesser General Public License as published by
42 the Free Software Foundation; either version 2 of the License, or
43 (at your option) any later version.
44
45 This program is distributed in the hope that it will be useful,
46 but WITHOUT ANY WARRANTY; without even the implied warranty of
47 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
48 GNU Lesser General Public License for more details.
49
50 You should have received a copy of the GNU Lesser General Public License
51 along with this program; if not, write to the Free Software
52 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
53 02110-1301 USA
54
55 */
56
57 /*
58
59 These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk
60
61 */
62
63 #ifdef HAVE_CONFIG_H
64 #include <config.h>
65 #endif /*HAVE_CONFIG_H*/
66 #include <vips/intl.h>
67
68 #include <stdio.h>
69 #include <stdlib.h>
70 #include <math.h>
71
72 #include <vips/vips.h>
73 #include <vips/internal.h>
74
75 #include "statistic.h"
76
77 typedef struct _VipsDeviate {
78 VipsStatistic parent_instance;
79
80 double sum;
81 double sum2;
82 double out;
83 } VipsDeviate;
84
85 typedef VipsStatisticClass VipsDeviateClass;
86
87 G_DEFINE_TYPE( VipsDeviate, vips_deviate, VIPS_TYPE_STATISTIC );
88
89 static int
vips_deviate_build(VipsObject * object)90 vips_deviate_build( VipsObject *object )
91 {
92 VipsObjectClass *class = VIPS_OBJECT_GET_CLASS( object );
93 VipsStatistic *statistic = VIPS_STATISTIC( object );
94 VipsDeviate *deviate = (VipsDeviate *) object;
95
96 gint64 vals;
97 double s, s2;
98
99 if( statistic->in &&
100 vips_check_noncomplex( class->nickname, statistic->in ) )
101 return( -1 );
102
103 if( VIPS_OBJECT_CLASS( vips_deviate_parent_class )->build( object ) )
104 return( -1 );
105
106 /*
107
108 NOTE: NR suggests a two-pass algorithm to minimise roundoff.
109 But that's too expensive for us :-( so do it the old one-pass
110 way.
111
112 */
113
114 /* Calculate and return deviation. Add a fabs to stop sqrt(<=0).
115 */
116 vals = (gint64)
117 vips_image_get_width( statistic->in ) *
118 vips_image_get_height( statistic->in ) *
119 vips_image_get_bands( statistic->in );
120 s = deviate->sum;
121 s2 = deviate->sum2;
122
123 g_object_set( object,
124 "out", sqrt( VIPS_FABS( s2 - (s * s / vals) ) / (vals - 1) ),
125 NULL );
126
127 return( 0 );
128 }
129
130 /* Start function: allocate space for an array in which we can accumulate the
131 * sum and sum of squares for this thread.
132 */
133 static void *
vips_deviate_start(VipsStatistic * statistic)134 vips_deviate_start( VipsStatistic *statistic )
135 {
136 return( (void *) g_new0( double, 2 ) );
137 }
138
139 /* Stop function. Add this little sum to the main sum.
140 */
141 static int
vips_deviate_stop(VipsStatistic * statistic,void * seq)142 vips_deviate_stop( VipsStatistic *statistic, void *seq )
143 {
144 VipsDeviate *deviate = (VipsDeviate *) statistic;
145 double *ss2 = (double *) seq;
146
147 deviate->sum += ss2[0];
148 deviate->sum2 += ss2[1];
149
150 g_free( ss2 );
151
152 return( 0 );
153 }
154
155 #define LOOP( TYPE ) { \
156 TYPE *p = (TYPE *) in; \
157 \
158 for( x = 0; x < sz; x++ ) { \
159 TYPE v = p[x]; \
160 \
161 sum += v; \
162 sum2 += (double) v * (double) v; \
163 } \
164 }
165
166 static int
vips_deviate_scan(VipsStatistic * statistic,void * seq,int x,int y,void * in,int n)167 vips_deviate_scan( VipsStatistic *statistic, void *seq,
168 int x, int y, void *in, int n )
169 {
170 const int sz = n * vips_image_get_bands( statistic->in );
171
172 double *ss2 = (double *) seq;
173
174 double sum;
175 double sum2;
176
177 sum = ss2[0];
178 sum2 = ss2[1];
179
180 /* Now generate code for all types.
181 */
182 switch( vips_image_get_format( statistic->in ) ) {
183 case VIPS_FORMAT_UCHAR: LOOP( unsigned char ); break;
184 case VIPS_FORMAT_CHAR: LOOP( signed char ); break;
185 case VIPS_FORMAT_USHORT: LOOP( unsigned short ); break;
186 case VIPS_FORMAT_SHORT: LOOP( signed short ); break;
187 case VIPS_FORMAT_UINT: LOOP( unsigned int ); break;
188 case VIPS_FORMAT_INT: LOOP( signed int ); break;
189 case VIPS_FORMAT_FLOAT: LOOP( float ); break;
190 case VIPS_FORMAT_DOUBLE: LOOP( double ); break;
191
192 default:
193 g_assert_not_reached();
194 }
195
196 ss2[0] = sum;
197 ss2[1] = sum2;
198
199 return( 0 );
200 }
201
202 static void
vips_deviate_class_init(VipsDeviateClass * class)203 vips_deviate_class_init( VipsDeviateClass *class )
204 {
205 GObjectClass *gobject_class = (GObjectClass *) class;
206 VipsObjectClass *object_class = (VipsObjectClass *) class;
207 VipsStatisticClass *sclass = VIPS_STATISTIC_CLASS( class );
208
209 gobject_class->set_property = vips_object_set_property;
210 gobject_class->get_property = vips_object_get_property;
211
212 object_class->nickname = "deviate";
213 object_class->description = _( "find image standard deviation" );
214 object_class->build = vips_deviate_build;
215
216 sclass->start = vips_deviate_start;
217 sclass->scan = vips_deviate_scan;
218 sclass->stop = vips_deviate_stop;
219
220 VIPS_ARG_DOUBLE( class, "out", 2,
221 _( "Output" ),
222 _( "Output value" ),
223 VIPS_ARGUMENT_REQUIRED_OUTPUT,
224 G_STRUCT_OFFSET( VipsDeviate, out ),
225 -INFINITY, INFINITY, 0.0 );
226 }
227
228 static void
vips_deviate_init(VipsDeviate * deviate)229 vips_deviate_init( VipsDeviate *deviate )
230 {
231 }
232
233 /**
234 * vips_deviate: (method)
235 * @in: input #VipsImage
236 * @out: (out): output pixel standard deviation
237 * @...: %NULL-terminated list of optional named arguments
238 *
239 * This operation finds the standard deviation of all pixels in @in. It
240 * operates on all bands of the input image: use vips_stats() if you need
241 * to calculate an average for each band.
242 *
243 * Non-complex images only.
244 *
245 * See also: vips_avg(), vips_stats()..
246 *
247 * Returns: 0 on success, -1 on error
248 */
249 int
vips_deviate(VipsImage * in,double * out,...)250 vips_deviate( VipsImage *in, double *out, ... )
251 {
252 va_list ap;
253 int result;
254
255 va_start( ap, out );
256 result = vips_call_split( "deviate", ap, in, out );
257 va_end( ap );
258
259 return( result );
260 }
261