1/*
2 * Copyright (C) 2011 The Android Open Source Project
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 *      http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16
17///////////////////////////////////////////////////
18// Blend.cpp
19// $Id: Blend.cpp,v 1.22 2011/06/24 04:22:14 mbansal Exp $
20
21#include <string.h>
22
23#include "Interp.h"
24#include "Blend.h"
25
26#include "Geometry.h"
27#include "trsMatrix.h"
28
29#include "Log.h"
30#define LOG_TAG "BLEND"
31
32Blend::Blend()
33{
34  m_wb.blendingType = BLEND_TYPE_NONE;
35}
36
37Blend::~Blend()
38{
39    if (m_pFrameVPyr) free(m_pFrameVPyr);
40    if (m_pFrameUPyr) free(m_pFrameUPyr);
41    if (m_pFrameYPyr) free(m_pFrameYPyr);
42}
43
44int Blend::initialize(int blendingType, int stripType, int frame_width, int frame_height)
45{
46    this->width = frame_width;
47    this->height = frame_height;
48    this->m_wb.blendingType = blendingType;
49    this->m_wb.stripType = stripType;
50
51    m_wb.blendRange = m_wb.blendRangeUV = BLEND_RANGE_DEFAULT;
52    m_wb.nlevs = m_wb.blendRange;
53    m_wb.nlevsC = m_wb.blendRangeUV;
54
55    if (m_wb.nlevs <= 0) m_wb.nlevs = 1; // Need levels for YUV processing
56    if (m_wb.nlevsC > m_wb.nlevs) m_wb.nlevsC = m_wb.nlevs;
57
58    m_wb.roundoffOverlap = 1.5;
59
60    m_pFrameYPyr = NULL;
61    m_pFrameUPyr = NULL;
62    m_pFrameVPyr = NULL;
63
64    m_pFrameYPyr = PyramidShort::allocatePyramidPacked(m_wb.nlevs, (unsigned short) width, (unsigned short) height, BORDER);
65    m_pFrameUPyr = PyramidShort::allocatePyramidPacked(m_wb.nlevsC, (unsigned short) (width), (unsigned short) (height), BORDER);
66    m_pFrameVPyr = PyramidShort::allocatePyramidPacked(m_wb.nlevsC, (unsigned short) (width), (unsigned short) (height), BORDER);
67
68    if (!m_pFrameYPyr || !m_pFrameUPyr || !m_pFrameVPyr)
69    {
70        LOGE("Error: Could not allocate pyramids for blending");
71        return BLEND_RET_ERROR_MEMORY;
72    }
73
74    return BLEND_RET_OK;
75}
76
77inline double max(double a, double b) { return a > b ? a : b; }
78inline double min(double a, double b) { return a < b ? a : b; }
79
80void Blend::AlignToMiddleFrame(MosaicFrame **frames, int frames_size)
81{
82    // Unwarp this frame and Warp the others to match
83    MosaicFrame *mb = NULL;
84    MosaicFrame *ref = frames[int(frames_size/2)];    // Middle frame
85
86    double invtrs[3][3];
87    inv33d(ref->trs, invtrs);
88
89    for(int mfit = 0; mfit < frames_size; mfit++)
90    {
91        mb = frames[mfit];
92        double temp[3][3];
93        mult33d(temp, invtrs, mb->trs);
94        memcpy(mb->trs, temp, sizeof(temp));
95        normProjMat33d(mb->trs);
96    }
97}
98
99int Blend::runBlend(MosaicFrame **oframes, MosaicFrame **rframes,
100        int frames_size,
101        ImageType &imageMosaicYVU, int &mosaicWidth, int &mosaicHeight,
102        float &progress, bool &cancelComputation)
103{
104    int ret;
105    int numCenters;
106
107    MosaicFrame **frames;
108
109    // For THIN strip mode, accept all frames for blending
110    if (m_wb.stripType == STRIP_TYPE_THIN)
111    {
112        frames = oframes;
113    }
114    else // For WIDE strip mode, first select the relevant frames to blend.
115    {
116        SelectRelevantFrames(oframes, frames_size, rframes, frames_size);
117        frames = rframes;
118    }
119
120    ComputeBlendParameters(frames, frames_size, true);
121    numCenters = frames_size;
122
123    if (numCenters == 0)
124    {
125        LOGE("Error: No frames to blend");
126        return BLEND_RET_ERROR;
127    }
128
129    if (!(m_AllSites = m_Triangulator.allocMemory(numCenters)))
130    {
131        return BLEND_RET_ERROR_MEMORY;
132    }
133
134    // Bounding rectangle (real numbers) of the final mosaic computed by projecting
135    // each input frame into the mosaic coordinate system.
136    BlendRect global_rect;
137
138    global_rect.lft = global_rect.bot = 2e30; // min values
139    global_rect.rgt = global_rect.top = -2e30; // max values
140    MosaicFrame *mb = NULL;
141    double halfwidth = width / 2.0;
142    double halfheight = height / 2.0;
143
144    double z, x0, y0, x1, y1, x2, y2, x3, y3;
145
146    // Corners of the left-most and right-most frames respectively in the
147    // mosaic coordinate system.
148    double xLeftCorners[2] = {2e30, 2e30};
149    double xRightCorners[2] = {-2e30, -2e30};
150
151    // Corners of the top-most and bottom-most frames respectively in the
152    // mosaic coordinate system.
153    double yTopCorners[2] = {2e30, 2e30};
154    double yBottomCorners[2] = {-2e30, -2e30};
155
156
157    // Determine the extents of the final mosaic
158    CSite *csite = m_AllSites ;
159    for(int mfit = 0; mfit < frames_size; mfit++)
160    {
161        mb = frames[mfit];
162
163        // Compute clipping for this frame's rect
164        FrameToMosaicRect(mb->width, mb->height, mb->trs, mb->brect);
165        // Clip global rect using this frame's rect
166        ClipRect(mb->brect, global_rect);
167
168        // Calculate the corner points
169        FrameToMosaic(mb->trs, 0.0,             0.0,            x0, y0);
170        FrameToMosaic(mb->trs, 0.0,             mb->height-1.0, x1, y1);
171        FrameToMosaic(mb->trs, mb->width-1.0,   mb->height-1.0, x2, y2);
172        FrameToMosaic(mb->trs, mb->width-1.0,   0.0,            x3, y3);
173
174        if(x0 < xLeftCorners[0] || x1 < xLeftCorners[1])    // If either of the left corners is lower
175        {
176            xLeftCorners[0] = x0;
177            xLeftCorners[1] = x1;
178        }
179
180        if(x3 > xRightCorners[0] || x2 > xRightCorners[1])    // If either of the right corners is higher
181        {
182            xRightCorners[0] = x3;
183            xRightCorners[1] = x2;
184        }
185
186        if(y0 < yTopCorners[0] || y3 < yTopCorners[1])    // If either of the top corners is lower
187        {
188            yTopCorners[0] = y0;
189            yTopCorners[1] = y3;
190        }
191
192        if(y1 > yBottomCorners[0] || y2 > yBottomCorners[1])    // If either of the bottom corners is higher
193        {
194            yBottomCorners[0] = y1;
195            yBottomCorners[1] = y2;
196        }
197
198
199        // Compute the centroid of the warped region
200        FindQuadCentroid(x0, y0, x1, y1, x2, y2, x3, y3, csite->getVCenter().x, csite->getVCenter().y);
201
202        csite->setMb(mb);
203        csite++;
204    }
205
206    // Get origin and sizes
207
208    // Bounding rectangle (int numbers) of the final mosaic computed by projecting
209    // each input frame into the mosaic coordinate system.
210    MosaicRect fullRect;
211
212    fullRect.left = (int) floor(global_rect.lft); // min-x
213    fullRect.top = (int) floor(global_rect.bot);  // min-y
214    fullRect.right = (int) ceil(global_rect.rgt); // max-x
215    fullRect.bottom = (int) ceil(global_rect.top);// max-y
216    Mwidth = (unsigned short) (fullRect.right - fullRect.left + 1);
217    Mheight = (unsigned short) (fullRect.bottom - fullRect.top + 1);
218
219    int xLeftMost, xRightMost;
220    int yTopMost, yBottomMost;
221
222    // Rounding up, so that we don't include the gray border.
223    xLeftMost = max(0, max(xLeftCorners[0], xLeftCorners[1]) - fullRect.left + 1);
224    xRightMost = min(Mwidth - 1, min(xRightCorners[0], xRightCorners[1]) - fullRect.left - 1);
225
226    yTopMost = max(0, max(yTopCorners[0], yTopCorners[1]) - fullRect.top + 1);
227    yBottomMost = min(Mheight - 1, min(yBottomCorners[0], yBottomCorners[1]) - fullRect.top - 1);
228
229    if (xRightMost <= xLeftMost || yBottomMost <= yTopMost)
230    {
231        LOGE("RunBlend: aborting -consistency check failed,"
232             "(xLeftMost, xRightMost, yTopMost, yBottomMost): (%d, %d, %d, %d)",
233             xLeftMost, xRightMost, yTopMost, yBottomMost);
234        return BLEND_RET_ERROR;
235    }
236
237    // Make sure image width is multiple of 4
238    Mwidth = (unsigned short) ((Mwidth + 3) & ~3);
239    Mheight = (unsigned short) ((Mheight + 3) & ~3);    // Round up.
240
241    ret = MosaicSizeCheck(LIMIT_SIZE_MULTIPLIER, LIMIT_HEIGHT_MULTIPLIER);
242    if (ret != BLEND_RET_OK)
243    {
244       LOGE("RunBlend: aborting - mosaic size check failed, "
245            "(frame_width, frame_height) vs (mosaic_width, mosaic_height): "
246            "(%d, %d) vs (%d, %d)", width, height, Mwidth, Mheight);
247       return ret;
248    }
249
250    LOGI("Allocate mosaic image for blending - size: %d x %d", Mwidth, Mheight);
251    YUVinfo *imgMos = YUVinfo::allocateImage(Mwidth, Mheight);
252    if (imgMos == NULL)
253    {
254        LOGE("RunBlend: aborting - couldn't alloc %d x %d mosaic image", Mwidth, Mheight);
255        return BLEND_RET_ERROR_MEMORY;
256    }
257
258    // Set the Y image to 255 so we can distinguish when frame idx are written to it
259    memset(imgMos->Y.ptr[0], 255, (imgMos->Y.width * imgMos->Y.height));
260    // Set the v and u images to black
261    memset(imgMos->V.ptr[0], 128, (imgMos->V.width * imgMos->V.height) << 1);
262
263    // Do the triangulation.  It returns a sorted list of edges
264    SEdgeVector *edge;
265    int n = m_Triangulator.triangulate(&edge, numCenters, width, height);
266    m_Triangulator.linkNeighbors(edge, n, numCenters);
267
268    // Bounding rectangle that determines the positioning of the rectangle that is
269    // cropped out of the computed mosaic to get rid of the gray borders.
270    MosaicRect cropping_rect;
271
272    if (m_wb.horizontal)
273    {
274        cropping_rect.left = xLeftMost;
275        cropping_rect.right = xRightMost;
276    }
277    else
278    {
279        cropping_rect.top = yTopMost;
280        cropping_rect.bottom = yBottomMost;
281    }
282
283    // Do merging and blending :
284    ret = DoMergeAndBlend(frames, numCenters, width, height, *imgMos, fullRect,
285            cropping_rect, progress, cancelComputation);
286
287    if (m_wb.blendingType == BLEND_TYPE_HORZ)
288        CropFinalMosaic(*imgMos, cropping_rect);
289
290
291    m_Triangulator.freeMemory();    // note: can be called even if delaunay_alloc() wasn't successful
292
293    imageMosaicYVU = imgMos->Y.ptr[0];
294
295
296    if (m_wb.blendingType == BLEND_TYPE_HORZ)
297    {
298        mosaicWidth = cropping_rect.right - cropping_rect.left + 1;
299        mosaicHeight = cropping_rect.bottom - cropping_rect.top + 1;
300    }
301    else
302    {
303        mosaicWidth = Mwidth;
304        mosaicHeight = Mheight;
305    }
306
307    return ret;
308}
309
310int Blend::MosaicSizeCheck(float sizeMultiplier, float heightMultiplier) {
311   if (Mwidth < width || Mheight < height) {
312        return BLEND_RET_ERROR;
313    }
314
315   if ((Mwidth * Mheight) > (width * height * sizeMultiplier)) {
316         return BLEND_RET_ERROR;
317   }
318
319   // We won't do blending for the cases where users swing the device too much
320   // in the secondary direction. We use a short side to determine the
321   // secondary direction because users may hold the device in landsape
322   // or portrait.
323   int shortSide = min(Mwidth, Mheight);
324   if (shortSide > height * heightMultiplier) {
325       return BLEND_RET_ERROR;
326   }
327
328   return BLEND_RET_OK;
329}
330
331int Blend::FillFramePyramid(MosaicFrame *mb)
332{
333    ImageType mbY, mbU, mbV;
334    // Lay this image, centered into the temporary buffer
335    mbY = mb->image;
336    mbU = mb->getU();
337    mbV = mb->getV();
338
339    int h, w;
340
341    for(h=0; h<height; h++)
342    {
343        ImageTypeShort yptr = m_pFrameYPyr->ptr[h];
344        ImageTypeShort uptr = m_pFrameUPyr->ptr[h];
345        ImageTypeShort vptr = m_pFrameVPyr->ptr[h];
346
347        for(w=0; w<width; w++)
348        {
349            yptr[w] = (short) ((*(mbY++)) << 3);
350            uptr[w] = (short) ((*(mbU++)) << 3);
351            vptr[w] = (short) ((*(mbV++)) << 3);
352        }
353    }
354
355    // Spread the image through the border
356    PyramidShort::BorderSpread(m_pFrameYPyr, BORDER, BORDER, BORDER, BORDER);
357    PyramidShort::BorderSpread(m_pFrameUPyr, BORDER, BORDER, BORDER, BORDER);
358    PyramidShort::BorderSpread(m_pFrameVPyr, BORDER, BORDER, BORDER, BORDER);
359
360    // Generate Laplacian pyramids
361    if (!PyramidShort::BorderReduce(m_pFrameYPyr, m_wb.nlevs) || !PyramidShort::BorderExpand(m_pFrameYPyr, m_wb.nlevs, -1) ||
362            !PyramidShort::BorderReduce(m_pFrameUPyr, m_wb.nlevsC) || !PyramidShort::BorderExpand(m_pFrameUPyr, m_wb.nlevsC, -1) ||
363            !PyramidShort::BorderReduce(m_pFrameVPyr, m_wb.nlevsC) || !PyramidShort::BorderExpand(m_pFrameVPyr, m_wb.nlevsC, -1))
364    {
365        LOGE("Error: Could not generate Laplacian pyramids");
366        return BLEND_RET_ERROR;
367    }
368    else
369    {
370        return BLEND_RET_OK;
371    }
372}
373
374int Blend::DoMergeAndBlend(MosaicFrame **frames, int nsite,
375             int width, int height, YUVinfo &imgMos, MosaicRect &rect,
376             MosaicRect &cropping_rect, float &progress, bool &cancelComputation)
377{
378    m_pMosaicYPyr = NULL;
379    m_pMosaicUPyr = NULL;
380    m_pMosaicVPyr = NULL;
381
382    m_pMosaicYPyr = PyramidShort::allocatePyramidPacked(m_wb.nlevs,(unsigned short)rect.Width(),(unsigned short)rect.Height(),BORDER);
383    m_pMosaicUPyr = PyramidShort::allocatePyramidPacked(m_wb.nlevsC,(unsigned short)rect.Width(),(unsigned short)rect.Height(),BORDER);
384    m_pMosaicVPyr = PyramidShort::allocatePyramidPacked(m_wb.nlevsC,(unsigned short)rect.Width(),(unsigned short)rect.Height(),BORDER);
385    if (!m_pMosaicYPyr || !m_pMosaicUPyr || !m_pMosaicVPyr)
386    {
387      LOGE("Error: Could not allocate pyramids for blending");
388      return BLEND_RET_ERROR_MEMORY;
389    }
390
391    MosaicFrame *mb;
392
393    CSite *esite = m_AllSites + nsite;
394    int site_idx;
395
396    // First go through each frame and for each mosaic pixel determine which frame it should come from
397    site_idx = 0;
398    for(CSite *csite = m_AllSites; csite < esite; csite++)
399    {
400        if(cancelComputation)
401        {
402            if (m_pMosaicVPyr) free(m_pMosaicVPyr);
403            if (m_pMosaicUPyr) free(m_pMosaicUPyr);
404            if (m_pMosaicYPyr) free(m_pMosaicYPyr);
405            return BLEND_RET_CANCELLED;
406        }
407
408        mb = csite->getMb();
409
410        mb->vcrect = mb->brect;
411        ClipBlendRect(csite, mb->vcrect);
412
413        ComputeMask(csite, mb->vcrect, mb->brect, rect, imgMos, site_idx);
414
415        site_idx++;
416    }
417
418    ////////// imgMos.Y, imgMos.V, imgMos.U are used as follows //////////////
419    ////////////////////// THIN STRIP MODE ///////////////////////////////////
420
421    // imgMos.Y is used to store the index of the image from which each pixel
422    // in the output mosaic can be read out for the thin-strip mode. Thus,
423    // there is no special handling for pixels around the seam. Also, imgMos.Y
424    // is set to 255 wherever we can't get its value from any input image e.g.
425    // in the gray border areas. imgMos.V and imgMos.U are set to 128 for the
426    // thin-strip mode.
427
428    ////////////////////// WIDE STRIP MODE ///////////////////////////////////
429
430    // imgMos.Y is used the same way as the thin-strip mode.
431    // imgMos.V is used to store the index of the neighboring image which
432    // should contribute to the color of an output pixel in a band around
433    // the seam. Thus, in this band, we will crossfade between the color values
434    // from the image index imgMos.Y and image index imgMos.V. imgMos.U is
435    // used to store the weight (multiplied by 100) that each image will
436    // contribute to the blending process. Thus, we start at 99% contribution
437    // from the first image, then go to 50% contribution from each image at
438    // the seam. Then, the contribution from the second image goes up to 99%.
439
440    // For WIDE mode, set the pixel masks to guide the blender to cross-fade
441    // between the images on either side of each seam:
442    if (m_wb.stripType == STRIP_TYPE_WIDE)
443    {
444        if(m_wb.horizontal)
445        {
446            // Set the number of pixels around the seam to cross-fade between
447            // the two component images,
448            int tw = STRIP_CROSS_FADE_WIDTH_PXLS;
449
450            // Proceed with the image index calculation for cross-fading
451            // only if the cross-fading width is larger than 0
452            if (tw > 0)
453            {
454                for(int y = 0; y < imgMos.Y.height; y++)
455                {
456                    // Since we compare two adjecant pixels to determine
457                    // whether there is a seam, the termination condition of x
458                    // is set to imgMos.Y.width - tw, so that x+1 below
459                    // won't exceed the imgMos' boundary.
460                    for(int x = tw; x < imgMos.Y.width - tw; )
461                    {
462                        // Determine where the seam is...
463                        if (imgMos.Y.ptr[y][x] != imgMos.Y.ptr[y][x+1] &&
464                                imgMos.Y.ptr[y][x] != 255 &&
465                                imgMos.Y.ptr[y][x+1] != 255)
466                        {
467                            // Find the image indices on both sides of the seam
468                            unsigned char idx1 = imgMos.Y.ptr[y][x];
469                            unsigned char idx2 = imgMos.Y.ptr[y][x+1];
470
471                            for (int o = tw; o >= 0; o--)
472                            {
473                                // Set the image index to use for cross-fading
474                                imgMos.V.ptr[y][x - o] = idx2;
475                                // Set the intensity weights to use for cross-fading
476                                imgMos.U.ptr[y][x - o] = 50 + (99 - 50) * o / tw;
477                            }
478
479                            for (int o = 1; o <= tw; o++)
480                            {
481                                // Set the image index to use for cross-fading
482                                imgMos.V.ptr[y][x + o] = idx1;
483                                // Set the intensity weights to use for cross-fading
484                                imgMos.U.ptr[y][x + o] = imgMos.U.ptr[y][x - o];
485                            }
486
487                            x += (tw + 1);
488                        }
489                        else
490                        {
491                            x++;
492                        }
493                    }
494                }
495            }
496        }
497        else
498        {
499            // Set the number of pixels around the seam to cross-fade between
500            // the two component images,
501            int tw = STRIP_CROSS_FADE_WIDTH_PXLS;
502
503            // Proceed with the image index calculation for cross-fading
504            // only if the cross-fading width is larger than 0
505            if (tw > 0)
506            {
507                for(int x = 0; x < imgMos.Y.width; x++)
508                {
509                    // Since we compare two adjecant pixels to determine
510                    // whether there is a seam, the termination condition of y
511                    // is set to imgMos.Y.height - tw, so that y+1 below
512                    // won't exceed the imgMos' boundary.
513                    for(int y = tw; y < imgMos.Y.height - tw; )
514                    {
515                        // Determine where the seam is...
516                        if (imgMos.Y.ptr[y][x] != imgMos.Y.ptr[y+1][x] &&
517                                imgMos.Y.ptr[y][x] != 255 &&
518                                imgMos.Y.ptr[y+1][x] != 255)
519                        {
520                            // Find the image indices on both sides of the seam
521                            unsigned char idx1 = imgMos.Y.ptr[y][x];
522                            unsigned char idx2 = imgMos.Y.ptr[y+1][x];
523
524                            for (int o = tw; o >= 0; o--)
525                            {
526                                // Set the image index to use for cross-fading
527                                imgMos.V.ptr[y - o][x] = idx2;
528                                // Set the intensity weights to use for cross-fading
529                                imgMos.U.ptr[y - o][x] = 50 + (99 - 50) * o / tw;
530                            }
531
532                            for (int o = 1; o <= tw; o++)
533                            {
534                                // Set the image index to use for cross-fading
535                                imgMos.V.ptr[y + o][x] = idx1;
536                                // Set the intensity weights to use for cross-fading
537                                imgMos.U.ptr[y + o][x] = imgMos.U.ptr[y - o][x];
538                            }
539
540                            y += (tw + 1);
541                        }
542                        else
543                        {
544                            y++;
545                        }
546                    }
547                }
548            }
549        }
550
551    }
552
553    // Now perform the actual blending using the frame assignment determined above
554    site_idx = 0;
555    for(CSite *csite = m_AllSites; csite < esite; csite++)
556    {
557        if(cancelComputation)
558        {
559            if (m_pMosaicVPyr) free(m_pMosaicVPyr);
560            if (m_pMosaicUPyr) free(m_pMosaicUPyr);
561            if (m_pMosaicYPyr) free(m_pMosaicYPyr);
562            return BLEND_RET_CANCELLED;
563        }
564
565        mb = csite->getMb();
566
567
568        if(FillFramePyramid(mb)!=BLEND_RET_OK)
569            return BLEND_RET_ERROR;
570
571        ProcessPyramidForThisFrame(csite, mb->vcrect, mb->brect, rect, imgMos, mb->trs, site_idx);
572
573        progress += TIME_PERCENT_BLEND/nsite;
574
575        site_idx++;
576    }
577
578
579    // Blend
580    PerformFinalBlending(imgMos, cropping_rect);
581
582    if (cropping_rect.Width() <= 0 || cropping_rect.Height() <= 0)
583    {
584        LOGE("Size of the cropping_rect is invalid - (width, height): (%d, %d)",
585                cropping_rect.Width(), cropping_rect.Height());
586        return BLEND_RET_ERROR;
587    }
588
589    if (m_pMosaicVPyr) free(m_pMosaicVPyr);
590    if (m_pMosaicUPyr) free(m_pMosaicUPyr);
591    if (m_pMosaicYPyr) free(m_pMosaicYPyr);
592
593    progress += TIME_PERCENT_FINAL;
594
595    return BLEND_RET_OK;
596}
597
598void Blend::CropFinalMosaic(YUVinfo &imgMos, MosaicRect &cropping_rect)
599{
600    int i, j, k;
601    ImageType yimg;
602    ImageType uimg;
603    ImageType vimg;
604
605
606    yimg = imgMos.Y.ptr[0];
607    uimg = imgMos.U.ptr[0];
608    vimg = imgMos.V.ptr[0];
609
610    k = 0;
611    for (j = cropping_rect.top; j <= cropping_rect.bottom; j++)
612    {
613        for (i = cropping_rect.left; i <= cropping_rect.right; i++)
614        {
615            yimg[k] = yimg[j*imgMos.Y.width+i];
616            k++;
617        }
618    }
619    for (j = cropping_rect.top; j <= cropping_rect.bottom; j++)
620    {
621       for (i = cropping_rect.left; i <= cropping_rect.right; i++)
622        {
623            yimg[k] = vimg[j*imgMos.Y.width+i];
624            k++;
625        }
626    }
627    for (j = cropping_rect.top; j <= cropping_rect.bottom; j++)
628    {
629       for (i = cropping_rect.left; i <= cropping_rect.right; i++)
630        {
631            yimg[k] = uimg[j*imgMos.Y.width+i];
632            k++;
633        }
634    }
635}
636
637int Blend::PerformFinalBlending(YUVinfo &imgMos, MosaicRect &cropping_rect)
638{
639    if (!PyramidShort::BorderExpand(m_pMosaicYPyr, m_wb.nlevs, 1) || !PyramidShort::BorderExpand(m_pMosaicUPyr, m_wb.nlevsC, 1) ||
640        !PyramidShort::BorderExpand(m_pMosaicVPyr, m_wb.nlevsC, 1))
641    {
642      LOGE("Error: Could not BorderExpand!");
643      return BLEND_RET_ERROR;
644    }
645
646    ImageTypeShort myimg;
647    ImageTypeShort muimg;
648    ImageTypeShort mvimg;
649    ImageType yimg;
650    ImageType uimg;
651    ImageType vimg;
652
653    int cx = (int)imgMos.Y.width/2;
654    int cy = (int)imgMos.Y.height/2;
655
656    // 2D boolean array that contains true wherever the mosaic image data is
657    // invalid (i.e. in the gray border).
658    bool **b = new bool*[imgMos.Y.height];
659
660    for(int j=0; j<imgMos.Y.height; j++)
661    {
662        b[j] = new bool[imgMos.Y.width];
663    }
664
665    // Copy the resulting image into the full image using the mask
666    int i, j;
667
668    yimg = imgMos.Y.ptr[0];
669    uimg = imgMos.U.ptr[0];
670    vimg = imgMos.V.ptr[0];
671
672    for (j = 0; j < imgMos.Y.height; j++)
673    {
674        myimg = m_pMosaicYPyr->ptr[j];
675        muimg = m_pMosaicUPyr->ptr[j];
676        mvimg = m_pMosaicVPyr->ptr[j];
677
678        for (i = 0; i<imgMos.Y.width; i++)
679        {
680            // A final mask was set up previously,
681            // if the value is zero skip it, otherwise replace it.
682            if (*yimg <255)
683            {
684                short value = (short) ((*myimg) >> 3);
685                if (value < 0) value = 0;
686                else if (value > 255) value = 255;
687                *yimg = (unsigned char) value;
688
689                value = (short) ((*muimg) >> 3);
690                if (value < 0) value = 0;
691                else if (value > 255) value = 255;
692                *uimg = (unsigned char) value;
693
694                value = (short) ((*mvimg) >> 3);
695                if (value < 0) value = 0;
696                else if (value > 255) value = 255;
697                *vimg = (unsigned char) value;
698
699                b[j][i] = false;
700
701            }
702            else
703            {   // set border color in here
704                *yimg = (unsigned char) 96;
705                *uimg = (unsigned char) 128;
706                *vimg = (unsigned char) 128;
707
708                b[j][i] = true;
709            }
710
711            yimg++;
712            uimg++;
713            vimg++;
714            myimg++;
715            muimg++;
716            mvimg++;
717        }
718    }
719
720    if(m_wb.horizontal)
721    {
722        //Scan through each row and increment top if the row contains any gray
723        for (j = 0; j < imgMos.Y.height; j++)
724        {
725            for (i = cropping_rect.left; i < cropping_rect.right; i++)
726            {
727                if (b[j][i])
728                {
729                    break; // to next row
730                }
731            }
732
733            if (i == cropping_rect.right)   //no gray pixel in this row!
734            {
735                cropping_rect.top = j;
736                break;
737            }
738        }
739
740        //Scan through each row and decrement bottom if the row contains any gray
741        for (j = imgMos.Y.height-1; j >= 0; j--)
742        {
743            for (i = cropping_rect.left; i < cropping_rect.right; i++)
744            {
745                if (b[j][i])
746                {
747                    break; // to next row
748                }
749            }
750
751            if (i == cropping_rect.right)   //no gray pixel in this row!
752            {
753                cropping_rect.bottom = j;
754                break;
755            }
756        }
757    }
758    else // Vertical Mosaic
759    {
760        //Scan through each column and increment left if the column contains any gray
761        for (i = 0; i < imgMos.Y.width; i++)
762        {
763            for (j = cropping_rect.top; j < cropping_rect.bottom; j++)
764            {
765                if (b[j][i])
766                {
767                    break; // to next column
768                }
769            }
770
771            if (j == cropping_rect.bottom)   //no gray pixel in this column!
772            {
773                cropping_rect.left = i;
774                break;
775            }
776        }
777
778        //Scan through each column and decrement right if the column contains any gray
779        for (i = imgMos.Y.width-1; i >= 0; i--)
780        {
781            for (j = cropping_rect.top; j < cropping_rect.bottom; j++)
782            {
783                if (b[j][i])
784                {
785                    break; // to next column
786                }
787            }
788
789            if (j == cropping_rect.bottom)   //no gray pixel in this column!
790            {
791                cropping_rect.right = i;
792                break;
793            }
794        }
795
796    }
797
798    RoundingCroppingSizeToMultipleOf8(cropping_rect);
799
800    for(int j=0; j<imgMos.Y.height; j++)
801    {
802        delete b[j];
803    }
804
805    delete b;
806
807    return BLEND_RET_OK;
808}
809
810void Blend::RoundingCroppingSizeToMultipleOf8(MosaicRect &rect) {
811    int height = rect.bottom - rect.top + 1;
812    int residue = height & 7;
813    rect.bottom -= residue;
814
815    int width = rect.right - rect.left + 1;
816    residue = width & 7;
817    rect.right -= residue;
818}
819
820void Blend::ComputeMask(CSite *csite, BlendRect &vcrect, BlendRect &brect, MosaicRect &rect, YUVinfo &imgMos, int site_idx)
821{
822    PyramidShort *dptr = m_pMosaicYPyr;
823
824    int nC = m_wb.nlevsC;
825    int l = (int) ((vcrect.lft - rect.left));
826    int b = (int) ((vcrect.bot - rect.top));
827    int r = (int) ((vcrect.rgt - rect.left));
828    int t = (int) ((vcrect.top - rect.top));
829
830    if (vcrect.lft == brect.lft)
831        l = (l <= 0) ? -BORDER : l - BORDER;
832    else if (l < -BORDER)
833        l = -BORDER;
834
835    if (vcrect.bot == brect.bot)
836        b = (b <= 0) ? -BORDER : b - BORDER;
837    else if (b < -BORDER)
838        b = -BORDER;
839
840    if (vcrect.rgt == brect.rgt)
841        r = (r >= dptr->width) ? dptr->width + BORDER - 1 : r + BORDER;
842    else if (r >= dptr->width + BORDER)
843        r = dptr->width + BORDER - 1;
844
845    if (vcrect.top == brect.top)
846        t = (t >= dptr->height) ? dptr->height + BORDER - 1 : t + BORDER;
847    else if (t >= dptr->height + BORDER)
848        t = dptr->height + BORDER - 1;
849
850    // Walk the Region of interest and populate the pyramid
851    for (int j = b; j <= t; j++)
852    {
853        int jj = j;
854        double sj = jj + rect.top;
855
856        for (int i = l; i <= r; i++)
857        {
858            int ii = i;
859            // project point and then triangulate to neighbors
860            double si = ii + rect.left;
861
862            double dself = hypotSq(csite->getVCenter().x - si, csite->getVCenter().y - sj);
863            int inMask = ((unsigned) ii < imgMos.Y.width &&
864                    (unsigned) jj < imgMos.Y.height) ? 1 : 0;
865
866            if(!inMask)
867                continue;
868
869            // scan the neighbors to see if this is a valid position
870            unsigned char mask = (unsigned char) 255;
871            SEdgeVector *ce;
872            int ecnt;
873            for (ce = csite->getNeighbor(), ecnt = csite->getNumNeighbors(); ecnt--; ce++)
874            {
875                double d1 = hypotSq(m_AllSites[ce->second].getVCenter().x - si,
876                        m_AllSites[ce->second].getVCenter().y - sj);
877                if (d1 < dself)
878                {
879                    break;
880                }
881            }
882
883            if (ecnt >= 0) continue;
884
885            imgMos.Y.ptr[jj][ii] = (unsigned char)site_idx;
886        }
887    }
888}
889
890void Blend::ProcessPyramidForThisFrame(CSite *csite, BlendRect &vcrect, BlendRect &brect, MosaicRect &rect, YUVinfo &imgMos, double trs[3][3], int site_idx)
891{
892    // Put the Region of interest (for all levels) into m_pMosaicYPyr
893    double inv_trs[3][3];
894    inv33d(trs, inv_trs);
895
896    // Process each pyramid level
897    PyramidShort *sptr = m_pFrameYPyr;
898    PyramidShort *suptr = m_pFrameUPyr;
899    PyramidShort *svptr = m_pFrameVPyr;
900
901    PyramidShort *dptr = m_pMosaicYPyr;
902    PyramidShort *duptr = m_pMosaicUPyr;
903    PyramidShort *dvptr = m_pMosaicVPyr;
904
905    int dscale = 0; // distance scale for the current level
906    int nC = m_wb.nlevsC;
907    for (int n = m_wb.nlevs; n--; dscale++, dptr++, sptr++, dvptr++, duptr++, svptr++, suptr++, nC--)
908    {
909        int l = (int) ((vcrect.lft - rect.left) / (1 << dscale));
910        int b = (int) ((vcrect.bot - rect.top) / (1 << dscale));
911        int r = (int) ((vcrect.rgt - rect.left) / (1 << dscale) + .5);
912        int t = (int) ((vcrect.top - rect.top) / (1 << dscale) + .5);
913
914        if (vcrect.lft == brect.lft)
915            l = (l <= 0) ? -BORDER : l - BORDER;
916        else if (l < -BORDER)
917            l = -BORDER;
918
919        if (vcrect.bot == brect.bot)
920            b = (b <= 0) ? -BORDER : b - BORDER;
921        else if (b < -BORDER)
922            b = -BORDER;
923
924        if (vcrect.rgt == brect.rgt)
925            r = (r >= dptr->width) ? dptr->width + BORDER - 1 : r + BORDER;
926        else if (r >= dptr->width + BORDER)
927            r = dptr->width + BORDER - 1;
928
929        if (vcrect.top == brect.top)
930            t = (t >= dptr->height) ? dptr->height + BORDER - 1 : t + BORDER;
931        else if (t >= dptr->height + BORDER)
932            t = dptr->height + BORDER - 1;
933
934        // Walk the Region of interest and populate the pyramid
935        for (int j = b; j <= t; j++)
936        {
937            int jj = (j << dscale);
938            double sj = jj + rect.top;
939
940            for (int i = l; i <= r; i++)
941            {
942                int ii = (i << dscale);
943                // project point and then triangulate to neighbors
944                double si = ii + rect.left;
945
946                int inMask = ((unsigned) ii < imgMos.Y.width &&
947                        (unsigned) jj < imgMos.Y.height) ? 1 : 0;
948
949                if(inMask && imgMos.Y.ptr[jj][ii] != site_idx &&
950                        imgMos.V.ptr[jj][ii] != site_idx &&
951                        imgMos.Y.ptr[jj][ii] != 255)
952                    continue;
953
954                // Setup weights for cross-fading
955                // Weight of the intensity already in the output pixel
956                double wt0 = 0.0;
957                // Weight of the intensity from the input pixel (current frame)
958                double wt1 = 1.0;
959
960                if (m_wb.stripType == STRIP_TYPE_WIDE)
961                {
962                    if(inMask && imgMos.Y.ptr[jj][ii] != 255)
963                    {
964                        // If not on a seam OR pyramid level exceeds
965                        // maximum level for cross-fading.
966                        if((imgMos.V.ptr[jj][ii] == 128) ||
967                            (dscale > STRIP_CROSS_FADE_MAX_PYR_LEVEL))
968                        {
969                            wt0 = 0.0;
970                            wt1 = 1.0;
971                        }
972                        else
973                        {
974                            wt0 = 1.0;
975                            wt1 = ((imgMos.Y.ptr[jj][ii] == site_idx) ?
976                                    (double)imgMos.U.ptr[jj][ii] / 100.0 :
977                                    1.0 - (double)imgMos.U.ptr[jj][ii] / 100.0);
978                        }
979                    }
980                }
981
982                // Project this mosaic point into the original frame coordinate space
983                double xx, yy;
984
985                MosaicToFrame(inv_trs, si, sj, xx, yy);
986
987                if (xx < 0.0 || yy < 0.0 || xx > width - 1.0 || yy > height - 1.0)
988                {
989                    if(inMask)
990                    {
991                        imgMos.Y.ptr[jj][ii] = 255;
992                        wt0 = 0.0f;
993                        wt1 = 1.0f;
994                    }
995                }
996
997                xx /= (1 << dscale);
998                yy /= (1 << dscale);
999
1000
1001                int x1 = (xx >= 0.0) ? (int) xx : (int) floor(xx);
1002                int y1 = (yy >= 0.0) ? (int) yy : (int) floor(yy);
1003
1004                // Final destination in extended pyramid
1005#ifndef LINEAR_INTERP
1006                if(inSegment(x1, sptr->width, BORDER-1) &&
1007                        inSegment(y1, sptr->height, BORDER-1))
1008                {
1009                    double xfrac = xx - x1;
1010                    double yfrac = yy - y1;
1011                    dptr->ptr[j][i] = (short) (wt0 * dptr->ptr[j][i] + .5 +
1012                            wt1 * ciCalc(sptr, x1, y1, xfrac, yfrac));
1013                    if (dvptr >= m_pMosaicVPyr && nC > 0)
1014                    {
1015                        duptr->ptr[j][i] = (short) (wt0 * duptr->ptr[j][i] + .5 +
1016                                wt1 * ciCalc(suptr, x1, y1, xfrac, yfrac));
1017                        dvptr->ptr[j][i] = (short) (wt0 * dvptr->ptr[j][i] + .5 +
1018                                wt1 * ciCalc(svptr, x1, y1, xfrac, yfrac));
1019                    }
1020                }
1021#else
1022                if(inSegment(x1, sptr->width, BORDER) && inSegment(y1, sptr->height, BORDER))
1023                {
1024                    int x2 = x1 + 1;
1025                    int y2 = y1 + 1;
1026                    double xfrac = xx - x1;
1027                    double yfrac = yy - y1;
1028                    double y1val = sptr->ptr[y1][x1] +
1029                        (sptr->ptr[y1][x2] - sptr->ptr[y1][x1]) * xfrac;
1030                    double y2val = sptr->ptr[y2][x1] +
1031                        (sptr->ptr[y2][x2] - sptr->ptr[y2][x1]) * xfrac;
1032                    dptr->ptr[j][i] = (short) (y1val + yfrac * (y2val - y1val));
1033
1034                    if (dvptr >= m_pMosaicVPyr && nC > 0)
1035                    {
1036                        y1val = suptr->ptr[y1][x1] +
1037                            (suptr->ptr[y1][x2] - suptr->ptr[y1][x1]) * xfrac;
1038                        y2val = suptr->ptr[y2][x1] +
1039                            (suptr->ptr[y2][x2] - suptr->ptr[y2][x1]) * xfrac;
1040
1041                        duptr->ptr[j][i] = (short) (y1val + yfrac * (y2val - y1val));
1042
1043                        y1val = svptr->ptr[y1][x1] +
1044                            (svptr->ptr[y1][x2] - svptr->ptr[y1][x1]) * xfrac;
1045                        y2val = svptr->ptr[y2][x1] +
1046                            (svptr->ptr[y2][x2] - svptr->ptr[y2][x1]) * xfrac;
1047
1048                        dvptr->ptr[j][i] = (short) (y1val + yfrac * (y2val - y1val));
1049                    }
1050                }
1051#endif
1052                else
1053                {
1054                    clipToSegment(x1, sptr->width, BORDER);
1055                    clipToSegment(y1, sptr->height, BORDER);
1056
1057                    dptr->ptr[j][i] = (short) (wt0 * dptr->ptr[j][i] + 0.5 +
1058                            wt1 * sptr->ptr[y1][x1] );
1059                    if (dvptr >= m_pMosaicVPyr && nC > 0)
1060                    {
1061                        dvptr->ptr[j][i] = (short) (wt0 * dvptr->ptr[j][i] +
1062                                0.5 + wt1 * svptr->ptr[y1][x1] );
1063                        duptr->ptr[j][i] = (short) (wt0 * duptr->ptr[j][i] +
1064                                0.5 + wt1 * suptr->ptr[y1][x1] );
1065                    }
1066                }
1067            }
1068        }
1069    }
1070}
1071
1072void Blend::MosaicToFrame(double trs[3][3], double x, double y, double &wx, double &wy)
1073{
1074    double X, Y, z;
1075    if (m_wb.theta == 0.0)
1076    {
1077        X = x;
1078        Y = y;
1079    }
1080    else if (m_wb.horizontal)
1081    {
1082        double alpha = x * m_wb.direction / m_wb.width;
1083        double length = (y - alpha * m_wb.correction) * m_wb.direction + m_wb.radius;
1084        double deltaTheta = m_wb.theta * alpha;
1085        double sinTheta = sin(deltaTheta);
1086        double cosTheta = sqrt(1.0 - sinTheta * sinTheta) * m_wb.direction;
1087        X = length * sinTheta + m_wb.x;
1088        Y = length * cosTheta + m_wb.y;
1089    }
1090    else
1091    {
1092        double alpha = y * m_wb.direction / m_wb.width;
1093        double length = (x - alpha * m_wb.correction) * m_wb.direction + m_wb.radius;
1094        double deltaTheta = m_wb.theta * alpha;
1095        double sinTheta = sin(deltaTheta);
1096        double cosTheta = sqrt(1.0 - sinTheta * sinTheta) * m_wb.direction;
1097        Y = length * sinTheta + m_wb.y;
1098        X = length * cosTheta + m_wb.x;
1099    }
1100    z = ProjZ(trs, X, Y, 1.0);
1101    wx = ProjX(trs, X, Y, z, 1.0);
1102    wy = ProjY(trs, X, Y, z, 1.0);
1103}
1104
1105void Blend::FrameToMosaic(double trs[3][3], double x, double y, double &wx, double &wy)
1106{
1107    // Project into the intermediate Mosaic coordinate system
1108    double z = ProjZ(trs, x, y, 1.0);
1109    double X = ProjX(trs, x, y, z, 1.0);
1110    double Y = ProjY(trs, x, y, z, 1.0);
1111
1112    if (m_wb.theta == 0.0)
1113    {
1114        // No rotation, then this is all we need to do.
1115        wx = X;
1116        wy = Y;
1117    }
1118    else if (m_wb.horizontal)
1119    {
1120        double deltaX = X - m_wb.x;
1121        double deltaY = Y - m_wb.y;
1122        double length = sqrt(deltaX * deltaX + deltaY * deltaY);
1123        double deltaTheta = asin(deltaX / length);
1124        double alpha = deltaTheta / m_wb.theta;
1125        wx = alpha * m_wb.width * m_wb.direction;
1126        wy = (length - m_wb.radius) * m_wb.direction + alpha * m_wb.correction;
1127    }
1128    else
1129    {
1130        double deltaX = X - m_wb.x;
1131        double deltaY = Y - m_wb.y;
1132        double length = sqrt(deltaX * deltaX + deltaY * deltaY);
1133        double deltaTheta = asin(deltaY / length);
1134        double alpha = deltaTheta / m_wb.theta;
1135        wy = alpha * m_wb.width * m_wb.direction;
1136        wx = (length - m_wb.radius) * m_wb.direction + alpha * m_wb.correction;
1137    }
1138}
1139
1140
1141
1142// Clip the region of interest as small as possible by using the Voronoi edges of
1143// the neighbors
1144void Blend::ClipBlendRect(CSite *csite, BlendRect &brect)
1145{
1146      SEdgeVector *ce;
1147      int ecnt;
1148      for (ce = csite->getNeighbor(), ecnt = csite->getNumNeighbors(); ecnt--; ce++)
1149      {
1150        // calculate the Voronoi bisector intersection
1151        const double epsilon = 1e-5;
1152        double dx = (m_AllSites[ce->second].getVCenter().x - m_AllSites[ce->first].getVCenter().x);
1153        double dy = (m_AllSites[ce->second].getVCenter().y - m_AllSites[ce->first].getVCenter().y);
1154        double xmid = m_AllSites[ce->first].getVCenter().x + dx/2.0;
1155        double ymid = m_AllSites[ce->first].getVCenter().y + dy/2.0;
1156        double inter;
1157
1158        if (dx > epsilon)
1159        {
1160          // neighbor is on right
1161          if ((inter = m_wb.roundoffOverlap + xmid - dy * (((dy >= 0.0) ? brect.bot : brect.top) - ymid) / dx) < brect.rgt)
1162            brect.rgt = inter;
1163        }
1164        else if (dx < -epsilon)
1165        {
1166          // neighbor is on left
1167          if ((inter = -m_wb.roundoffOverlap + xmid - dy * (((dy >= 0.0) ? brect.bot : brect.top) - ymid) / dx) > brect.lft)
1168            brect.lft = inter;
1169        }
1170        if (dy > epsilon)
1171        {
1172          // neighbor is above
1173          if ((inter = m_wb.roundoffOverlap + ymid - dx * (((dx >= 0.0) ? brect.lft : brect.rgt) - xmid) / dy) < brect.top)
1174            brect.top = inter;
1175        }
1176        else if (dy < -epsilon)
1177        {
1178          // neighbor is below
1179          if ((inter = -m_wb.roundoffOverlap + ymid - dx * (((dx >= 0.0) ? brect.lft : brect.rgt) - xmid) / dy) > brect.bot)
1180            brect.bot = inter;
1181        }
1182      }
1183}
1184
1185void Blend::FrameToMosaicRect(int width, int height, double trs[3][3], BlendRect &brect)
1186{
1187    // We need to walk the perimeter since the borders can be bent.
1188    brect.lft = brect.bot = 2e30;
1189    brect.rgt = brect.top = -2e30;
1190    double xpos, ypos;
1191    double lasty = height - 1.0;
1192    double lastx = width - 1.0;
1193    int i;
1194
1195    for (i = width; i--;)
1196    {
1197
1198        FrameToMosaic(trs, (double) i, 0.0, xpos, ypos);
1199        ClipRect(xpos, ypos, brect);
1200        FrameToMosaic(trs, (double) i, lasty, xpos, ypos);
1201        ClipRect(xpos, ypos, brect);
1202    }
1203    for (i = height; i--;)
1204    {
1205        FrameToMosaic(trs, 0.0, (double) i, xpos, ypos);
1206        ClipRect(xpos, ypos, brect);
1207        FrameToMosaic(trs, lastx, (double) i, xpos, ypos);
1208        ClipRect(xpos, ypos, brect);
1209    }
1210}
1211
1212void Blend::SelectRelevantFrames(MosaicFrame **frames, int frames_size,
1213        MosaicFrame **relevant_frames, int &relevant_frames_size)
1214{
1215    MosaicFrame *first = frames[0];
1216    MosaicFrame *last = frames[frames_size-1];
1217    MosaicFrame *mb;
1218
1219    double fxpos = first->trs[0][2], fypos = first->trs[1][2];
1220
1221    double midX = last->width / 2.0;
1222    double midY = last->height / 2.0;
1223    double z = ProjZ(first->trs, midX, midY, 1.0);
1224    double firstX, firstY;
1225    double prevX = firstX = ProjX(first->trs, midX, midY, z, 1.0);
1226    double prevY = firstY = ProjY(first->trs, midX, midY, z, 1.0);
1227
1228    relevant_frames[0] = first; // Add first frame by default
1229    relevant_frames_size = 1;
1230
1231    for (int i = 0; i < frames_size - 1; i++)
1232    {
1233        mb = frames[i];
1234        double currX, currY;
1235        z = ProjZ(mb->trs, midX, midY, 1.0);
1236        currX = ProjX(mb->trs, midX, midY, z, 1.0);
1237        currY = ProjY(mb->trs, midX, midY, z, 1.0);
1238        double deltaX = currX - prevX;
1239        double deltaY = currY - prevY;
1240        double center2centerDist = sqrt(deltaY * deltaY + deltaX * deltaX);
1241
1242        if (fabs(deltaX) > STRIP_SEPARATION_THRESHOLD_PXLS ||
1243                fabs(deltaY) > STRIP_SEPARATION_THRESHOLD_PXLS)
1244        {
1245            relevant_frames[relevant_frames_size] = mb;
1246            relevant_frames_size++;
1247
1248            prevX = currX;
1249            prevY = currY;
1250        }
1251    }
1252
1253    // Add last frame by default
1254    relevant_frames[relevant_frames_size] = last;
1255    relevant_frames_size++;
1256}
1257
1258void Blend::ComputeBlendParameters(MosaicFrame **frames, int frames_size, int is360)
1259{
1260    // For FULL and PAN modes, we do not unwarp the mosaic into a rectangular coordinate system
1261    // and so we set the theta to 0 and return.
1262    if (m_wb.blendingType != BLEND_TYPE_CYLPAN && m_wb.blendingType != BLEND_TYPE_HORZ)
1263    {
1264        m_wb.theta = 0.0;
1265        return;
1266    }
1267
1268    MosaicFrame *first = frames[0];
1269    MosaicFrame *last = frames[frames_size-1];
1270    MosaicFrame *mb;
1271
1272    double lxpos = last->trs[0][2], lypos = last->trs[1][2];
1273    double fxpos = first->trs[0][2], fypos = first->trs[1][2];
1274
1275    // Calculate warp to produce proper stitching.
1276    // get x, y displacement
1277    double midX = last->width / 2.0;
1278    double midY = last->height / 2.0;
1279    double z = ProjZ(first->trs, midX, midY, 1.0);
1280    double firstX, firstY;
1281    double prevX = firstX = ProjX(first->trs, midX, midY, z, 1.0);
1282    double prevY = firstY = ProjY(first->trs, midX, midY, z, 1.0);
1283
1284    double arcLength, lastTheta;
1285    m_wb.theta = lastTheta = arcLength = 0.0;
1286
1287    // Step through all the frames to compute the total arc-length of the cone
1288    // swept while capturing the mosaic (in the original conical coordinate system).
1289    for (int i = 0; i < frames_size; i++)
1290    {
1291        mb = frames[i];
1292        double currX, currY;
1293        z = ProjZ(mb->trs, midX, midY, 1.0);
1294        currX = ProjX(mb->trs, midX, midY, z, 1.0);
1295        currY = ProjY(mb->trs, midX, midY, z, 1.0);
1296        double deltaX = currX - prevX;
1297        double deltaY = currY - prevY;
1298
1299        // The arcLength is computed by summing the lengths of the chords
1300        // connecting the pairwise projected image centers of the input image frames.
1301        arcLength += sqrt(deltaY * deltaY + deltaX * deltaX);
1302
1303        if (!is360)
1304        {
1305            double thisTheta = asin(mb->trs[1][0]);
1306            m_wb.theta += thisTheta - lastTheta;
1307            lastTheta = thisTheta;
1308        }
1309
1310        prevX = currX;
1311        prevY = currY;
1312    }
1313
1314    // Stretch this to end at the proper alignment i.e. the width of the
1315    // rectangle is determined by the arcLength computed above and the cone
1316    // sector angle is determined using the rotation of the last frame.
1317    m_wb.width = arcLength;
1318    if (is360) m_wb.theta = asin(last->trs[1][0]);
1319
1320    // If there is no rotation, we're done.
1321    if (m_wb.theta != 0.0)
1322    {
1323        double dx = prevX - firstX;
1324        double dy = prevY - firstY;
1325
1326        // If the mosaic was captured by sweeping horizontally
1327        if (abs(lxpos - fxpos) > abs(lypos - fypos))
1328        {
1329            m_wb.horizontal = 1;
1330            // Calculate radius position to make ends exactly the same Y offset
1331            double radiusTheta = dx / cos(3.14159 / 2.0 - m_wb.theta);
1332            m_wb.radius = dy + radiusTheta * cos(m_wb.theta);
1333            if (m_wb.radius < 0.0) m_wb.radius = -m_wb.radius;
1334        }
1335        else
1336        {
1337            m_wb.horizontal = 0;
1338            // Calculate radius position to make ends exactly the same Y offset
1339            double radiusTheta = dy / cos(3.14159 / 2.0 - m_wb.theta);
1340            m_wb.radius = dx + radiusTheta * cos(m_wb.theta);
1341            if (m_wb.radius < 0.0) m_wb.radius = -m_wb.radius;
1342        }
1343
1344        // Determine major direction
1345        if (m_wb.horizontal)
1346        {
1347            // Horizontal strip
1348            // m_wb.x,y record the origin of the rectangle coordinate system.
1349            if (is360) m_wb.x = firstX;
1350            else
1351            {
1352                if (lxpos - fxpos < 0)
1353                {
1354                    m_wb.x = firstX + midX;
1355                    z = ProjZ(last->trs, 0.0, midY, 1.0);
1356                    prevX = ProjX(last->trs, 0.0, midY, z, 1.0);
1357                    prevY = ProjY(last->trs, 0.0, midY, z, 1.0);
1358                }
1359                else
1360                {
1361                    m_wb.x = firstX - midX;
1362                    z = ProjZ(last->trs, last->width - 1.0, midY, 1.0);
1363                    prevX = ProjX(last->trs, last->width - 1.0, midY, z, 1.0);
1364                    prevY = ProjY(last->trs, last->width - 1.0, midY, z, 1.0);
1365                }
1366            }
1367            dy = prevY - firstY;
1368            if (dy < 0.0) m_wb.direction = 1.0;
1369            else m_wb.direction = -1.0;
1370            m_wb.y = firstY - m_wb.radius * m_wb.direction;
1371            if (dy * m_wb.theta > 0.0) m_wb.width = -m_wb.width;
1372        }
1373        else
1374        {
1375            // Vertical strip
1376            if (is360) m_wb.y = firstY;
1377            else
1378            {
1379                if (lypos - fypos < 0)
1380                {
1381                    m_wb.x = firstY + midY;
1382                    z = ProjZ(last->trs, midX, 0.0, 1.0);
1383                    prevX = ProjX(last->trs, midX, 0.0, z, 1.0);
1384                    prevY = ProjY(last->trs, midX, 0.0, z, 1.0);
1385                }
1386                else
1387                {
1388                    m_wb.x = firstX - midX;
1389                    z = ProjZ(last->trs, midX, last->height - 1.0, 1.0);
1390                    prevX = ProjX(last->trs, midX, last->height - 1.0, z, 1.0);
1391                    prevY = ProjY(last->trs, midX, last->height - 1.0, z, 1.0);
1392                }
1393            }
1394            dx = prevX - firstX;
1395            if (dx < 0.0) m_wb.direction = 1.0;
1396            else m_wb.direction = -1.0;
1397            m_wb.x = firstX - m_wb.radius * m_wb.direction;
1398            if (dx * m_wb.theta > 0.0) m_wb.width = -m_wb.width;
1399        }
1400
1401        // Calculate the correct correction factor
1402        double deltaX = prevX - m_wb.x;
1403        double deltaY = prevY - m_wb.y;
1404        double length = sqrt(deltaX * deltaX + deltaY * deltaY);
1405        double deltaTheta = (m_wb.horizontal) ? deltaX : deltaY;
1406        deltaTheta = asin(deltaTheta / length);
1407        m_wb.correction = ((m_wb.radius - length) * m_wb.direction) /
1408            (deltaTheta / m_wb.theta);
1409    }
1410}
1411