1 /* Merge two images left-right.
2  *
3  * Copyright: 1990, 1991 N. Dessipris
4  * Author: N. Dessipris
5  * Written on: 20/09/1990
6  * Updated on: 17/04/1991
7  * 1/6/92: JC
8  * 	- check for difference bug fixed
9  *	- geometry calculations improved and simplified
10  *	- small speedups
11  Kirk Martinez for Sys5 29/4/93
12  * 7/8/93 JC
13  *	- ANSIfied
14  *	- memory leaks fixed, ready for partial v2
15  *	- now does IM_CODING_LABQ too
16  * 8/11/93 JC
17  *	- now propogates both input histories
18  *	- adds magic lines for global mosaic optimisation
19  *
20  *
21  *
22  May/1994 Ahmed Abbood
23  *
24  *	- Modified to use partials on all IO
25  *
26  June/1995 Ahmed Abbood
27  *
28  *	- Modified to work with different types of images.
29  *
30  * 16/6/95 JC
31  *	- tidied up a little
32  *	- added to VIPS!
33  * 7/9/95 JC
34  *	- split into two parts: im_lrmerge() and im__lrmerge()
35  *	- latter called by im_lrmosaic()
36  *	- just the same as public im_lrmerge(), but adds no history
37  *	- necessary for im_global_balance()
38  *	- small bugs fixed
39  * 10/10/95 JC
40  *	- better checks that parameters are sensible
41  * 11/10/95 JC
42  *	- Kirk spotted what a load of rubbish Ahmed's code is
43  *	- rewritten - many, many bugs fixed
44  * 24/1/97 JC
45  *	- now outputs bounding area of input images, rather than clipping
46  *	- ignores 0 pixels in blend
47  *	- small tidies
48  * 7/2/97 JC
49  *	- new blend, caching
50  * 25/2/97 JC
51  *	- old blend back, much simpler
52  *	- speed this up at some point if you think of an easy way to do it
53  * 29/7/97 JC
54  *	- IM_CODING_LABQ blend now works, was bug in im_wrapone()
55  *	- small tidies
56  * 10/1/98 JC
57  *	- merge LUTs now shared between all running mergers
58  *	- frees memory explicitly in im__stop_merge, for much better memory
59  *	  use in large mosaics, huge improvement!
60  * 18/2/98 JC
61  *	- im_demand_hint() call added
62  * 19/2/98 JC
63  *	- now works for any dx/dy by calling im_insert() for bizarre cases
64  * 26/9/99 JC
65  *	- ooops, blend lut was wrong! wonder how long that's been broken,
66  *	  since feb97 I guess
67  * 2/2/01 JC
68  *	- added tunable max blend width
69  * 8/3/01 JC
70  *	- switched to integer arithmetic for integer blends
71  * 7/11/01 JC
72  *	- more sophisticated transparency handling
73  *	- tiny blend speed up
74  * 19/3/02 JC
75  * 	- move fl cache to main state for better sharing
76  * 15/8/02 JC
77  *	- records Xoffset/Yoffset
78  * 20/6/05
79  *	- now requires all bands == 0 for transparency (used to just check
80  *	  band 0)
81  * 24/1/11
82  * 	- gtk-doc
83  * 	- match formats and bands automatically
84  * 22/5/14
85  * 	- wrap as a class
86  * 18/6/20 kleisauke
87  * 	- convert to vips8
88  */
89 
90 /*
91 
92     This file is part of VIPS.
93 
94     VIPS is free software; you can redistribute it and/or modify
95     it under the terms of the GNU Lesser General Public License as published by
96     the Free Software Foundation; either version 2 of the License, or
97     (at your option) any later version.
98 
99     This program is distributed in the hope that it will be useful,
100     but WITHOUT ANY WARRANTY; without even the implied warranty of
101     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
102     GNU Lesser General Public License for more details.
103 
104     You should have received a copy of the GNU Lesser General Public License
105     along with this program; if not, write to the Free Software
106     Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
107     02110-1301  USA
108 
109  */
110 
111 /*
112 
113     These files are distributed with VIPS - http://www.vips.ecs.soton.ac.uk
114 
115  */
116 
117 #ifdef HAVE_CONFIG_H
118 #include <config.h>
119 #endif /*HAVE_CONFIG_H*/
120 #include <vips/intl.h>
121 
122 #include <stdio.h>
123 #include <stdlib.h>
124 #include <string.h>
125 #include <math.h>
126 #include <limits.h>
127 
128 /* Define for debug output.
129 #define DEBUG
130  */
131 
132 #include <vips/vips.h>
133 #include <vips/thread.h>
134 #include <vips/transform.h>
135 #include <vips/internal.h>
136 
137 #include "pmosaicing.h"
138 
139 /* Blend luts. Shared between all lr and tb blends.
140  */
141 double *vips__coef1 = NULL;
142 double *vips__coef2 = NULL;
143 int *vips__icoef1 = NULL;
144 int *vips__icoef2 = NULL;
145 
146 /* Create a lut for the merging area. Always BLEND_SIZE entries, we
147  * scale later when we index it.
148  */
149 int
vips__make_blend_luts(void)150 vips__make_blend_luts( void )
151 {
152 	int x;
153 
154 	/* Already done?
155 	 */
156 	if( vips__coef1 && vips__coef2 )
157 		return( 0 );
158 
159 	/* Allocate and fill.
160 	 */
161 	vips__coef1 = VIPS_ARRAY( NULL, BLEND_SIZE, double );
162 	vips__coef2 = VIPS_ARRAY( NULL, BLEND_SIZE, double );
163 	vips__icoef1 = VIPS_ARRAY( NULL, BLEND_SIZE, int );
164 	vips__icoef2 = VIPS_ARRAY( NULL, BLEND_SIZE, int );
165 	if( !vips__coef1 || !vips__coef2 || !vips__icoef1 || !vips__icoef2 )
166 		return( -1 );
167 
168 	for( x = 0; x < BLEND_SIZE; x++ ) {
169 		double a = VIPS_PI * x / (BLEND_SIZE - 1.0);
170 
171 		vips__coef1[x] = (cos( a ) + 1.0) / 2.0;
172 		vips__coef2[x] = 1.0 - vips__coef1[x];
173 		vips__icoef1[x] = vips__coef1[x] * BLEND_SCALE;
174 		vips__icoef2[x] = vips__coef2[x] * BLEND_SCALE;
175 	}
176 
177 	return( 0 );
178 }
179 
180 /* Return the position of the first non-zero pel from the left.
181  */
182 static int
find_first(VipsRegion * ir,int * pos,int x,int y,int w)183 find_first( VipsRegion *ir, int *pos, int x, int y, int w )
184 {
185 	VipsPel *pr = VIPS_REGION_ADDR( ir, x, y );
186 	VipsImage *im = ir->im;
187 	int ne = w * im->Bands;
188 	int i;
189 
190 	/* Double the number of bands in a complex.
191 	 */
192 	if( vips_band_format_iscomplex( im->BandFmt ) )
193 		ne *= 2;
194 
195 /* Search for the first non-zero band element from the left edge of the image.
196  */
197 #define lsearch( TYPE ) { \
198 	TYPE *p = (TYPE *) pr; \
199 	\
200 	for( i = 0; i < ne; i++ ) \
201 		if( p[i] )\
202 			break;\
203 }
204 
205 	switch( im->BandFmt ) {
206 	case VIPS_FORMAT_UCHAR:		lsearch( unsigned char ); break;
207 	case VIPS_FORMAT_CHAR:		lsearch( signed char ); break;
208 	case VIPS_FORMAT_USHORT:	lsearch( unsigned short ); break;
209 	case VIPS_FORMAT_SHORT:		lsearch( signed short ); break;
210 	case VIPS_FORMAT_UINT:		lsearch( unsigned int ); break;
211 	case VIPS_FORMAT_INT:		lsearch( signed int );  break;
212 	case VIPS_FORMAT_FLOAT:		lsearch( float ); break;
213 	case VIPS_FORMAT_DOUBLE:	lsearch( double ); break;
214 	case VIPS_FORMAT_COMPLEX:	lsearch( float ); break;
215 	case VIPS_FORMAT_DPCOMPLEX:	lsearch( double ); break;
216 
217 	default:
218 		g_assert_not_reached();
219 		return( -1 );
220 	}
221 
222 	/* i is first non-zero band element, we want first non-zero pixel.
223 	 */
224 	*pos = x + i / im->Bands;
225 
226 	return( 0 );
227 }
228 
229 /* Return the position of the first non-zero pel from the right.
230  */
231 static int
find_last(VipsRegion * ir,int * pos,int x,int y,int w)232 find_last( VipsRegion *ir, int *pos, int x, int y, int w )
233 {
234 	VipsPel *pr = VIPS_REGION_ADDR( ir, x, y );
235 	VipsImage *im = ir->im;
236 	int ne = w * im->Bands;
237 	int i;
238 
239 	/* Double the number of bands in a complex.
240 	 */
241 	if( vips_band_format_iscomplex( im->BandFmt ) )
242 		ne *= 2;
243 
244 /* Search for the first non-zero band element from the right.
245  */
246 #define rsearch( TYPE ) { \
247 	TYPE *p = (TYPE *) pr; \
248 	\
249 	for( i = ne - 1; i >= 0; i-- )\
250 		if( p[i] )\
251 			break;\
252 }
253 
254 	switch( im->BandFmt ) {
255 	case VIPS_FORMAT_UCHAR:		rsearch( unsigned char ); break;
256 	case VIPS_FORMAT_CHAR:		rsearch( signed char ); break;
257 	case VIPS_FORMAT_USHORT:	rsearch( unsigned short ); break;
258 	case VIPS_FORMAT_SHORT:		rsearch( signed short ); break;
259 	case VIPS_FORMAT_UINT:		rsearch( unsigned int ); break;
260 	case VIPS_FORMAT_INT:		rsearch( signed int );  break;
261 	case VIPS_FORMAT_FLOAT:		rsearch( float ); break;
262 	case VIPS_FORMAT_DOUBLE:	rsearch( double ); break;
263 	case VIPS_FORMAT_COMPLEX:	rsearch( float ); break;
264 	case VIPS_FORMAT_DPCOMPLEX:	rsearch( double ); break;
265 
266 	default:
267 		vips_error( "lrmerge", "%s", _( "internal error" ) );
268 		return( -1 );
269 	}
270 
271 	/* i is first non-zero band element, we want first non-zero pixel.
272 	 */
273 	*pos = x + i / im->Bands;
274 
275 	return( 0 );
276 }
277 
278 /* Make sure we have first/last for this area.
279  */
280 static int
make_firstlast(MergeInfo * inf,Overlapping * ovlap,VipsRect * oreg)281 make_firstlast( MergeInfo *inf, Overlapping *ovlap, VipsRect *oreg )
282 {
283 	VipsRegion *rir = inf->rir;
284 	VipsRegion *sir = inf->sir;
285 	VipsRect rr, sr;
286 	int y, yr, ys;
287 	int missing;
288 
289 	/* We're going to build first/last ... lock it from other generate
290 	 * threads. In fact it's harmless if we do get two writers, but we may
291 	 * avoid duplicating work.
292 	 */
293 	g_mutex_lock( ovlap->fl_lock );
294 
295 	/* Do we already have first/last for this area? Bail out if we do.
296 	 */
297 	missing = 0;
298 	for( y = oreg->top; y < VIPS_RECT_BOTTOM( oreg ); y++ ) {
299 		const int j = y - ovlap->overlap.top;
300 		const int first = ovlap->first[j];
301 
302 		if( first < 0 ) {
303 			missing = 1;
304 			break;
305 		}
306 	}
307 	if( !missing ) {
308 		/* No work to do!
309 		 */
310 		g_mutex_unlock( ovlap->fl_lock );
311 		return( 0 );
312 	}
313 
314 	/* Entire width of overlap in ref for scan-lines we want.
315 	 */
316 	rr.left = ovlap->overlap.left;
317 	rr.top = oreg->top;
318 	rr.width = ovlap->overlap.width;
319 	rr.height = oreg->height;
320 	rr.left -= ovlap->rarea.left;
321 	rr.top -= ovlap->rarea.top;
322 
323 	/* Entire width of overlap in sec for scan-lines we want.
324 	 */
325 	sr.left = ovlap->overlap.left;
326 	sr.top = oreg->top;
327 	sr.width = ovlap->overlap.width;
328 	sr.height = oreg->height;
329 	sr.left -= ovlap->sarea.left;
330 	sr.top -= ovlap->sarea.top;
331 
332 #ifdef DEBUG
333 	printf( "lrmerge: making first/last for areas:\n" );
334 	printf( "ref: left = %d, top = %d, width = %d, height = %d\n",
335 		rr.left, rr.top, rr.width, rr.height );
336 	printf( "sec: left = %d, top = %d, width = %d, height = %d\n",
337 		sr.left, sr.top, sr.width, sr.height );
338 #endif
339 
340 	/* Make pixels.
341 	 */
342 	if( vips_region_prepare( rir, &rr ) ||
343 		vips_region_prepare( sir, &sr ) ) {
344 		g_mutex_unlock( ovlap->fl_lock );
345 		return( -1 );
346 	}
347 
348 	/* Make first/last cache.
349 	 */
350 	for( y = oreg->top, yr = rr.top, ys = sr.top;
351 		y < VIPS_RECT_BOTTOM( oreg ); y++, yr++, ys++ ) {
352 		const int j = y - ovlap->overlap.top;
353 		int *first = &ovlap->first[j];
354 		int *last = &ovlap->last[j];
355 
356 		/* Done this line already?
357 		 */
358 		if( *first < 0 ) {
359 			/* Search for start/end of overlap on this scan-line.
360 			 */
361 			if( find_first( sir, first,
362 				sr.left, ys, sr.width ) ||
363 				find_last( rir, last,
364 					rr.left, yr, rr.width ) ) {
365 				g_mutex_unlock( ovlap->fl_lock );
366 				return( -1 );
367 			}
368 
369 			/* Translate to output space.
370 			 */
371 			*first += ovlap->sarea.left;
372 			*last += ovlap->rarea.left;
373 
374 			/* Clip to maximum blend width, if necessary.
375 			 */
376 			if( ovlap->mwidth >= 0 &&
377 				*last - *first > ovlap->mwidth ) {
378 				int shrinkby = (*last - *first) - ovlap->mwidth;
379 
380 				*first += shrinkby / 2;
381 				*last -= shrinkby / 2;
382 			}
383 		}
384 	}
385 
386 	g_mutex_unlock( ovlap->fl_lock );
387 
388 	return( 0 );
389 }
390 
391 /* Test pixel == 0.
392  */
393 #define TEST_ZERO( TYPE, T, RESULT ) { \
394 	TYPE *tt = (T); \
395 	int ii; \
396 	\
397 	for( ii = 0; ii < cb; ii++ ) \
398 		if( tt[i + ii] ) \
399 			break; \
400 	if( ii == cb )  \
401 		(RESULT) = 1; \
402 }
403 
404 /* Blend two integer images.
405  */
406 #define iblend( TYPE, B, IN1, IN2, OUT ) { \
407 	TYPE *tr = (TYPE *) (IN1); \
408 	TYPE *ts = (TYPE *) (IN2); \
409 	TYPE *tq = (TYPE *) (OUT); \
410 	const int cb = (B); \
411 	const int left = VIPS_CLIP( 0, first - oreg->left, oreg->width ); \
412 	const int right = VIPS_CLIP( left, last - oreg->left, oreg->width ); \
413 	int ref_zero; \
414 	int sec_zero; \
415 	int x, b; \
416 	int i; \
417 	\
418 	/* Left of the blend area. \
419 	 */ \
420 	for( i = 0, x = 0; x < left; x++ ) { \
421 		ref_zero = 0; \
422 		TEST_ZERO( TYPE, tr, ref_zero ); \
423 		if( !ref_zero ) \
424 			for( b = 0; b < cb; b++, i++ ) \
425 				tq[i] = tr[i]; \
426 		else \
427 			for( b = 0; b < cb; b++, i++ ) \
428 				tq[i] = ts[i]; \
429 	} \
430 	\
431 	/* In blend area. \
432 	 */ \
433 	for( x = left; x < right; x++ ) { \
434 		ref_zero = 0; \
435 		sec_zero = 0; \
436 		TEST_ZERO( TYPE, tr, ref_zero ); \
437 		TEST_ZERO( TYPE, ts, sec_zero ); \
438 		\
439 		if( !ref_zero && !sec_zero ) { \
440 			int inx = ((x + oreg->left - first) <<  \
441 				BLEND_SHIFT) / bwidth; \
442 			int c1 = vips__icoef1[inx]; \
443 			int c2 = vips__icoef2[inx]; \
444 			\
445 			for( b = 0; b < cb; b++, i++ ) \
446 				tq[i] = c1 * tr[i] / BLEND_SCALE + \
447 					c2 * ts[i] / BLEND_SCALE; \
448 		} \
449 		else if( !ref_zero ) \
450 			for( b = 0; b < cb; b++, i++ ) \
451 				tq[i] = tr[i]; \
452 		else \
453 			for( b = 0; b < cb; b++, i++ ) \
454 				tq[i] = ts[i]; \
455 	} \
456 	\
457 	/* Right of blend.
458 	 */ \
459 	for( x = right; x < oreg->width; x++ ) { \
460 		sec_zero = 0; \
461 		TEST_ZERO( TYPE, ts, sec_zero ); \
462 		if( !sec_zero ) \
463 			for( b = 0; b < cb; b++, i++ )  \
464 				tq[i] = ts[i]; \
465 		else \
466 			for( b = 0; b < cb; b++, i++ )  \
467 				tq[i] = tr[i]; \
468 	} \
469 }
470 
471 /* Blend two float images.
472  */
473 #define fblend( TYPE, B, IN1, IN2, OUT ) { \
474 	TYPE *tr = (TYPE *) (IN1); \
475 	TYPE *ts = (TYPE *) (IN2); \
476 	TYPE *tq = (TYPE *) (OUT); \
477 	const int cb = (B); \
478 	const int left = VIPS_CLIP( 0, first - oreg->left, oreg->width ); \
479 	const int right = VIPS_CLIP( left, last - oreg->left, oreg->width ); \
480 	int ref_zero; \
481 	int sec_zero; \
482 	int x, b; \
483 	int i; \
484 	\
485 	/* Left of the blend area. \
486 	 */ \
487 	for( i = 0, x = 0; x < left; x++ ) { \
488 		ref_zero = 0; \
489 		TEST_ZERO( TYPE, tr, ref_zero ); \
490 		if( !ref_zero ) \
491 			for( b = 0; b < cb; b++, i++ ) \
492 				tq[i] = tr[i]; \
493 		else \
494 			for( b = 0; b < cb; b++, i++ ) \
495 				tq[i] = ts[i]; \
496 	} \
497 	\
498 	/* In blend area. \
499 	 */ \
500 	for( x = left; x < right; x++ ) { \
501 		ref_zero = 0; \
502 		sec_zero = 0; \
503 		TEST_ZERO( TYPE, tr, ref_zero ); \
504 		TEST_ZERO( TYPE, ts, sec_zero ); \
505 		\
506 		if( !ref_zero && !sec_zero ) { \
507 			int inx = ((x + oreg->left - first) <<  \
508 				BLEND_SHIFT) / bwidth; \
509 			double c1 = vips__coef1[inx];  \
510 			double c2 = vips__coef2[inx];  \
511 			\
512 			for( b = 0; b < cb; b++, i++ ) \
513 				tq[i] = c1 * tr[i] + c2 * ts[i]; \
514 		} \
515 		else if( !ref_zero ) \
516 			for( b = 0; b < cb; b++, i++ ) \
517 				tq[i] = tr[i]; \
518 		else \
519 			for( b = 0; b < cb; b++, i++ ) \
520 				tq[i] = ts[i]; \
521 	} \
522 	\
523 	/* Right of blend.
524 	 */ \
525 	for( x = right; x < oreg->width; x++ ) { \
526 		sec_zero = 0; \
527 		TEST_ZERO( TYPE, ts, sec_zero ); \
528 		if( !sec_zero ) \
529 			for( b = 0; b < cb; b++, i++ )  \
530 				tq[i] = ts[i]; \
531 		else \
532 			for( b = 0; b < cb; b++, i++ )  \
533 				tq[i] = tr[i]; \
534 	} \
535 }
536 
537 /* Left-right blend function for non-labpack images.
538  */
539 static int
lr_blend(VipsRegion * or,MergeInfo * inf,Overlapping * ovlap,VipsRect * oreg)540 lr_blend( VipsRegion *or, MergeInfo *inf, Overlapping *ovlap, VipsRect *oreg )
541 {
542 	VipsRegion *rir = inf->rir;
543 	VipsRegion *sir = inf->sir;
544 	VipsImage *im = or->im;
545 
546 	VipsRect prr, psr;
547 	int y, yr, ys;
548 
549 	/* Make sure we have a complete first/last set for this area.
550 	 */
551 	if( make_firstlast( inf, ovlap, oreg ) )
552 		return( -1 );
553 
554 	/* Part of rr which we will output.
555 	 */
556 	prr = *oreg;
557 	prr.left -= ovlap->rarea.left;
558 	prr.top -= ovlap->rarea.top;
559 
560 	/* Part of sr which we will output.
561 	 */
562 	psr = *oreg;
563 	psr.left -= ovlap->sarea.left;
564 	psr.top -= ovlap->sarea.top;
565 
566 	/* Make pixels.
567 	 */
568 	if( vips_region_prepare( rir, &prr ) ||
569 		vips_region_prepare( sir, &psr ) )
570 		return( -1 );
571 
572 	/* Loop down overlap area.
573 	 */
574 	for( y = oreg->top, yr = prr.top, ys = psr.top;
575 		y < VIPS_RECT_BOTTOM( oreg ); y++, yr++, ys++ ) {
576 		VipsPel *pr = VIPS_REGION_ADDR( rir, prr.left, yr );
577 		VipsPel *ps = VIPS_REGION_ADDR( sir, psr.left, ys );
578 		VipsPel *q = VIPS_REGION_ADDR( or, oreg->left, y );
579 
580 		const int j = y - ovlap->overlap.top;
581 		const int first = ovlap->first[j];
582 		const int last = ovlap->last[j];
583 		const int bwidth = last - first;
584 
585 		switch( im->BandFmt ) {
586 		case VIPS_FORMAT_UCHAR:
587 			iblend( unsigned char, im->Bands, pr, ps, q ); break;
588 		case VIPS_FORMAT_CHAR:
589 			iblend( signed char, im->Bands, pr, ps, q ); break;
590 		case VIPS_FORMAT_USHORT:
591 			iblend( unsigned short, im->Bands, pr, ps, q ); break;
592 		case VIPS_FORMAT_SHORT:
593 			iblend( signed short, im->Bands, pr, ps, q ); break;
594 		case VIPS_FORMAT_UINT:
595 			iblend( unsigned int, im->Bands, pr, ps, q ); break;
596 		case VIPS_FORMAT_INT:
597 			iblend( signed int, im->Bands, pr, ps, q );  break;
598 		case VIPS_FORMAT_FLOAT:
599 			fblend( float, im->Bands, pr, ps, q ); break;
600 		case VIPS_FORMAT_DOUBLE:
601 			fblend( double, im->Bands, pr, ps, q ); break;
602 		case VIPS_FORMAT_COMPLEX:
603 			fblend( float, im->Bands * 2, pr, ps, q ); break;
604 		case VIPS_FORMAT_DPCOMPLEX:
605 			fblend( double, im->Bands * 2, pr, ps, q ); break;
606 
607 		default:
608 			g_assert_not_reached();
609 			return( -1 );
610 		}
611 	}
612 
613 	return( 0 );
614 }
615 
616 /* Left-right blend function for VIPS_CODING_LABQ images.
617  */
618 static int
lr_blend_labpack(VipsRegion * or,MergeInfo * inf,Overlapping * ovlap,VipsRect * oreg)619 lr_blend_labpack( VipsRegion *or, MergeInfo *inf, Overlapping *ovlap,
620 	VipsRect *oreg )
621 {
622 	VipsRegion *rir = inf->rir;
623 	VipsRegion *sir = inf->sir;
624 	VipsRect prr, psr;
625 	int y, yr, ys;
626 
627 	/* Make sure we have a complete first/last set for this area. This
628 	 * will just look at the top 8 bits of L, not all 10, but should be OK.
629 	 */
630 	if( make_firstlast( inf, ovlap, oreg ) )
631 		return( -1 );
632 
633 	/* Part of rr which we will output.
634 	 */
635 	prr = *oreg;
636 	prr.left -= ovlap->rarea.left;
637 	prr.top -= ovlap->rarea.top;
638 
639 	/* Part of sr which we will output.
640 	 */
641 	psr = *oreg;
642 	psr.left -= ovlap->sarea.left;
643 	psr.top -= ovlap->sarea.top;
644 
645 	/* Make pixels.
646 	 */
647 	if( vips_region_prepare( rir, &prr ) ||
648 		vips_region_prepare( sir, &psr ) )
649 		return( -1 );
650 
651 	/* Loop down overlap area.
652 	 */
653 	for( y = oreg->top, yr = prr.top, ys = psr.top;
654 		y < VIPS_RECT_BOTTOM( oreg ); y++, yr++, ys++ ) {
655 		VipsPel *pr = VIPS_REGION_ADDR( rir, prr.left, yr );
656 		VipsPel *ps = VIPS_REGION_ADDR( sir, psr.left, ys );
657 		VipsPel *q = VIPS_REGION_ADDR( or, oreg->left, y );
658 
659 		const int j = y - ovlap->overlap.top;
660 		const int first = ovlap->first[j];
661 		const int last = ovlap->last[j];
662 		const int bwidth = last - first;
663 
664 		float *fq = inf->merge;
665 		float *r = inf->from1;
666 		float *s = inf->from2;
667 
668 		/* Unpack two bits we want.
669 		 */
670 		vips__LabQ2Lab_vec( r, pr, oreg->width );
671 		vips__LabQ2Lab_vec( s, ps, oreg->width );
672 
673 		/* Blend as floats.
674 		 */
675 		fblend( float, 3, r, s, fq );
676 
677 		/* Re-pack to output buffer.
678 		 */
679 		vips__Lab2LabQ_vec( q, inf->merge, oreg->width );
680 	}
681 
682 	return( 0 );
683 }
684 
685 static void
lock_free(VipsImage * image,GMutex * lock)686 lock_free( VipsImage *image, GMutex *lock )
687 {
688 	VIPS_FREEF( vips_g_mutex_free, lock );
689 }
690 
691 /* Build basic per-call state and do some geometry calculations. Shared with
692  * tbmerge, so not static.
693  */
694 Overlapping *
vips__build_mergestate(const char * domain,VipsImage * ref,VipsImage * sec,VipsImage * out,int dx,int dy,int mwidth)695 vips__build_mergestate( const char *domain,
696 	VipsImage *ref, VipsImage *sec, VipsImage *out,
697 	int dx, int dy, int mwidth )
698 {
699 	VipsImage **t = (VipsImage **)
700 		vips_object_local_array( VIPS_OBJECT( out ), 4 );
701 
702 	VipsImage **arry;
703    	Overlapping *ovlap;
704 	int x;
705 
706 	/* TODO(kleisauke): Copied from vips_insert, perhaps we
707 	 * need a separate function for this?
708 	 * (just like im__insert_base)
709 	 */
710 	if( vips_image_pio_input( ref ) ||
711 		vips_image_pio_input( sec ) ||
712 		vips_check_bands_1orn( domain, ref, sec ) ||
713 		vips_check_coding_known( domain, ref ) ||
714 		vips_check_coding_same( domain, ref, sec ) )
715 		return( NULL );
716 
717 	/* Cast our input images up to a common format and bands.
718 	 */
719 	if( vips__formatalike( ref, sec, &t[0], &t[1] ) ||
720 		vips__bandalike( domain, t[0], t[1], &t[2], &t[3] ) )
721 		return( NULL );
722 
723 	if( !(arry = vips_allocate_input_array( out,
724 		t[2], t[3], NULL )) )
725 		return( NULL );
726 
727 	if( vips_image_pipeline_array( out,
728 		VIPS_DEMAND_STYLE_SMALLTILE, arry ) )
729 		return( NULL );
730 
731 	if( mwidth < -1 ) {
732 		vips_error( domain, "%s", _( "mwidth must be -1 or >= 0" ) );
733 		return( NULL );
734 	}
735 
736 	if( !(ovlap = VIPS_NEW( out, Overlapping )) )
737 		return( NULL );
738 
739 	ovlap->ref = arry[0];
740 	ovlap->sec = arry[1];
741 	ovlap->out = out;
742 	ovlap->dx = dx;
743 	ovlap->dy = dy;
744 	ovlap->mwidth = mwidth;
745 
746 	/* Area occupied by ref image. Place at (0,0) to start with.
747 	 */
748    	ovlap->rarea.left = 0;
749    	ovlap->rarea.top = 0;
750    	ovlap->rarea.width = ovlap->ref->Xsize;
751    	ovlap->rarea.height = ovlap->ref->Ysize;
752 
753 	/* Area occupied by sec image.
754 	 */
755    	ovlap->sarea.left = -dx;
756    	ovlap->sarea.top = -dy;
757    	ovlap->sarea.width = ovlap->sec->Xsize;
758    	ovlap->sarea.height = ovlap->sec->Ysize;
759 
760 	/* Compute overlap.
761 	 */
762 	vips_rect_intersectrect( &ovlap->rarea, &ovlap->sarea, &ovlap->overlap );
763 	if( vips_rect_isempty( &ovlap->overlap ) ) {
764 		vips_error( domain, "%s", _( "no overlap" ) );
765 		return( NULL );
766 	}
767 
768 	/* Find position and size of output image.
769 	 */
770 	vips_rect_unionrect( &ovlap->rarea, &ovlap->sarea, &ovlap->oarea );
771 
772 	/* Now: translate everything, so that the output image, not the left
773 	 * image, is at (0,0).
774 	 */
775 	ovlap->rarea.left -= ovlap->oarea.left;
776 	ovlap->rarea.top -= ovlap->oarea.top;
777 	ovlap->sarea.left -= ovlap->oarea.left;
778 	ovlap->sarea.top -= ovlap->oarea.top;
779 	ovlap->overlap.left -= ovlap->oarea.left;
780 	ovlap->overlap.top -= ovlap->oarea.top;
781 	ovlap->oarea.left = 0;
782 	ovlap->oarea.top = 0;
783 
784 	/* Make sure blend luts are built.
785 	 */
786 	vips__make_blend_luts();
787 
788 	/* Size of first/last cache. Could be either of these ... just pick
789 	 * the larger.
790 	 */
791 	ovlap->flsize = VIPS_MAX( ovlap->overlap.width, ovlap->overlap.height );
792 
793 	/* Build first/last cache.
794 	 */
795 	ovlap->first = VIPS_ARRAY( out, ovlap->flsize, int );
796 	ovlap->last = VIPS_ARRAY( out, ovlap->flsize, int );
797 	if( !ovlap->first || !ovlap->last )
798 		return( NULL );
799 	for( x = 0; x < ovlap->flsize; x++ )
800 		ovlap->first[x] = -1;
801 
802 	ovlap->fl_lock = vips_g_mutex_new();
803 
804 	g_signal_connect( out, "close",
805 		G_CALLBACK( lock_free ), ovlap->fl_lock );
806 
807 	return( ovlap );
808 }
809 
810 /* Build per-call state.
811  */
812 static Overlapping *
build_lrstate(VipsImage * ref,VipsImage * sec,VipsImage * out,int dx,int dy,int mwidth)813 build_lrstate( VipsImage *ref, VipsImage *sec, VipsImage *out,
814 	int dx, int dy, int mwidth )
815 {
816    	Overlapping *ovlap;
817 
818 	if( !(ovlap = vips__build_mergestate( "lrmerge",
819 		ref, sec, out, dx, dy, mwidth )) )
820 		return( NULL );
821 
822 	/* Select blender.
823 	 */
824 	switch( ovlap->ref->Coding ) {
825 	case VIPS_CODING_LABQ:
826 		ovlap->blend = lr_blend_labpack;
827 		break;
828 
829 	case VIPS_CODING_NONE:
830 		ovlap->blend = lr_blend;
831 		break;
832 
833 	default:
834 		vips_error( "lrmerge", "%s", _( "unknown coding type" ) );
835 		return( NULL );
836 	}
837 
838 	/* Find the parts of output which come just from ref and just from sec.
839 	 */
840 	ovlap->rpart = ovlap->rarea;
841 	ovlap->spart = ovlap->sarea;
842 	ovlap->rpart.width -= ovlap->overlap.width;
843 	ovlap->spart.left += ovlap->overlap.width;
844 	ovlap->spart.width -= ovlap->overlap.width;
845 
846 	/* Is there too much overlap? ie. right edge of ref image is greater
847 	 * than right edge of sec image, or left > left.
848 	 */
849 	if( VIPS_RECT_RIGHT( &ovlap->rarea ) >
850 		VIPS_RECT_RIGHT( &ovlap->sarea ) ||
851 		ovlap->rarea.left > ovlap->sarea.left ) {
852 		vips_error( "lrmerge", "%s", _( "too much overlap" ) );
853 		return( NULL );
854 	}
855 
856 	/* Max number of pixels we may have to blend over.
857 	 */
858 	ovlap->blsize = ovlap->overlap.width;
859 
860 	return( ovlap );
861 }
862 
863 /* The area being demanded can be filled using only pels from either the ref
864  * or the sec images. Attach output to the appropriate part of the input image.
865  * area is the position that ir->im occupies in the output image.
866  *
867  * Shared with tbmerge, so not static.
868  */
869 int
vips__attach_input(VipsRegion * or,VipsRegion * ir,VipsRect * area)870 vips__attach_input( VipsRegion *or, VipsRegion *ir, VipsRect *area )
871 {
872 	VipsRect r = or->valid;
873 
874 	/* Translate to source coordinate space.
875 	 */
876 	r.left -= area->left;
877 	r.top -= area->top;
878 
879 	/* Demand input.
880 	 */
881 	if( vips_region_prepare( ir, &r ) )
882 		return( -1 );
883 
884 	/* Attach or to ir.
885 	 */
886 	if( vips_region_region( or, ir, &or->valid, r.left, r.top ) )
887 		 return( -1 );
888 
889 	return( 0 );
890 }
891 
892 /* The area being demanded requires pixels from the ref and sec images. As
893  * above, but just do a sub-area of the output, and make sure we copy rather
894  * than just pointer-fiddling. reg is the sub-area of or->valid we should do.
895  *
896  * Shared with tbmerge, so not static.
897  */
898 int
vips__copy_input(VipsRegion * or,VipsRegion * ir,VipsRect * area,VipsRect * reg)899 vips__copy_input( VipsRegion *or, VipsRegion *ir,
900 	VipsRect *area, VipsRect *reg )
901 {
902 	VipsRect r = *reg;
903 
904 	/* Translate to source coordinate space.
905 	 */
906 	r.left -= area->left;
907 	r.top -= area->top;
908 
909 	/* Paint this area of ir into or.
910 	 */
911 	if( vips_region_prepare_to( ir, or, &r, reg->left, reg->top ) )
912 		return( -1 );
913 
914 	return( 0 );
915 }
916 
917 /* Generate function for merge. This is shared between lrmerge and
918  * tbmerge.
919  */
920 int
vips__merge_gen(VipsRegion * or,void * seq,void * a,void * b,gboolean * stop)921 vips__merge_gen( VipsRegion *or, void *seq, void *a, void *b,
922 	gboolean *stop )
923 {
924 	MergeInfo *inf = (MergeInfo *) seq;
925 	Overlapping *ovlap = (Overlapping *) a;
926 	VipsRect *r = &or->valid;
927 	VipsRect rreg, sreg, oreg;
928 
929 	/* Find intersection with overlap, ref and sec parts.
930 	 */
931 	vips_rect_intersectrect( r, &ovlap->rpart, &rreg );
932 	vips_rect_intersectrect( r, &ovlap->spart, &sreg );
933 
934 	/* Do easy cases first: can we satisfy this demand with pixels just
935 	 * from ref, or just from sec.
936 	 */
937 	if( vips_rect_equalsrect( r, &rreg ) ) {
938 		if( vips__attach_input( or, inf->rir, &ovlap->rarea ) )
939 			return( -1 );
940 	}
941 	else if( vips_rect_equalsrect( r, &sreg ) ) {
942 		if( vips__attach_input( or, inf->sir, &ovlap->sarea ) )
943 			return( -1 );
944 	}
945 	else {
946 		/* Difficult case - do in three stages: black out whole area,
947 		 * copy in parts of ref and sec we touch, write blend area.
948 		 * This could be sped up somewhat ... we will usually black
949 		 * out far too much, and write to the blend area three times.
950 		 * Upgrade in the future!
951 		 */
952 
953 		/* Need intersections with whole of left & right, and overlap
954 		 * too.
955 		 */
956 		vips_rect_intersectrect( r, &ovlap->rarea, &rreg );
957 		vips_rect_intersectrect( r, &ovlap->sarea, &sreg );
958 		vips_rect_intersectrect( r, &ovlap->overlap, &oreg );
959 
960 		vips_region_black( or );
961 		if( !vips_rect_isempty( &rreg ) )
962 			if( vips__copy_input( or,
963 				inf->rir, &ovlap->rarea, &rreg ) )
964 				return( -1 );
965 		if( !vips_rect_isempty( &sreg ) )
966 			if( vips__copy_input( or,
967 				inf->sir, &ovlap->sarea, &sreg ) )
968 				return( -1 );
969 
970 		/* Nasty: inf->rir and inf->sir now point to the same bit of
971 		 * memory (part of or), and we've written twice. We need to
972 		 * make sure we get fresh pixels for the blend, so we must
973 		 * invalidate them both. Should maybe add a call to the API
974 		 * for this.
975 		 */
976 		inf->rir->valid.width = inf->sir->valid.width = 0;
977 
978 		/* Now blat in the blended area.
979 		 */
980 		if( !vips_rect_isempty( &oreg ) )
981 			if( ovlap->blend( or, inf, ovlap, &oreg ) )
982 				return( -1 );
983 	}
984 
985 	return( 0 );
986 }
987 
988 /* Stop function. Shared with tbmerge. Free explicitly to reduce mem
989  * requirements quickly for large mosaics.
990  */
991 int
vips__stop_merge(void * seq,void * a,void * b)992 vips__stop_merge( void *seq, void *a, void *b )
993 {
994 	MergeInfo *inf = (MergeInfo *) seq;
995 
996 	VIPS_UNREF( inf->rir );
997 	VIPS_UNREF( inf->sir );
998 	VIPS_FREE( inf->from1 );
999 	VIPS_FREE( inf->from2 );
1000 	VIPS_FREE( inf->merge );
1001 	g_free( inf );
1002 
1003 	return( 0 );
1004 }
1005 
1006 /* Start function. Shared with tbmerge.
1007  */
1008 void *
vips__start_merge(VipsImage * out,void * a,void * b)1009 vips__start_merge( VipsImage *out, void *a, void *b )
1010 {
1011 	Overlapping *ovlap = (Overlapping *) a;
1012 	MergeInfo *inf;
1013 
1014 	if( !(inf = VIPS_NEW( NULL, MergeInfo )) )
1015 		return( NULL );
1016 
1017 	inf->rir = NULL;
1018 	inf->sir = NULL;
1019 	inf->from1 = NULL;
1020 	inf->from2 = NULL;
1021 	inf->merge = NULL;
1022 
1023 	/* If this is going to be a VIPS_CODING_LABQ, we need VIPS_CODING_LABQ
1024 	 * blend buffers.
1025 	 */
1026 	if( out->Coding == VIPS_CODING_LABQ ) {
1027 		inf->from1 = VIPS_ARRAY( NULL, ovlap->blsize * 3, float );
1028 		inf->from2 = VIPS_ARRAY( NULL, ovlap->blsize * 3, float );
1029 		inf->merge = VIPS_ARRAY( NULL, ovlap->blsize * 3, float );
1030 		if( !inf->from1 || !inf->from2 || !inf->merge ) {
1031 			vips__stop_merge( inf, NULL, NULL );
1032 			return( NULL );
1033 		}
1034 	}
1035 
1036 	inf->rir = vips_region_new( ovlap->ref );
1037 	inf->sir = vips_region_new( ovlap->sec );
1038 
1039 	if( !inf->rir || !inf->sir ) {
1040 		vips__stop_merge( inf, NULL, NULL );
1041 		return( NULL );
1042 	}
1043 
1044 	return( inf );
1045 }
1046 
1047 int
vips__lrmerge(VipsImage * ref,VipsImage * sec,VipsImage * out,int dx,int dy,int mwidth)1048 vips__lrmerge( VipsImage *ref, VipsImage *sec, VipsImage *out,
1049 	int dx, int dy, int mwidth )
1050 {
1051 	Overlapping *ovlap;
1052 
1053 #ifdef DEBUG
1054 	printf( "lrmerge %s %s %s %d %d %d\n",
1055 		ref->filename, sec->filename, out->filename,
1056 		dx, dy, mwidth );
1057 	printf( "ref is %d x %d pixels\n", ref->Xsize, ref->Ysize );
1058 	printf( "sec is %d x %d pixels\n", sec->Xsize, sec->Ysize );
1059 #endif
1060 
1061 	if( dx > 0 || dx < 1 - ref->Xsize ) {
1062 		VipsImage *x;
1063 
1064 #ifdef DEBUG
1065 		printf( "lrmerge: no overlap, using insert\n" );
1066 #endif
1067 
1068 		/* No overlap, use insert instead.
1069 		 */
1070   		if( vips_insert( ref, sec, &x, -dx, -dy,
1071   			"expand", TRUE,
1072 			NULL ) )
1073 			return( -1 );
1074 		if( vips_image_write( x, out ) ) {
1075 			g_object_unref( x );
1076 			return( -1 );
1077 		}
1078 		g_object_unref( x );
1079 
1080 		out->Xoffset = -dx;
1081 		out->Yoffset = -dy;
1082 
1083 		return( 0 );
1084 	}
1085 
1086 	if( !(ovlap = build_lrstate( ref, sec, out, dx, dy, mwidth )) )
1087 		return( -1 );
1088 
1089 	if( vips_image_pipelinev( out,
1090 		VIPS_DEMAND_STYLE_THINSTRIP, ovlap->ref, ovlap->sec, NULL ) )
1091 		return( -1 );
1092 
1093 	out->Xsize = ovlap->oarea.width;
1094 	out->Ysize = ovlap->oarea.height;
1095 	out->Xoffset = -dx;
1096 	out->Yoffset = -dy;
1097 
1098 	if( vips_image_generate( out,
1099 		vips__start_merge, vips__merge_gen, vips__stop_merge,
1100 		ovlap, NULL ) )
1101 		return( -1 );
1102 
1103 	return ( 0 );
1104 }
1105 
1106 const char *
vips__get_mosaic_name(VipsImage * image)1107 vips__get_mosaic_name( VipsImage *image )
1108 {
1109 	const char *name;
1110 
1111 	if( vips_image_get_typeof( image, "mosaic-name" ) ) {
1112 		if( vips_image_get_string( image, "mosaic-name", &name ) )
1113 			return( NULL );
1114 	}
1115 	else
1116 		name = image->filename;
1117 
1118 	return( name );
1119 }
1120 
1121 void
vips__add_mosaic_name(VipsImage * image)1122 vips__add_mosaic_name( VipsImage *image )
1123 {
1124 	static int global_serial = 0;
1125 
1126 	/* TODO(kleisauke): Could we call vips_image_temp_name instead?
1127 	 */
1128 	int serial = g_atomic_int_add( &global_serial, 1 );
1129 
1130 	char name[256];
1131 
1132 	/* We must override any inherited name, so don't test for doesn't
1133 	 * exist before setting.
1134 	 */
1135 	vips_snprintf( name, 256, "mosaic-temp-%d", serial );
1136 	vips_image_set_string( image, "mosaic-name", name );
1137 }
1138 
1139