Blend.cpp revision e1178a73fd5756771d25d0b8375452450f509e99
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 29Blend::Blend() 30{ 31 m_wb.blendingType = BLEND_TYPE_NONE; 32} 33 34Blend::~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 41int Blend::initialize(int blendingType, int frame_width, int frame_height) 42{ 43 this->width = frame_width; 44 this->height = frame_height; 45 this->m_wb.blendingType = blendingType; 46 47 m_wb.blendRange = m_wb.blendRangeUV = BLEND_RANGE_DEFAULT; 48 m_wb.nlevs = m_wb.blendRange; 49 m_wb.nlevsC = m_wb.blendRangeUV; 50 51 if (m_wb.nlevs <= 0) m_wb.nlevs = 1; // Need levels for YUV processing 52 if (m_wb.nlevsC > m_wb.nlevs) m_wb.nlevsC = m_wb.nlevs; 53 54 m_wb.roundoffOverlap = 1.5; 55 56 m_pFrameYPyr = NULL; 57 m_pFrameUPyr = NULL; 58 m_pFrameVPyr = NULL; 59 60 m_pFrameYPyr = PyramidShort::allocatePyramidPacked(m_wb.nlevs, (unsigned short) width, (unsigned short) height, BORDER); 61 m_pFrameUPyr = PyramidShort::allocatePyramidPacked(m_wb.nlevsC, (unsigned short) (width), (unsigned short) (height), BORDER); 62 m_pFrameVPyr = PyramidShort::allocatePyramidPacked(m_wb.nlevsC, (unsigned short) (width), (unsigned short) (height), BORDER); 63 64 if (!m_pFrameYPyr || !m_pFrameUPyr || !m_pFrameVPyr) 65 { 66 LOGIE("Error: Could not allocate pyramids for blending\n"); 67 return BLEND_RET_ERROR_MEMORY; 68 } 69 70 return BLEND_RET_OK; 71} 72 73void Blend::AlignToMiddleFrame(MosaicFrame **frames, int frames_size) 74{ 75 // Unwarp this frame and Warp the others to match 76 MosaicFrame *mb = NULL; 77 MosaicFrame *ref = frames[int(frames_size/2)]; // Middle frame 78 79 double invtrs[3][3]; 80 inv33d(ref->trs, invtrs); 81 82 for(int mfit = 0; mfit < frames_size; mfit++) 83 { 84 mb = frames[mfit]; 85 double temp[3][3]; 86 mult33d(temp, invtrs, mb->trs); 87 memcpy(mb->trs, temp, sizeof(temp)); 88 normProjMat33d(mb->trs); 89 } 90} 91 92int Blend::runBlend(MosaicFrame **frames, int frames_size, 93 ImageType &imageMosaicYVU, int &mosaicWidth, int &mosaicHeight, 94 float &progress, bool &cancelComputation) 95{ 96 int ret; 97 int numCenters; 98 99 ComputeBlendParameters(frames, frames_size, true); 100 numCenters = frames_size; 101 102 if (numCenters == 0) 103 { 104 LOGIE("Error: No frames to blend\n"); 105 return BLEND_RET_ERROR; 106 } 107 108 if (!(m_AllSites = m_Triangulator.allocMemory(numCenters))) 109 { 110 return BLEND_RET_ERROR_MEMORY; 111 } 112 113 BlendRect global_rect; 114 115 global_rect.lft = global_rect.bot = 2e30; // min values 116 global_rect.rgt = global_rect.top = -2e30; // max values 117 MosaicFrame *mb = NULL; 118 double halfwidth = width / 2.0; 119 double halfheight = height / 2.0; 120 121 double z, x0, y0, x1, y1, x2, y2, x3, y3; 122 123 // Determine the extents of the final mosaic 124 CSite *csite = m_AllSites ; 125 for(int mfit = 0; mfit < frames_size; mfit++) 126 { 127 mb = frames[mfit]; 128 129 // Compute clipping for this frame's rect 130 FrameToMosaicRect(mb->width, mb->height, mb->trs, mb->brect); 131 // Clip global rect using this frame's rect 132 ClipRect(mb->brect, global_rect); 133 134 // Calculate the corner points 135 FrameToMosaic(mb->trs, 0.0, 0.0, x0, y0); 136 FrameToMosaic(mb->trs, 0.0, mb->height-1.0, x1, y1); 137 FrameToMosaic(mb->trs, mb->width-1.0, mb->height-1.0, x2, y2); 138 FrameToMosaic(mb->trs, mb->width-1.0, 0.0, x3, y3); 139 140 // Compute the centroid of the warped region 141 FindQuadCentroid(x0, y0, x1, y1, x2, y2, x3, y3, csite->getVCenter().x, csite->getVCenter().y); 142 143 csite->setMb(mb); 144 csite++; 145 } 146 147 // Get origin and sizes 148 MosaicRect fullRect; 149 150 fullRect.left = (int) floor(global_rect.lft); // min-x 151 fullRect.top = (int) floor(global_rect.bot); // min-y 152 fullRect.right = (int) ceil(global_rect.rgt); // max-x 153 fullRect.bottom = (int) ceil(global_rect.top);// max-y 154 Mwidth = (unsigned short) (fullRect.right - fullRect.left + 1); 155 Mheight = (unsigned short) (fullRect.bottom - fullRect.top + 1); 156 157 // Make sure image width is multiple of 4 158 Mwidth = (unsigned short) ((Mwidth + 3) & ~3); 159 Mheight = (unsigned short) ((Mheight + 3) & ~3); // Round up. 160 161 if (Mwidth < width || Mheight < height) 162 { 163 LOGIE("RunBlend: aborting - consistency check failed, w=%d, h=%d\n", Mwidth, Mheight); 164 return BLEND_RET_ERROR; 165 } 166 167 YUVinfo *imgMos = YUVinfo::allocateImage(Mwidth, Mheight); 168 if (imgMos == NULL) 169 { 170 LOGIE("RunBlend: aborting - couldn't alloc %d x %d mosaic image\n", Mwidth, Mheight); 171 return BLEND_RET_ERROR_MEMORY; 172 } 173 174 // Set the Y image to 255 so we can distinguish when frame idx are written to it 175 memset(imgMos->Y.ptr[0], 255, (imgMos->Y.width * imgMos->Y.height)); 176 // Set the v and u images to black 177 memset(imgMos->V.ptr[0], 128, (imgMos->V.width * imgMos->V.height) << 1); 178 179 // Do the triangulation. It returns a sorted list of edges 180 SEdgeVector *edge; 181 int n = m_Triangulator.triangulate(&edge, numCenters, width, height); 182 m_Triangulator.linkNeighbors(edge, n, numCenters); 183 184 MosaicRect cropping_rect; 185 186 // Do merging and blending : 187 ret = DoMergeAndBlend(frames, numCenters, width, height, *imgMos, fullRect, 188 cropping_rect, progress, cancelComputation); 189 190 if (m_wb.blendingType == BLEND_TYPE_HORZ) 191 CropFinalMosaic(*imgMos, cropping_rect); 192 193 194 m_Triangulator.freeMemory(); // note: can be called even if delaunay_alloc() wasn't successful 195 196 imageMosaicYVU = imgMos->Y.ptr[0]; 197 198 199 if (m_wb.blendingType == BLEND_TYPE_HORZ) 200 { 201 mosaicWidth = cropping_rect.right - cropping_rect.left + 1; 202 mosaicHeight = cropping_rect.bottom - cropping_rect.top + 1; 203 } 204 else 205 { 206 mosaicWidth = Mwidth; 207 mosaicHeight = Mheight; 208 } 209 210 return ret; 211} 212 213 214int Blend::FillFramePyramid(MosaicFrame *mb) 215{ 216 ImageType mbY, mbU, mbV; 217 // Lay this image, centered into the temporary buffer 218 mbY = mb->image; 219 mbU = mb->getU(); 220 mbV = mb->getV(); 221 222 int h, w; 223 224 for(h=0; h<height; h++) 225 { 226 ImageTypeShort yptr = m_pFrameYPyr->ptr[h]; 227 ImageTypeShort uptr = m_pFrameUPyr->ptr[h]; 228 ImageTypeShort vptr = m_pFrameVPyr->ptr[h]; 229 230 for(w=0; w<width; w++) 231 { 232 yptr[w] = (short) ((*(mbY++)) << 3); 233 uptr[w] = (short) ((*(mbU++)) << 3); 234 vptr[w] = (short) ((*(mbV++)) << 3); 235 } 236 } 237 238 // Spread the image through the border 239 PyramidShort::BorderSpread(m_pFrameYPyr, BORDER, BORDER, BORDER, BORDER); 240 PyramidShort::BorderSpread(m_pFrameUPyr, BORDER, BORDER, BORDER, BORDER); 241 PyramidShort::BorderSpread(m_pFrameVPyr, BORDER, BORDER, BORDER, BORDER); 242 243 // Generate Laplacian pyramids 244 if (!PyramidShort::BorderReduce(m_pFrameYPyr, m_wb.nlevs) || !PyramidShort::BorderExpand(m_pFrameYPyr, m_wb.nlevs, -1) || 245 !PyramidShort::BorderReduce(m_pFrameUPyr, m_wb.nlevsC) || !PyramidShort::BorderExpand(m_pFrameUPyr, m_wb.nlevsC, -1) || 246 !PyramidShort::BorderReduce(m_pFrameVPyr, m_wb.nlevsC) || !PyramidShort::BorderExpand(m_pFrameVPyr, m_wb.nlevsC, -1)) 247 { 248 LOGIE("Error: Could not generate Laplacian pyramids\n"); 249 return BLEND_RET_ERROR; 250 } 251 else 252 { 253 return BLEND_RET_OK; 254 } 255} 256 257int Blend::DoMergeAndBlend(MosaicFrame **frames, int nsite, 258 int width, int height, YUVinfo &imgMos, MosaicRect &rect, 259 MosaicRect &cropping_rect, float &progress, bool &cancelComputation) 260{ 261 m_pMosaicYPyr = NULL; 262 m_pMosaicUPyr = NULL; 263 m_pMosaicVPyr = NULL; 264 265 m_pMosaicYPyr = PyramidShort::allocatePyramidPacked(m_wb.nlevs,(unsigned short)rect.Width(),(unsigned short)rect.Height(),BORDER); 266 m_pMosaicUPyr = PyramidShort::allocatePyramidPacked(m_wb.nlevsC,(unsigned short)rect.Width(),(unsigned short)rect.Height(),BORDER); 267 m_pMosaicVPyr = PyramidShort::allocatePyramidPacked(m_wb.nlevsC,(unsigned short)rect.Width(),(unsigned short)rect.Height(),BORDER); 268 if (!m_pMosaicYPyr || !m_pMosaicUPyr || !m_pMosaicVPyr) 269 { 270 LOGIE("Error: Could not allocate pyramids for blending\n"); 271 return BLEND_RET_ERROR_MEMORY; 272 } 273 274 MosaicFrame *mb; 275 276 CSite *esite = m_AllSites + nsite; 277 int site_idx; 278 279 // First go through each frame and for each mosaic pixel determine which frame it should come from 280 site_idx = 0; 281 for(CSite *csite = m_AllSites; csite < esite; csite++) 282 { 283 if(cancelComputation) 284 { 285 if (m_pMosaicVPyr) free(m_pMosaicVPyr); 286 if (m_pMosaicUPyr) free(m_pMosaicUPyr); 287 if (m_pMosaicYPyr) free(m_pMosaicYPyr); 288 return BLEND_RET_CANCELLED; 289 } 290 291 mb = csite->getMb(); 292 293 mb->vcrect = mb->brect; 294 ClipBlendRect(csite, mb->vcrect); 295 296 ComputeMask(csite, mb->vcrect, mb->brect, rect, imgMos, site_idx); 297 298 site_idx++; 299 } 300 301 // Now perform the actual blending using the frame assignment determined above 302 site_idx = 0; 303 for(CSite *csite = m_AllSites; csite < esite; csite++) 304 { 305 if(cancelComputation) 306 { 307 if (m_pMosaicVPyr) free(m_pMosaicVPyr); 308 if (m_pMosaicUPyr) free(m_pMosaicUPyr); 309 if (m_pMosaicYPyr) free(m_pMosaicYPyr); 310 return BLEND_RET_CANCELLED; 311 } 312 313 mb = csite->getMb(); 314 315 316 if(FillFramePyramid(mb)!=BLEND_RET_OK) 317 return BLEND_RET_ERROR; 318 319 ProcessPyramidForThisFrame(csite, mb->vcrect, mb->brect, rect, imgMos, mb->trs, site_idx); 320 321 progress += TIME_PERCENT_BLEND/nsite; 322 323 site_idx++; 324 } 325 326 327 // Blend 328 PerformFinalBlending(imgMos, cropping_rect); 329 330 if (m_pMosaicVPyr) free(m_pMosaicVPyr); 331 if (m_pMosaicUPyr) free(m_pMosaicUPyr); 332 if (m_pMosaicYPyr) free(m_pMosaicYPyr); 333 334 progress += TIME_PERCENT_FINAL; 335 336 return BLEND_RET_OK; 337} 338 339void Blend::CropFinalMosaic(YUVinfo &imgMos, MosaicRect &cropping_rect) 340{ 341 int i, j, k; 342 ImageType yimg; 343 ImageType uimg; 344 ImageType vimg; 345 346 347 yimg = imgMos.Y.ptr[0]; 348 uimg = imgMos.U.ptr[0]; 349 vimg = imgMos.V.ptr[0]; 350 351 k = 0; 352 for (j = cropping_rect.top; j <= cropping_rect.bottom; j++) 353 { 354 for (i = cropping_rect.left; i <= cropping_rect.right; i++) 355 { 356 yimg[k] = yimg[j*imgMos.Y.width+i]; 357 k++; 358 } 359 } 360 for (j = cropping_rect.top; j <= cropping_rect.bottom; j++) 361 { 362 for (i = cropping_rect.left; i <= cropping_rect.right; i++) 363 { 364 yimg[k] = vimg[j*imgMos.Y.width+i]; 365 k++; 366 } 367 } 368 for (j = cropping_rect.top; j <= cropping_rect.bottom; j++) 369 { 370 for (i = cropping_rect.left; i <= cropping_rect.right; i++) 371 { 372 yimg[k] = uimg[j*imgMos.Y.width+i]; 373 k++; 374 } 375 } 376} 377 378int Blend::PerformFinalBlending(YUVinfo &imgMos, MosaicRect &cropping_rect) 379{ 380 if (!PyramidShort::BorderExpand(m_pMosaicYPyr, m_wb.nlevs, 1) || !PyramidShort::BorderExpand(m_pMosaicUPyr, m_wb.nlevsC, 1) || 381 !PyramidShort::BorderExpand(m_pMosaicVPyr, m_wb.nlevsC, 1)) 382 { 383 LOGIE("Error: Could not BorderExpand!\n"); 384 return BLEND_RET_ERROR; 385 } 386 387 ImageTypeShort myimg; 388 ImageTypeShort muimg; 389 ImageTypeShort mvimg; 390 ImageType yimg; 391 ImageType uimg; 392 ImageType vimg; 393 394 int cx = (int)imgMos.Y.width/2; 395 int cy = (int)imgMos.Y.height/2; 396 397 cropping_rect.left = 0; 398 cropping_rect.right = imgMos.Y.width-1; 399 cropping_rect.top = 0; 400 cropping_rect.bottom = imgMos.Y.height-1; 401 402 // Copy the resulting image into the full image using the mask 403 int i, j; 404 405 yimg = imgMos.Y.ptr[0]; 406 uimg = imgMos.U.ptr[0]; 407 vimg = imgMos.V.ptr[0]; 408 409 for (j = 0; j < imgMos.Y.height; j++) 410 { 411 myimg = m_pMosaicYPyr->ptr[j]; 412 muimg = m_pMosaicUPyr->ptr[j]; 413 mvimg = m_pMosaicVPyr->ptr[j]; 414 415 for (i = imgMos.Y.width; i--;) 416 { 417 // A final mask was set up previously, 418 // if the value is zero skip it, otherwise replace it. 419 if (*yimg <255) 420 { 421 short value = (short) ((*myimg) >> 3); 422 if (value < 0) value = 0; 423 else if (value > 255) value = 255; 424 *yimg = (unsigned char) value; 425 426 value = (short) ((*muimg) >> 3); 427 if (value < 0) value = 0; 428 else if (value > 255) value = 255; 429 *uimg = (unsigned char) value; 430 431 value = (short) ((*mvimg) >> 3); 432 if (value < 0) value = 0; 433 else if (value > 255) value = 255; 434 *vimg = (unsigned char) value; 435 436 } 437 else 438 { // set border color in here 439 *yimg = (unsigned char) 96; 440 } 441 442 yimg++; 443 uimg++; 444 vimg++; 445 myimg++; 446 muimg++; 447 mvimg++; 448 } 449 } 450 451 yimg = imgMos.Y.ptr[0]; 452 uimg = imgMos.U.ptr[0]; 453 vimg = imgMos.V.ptr[0]; 454 455 int mincol = 0; 456 int maxcol = imgMos.Y.width-1; 457 458 const int MIN_X_FS = width*1/8; 459 const int MAX_X_FS = imgMos.Y.width - width*1/8; 460 461 for (j = 0; j < imgMos.Y.height; j++) 462 { 463 int minx = MIN_X_FS; 464 int maxx = MAX_X_FS; 465 466 for (i = 0; i < imgMos.Y.width; i++) 467 { 468 if(*yimg == 96 && *uimg == 128 && *vimg == 128) 469 { 470 } 471 else 472 { 473 if(i<minx) 474 minx = i; 475 if(i>maxx) 476 maxx = i; 477 } 478 479 yimg++; 480 uimg++; 481 vimg++; 482 } 483 484 // If we never touched the values, reset them to the image limits 485 if(minx == MIN_X_FS) 486 minx = 0; 487 if(maxx == MAX_X_FS) 488 maxx = imgMos.Y.width-1; 489 490 if(minx>mincol) 491 mincol = minx; 492 if(maxx<maxcol) 493 maxcol = maxx; 494 } 495 496 cropping_rect.left = mincol; 497 cropping_rect.right = maxcol; 498 499 // Crop rows 500 yimg = imgMos.Y.ptr[0]; 501 uimg = imgMos.U.ptr[0]; 502 vimg = imgMos.V.ptr[0]; 503 504 int minrow = 0; 505 int maxrow = imgMos.Y.height-1; 506 507 const int MIN_Y_FS = height*1/8; 508 const int MAX_Y_FS = imgMos.Y.height - height*1/8; 509 510 for (i = 0; i < imgMos.Y.width; i++) 511 { 512 int miny = MIN_Y_FS; 513 int maxy = MAX_Y_FS; 514 515 for (j = 0; j < imgMos.Y.height; j++) 516 { 517 if(*yimg == 96 && *uimg == 128 && *vimg == 128) 518 { 519 } 520 else 521 { 522 if(j<miny) 523 miny = j; 524 if(j>maxy) 525 maxy = j; 526 } 527 528 yimg++; 529 uimg++; 530 vimg++; 531 } 532 533 // If we never touched the values, reset them to the image limits 534 if(miny == MIN_Y_FS) 535 miny = 0; 536 if(maxy == MAX_Y_FS) 537 maxy = imgMos.Y.height-1; 538 539 if(miny>minrow) 540 minrow = miny; 541 if(maxy<maxrow) 542 maxrow = maxy; 543 } 544 545 cropping_rect.top = minrow; 546 cropping_rect.bottom = maxrow; 547 548 return BLEND_RET_OK; 549} 550 551void Blend::ComputeMask(CSite *csite, BlendRect &vcrect, BlendRect &brect, MosaicRect &rect, YUVinfo &imgMos, int site_idx) 552{ 553 PyramidShort *dptr = m_pMosaicYPyr; 554 555 int nC = m_wb.nlevsC; 556 int l = (int) ((vcrect.lft - rect.left)); 557 int b = (int) ((vcrect.bot - rect.top)); 558 int r = (int) ((vcrect.rgt - rect.left)); 559 int t = (int) ((vcrect.top - rect.top)); 560 561 if (vcrect.lft == brect.lft) 562 l = (l <= 0) ? -BORDER : l - BORDER; 563 else if (l < -BORDER) 564 l = -BORDER; 565 566 if (vcrect.bot == brect.bot) 567 b = (b <= 0) ? -BORDER : b - BORDER; 568 else if (b < -BORDER) 569 b = -BORDER; 570 571 if (vcrect.rgt == brect.rgt) 572 r = (r >= dptr->width) ? dptr->width + BORDER - 1 : r + BORDER; 573 else if (r >= dptr->width + BORDER) 574 r = dptr->width + BORDER - 1; 575 576 if (vcrect.top == brect.top) 577 t = (t >= dptr->height) ? dptr->height + BORDER - 1 : t + BORDER; 578 else if (t >= dptr->height + BORDER) 579 t = dptr->height + BORDER - 1; 580 581 // Walk the Region of interest and populate the pyramid 582 for (int j = b; j <= t; j++) 583 { 584 int jj = j; 585 double sj = jj + rect.top; 586 587 for (int i = l; i <= r; i++) 588 { 589 int ii = i; 590 // project point and then triangulate to neighbors 591 double si = ii + rect.left; 592 593 double dself = hypotSq(csite->getVCenter().x - si, csite->getVCenter().y - sj); 594 int inMask = ((unsigned) ii < imgMos.Y.width && 595 (unsigned) jj < imgMos.Y.height) ? 1 : 0; 596 597 if(!inMask) 598 continue; 599 600 // scan the neighbors to see if this is a valid position 601 unsigned char mask = (unsigned char) 255; 602 SEdgeVector *ce; 603 int ecnt; 604 for (ce = csite->getNeighbor(), ecnt = csite->getNumNeighbors(); ecnt--; ce++) 605 { 606 double d1 = hypotSq(m_AllSites[ce->second].getVCenter().x - si, 607 m_AllSites[ce->second].getVCenter().y - sj); 608 if (d1 < dself) 609 { 610 break; 611 } 612 } 613 614 if (ecnt >= 0) continue; 615 616 imgMos.Y.ptr[jj][ii] = (unsigned char)site_idx; 617 } 618 } 619} 620 621void Blend::ProcessPyramidForThisFrame(CSite *csite, BlendRect &vcrect, BlendRect &brect, MosaicRect &rect, YUVinfo &imgMos, double trs[3][3], int site_idx) 622{ 623 // Put the Region of interest (for all levels) into m_pMosaicYPyr 624 double inv_trs[3][3]; 625 inv33d(trs, inv_trs); 626 627 // Process each pyramid level 628 PyramidShort *sptr = m_pFrameYPyr; 629 PyramidShort *suptr = m_pFrameUPyr; 630 PyramidShort *svptr = m_pFrameVPyr; 631 632 PyramidShort *dptr = m_pMosaicYPyr; 633 PyramidShort *duptr = m_pMosaicUPyr; 634 PyramidShort *dvptr = m_pMosaicVPyr; 635 636 int dscale = 0; // distance scale for the current level 637 int nC = m_wb.nlevsC; 638 for (int n = m_wb.nlevs; n--; dscale++, dptr++, sptr++, dvptr++, duptr++, svptr++, suptr++, nC--) 639 { 640 int l = (int) ((vcrect.lft - rect.left) / (1 << dscale)); 641 int b = (int) ((vcrect.bot - rect.top) / (1 << dscale)); 642 int r = (int) ((vcrect.rgt - rect.left) / (1 << dscale) + .5); 643 int t = (int) ((vcrect.top - rect.top) / (1 << dscale) + .5); 644 645 if (vcrect.lft == brect.lft) 646 l = (l <= 0) ? -BORDER : l - BORDER; 647 else if (l < -BORDER) 648 l = -BORDER; 649 650 if (vcrect.bot == brect.bot) 651 b = (b <= 0) ? -BORDER : b - BORDER; 652 else if (b < -BORDER) 653 b = -BORDER; 654 655 if (vcrect.rgt == brect.rgt) 656 r = (r >= dptr->width) ? dptr->width + BORDER - 1 : r + BORDER; 657 else if (r >= dptr->width + BORDER) 658 r = dptr->width + BORDER - 1; 659 660 if (vcrect.top == brect.top) 661 t = (t >= dptr->height) ? dptr->height + BORDER - 1 : t + BORDER; 662 else if (t >= dptr->height + BORDER) 663 t = dptr->height + BORDER - 1; 664 665 // Walk the Region of interest and populate the pyramid 666 for (int j = b; j <= t; j++) 667 { 668 int jj = (j << dscale); 669 double sj = jj + rect.top; 670 671 for (int i = l; i <= r; i++) 672 { 673 int ii = (i << dscale); 674 // project point and then triangulate to neighbors 675 double si = ii + rect.left; 676 677 int inMask = ((unsigned) ii < imgMos.Y.width && 678 (unsigned) jj < imgMos.Y.height) ? 1 : 0; 679 680 if(inMask && imgMos.Y.ptr[jj][ii]!=site_idx && imgMos.Y.ptr[jj][ii]!=255) 681 continue; 682 683 // Project this mosaic point into the original frame coordinate space 684 double xx, yy; 685 686 MosaicToFrame(inv_trs, si, sj, xx, yy); 687 688 if (xx < 0.0 || yy < 0.0 || xx > width - 1.0 || yy > height - 1.0) 689 { 690 if(inMask) 691 { 692 imgMos.Y.ptr[jj][ii] = 255; 693 } 694 } 695 696 xx /= (1 << dscale); 697 yy /= (1 << dscale); 698 699 700 int x1 = (xx >= 0.0) ? (int) xx : (int) floor(xx); 701 int y1 = (yy >= 0.0) ? (int) yy : (int) floor(yy); 702 703 // Final destination in extended pyramid 704#ifndef LINEAR_INTERP 705 if(inSegment(x1, sptr->width, BORDER-1) && inSegment(y1, sptr->height, BORDER-1)) 706 { 707 double xfrac = xx - x1; 708 double yfrac = yy - y1; 709 dptr->ptr[j][i] = (short) (.5 + ciCalc(sptr, x1, y1, xfrac, yfrac)); 710 if (dvptr >= m_pMosaicVPyr && nC > 0) 711 { 712 duptr->ptr[j][i] = (short) (.5 + ciCalc(suptr, x1, y1, xfrac, yfrac)); 713 dvptr->ptr[j][i] = (short) (.5 + ciCalc(svptr, x1, y1, xfrac, yfrac)); 714 } 715 } 716#else 717 if(inSegment(x1, sptr->width, BORDER) && inSegment(y1, sptr->height, BORDER)) 718 { 719 int x2 = x1 + 1; 720 int y2 = y1 + 1; 721 double xfrac = xx - x1; 722 double yfrac = yy - y1; 723 double y1val = sptr->ptr[y1][x1] + 724 (sptr->ptr[y1][x2] - sptr->ptr[y1][x1]) * xfrac; 725 double y2val = sptr->ptr[y2][x1] + 726 (sptr->ptr[y2][x2] - sptr->ptr[y2][x1]) * xfrac; 727 dptr->ptr[j][i] = (short) (y1val + yfrac * (y2val - y1val)); 728 729 if (dvptr >= m_pMosaicVPyr && nC > 0) 730 { 731 y1val = suptr->ptr[y1][x1] + 732 (suptr->ptr[y1][x2] - suptr->ptr[y1][x1]) * xfrac; 733 y2val = suptr->ptr[y2][x1] + 734 (suptr->ptr[y2][x2] - suptr->ptr[y2][x1]) * xfrac; 735 736 duptr->ptr[j][i] = (short) (y1val + yfrac * (y2val - y1val)); 737 738 y1val = svptr->ptr[y1][x1] + 739 (svptr->ptr[y1][x2] - svptr->ptr[y1][x1]) * xfrac; 740 y2val = svptr->ptr[y2][x1] + 741 (svptr->ptr[y2][x2] - svptr->ptr[y2][x1]) * xfrac; 742 743 dvptr->ptr[j][i] = (short) (y1val + yfrac * (y2val - y1val)); 744 } 745 } 746#endif 747 else 748 { 749 clipToSegment(x1, sptr->width, BORDER); 750 clipToSegment(y1, sptr->height, BORDER); 751 752 dptr->ptr[j][i] = sptr->ptr[y1][x1]; 753 if (dvptr >= m_pMosaicVPyr && nC > 0) 754 { 755 dvptr->ptr[j][i] = svptr->ptr[y1][x1]; 756 duptr->ptr[j][i] = suptr->ptr[y1][x1]; 757 } 758 } 759 } 760 } 761 } 762} 763 764void Blend::MosaicToFrame(double trs[3][3], double x, double y, double &wx, double &wy) 765{ 766 double X, Y, z; 767 if (m_wb.theta == 0.0) 768 { 769 X = x; 770 Y = y; 771 } 772 else if (m_wb.horizontal) 773 { 774 double alpha = x * m_wb.direction / m_wb.width; 775 double length = (y - alpha * m_wb.correction) * m_wb.direction + m_wb.radius; 776 double deltaTheta = m_wb.theta * alpha; 777 double sinTheta = sin(deltaTheta); 778 double cosTheta = sqrt(1.0 - sinTheta * sinTheta) * m_wb.direction; 779 X = length * sinTheta + m_wb.x; 780 Y = length * cosTheta + m_wb.y; 781 } 782 else 783 { 784 double alpha = y * m_wb.direction / m_wb.width; 785 double length = (x - alpha * m_wb.correction) * m_wb.direction + m_wb.radius; 786 double deltaTheta = m_wb.theta * alpha; 787 double sinTheta = sin(deltaTheta); 788 double cosTheta = sqrt(1.0 - sinTheta * sinTheta) * m_wb.direction; 789 Y = length * sinTheta + m_wb.y; 790 X = length * cosTheta + m_wb.x; 791 } 792 z = ProjZ(trs, X, Y, 1.0); 793 wx = ProjX(trs, X, Y, z, 1.0); 794 wy = ProjY(trs, X, Y, z, 1.0); 795} 796 797void Blend::FrameToMosaic(double trs[3][3], double x, double y, double &wx, double &wy) 798{ 799 // Project into the intermediate Mosaic coordinate system 800 double z = ProjZ(trs, x, y, 1.0); 801 double X = ProjX(trs, x, y, z, 1.0); 802 double Y = ProjY(trs, x, y, z, 1.0); 803 804 if (m_wb.theta == 0.0) 805 { 806 // No rotation, then this is all we need to do. 807 wx = X; 808 wy = Y; 809 } 810 else if (m_wb.horizontal) 811 { 812 double deltaX = X - m_wb.x; 813 double deltaY = Y - m_wb.y; 814 double length = sqrt(deltaX * deltaX + deltaY * deltaY); 815 double deltaTheta = asin(deltaX / length); 816 double alpha = deltaTheta / m_wb.theta; 817 wx = alpha * m_wb.width * m_wb.direction; 818 wy = (length - m_wb.radius) * m_wb.direction + alpha * m_wb.correction; 819 } 820 else 821 { 822 double deltaX = X - m_wb.x; 823 double deltaY = Y - m_wb.y; 824 double length = sqrt(deltaX * deltaX + deltaY * deltaY); 825 double deltaTheta = asin(deltaY / length); 826 double alpha = deltaTheta / m_wb.theta; 827 wy = alpha * m_wb.width * m_wb.direction; 828 wx = (length - m_wb.radius) * m_wb.direction + alpha * m_wb.correction; 829 } 830} 831 832 833 834// Clip the region of interest as small as possible by using the Voronoi edges of 835// the neighbors 836void Blend::ClipBlendRect(CSite *csite, BlendRect &brect) 837{ 838 SEdgeVector *ce; 839 int ecnt; 840 for (ce = csite->getNeighbor(), ecnt = csite->getNumNeighbors(); ecnt--; ce++) 841 { 842 // calculate the Voronoi bisector intersection 843 const double epsilon = 1e-5; 844 double dx = (m_AllSites[ce->second].getVCenter().x - m_AllSites[ce->first].getVCenter().x); 845 double dy = (m_AllSites[ce->second].getVCenter().y - m_AllSites[ce->first].getVCenter().y); 846 double xmid = m_AllSites[ce->first].getVCenter().x + dx/2.0; 847 double ymid = m_AllSites[ce->first].getVCenter().y + dy/2.0; 848 double inter; 849 850 if (dx > epsilon) 851 { 852 // neighbor is on right 853 if ((inter = m_wb.roundoffOverlap + xmid - dy * (((dy >= 0.0) ? brect.bot : brect.top) - ymid) / dx) < brect.rgt) 854 brect.rgt = inter; 855 } 856 else if (dx < -epsilon) 857 { 858 // neighbor is on left 859 if ((inter = -m_wb.roundoffOverlap + xmid - dy * (((dy >= 0.0) ? brect.bot : brect.top) - ymid) / dx) > brect.lft) 860 brect.lft = inter; 861 } 862 if (dy > epsilon) 863 { 864 // neighbor is above 865 if ((inter = m_wb.roundoffOverlap + ymid - dx * (((dx >= 0.0) ? brect.lft : brect.rgt) - xmid) / dy) < brect.top) 866 brect.top = inter; 867 } 868 else if (dy < -epsilon) 869 { 870 // neighbor is below 871 if ((inter = -m_wb.roundoffOverlap + ymid - dx * (((dx >= 0.0) ? brect.lft : brect.rgt) - xmid) / dy) > brect.bot) 872 brect.bot = inter; 873 } 874 } 875} 876 877void Blend::FrameToMosaicRect(int width, int height, double trs[3][3], BlendRect &brect) 878{ 879 // We need to walk the perimeter since the borders can be bent. 880 brect.lft = brect.bot = 2e30; 881 brect.rgt = brect.top = -2e30; 882 double xpos, ypos; 883 double lasty = height - 1.0; 884 double lastx = width - 1.0; 885 int i; 886 887 for (i = width; i--;) 888 { 889 890 FrameToMosaic(trs, (double) i, 0.0, xpos, ypos); 891 ClipRect(xpos, ypos, brect); 892 FrameToMosaic(trs, (double) i, lasty, xpos, ypos); 893 ClipRect(xpos, ypos, brect); 894 } 895 for (i = height; i--;) 896 { 897 FrameToMosaic(trs, 0.0, (double) i, xpos, ypos); 898 ClipRect(xpos, ypos, brect); 899 FrameToMosaic(trs, lastx, (double) i, xpos, ypos); 900 ClipRect(xpos, ypos, brect); 901 } 902} 903 904 905 906void Blend::ComputeBlendParameters(MosaicFrame **frames, int frames_size, int is360) 907{ 908 if (m_wb.blendingType != BLEND_TYPE_CYLPAN && m_wb.blendingType != BLEND_TYPE_HORZ) 909 { 910 m_wb.theta = 0.0; 911 return; 912 } 913 914 MosaicFrame *first = frames[0]; 915 MosaicFrame *last = frames[frames_size-1]; 916 MosaicFrame *mb; 917 918 double lxpos = last->trs[0][2], lypos = last->trs[1][2]; 919 double fxpos = first->trs[0][2], fypos = first->trs[1][2]; 920 921 // Calculate warp to produce proper stitching. 922 // get x, y displacement 923 double midX = last->width / 2.0; 924 double midY = last->height / 2.0; 925 double z = ProjZ(first->trs, midX, midY, 1.0); 926 double firstX, firstY; 927 double prevX = firstX = ProjX(first->trs, midX, midY, z, 1.0); 928 double prevY = firstY = ProjY(first->trs, midX, midY, z, 1.0); 929 930 double arcLength, lastTheta; 931 m_wb.theta = lastTheta = arcLength = 0.0; 932 for (int i = 0; i < frames_size; i++) 933 { 934 mb = frames[i]; 935 double currX, currY; 936 z = ProjZ(mb->trs, midX, midY, 1.0); 937 currX = ProjX(mb->trs, midX, midY, z, 1.0); 938 currY = ProjY(mb->trs, midX, midY, z, 1.0); 939 double deltaX = currX - prevX; 940 double deltaY = currY - prevY; 941 arcLength += sqrt(deltaY * deltaY + deltaX * deltaX); 942 if (!is360) 943 { 944 double thisTheta = asin(mb->trs[1][0]); 945 m_wb.theta += thisTheta - lastTheta; 946 lastTheta = thisTheta; 947 } 948 prevX = currX; 949 prevY = currY; 950 } 951 952 // In case of BMP output, stretch this to end at the proper alignment 953 m_wb.width = arcLength; 954 if (is360) m_wb.theta = asin(last->trs[1][0]); 955 956 // If there is no rotation, we're done. 957 if (m_wb.theta != 0.0) 958 { 959 double dx = prevX - firstX; 960 double dy = prevY - firstY; 961 if (abs(lxpos - fxpos) > abs(lypos - fypos)) 962 { 963 m_wb.horizontal = 1; 964 // Calculate radius position to make ends exactly the same Y offset 965 double radiusTheta = dx / cos(3.14159 / 2.0 - m_wb.theta); 966 m_wb.radius = dy + radiusTheta * cos(m_wb.theta); 967 if (m_wb.radius < 0.0) m_wb.radius = -m_wb.radius; 968 } 969 else 970 { 971 m_wb.horizontal = 0; 972 // Calculate radius position to make ends exactly the same Y offset 973 double radiusTheta = dy / cos(3.14159 / 2.0 - m_wb.theta); 974 m_wb.radius = dx + radiusTheta * cos(m_wb.theta); 975 if (m_wb.radius < 0.0) m_wb.radius = -m_wb.radius; 976 } 977 978 // Determine major direction 979 if (m_wb.horizontal) 980 { 981 // Horizontal strip 982 if (is360) m_wb.x = firstX; 983 else 984 { 985 if (lxpos - fxpos < 0) 986 { 987 m_wb.x = firstX + midX; 988 z = ProjZ(last->trs, 0.0, midY, 1.0); 989 prevX = ProjX(last->trs, 0.0, midY, z, 1.0); 990 prevY = ProjY(last->trs, 0.0, midY, z, 1.0); 991 } 992 else 993 { 994 m_wb.x = firstX - midX; 995 z = ProjZ(last->trs, last->width - 1.0, midY, 1.0); 996 prevX = ProjX(last->trs, last->width - 1.0, midY, z, 1.0); 997 prevY = ProjY(last->trs, last->width - 1.0, midY, z, 1.0); 998 } 999 } 1000 dy = prevY - firstY; 1001 if (dy < 0.0) m_wb.direction = 1.0; 1002 else m_wb.direction = -1.0; 1003 m_wb.y = firstY - m_wb.radius * m_wb.direction; 1004 if (dy * m_wb.theta > 0.0) m_wb.width = -m_wb.width; 1005 } 1006 else 1007 { 1008 // Vertical strip 1009 if (is360) m_wb.y = firstY; 1010 else 1011 { 1012 if (lypos - fypos < 0) 1013 { 1014 m_wb.x = firstY + midY; 1015 z = ProjZ(last->trs, midX, 0.0, 1.0); 1016 prevX = ProjX(last->trs, midX, 0.0, z, 1.0); 1017 prevY = ProjY(last->trs, midX, 0.0, z, 1.0); 1018 } 1019 else 1020 { 1021 m_wb.x = firstX - midX; 1022 z = ProjZ(last->trs, midX, last->height - 1.0, 1.0); 1023 prevX = ProjX(last->trs, midX, last->height - 1.0, z, 1.0); 1024 prevY = ProjY(last->trs, midX, last->height - 1.0, z, 1.0); 1025 } 1026 } 1027 dx = prevX - firstX; 1028 if (dx < 0.0) m_wb.direction = 1.0; 1029 else m_wb.direction = -1.0; 1030 m_wb.x = firstX - m_wb.radius * m_wb.direction; 1031 if (dx * m_wb.theta > 0.0) m_wb.width = -m_wb.width; 1032 } 1033 1034 // Calculate the correct correction factor 1035 double deltaX = prevX - m_wb.x; 1036 double deltaY = prevY - m_wb.y; 1037 double length = sqrt(deltaX * deltaX + deltaY * deltaY); 1038 double deltaTheta = (m_wb.horizontal) ? deltaX : deltaY; 1039 deltaTheta = asin(deltaTheta / length); 1040 m_wb.correction = ((m_wb.radius - length) * m_wb.direction) / 1041 (deltaTheta / m_wb.theta); 1042 } 1043} 1044