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