1 /* im_gradcor
2 *
3 * Copyright: 2007 Nottingham Trent University
4 *
5 * Author: Tom Vajzovic
6 * Written on: 2007-06-07
7 * 3/2/10
8 * - gtkdoc
9 * - cleanups
10 */
11
12 /*
13
14 This file is part of VIPS.
15
16 VIPS is free software; you can redistribute it and/or modify
17 it under the terms of the GNU Lesser General Public License as published by
18 the Free Software Foundation; either version 2 of the License, or
19 (at your option) any later version.
20
21 This program is distributed in the hope that it will be useful,
22 but WITHOUT ANY WARRANTY; without even the implied warranty of
23 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 GNU Lesser General Public License for more details.
25
26 You should have received a copy of the GNU Lesser General Public License
27 along with this program; if not, write to the Free Software
28 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
29 02110-1301 USA
30
31 */
32
33 /*
34
35 These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk
36
37 */
38
39
40 /** HEADERS **/
41
42 #ifdef HAVE_CONFIG_H
43 #include <config.h>
44 #endif /*HAVE_CONFIG_H*/
45 #include <vips/intl.h>
46
47 #include <stdlib.h>
48
49 #include <vips/vips.h>
50 #include <vips/vips7compat.h>
51
52
53 /** LOCAL TYPES **/
54
55 typedef struct {
56 REGION *reg;
57 int *region_xgrad;
58 int *region_ygrad;
59 size_t region_xgrad_area;
60 size_t region_ygrad_area;
61 }
62 gradcor_seq_t;
63
64
65 /** LOCAL FUNCTION DECLARATIONS **/
66
67 static void *gradcor_start( IMAGE *out, void *vptr_large, void *unrequired );
68 static int gradcor_stop( void *vptr_seq, void *unrequired, void *unreq2 );
69 static int gradcor_gen( REGION *to_make, void *vptr_seq, void *unrequired, void *vptr_grads );
70
71 #define XGRAD_GEN_DECLARATION( TYPE ) static int xgrad_gen_ ## TYPE( REGION *to_make, void *vptr_make_from, void *unrequired, void *unreq2 )
72 #define YGRAD_GEN_DECLARATION( TYPE ) static int ygrad_gen_ ## TYPE( REGION *to_make, void *vptr_make_from, void *unrequired, void *unreq2 )
73
74 XGRAD_GEN_DECLARATION( guint8 );
75 YGRAD_GEN_DECLARATION( guint8 );
76 XGRAD_GEN_DECLARATION( gint8 );
77 YGRAD_GEN_DECLARATION( gint8 );
78 XGRAD_GEN_DECLARATION( guint16 );
79 YGRAD_GEN_DECLARATION( guint16 );
80 XGRAD_GEN_DECLARATION( gint16 );
81 YGRAD_GEN_DECLARATION( gint16 );
82 XGRAD_GEN_DECLARATION( guint32 );
83 YGRAD_GEN_DECLARATION( guint32 );
84 XGRAD_GEN_DECLARATION( gint32 );
85 YGRAD_GEN_DECLARATION( gint32 );
86 #if 0
87 XGRAD_GEN_DECLARATION( float );
88 YGRAD_GEN_DECLARATION( float );
89 XGRAD_GEN_DECLARATION( double );
90 YGRAD_GEN_DECLARATION( double );
91 #endif
92
93
94 /** EXPORTED FUNCTION DEFINITIONS **/
95
im_gradcor_raw(IMAGE * large,IMAGE * small,IMAGE * out)96 int im_gradcor_raw( IMAGE *large, IMAGE *small, IMAGE *out ){
97 #define FUNCTION_NAME "im_gradcor_raw"
98
99 if( im_piocheck( large, out ) || im_pincheck( small ) )
100 return -1;
101
102 if( im_check_uncoded( "im_gradcor", large ) ||
103 im_check_mono( "im_gradcor", large ) ||
104 im_check_uncoded( "im_gradcor", small ) ||
105 im_check_mono( "im_gradcor", small ) ||
106 im_check_format_same( "im_gradcor", large, small ) ||
107 im_check_int( "im_gradcor", large ) )
108 return( -1 );
109
110 if( large-> Xsize < small-> Xsize || large-> Ysize < small-> Ysize ){
111 im_error( FUNCTION_NAME, "second image must be smaller than first" );
112 return -1;
113 }
114 if( im_cp_desc( out, large ) )
115 return -1;
116
117 out-> Xsize= 1 + large-> Xsize - small-> Xsize;
118 out-> Ysize= 1 + large-> Ysize - small-> Ysize;
119 out-> BandFmt= IM_BANDFMT_FLOAT;
120
121 if( im_demand_hint( out, IM_FATSTRIP, large, NULL ) )
122 return -1;
123
124 {
125 IMAGE *xgrad= im_open_local( out, FUNCTION_NAME ": xgrad", "t" );
126 IMAGE *ygrad= im_open_local( out, FUNCTION_NAME ": ygrad", "t" );
127 IMAGE **grads= im_allocate_input_array( out, xgrad, ygrad, NULL );
128
129 return ! xgrad || ! ygrad || ! grads
130 || im_grad_x( small, xgrad )
131 || im_grad_y( small, ygrad )
132 || im_generate( out, gradcor_start, gradcor_gen, gradcor_stop, (void*) large, (void*) grads );
133 }
134 #undef FUNCTION_NAME
135 }
136
137 /**
138 * im_gradcor:
139 * @in: input image
140 * @ref: reference image
141 * @out: output image
142 *
143 * Calculate a correlation surface.
144 *
145 * @ref is placed at every position in @in and a correlation coefficient
146 * calculated. One-band, integer images only. @in and @ref must have the
147 * same #VipsBandFmt. The output
148 * image is always %IM_BANDFMT_FLOAT. @ref must be smaller than @in. The output
149 * image is the same size as the input.
150 *
151 * The method takes the gradient images of the two images then takes the
152 * dot-product correlation of the two vector images.
153 * The vector expression of this method is my (tcv) own creation. It is
154 * equivalent to the complex-number method of:
155 *
156 * ARGYRIOU, V. et al. 2003. Estimation of sub-pixel motion using
157 * gradient cross correlation. Electronics Letters, 39 (13).
158 *
159 * See also: im_spcor().
160 *
161 * Returns: 0 on success, -1 on error
162 */
163 int
im_gradcor(IMAGE * in,IMAGE * ref,IMAGE * out)164 im_gradcor( IMAGE *in, IMAGE *ref, IMAGE *out )
165 {
166 #define FUNCTION_NAME "im_gradcor"
167 IMAGE *t1 = im_open_local( out, FUNCTION_NAME " intermediate", "p" );
168
169 if( !t1 ||
170 im_embed( in, t1, 1,
171 ref->Xsize / 2, ref->Ysize / 2,
172 in->Xsize + ref->Xsize - 1,
173 in->Ysize + ref->Ysize - 1 ) ||
174 im_gradcor_raw( t1, ref, out ) )
175 return( -1 );
176
177 out->Xoffset = 0;
178 out->Yoffset = 0;
179
180 return( 0 );
181 #undef FUNCTION_NAME
182 }
183
184 /**
185 * im_grad_x:
186 * @in: input image
187 * @out: output image
188 *
189 * Find horizontal differences between adjacent pixels.
190 *
191 * Generates an image where the value of each pixel is the difference between
192 * it and the pixel to its right. The output has the same height as the input
193 * and one pixel less width. One-band integer formats only. The result is
194 * always %IM_BANDFMT_INT.
195 *
196 * This operation is much faster than (though equivalent to) im_conv() with the
197 * mask [[-1, 1]].
198 *
199 * See also: im_grad_y(), im_conv().
200 *
201 * Returns: 0 on success, -1 on error
202 */
im_grad_x(IMAGE * in,IMAGE * out)203 int im_grad_x( IMAGE *in, IMAGE *out ){
204 #define FUNCTION_NAME "im_grad_x"
205
206 if( im_piocheck( in, out ) )
207 return -1;
208
209 if( im_check_uncoded( "im_grad_x", in ) ||
210 im_check_mono( "im_grad_x", in ) ||
211 im_check_int( "im_grad_x", in ) )
212 return( -1 );
213 if( im_cp_desc( out, in ) )
214 return -1;
215
216 -- out-> Xsize;
217 out-> BandFmt= IM_BANDFMT_INT; /* do not change without updating im_gradcor() */
218
219 if( im_demand_hint( out, IM_THINSTRIP, in, NULL ) )
220 return -1;
221
222 #define RETURN_GENERATE( TYPE ) return im_generate( out, im_start_one, xgrad_gen_ ## TYPE, im_stop_one, (void*) in, NULL )
223
224 switch( in-> BandFmt ){
225
226 case IM_BANDFMT_UCHAR:
227 RETURN_GENERATE( guint8 );
228
229 case IM_BANDFMT_CHAR:
230 RETURN_GENERATE( gint8 );
231
232 case IM_BANDFMT_USHORT:
233 RETURN_GENERATE( guint16 );
234
235 case IM_BANDFMT_SHORT:
236 RETURN_GENERATE( gint16 );
237
238 case IM_BANDFMT_UINT:
239 RETURN_GENERATE( guint32 );
240
241 case IM_BANDFMT_INT:
242 RETURN_GENERATE( gint32 );
243 #if 0
244 case IM_BANDFMT_FLOAT:
245 RETURN_GENERATE( float );
246 case IM_BANDFMT_DOUBLE:
247 RETURN_GENERATE( double );
248 #endif
249 #undef RETURN_GENERATE
250 default:
251 g_assert( 0 );
252 }
253
254 /* Keep gcc happy.
255 */
256 return 0;
257 #undef FUNCTION_NAME
258 }
259
260
261 /**
262 * im_grad_y:
263 * @in: input image
264 * @out: output image
265 *
266 * Find vertical differences between adjacent pixels.
267 *
268 * Generates an image where the value of each pixel is the difference between
269 * it and the pixel below it. The output has the same width as the input
270 * and one pixel less height. One-band integer formats only. The result is
271 * always %IM_BANDFMT_INT.
272 *
273 * This operation is much faster than (though equivalent to) im_conv() with the
274 * mask [[-1], [1]].
275 *
276 * See also: im_grad_x(), im_conv().
277 *
278 * Returns: 0 on success, -1 on error
279 */
im_grad_y(IMAGE * in,IMAGE * out)280 int im_grad_y( IMAGE *in, IMAGE *out ){
281 #define FUNCTION_NAME "im_grad_y"
282
283 if( im_piocheck( in, out ) )
284 return -1;
285
286 if( im_check_uncoded( "im_grad_y", in ) ||
287 im_check_mono( "im_grad_y", in ) ||
288 im_check_int( "im_grad_y", in ) )
289 return( -1 );
290
291 if( im_cp_desc( out, in ) )
292 return -1;
293
294 -- out-> Ysize;
295 out-> BandFmt= IM_BANDFMT_INT; /* do not change without updating im_gradcor() */
296
297 if( im_demand_hint( out, IM_FATSTRIP, in, NULL ) )
298 return -1;
299
300 #define RETURN_GENERATE( TYPE ) return im_generate( out, im_start_one, ygrad_gen_ ## TYPE, im_stop_one, (void*) in, NULL )
301
302 switch( in-> BandFmt ){
303
304 case IM_BANDFMT_UCHAR:
305 RETURN_GENERATE( guint8 );
306
307 case IM_BANDFMT_CHAR:
308 RETURN_GENERATE( gint8 );
309
310 case IM_BANDFMT_USHORT:
311 RETURN_GENERATE( guint16 );
312
313 case IM_BANDFMT_SHORT:
314 RETURN_GENERATE( gint16 );
315
316 case IM_BANDFMT_UINT:
317 RETURN_GENERATE( guint32 );
318
319 case IM_BANDFMT_INT:
320 RETURN_GENERATE( gint32 );
321 #if 0
322 case IM_BANDFMT_FLOAT:
323 RETURN_GENERATE( float );
324 case IM_BANDFMT_DOUBLE:
325 RETURN_GENERATE( double );
326 #endif
327 #undef RETURN_GENERATE
328 default:
329 g_assert( 0 );
330 }
331
332 /* Keep gcc happy.
333 */
334 return 0;
335 #undef FUNCTION_NAME
336 }
337
338
339 /** LOCAL FUNCTION DEFINITIONS **/
340
gradcor_start(IMAGE * out,void * vptr_large,void * unrequired)341 static void *gradcor_start( IMAGE *out, void *vptr_large, void *unrequired ){
342
343 gradcor_seq_t *seq= IM_NEW( NULL, gradcor_seq_t );
344 if( ! seq )
345 return NULL;
346
347 seq-> region_xgrad= (int*) NULL;
348 seq-> region_ygrad= (int*) NULL;
349 seq-> region_xgrad_area= 0;
350 seq-> region_ygrad_area= 0;
351
352 seq-> reg= im_region_create( (IMAGE*) vptr_large );
353 if( ! seq-> reg ){
354 im_free( (void*) seq );
355 return NULL;
356 }
357 return (void*) seq;
358 }
359
gradcor_stop(void * vptr_seq,void * unrequired,void * unreq2)360 static int gradcor_stop( void *vptr_seq, void *unrequired, void *unreq2 ){
361
362 gradcor_seq_t *seq= (gradcor_seq_t*) vptr_seq;
363 if( seq ){
364 im_free( (void*) seq-> region_xgrad );
365 im_free( (void*) seq-> region_ygrad );
366 im_region_free( seq-> reg );
367 seq-> region_xgrad= (int*) NULL;
368 seq-> region_ygrad= (int*) NULL;
369 seq-> reg= (REGION*) NULL;
370 im_free( (void*) seq );
371 }
372 return 0;
373 }
374
gradcor_gen(REGION * to_make,void * vptr_seq,void * unrequired,void * vptr_grads)375 static int gradcor_gen( REGION *to_make, void *vptr_seq, void *unrequired, void *vptr_grads ){
376
377 gradcor_seq_t *seq= (gradcor_seq_t*) vptr_seq;
378 REGION *make_from= seq-> reg;
379
380 IMAGE **grads= (IMAGE**) vptr_grads;
381 IMAGE *small_xgrad= grads[0];
382 IMAGE *small_ygrad= grads[1];
383
384 Rect require= {
385 to_make-> valid. left,
386 to_make-> valid. top,
387 to_make-> valid. width + small_xgrad-> Xsize,
388 to_make-> valid. height + small_ygrad-> Ysize
389 };
390 size_t region_xgrad_width= require. width - 1;
391 size_t region_ygrad_height= require. height - 1;
392
393 if( im_prepare( make_from, &require ) )
394 return -1;
395
396 #define FILL_BUFFERS( TYPE ) /* fill region_xgrad */ \
397 { \
398 TYPE *reading= (TYPE*) IM_REGION_ADDR( make_from, require. left, require. top ); \
399 size_t read_skip= ( IM_REGION_LSKIP( make_from ) / sizeof(TYPE) ) - region_xgrad_width; \
400 size_t area_need= region_xgrad_width * require. height; \
401 \
402 if( seq-> region_xgrad_area < area_need ){ \
403 free( seq-> region_xgrad ); \
404 seq-> region_xgrad= malloc( area_need * sizeof(int) ); \
405 if( ! seq-> region_xgrad ) \
406 return -1; \
407 seq-> region_xgrad_area= area_need; \
408 } \
409 { \
410 int *writing= seq-> region_xgrad; \
411 int *write_end= writing + area_need; \
412 int *write_stop; \
413 for( ; writing < write_end; reading+= read_skip ) \
414 for( write_stop= writing + region_xgrad_width; writing < write_stop; ++reading, ++writing ) \
415 *writing= reading[1] - reading[0]; \
416 } \
417 } \
418 { /* fill region_ygrad */ \
419 TYPE *reading= (TYPE*) IM_REGION_ADDR( make_from, require. left, require. top ); \
420 size_t read_line= IM_REGION_LSKIP( make_from ) / sizeof(TYPE); \
421 size_t read_skip= read_line - require. width; \
422 size_t area_need= require. width * region_ygrad_height; \
423 \
424 if( seq-> region_ygrad_area < area_need ){ \
425 free( seq-> region_ygrad ); \
426 seq-> region_ygrad= malloc( area_need * sizeof(int) ); \
427 if( ! seq-> region_ygrad ) \
428 return -1; \
429 seq-> region_ygrad_area= area_need; \
430 } \
431 { \
432 int *writing= seq-> region_ygrad; \
433 int *write_end= writing + area_need; \
434 int *write_stop; \
435 for( ; writing < write_end; reading+= read_skip ) \
436 for( write_stop= writing + require. width; writing < write_stop; ++reading, ++writing ) \
437 *writing= reading[ read_line ] - reading[0]; \
438 } \
439 }
440 switch( make_from-> im-> BandFmt ){
441 case IM_BANDFMT_UCHAR:
442 FILL_BUFFERS( unsigned char )
443 break;
444 case IM_BANDFMT_CHAR:
445 FILL_BUFFERS( signed char )
446 break;
447 case IM_BANDFMT_USHORT:
448 FILL_BUFFERS( unsigned short int )
449 break;
450 case IM_BANDFMT_SHORT:
451 FILL_BUFFERS( signed short int )
452 break;
453 case IM_BANDFMT_UINT:
454 FILL_BUFFERS( unsigned int )
455 break;
456 case IM_BANDFMT_INT:
457 FILL_BUFFERS( signed int )
458 break;
459 default:
460 g_assert( 0 );
461 }
462 { /* write to output */
463 size_t write_skip= IM_REGION_LSKIP( to_make ) / sizeof( float );
464 float *writing= (float*) IM_REGION_ADDR_TOPLEFT( to_make );
465 float *write_end= writing + write_skip * to_make-> valid. height;
466 float *write_stop;
467 size_t write_width= to_make-> valid. width;
468
469 size_t small_xgrad_width= small_xgrad-> Xsize;
470 size_t small_ygrad_width= small_ygrad-> Xsize;
471 int *small_xgrad_end= (int*) small_xgrad-> data + small_xgrad_width * small_xgrad-> Ysize;
472 int *small_ygrad_end= (int*) small_ygrad-> data + small_ygrad_width * small_ygrad-> Ysize;
473
474 int *region_xgrad_start= seq-> region_xgrad;
475 int *region_ygrad_start= seq-> region_ygrad;
476 size_t region_xgrad_start_skip= region_xgrad_width - write_width;
477 size_t region_ygrad_start_skip= require. width - write_width;
478
479 size_t region_xgrad_read_skip= region_xgrad_width - small_xgrad_width;
480 size_t region_ygrad_read_skip= require. width - small_ygrad_width;
481
482 write_skip-= write_width;
483
484 for( ; writing < write_end; writing+= write_skip, region_xgrad_start+= region_xgrad_start_skip, region_ygrad_start+= region_ygrad_start_skip )
485 for( write_stop= writing + write_width; writing < write_stop; ++writing, ++region_xgrad_start, ++region_ygrad_start ){
486 gint64 sum= 0;
487 {
488 int *small_xgrad_read= (int*) small_xgrad-> data;
489 int *small_xgrad_stop;
490 int *region_xgrad_read= region_xgrad_start;
491
492 for( ; small_xgrad_read < small_xgrad_end; region_xgrad_read+= region_xgrad_read_skip )
493 for( small_xgrad_stop= small_xgrad_read + small_xgrad_width; small_xgrad_read < small_xgrad_stop; ++small_xgrad_read, ++region_xgrad_read )
494 sum+= *small_xgrad_read * *region_xgrad_read;
495 }
496 {
497 int *small_ygrad_read= (int*) small_ygrad-> data;
498 int *small_ygrad_stop;
499 int *region_ygrad_read= region_ygrad_start;
500
501 for( ; small_ygrad_read < small_ygrad_end; region_ygrad_read+= region_ygrad_read_skip )
502 for( small_ygrad_stop= small_ygrad_read + small_ygrad_width; small_ygrad_read < small_ygrad_stop; ++small_ygrad_read, ++region_ygrad_read )
503 sum+= *small_ygrad_read * *region_ygrad_read;
504 }
505 *writing= sum;
506 }
507 }
508 return 0;
509 }
510
511 #define XGRAD_GEN_DEFINITION( TYPE ) \
512 static int xgrad_gen_ ## TYPE( REGION *to_make, void *vptr_make_from, void *unrequired, void *unreq2 ){ \
513 \
514 REGION *make_from= (REGION*) vptr_make_from; \
515 Rect require= { \
516 to_make-> valid. left, \
517 to_make-> valid. top, \
518 to_make-> valid. width + 1, \
519 to_make-> valid. height \
520 }; \
521 if( im_prepare( make_from, &require ) ) \
522 return -1; \
523 \
524 { \
525 int *writing= (int*) IM_REGION_ADDR_TOPLEFT( to_make ); \
526 size_t write_skip= IM_REGION_LSKIP( to_make ) / sizeof(int); \
527 int *write_end= writing + write_skip * to_make-> valid. height; \
528 size_t write_width= to_make-> valid. width; \
529 int *write_stop; \
530 \
531 TYPE *reading= (TYPE*) IM_REGION_ADDR( make_from, require. left, require. top ); \
532 size_t read_skip= ( IM_REGION_LSKIP( make_from ) / sizeof(TYPE) ) - write_width; \
533 \
534 write_skip-= write_width; \
535 \
536 for( ; writing < write_end; writing+= write_skip, reading+= read_skip ) \
537 for( write_stop= writing + write_width; writing < write_stop; ++writing, ++reading ) \
538 *writing= (int)( reading[1] - reading[0] ); \
539 } \
540 return 0; \
541 }
542
543 #define YGRAD_GEN_DEFINITION( TYPE ) \
544 static int ygrad_gen_ ## TYPE( REGION *to_make, void *vptr_make_from, void *unrequired, void *unreq2 ){ \
545 \
546 REGION *make_from= (REGION*) vptr_make_from; \
547 Rect require= { \
548 to_make-> valid. left, \
549 to_make-> valid. top, \
550 to_make-> valid. width, \
551 to_make-> valid. height + 1 \
552 }; \
553 if( im_prepare( make_from, &require ) ) \
554 return -1; \
555 \
556 { \
557 int *writing= (int*) IM_REGION_ADDR_TOPLEFT( to_make ); \
558 size_t write_skip= IM_REGION_LSKIP( to_make ) / sizeof(int); \
559 int *write_end= writing + write_skip * to_make-> valid. height; \
560 size_t write_width= to_make-> valid. width; \
561 int *write_stop; \
562 \
563 TYPE *reading= (TYPE*) IM_REGION_ADDR( make_from, require. left, require. top ); \
564 size_t read_line= IM_REGION_LSKIP( make_from ) / sizeof(TYPE); \
565 size_t read_skip= read_line - write_width; \
566 \
567 write_skip-= write_width; \
568 \
569 for( ; writing < write_end; writing+= write_skip, reading+= read_skip ) \
570 for( write_stop= writing + write_width; writing < write_stop; ++writing, ++reading ) \
571 *writing= (int)( reading[ read_line ] - reading[0] ); \
572 } \
573 return 0; \
574 }
575
576 XGRAD_GEN_DEFINITION( guint8 )
577 YGRAD_GEN_DEFINITION( guint8 )
578 XGRAD_GEN_DEFINITION( gint8 )
579 YGRAD_GEN_DEFINITION( gint8 )
580 XGRAD_GEN_DEFINITION( guint16 )
581 YGRAD_GEN_DEFINITION( guint16 )
582 XGRAD_GEN_DEFINITION( gint16 )
583 YGRAD_GEN_DEFINITION( gint16 )
584 XGRAD_GEN_DEFINITION( guint32 )
585 YGRAD_GEN_DEFINITION( guint32 )
586 XGRAD_GEN_DEFINITION( gint32 )
587 YGRAD_GEN_DEFINITION( gint32 )
588 #if 0
589 XGRAD_GEN_DEFINITION( float )
590 YGRAD_GEN_DEFINITION( float )
591 XGRAD_GEN_DEFINITION( double )
592 YGRAD_GEN_DEFINITION( double )
593 #endif
594