Blend.cpp revision 50b3c890986aadb3780b4da8c0b8dbb0f1422eba
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) 95{ 96 int ret; 97 int numCenters; 98 99 ComputeBlendParameters(frames, frames_size, false); 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); 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 BLEND_RET_OK; 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) 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 mb = csite->getMb(); 284 285 mb->vcrect = mb->brect; 286 ClipBlendRect(csite, mb->vcrect); 287 288 ComputeMask(csite, mb->vcrect, mb->brect, rect, imgMos, site_idx); 289 290 site_idx++; 291 } 292 293 // Now perform the actual blending using the frame assignment determined above 294 site_idx = 0; 295 for(CSite *csite = m_AllSites; csite < esite; csite++) 296 { 297 mb = csite->getMb(); 298 299 300 if(FillFramePyramid(mb)!=BLEND_RET_OK) 301 return BLEND_RET_ERROR; 302 303 ProcessPyramidForThisFrame(csite, mb->vcrect, mb->brect, rect, imgMos, mb->trs, site_idx); 304 305 progress += TIME_PERCENT_BLEND/nsite; 306 307 site_idx++; 308 } 309 310 311 // Blend 312 PerformFinalBlending(imgMos, cropping_rect); 313 314 if (m_pMosaicVPyr) free(m_pMosaicVPyr); 315 if (m_pMosaicUPyr) free(m_pMosaicUPyr); 316 if (m_pMosaicYPyr) free(m_pMosaicYPyr); 317 318 progress += TIME_PERCENT_FINAL; 319 320 return BLEND_RET_OK; 321} 322 323void Blend::CropFinalMosaic(YUVinfo &imgMos, MosaicRect &cropping_rect) 324{ 325 int i, j, k; 326 ImageType yimg; 327 ImageType uimg; 328 ImageType vimg; 329 330 331 yimg = imgMos.Y.ptr[0]; 332 uimg = imgMos.U.ptr[0]; 333 vimg = imgMos.V.ptr[0]; 334 335 k = 0; 336 for (j = cropping_rect.top; j <= cropping_rect.bottom; j++) 337 { 338 for (i = cropping_rect.left; i <= cropping_rect.right; i++) 339 { 340 yimg[k] = yimg[j*imgMos.Y.width+i]; 341 k++; 342 } 343 } 344 for (j = cropping_rect.top; j <= cropping_rect.bottom; j++) 345 { 346 for (i = cropping_rect.left; i <= cropping_rect.right; i++) 347 { 348 yimg[k] = vimg[j*imgMos.Y.width+i]; 349 k++; 350 } 351 } 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] = uimg[j*imgMos.Y.width+i]; 357 k++; 358 } 359 } 360} 361 362int Blend::PerformFinalBlending(YUVinfo &imgMos, MosaicRect &cropping_rect) 363{ 364 if (!PyramidShort::BorderExpand(m_pMosaicYPyr, m_wb.nlevs, 1) || !PyramidShort::BorderExpand(m_pMosaicUPyr, m_wb.nlevsC, 1) || 365 !PyramidShort::BorderExpand(m_pMosaicVPyr, m_wb.nlevsC, 1)) 366 { 367 LOGIE("Error: Could not BorderExpand!\n"); 368 return BLEND_RET_ERROR; 369 } 370 371 ImageTypeShort myimg; 372 ImageTypeShort muimg; 373 ImageTypeShort mvimg; 374 ImageType yimg; 375 ImageType uimg; 376 ImageType vimg; 377 378 int cx = (int)imgMos.Y.width/2; 379 int cy = (int)imgMos.Y.height/2; 380 381 cropping_rect.left = 0; 382 cropping_rect.right = imgMos.Y.width-1; 383 cropping_rect.top = 0; 384 cropping_rect.bottom = imgMos.Y.height-1; 385 386 // Copy the resulting image into the full image using the mask 387 int i, j; 388 389 yimg = imgMos.Y.ptr[0]; 390 uimg = imgMos.U.ptr[0]; 391 vimg = imgMos.V.ptr[0]; 392 393 for (j = 0; j < imgMos.Y.height; j++) 394 { 395 myimg = m_pMosaicYPyr->ptr[j]; 396 muimg = m_pMosaicUPyr->ptr[j]; 397 mvimg = m_pMosaicVPyr->ptr[j]; 398 399 for (i = imgMos.Y.width; i--;) 400 { 401 // A final mask was set up previously, 402 // if the value is zero skip it, otherwise replace it. 403 if (*yimg <255) 404 { 405 short value = (short) ((*myimg) >> 3); 406 if (value < 0) value = 0; 407 else if (value > 255) value = 255; 408 *yimg = (unsigned char) value; 409 410 value = (short) ((*muimg) >> 3); 411 if (value < 0) value = 0; 412 else if (value > 255) value = 255; 413 *uimg = (unsigned char) value; 414 415 value = (short) ((*mvimg) >> 3); 416 if (value < 0) value = 0; 417 else if (value > 255) value = 255; 418 *vimg = (unsigned char) value; 419 420 } 421 else 422 { // set border color in here 423 *yimg = (unsigned char) 96; 424 } 425 426 yimg++; 427 uimg++; 428 vimg++; 429 myimg++; 430 muimg++; 431 mvimg++; 432 } 433 } 434 435 yimg = imgMos.Y.ptr[0]; 436 uimg = imgMos.U.ptr[0]; 437 vimg = imgMos.V.ptr[0]; 438 439 int mincol = 0; 440 int maxcol = imgMos.Y.width-1; 441 442 const int MIN_X_FS = width*1/8; 443 const int MAX_X_FS = imgMos.Y.width - width*1/8; 444 445 for (j = 0; j < imgMos.Y.height; j++) 446 { 447 int minx = MIN_X_FS; 448 int maxx = MAX_X_FS; 449 450 for (i = 0; i < imgMos.Y.width; i++) 451 { 452 if(*yimg == 96 && *uimg == 128 && *vimg == 128) 453 { 454 } 455 else 456 { 457 if(i<minx) 458 minx = i; 459 if(i>maxx) 460 maxx = i; 461 } 462 463 yimg++; 464 uimg++; 465 vimg++; 466 } 467 468 // If we never touched the values, reset them to the image limits 469 if(minx == MIN_X_FS) 470 minx = 0; 471 if(maxx == MAX_X_FS) 472 maxx = imgMos.Y.width-1; 473 474 if(minx>mincol) 475 mincol = minx; 476 if(maxx<maxcol) 477 maxcol = maxx; 478 } 479 480 cropping_rect.left = mincol; 481 cropping_rect.right = maxcol; 482 483 // Crop rows 484 yimg = imgMos.Y.ptr[0]; 485 uimg = imgMos.U.ptr[0]; 486 vimg = imgMos.V.ptr[0]; 487 488 int minrow = 0; 489 int maxrow = imgMos.Y.height-1; 490 491 const int MIN_Y_FS = height*1/8; 492 const int MAX_Y_FS = imgMos.Y.height - height*1/8; 493 494 for (i = 0; i < imgMos.Y.width; i++) 495 { 496 int miny = MIN_Y_FS; 497 int maxy = MAX_Y_FS; 498 499 for (j = 0; j < imgMos.Y.height; j++) 500 { 501 if(*yimg == 96 && *uimg == 128 && *vimg == 128) 502 { 503 } 504 else 505 { 506 if(j<miny) 507 miny = j; 508 if(j>maxy) 509 maxy = j; 510 } 511 512 yimg++; 513 uimg++; 514 vimg++; 515 } 516 517 // If we never touched the values, reset them to the image limits 518 if(miny == MIN_Y_FS) 519 miny = 0; 520 if(maxy == MAX_Y_FS) 521 maxy = imgMos.Y.height-1; 522 523 if(miny>minrow) 524 minrow = miny; 525 if(maxy<maxrow) 526 maxrow = maxy; 527 } 528 529 cropping_rect.top = minrow; 530 cropping_rect.bottom = maxrow; 531 532 return BLEND_RET_OK; 533} 534 535void Blend::ComputeMask(CSite *csite, BlendRect &vcrect, BlendRect &brect, MosaicRect &rect, YUVinfo &imgMos, int site_idx) 536{ 537 PyramidShort *dptr = m_pMosaicYPyr; 538 539 int nC = m_wb.nlevsC; 540 int l = (int) ((vcrect.lft - rect.left)); 541 int b = (int) ((vcrect.bot - rect.top)); 542 int r = (int) ((vcrect.rgt - rect.left)); 543 int t = (int) ((vcrect.top - rect.top)); 544 545 if (vcrect.lft == brect.lft) 546 l = (l <= 0) ? -BORDER : l - BORDER; 547 else if (l < -BORDER) 548 l = -BORDER; 549 550 if (vcrect.bot == brect.bot) 551 b = (b <= 0) ? -BORDER : b - BORDER; 552 else if (b < -BORDER) 553 b = -BORDER; 554 555 if (vcrect.rgt == brect.rgt) 556 r = (r >= dptr->width) ? dptr->width + BORDER - 1 : r + BORDER; 557 else if (r >= dptr->width + BORDER) 558 r = dptr->width + BORDER - 1; 559 560 if (vcrect.top == brect.top) 561 t = (t >= dptr->height) ? dptr->height + BORDER - 1 : t + BORDER; 562 else if (t >= dptr->height + BORDER) 563 t = dptr->height + BORDER - 1; 564 565 // Walk the Region of interest and populate the pyramid 566 for (int j = b; j <= t; j++) 567 { 568 int jj = j; 569 double sj = jj + rect.top; 570 571 for (int i = l; i <= r; i++) 572 { 573 int ii = i; 574 // project point and then triangulate to neighbors 575 double si = ii + rect.left; 576 577 double dself = hypotSq(csite->getVCenter().x - si, csite->getVCenter().y - sj); 578 int inMask = ((unsigned) ii < imgMos.Y.width && 579 (unsigned) jj < imgMos.Y.height) ? 1 : 0; 580 581 if(!inMask) 582 continue; 583 584 // scan the neighbors to see if this is a valid position 585 unsigned char mask = (unsigned char) 255; 586 SEdgeVector *ce; 587 int ecnt; 588 for (ce = csite->getNeighbor(), ecnt = csite->getNumNeighbors(); ecnt--; ce++) 589 { 590 double d1 = hypotSq(m_AllSites[ce->second].getVCenter().x - si, 591 m_AllSites[ce->second].getVCenter().y - sj); 592 if (d1 < dself) 593 { 594 break; 595 } 596 } 597 598 if (ecnt >= 0) continue; 599 600 imgMos.Y.ptr[jj][ii] = (unsigned char)site_idx; 601 } 602 } 603} 604 605void Blend::ProcessPyramidForThisFrame(CSite *csite, BlendRect &vcrect, BlendRect &brect, MosaicRect &rect, YUVinfo &imgMos, double trs[3][3], int site_idx) 606{ 607 // Put the Region of interest (for all levels) into m_pMosaicYPyr 608 double inv_trs[3][3]; 609 inv33d(trs, inv_trs); 610 611 // Process each pyramid level 612 PyramidShort *sptr = m_pFrameYPyr; 613 PyramidShort *suptr = m_pFrameUPyr; 614 PyramidShort *svptr = m_pFrameVPyr; 615 616 PyramidShort *dptr = m_pMosaicYPyr; 617 PyramidShort *duptr = m_pMosaicUPyr; 618 PyramidShort *dvptr = m_pMosaicVPyr; 619 620 int dscale = 0; // distance scale for the current level 621 int nC = m_wb.nlevsC; 622 for (int n = m_wb.nlevs; n--; dscale++, dptr++, sptr++, dvptr++, duptr++, svptr++, suptr++, nC--) 623 { 624 int l = (int) ((vcrect.lft - rect.left) / (1 << dscale)); 625 int b = (int) ((vcrect.bot - rect.top) / (1 << dscale)); 626 int r = (int) ((vcrect.rgt - rect.left) / (1 << dscale) + .5); 627 int t = (int) ((vcrect.top - rect.top) / (1 << dscale) + .5); 628 629 if (vcrect.lft == brect.lft) 630 l = (l <= 0) ? -BORDER : l - BORDER; 631 else if (l < -BORDER) 632 l = -BORDER; 633 634 if (vcrect.bot == brect.bot) 635 b = (b <= 0) ? -BORDER : b - BORDER; 636 else if (b < -BORDER) 637 b = -BORDER; 638 639 if (vcrect.rgt == brect.rgt) 640 r = (r >= dptr->width) ? dptr->width + BORDER - 1 : r + BORDER; 641 else if (r >= dptr->width + BORDER) 642 r = dptr->width + BORDER - 1; 643 644 if (vcrect.top == brect.top) 645 t = (t >= dptr->height) ? dptr->height + BORDER - 1 : t + BORDER; 646 else if (t >= dptr->height + BORDER) 647 t = dptr->height + BORDER - 1; 648 649 // Walk the Region of interest and populate the pyramid 650 for (int j = b; j <= t; j++) 651 { 652 int jj = (j << dscale); 653 double sj = jj + rect.top; 654 655 for (int i = l; i <= r; i++) 656 { 657 int ii = (i << dscale); 658 // project point and then triangulate to neighbors 659 double si = ii + rect.left; 660 661 int inMask = ((unsigned) ii < imgMos.Y.width && 662 (unsigned) jj < imgMos.Y.height) ? 1 : 0; 663 664 if(inMask && imgMos.Y.ptr[jj][ii]!=site_idx && imgMos.Y.ptr[jj][ii]!=255) 665 continue; 666 667 // Project this mosaic point into the original frame coordinate space 668 double xx, yy; 669 670 MosaicToFrame(inv_trs, si, sj, xx, yy); 671 672 if (xx < 0.0 || yy < 0.0 || xx > width - 1.0 || yy > height - 1.0) 673 { 674 if(inMask) 675 { 676 imgMos.Y.ptr[jj][ii] = 255; 677 } 678 } 679 680 xx /= (1 << dscale); 681 yy /= (1 << dscale); 682 683 684 int x1 = (xx >= 0.0) ? (int) xx : (int) floor(xx); 685 int y1 = (yy >= 0.0) ? (int) yy : (int) floor(yy); 686 687 // Final destination in extended pyramid 688#ifndef LINEAR_INTERP 689 if(inSegment(x1, sptr->width, BORDER-1) && inSegment(y1, sptr->height, BORDER-1)) 690 { 691 double xfrac = xx - x1; 692 double yfrac = yy - y1; 693 dptr->ptr[j][i] = (short) (.5 + ciCalc(sptr, x1, y1, xfrac, yfrac)); 694 if (dvptr >= m_pMosaicVPyr && nC > 0) 695 { 696 duptr->ptr[j][i] = (short) (.5 + ciCalc(suptr, x1, y1, xfrac, yfrac)); 697 dvptr->ptr[j][i] = (short) (.5 + ciCalc(svptr, x1, y1, xfrac, yfrac)); 698 } 699 } 700#else 701 if(inSegment(x1, sptr->width, BORDER) && inSegment(y1, sptr->height, BORDER)) 702 { 703 int x2 = x1 + 1; 704 int y2 = y1 + 1; 705 double xfrac = xx - x1; 706 double yfrac = yy - y1; 707 double y1val = sptr->ptr[y1][x1] + 708 (sptr->ptr[y1][x2] - sptr->ptr[y1][x1]) * xfrac; 709 double y2val = sptr->ptr[y2][x1] + 710 (sptr->ptr[y2][x2] - sptr->ptr[y2][x1]) * xfrac; 711 dptr->ptr[j][i] = (short) (y1val + yfrac * (y2val - y1val)); 712 713 if (dvptr >= m_pMosaicVPyr && nC > 0) 714 { 715 y1val = suptr->ptr[y1][x1] + 716 (suptr->ptr[y1][x2] - suptr->ptr[y1][x1]) * xfrac; 717 y2val = suptr->ptr[y2][x1] + 718 (suptr->ptr[y2][x2] - suptr->ptr[y2][x1]) * xfrac; 719 720 duptr->ptr[j][i] = (short) (y1val + yfrac * (y2val - y1val)); 721 722 y1val = svptr->ptr[y1][x1] + 723 (svptr->ptr[y1][x2] - svptr->ptr[y1][x1]) * xfrac; 724 y2val = svptr->ptr[y2][x1] + 725 (svptr->ptr[y2][x2] - svptr->ptr[y2][x1]) * xfrac; 726 727 dvptr->ptr[j][i] = (short) (y1val + yfrac * (y2val - y1val)); 728 } 729 } 730#endif 731 else 732 { 733 clipToSegment(x1, sptr->width, BORDER); 734 clipToSegment(y1, sptr->height, BORDER); 735 736 dptr->ptr[j][i] = sptr->ptr[y1][x1]; 737 if (dvptr >= m_pMosaicVPyr && nC > 0) 738 { 739 dvptr->ptr[j][i] = svptr->ptr[y1][x1]; 740 duptr->ptr[j][i] = suptr->ptr[y1][x1]; 741 } 742 } 743 } 744 } 745 } 746} 747 748void Blend::MosaicToFrame(double trs[3][3], double x, double y, double &wx, double &wy) 749{ 750 double X, Y, z; 751 if (m_wb.theta == 0.0) 752 { 753 X = x; 754 Y = y; 755 } 756 else if (m_wb.horizontal) 757 { 758 double alpha = x * m_wb.direction / m_wb.width; 759 double length = (y - alpha * m_wb.correction) * m_wb.direction + m_wb.radius; 760 double deltaTheta = m_wb.theta * alpha; 761 double sinTheta = sin(deltaTheta); 762 double cosTheta = sqrt(1.0 - sinTheta * sinTheta) * m_wb.direction; 763 X = length * sinTheta + m_wb.x; 764 Y = length * cosTheta + m_wb.y; 765 } 766 else 767 { 768 double alpha = y * m_wb.direction / m_wb.width; 769 double length = (x - alpha * m_wb.correction) * m_wb.direction + m_wb.radius; 770 double deltaTheta = m_wb.theta * alpha; 771 double sinTheta = sin(deltaTheta); 772 double cosTheta = sqrt(1.0 - sinTheta * sinTheta) * m_wb.direction; 773 Y = length * sinTheta + m_wb.y; 774 X = length * cosTheta + m_wb.x; 775 } 776 z = ProjZ(trs, X, Y, 1.0); 777 wx = ProjX(trs, X, Y, z, 1.0); 778 wy = ProjY(trs, X, Y, z, 1.0); 779} 780 781void Blend::FrameToMosaic(double trs[3][3], double x, double y, double &wx, double &wy) 782{ 783 // Project into the intermediate Mosaic coordinate system 784 double z = ProjZ(trs, x, y, 1.0); 785 double X = ProjX(trs, x, y, z, 1.0); 786 double Y = ProjY(trs, x, y, z, 1.0); 787 788 if (m_wb.theta == 0.0) 789 { 790 // No rotation, then this is all we need to do. 791 wx = X; 792 wy = Y; 793 } 794 else if (m_wb.horizontal) 795 { 796 double deltaX = X - m_wb.x; 797 double deltaY = Y - m_wb.y; 798 double length = sqrt(deltaX * deltaX + deltaY * deltaY); 799 double deltaTheta = asin(deltaX / length); 800 double alpha = deltaTheta / m_wb.theta; 801 wx = alpha * m_wb.width * m_wb.direction; 802 wy = (length - m_wb.radius) * m_wb.direction + alpha * m_wb.correction; 803 } 804 else 805 { 806 double deltaX = X - m_wb.x; 807 double deltaY = Y - m_wb.y; 808 double length = sqrt(deltaX * deltaX + deltaY * deltaY); 809 double deltaTheta = asin(deltaY / length); 810 double alpha = deltaTheta / m_wb.theta; 811 wy = alpha * m_wb.width * m_wb.direction; 812 wx = (length - m_wb.radius) * m_wb.direction + alpha * m_wb.correction; 813 } 814} 815 816 817 818// Clip the region of interest as small as possible by using the Voronoi edges of 819// the neighbors 820void Blend::ClipBlendRect(CSite *csite, BlendRect &brect) 821{ 822 SEdgeVector *ce; 823 int ecnt; 824 for (ce = csite->getNeighbor(), ecnt = csite->getNumNeighbors(); ecnt--; ce++) 825 { 826 // calculate the Voronoi bisector intersection 827 const double epsilon = 1e-5; 828 double dx = (m_AllSites[ce->second].getVCenter().x - m_AllSites[ce->first].getVCenter().x); 829 double dy = (m_AllSites[ce->second].getVCenter().y - m_AllSites[ce->first].getVCenter().y); 830 double xmid = m_AllSites[ce->first].getVCenter().x + dx/2.0; 831 double ymid = m_AllSites[ce->first].getVCenter().y + dy/2.0; 832 double inter; 833 834 if (dx > epsilon) 835 { 836 // neighbor is on right 837 if ((inter = m_wb.roundoffOverlap + xmid - dy * (((dy >= 0.0) ? brect.bot : brect.top) - ymid) / dx) < brect.rgt) 838 brect.rgt = inter; 839 } 840 else if (dx < -epsilon) 841 { 842 // neighbor is on left 843 if ((inter = -m_wb.roundoffOverlap + xmid - dy * (((dy >= 0.0) ? brect.bot : brect.top) - ymid) / dx) > brect.lft) 844 brect.lft = inter; 845 } 846 if (dy > epsilon) 847 { 848 // neighbor is above 849 if ((inter = m_wb.roundoffOverlap + ymid - dx * (((dx >= 0.0) ? brect.lft : brect.rgt) - xmid) / dy) < brect.top) 850 brect.top = inter; 851 } 852 else if (dy < -epsilon) 853 { 854 // neighbor is below 855 if ((inter = -m_wb.roundoffOverlap + ymid - dx * (((dx >= 0.0) ? brect.lft : brect.rgt) - xmid) / dy) > brect.bot) 856 brect.bot = inter; 857 } 858 } 859} 860 861void Blend::FrameToMosaicRect(int width, int height, double trs[3][3], BlendRect &brect) 862{ 863 // We need to walk the perimeter since the borders can be bent. 864 brect.lft = brect.bot = 2e30; 865 brect.rgt = brect.top = -2e30; 866 double xpos, ypos; 867 double lasty = height - 1.0; 868 double lastx = width - 1.0; 869 int i; 870 871 for (i = width; i--;) 872 { 873 874 FrameToMosaic(trs, (double) i, 0.0, xpos, ypos); 875 ClipRect(xpos, ypos, brect); 876 FrameToMosaic(trs, (double) i, lasty, xpos, ypos); 877 ClipRect(xpos, ypos, brect); 878 } 879 for (i = height; i--;) 880 { 881 FrameToMosaic(trs, 0.0, (double) i, xpos, ypos); 882 ClipRect(xpos, ypos, brect); 883 FrameToMosaic(trs, lastx, (double) i, xpos, ypos); 884 ClipRect(xpos, ypos, brect); 885 } 886} 887 888 889 890void Blend::ComputeBlendParameters(MosaicFrame **frames, int frames_size, int is360) 891{ 892 if (m_wb.blendingType != BLEND_TYPE_CYLPAN && m_wb.blendingType != BLEND_TYPE_HORZ) 893 { 894 m_wb.theta = 0.0; 895 return; 896 } 897 898 MosaicFrame *first = frames[0]; 899 MosaicFrame *last = frames[frames_size-1]; 900 MosaicFrame *mb; 901 902 double lxpos = last->trs[0][2], lypos = last->trs[1][2]; 903 double fxpos = first->trs[0][2], fypos = first->trs[1][2]; 904 905 // Calculate warp to produce proper stitching. 906 // get x, y displacement 907 double midX = last->width / 2.0; 908 double midY = last->height / 2.0; 909 double z = ProjZ(first->trs, midX, midY, 1.0); 910 double firstX, firstY; 911 double prevX = firstX = ProjX(first->trs, midX, midY, z, 1.0); 912 double prevY = firstY = ProjY(first->trs, midX, midY, z, 1.0); 913 914 double arcLength, lastTheta; 915 m_wb.theta = lastTheta = arcLength = 0.0; 916 for (int i = 0; i < frames_size; i++) 917 { 918 mb = frames[i]; 919 double currX, currY; 920 z = ProjZ(mb->trs, midX, midY, 1.0); 921 currX = ProjX(mb->trs, midX, midY, z, 1.0); 922 currY = ProjY(mb->trs, midX, midY, z, 1.0); 923 double deltaX = currX - prevX; 924 double deltaY = currY - prevY; 925 arcLength += sqrt(deltaY * deltaY + deltaX * deltaX); 926 if (!is360) 927 { 928 double thisTheta = asin(mb->trs[1][0]); 929 m_wb.theta += thisTheta - lastTheta; 930 lastTheta = thisTheta; 931 } 932 prevX = currX; 933 prevY = currY; 934 } 935 936 // In case of BMP output, stretch this to end at the proper alignment 937 m_wb.width = arcLength; 938 if (is360) m_wb.theta = asin(last->trs[1][0]); 939 940 // If there is no rotation, we're done. 941 if (m_wb.theta != 0.0) 942 { 943 double dx = prevX - firstX; 944 double dy = prevY - firstY; 945 if (abs(lxpos - fxpos) > abs(lypos - fypos)) 946 { 947 m_wb.horizontal = 1; 948 // Calculate radius position to make ends exactly the same Y offset 949 double radiusTheta = dx / cos(3.14159 / 2.0 - m_wb.theta); 950 m_wb.radius = dy + radiusTheta * cos(m_wb.theta); 951 if (m_wb.radius < 0.0) m_wb.radius = -m_wb.radius; 952 } 953 else 954 { 955 m_wb.horizontal = 0; 956 // Calculate radius position to make ends exactly the same Y offset 957 double radiusTheta = dy / cos(3.14159 / 2.0 - m_wb.theta); 958 m_wb.radius = dx + radiusTheta * cos(m_wb.theta); 959 if (m_wb.radius < 0.0) m_wb.radius = -m_wb.radius; 960 } 961 962 // Determine major direction 963 if (m_wb.horizontal) 964 { 965 // Horizontal strip 966 if (is360) m_wb.x = firstX; 967 else 968 { 969 if (lxpos - fxpos < 0) 970 { 971 m_wb.x = firstX + midX; 972 z = ProjZ(last->trs, 0.0, midY, 1.0); 973 prevX = ProjX(last->trs, 0.0, midY, z, 1.0); 974 prevY = ProjY(last->trs, 0.0, midY, z, 1.0); 975 } 976 else 977 { 978 m_wb.x = firstX - midX; 979 z = ProjZ(last->trs, last->width - 1.0, midY, 1.0); 980 prevX = ProjX(last->trs, last->width - 1.0, midY, z, 1.0); 981 prevY = ProjY(last->trs, last->width - 1.0, midY, z, 1.0); 982 } 983 } 984 dy = prevY - firstY; 985 if (dy < 0.0) m_wb.direction = 1.0; 986 else m_wb.direction = -1.0; 987 m_wb.y = firstY - m_wb.radius * m_wb.direction; 988 if (dy * m_wb.theta > 0.0) m_wb.width = -m_wb.width; 989 } 990 else 991 { 992 // Vertical strip 993 if (is360) m_wb.y = firstY; 994 else 995 { 996 if (lypos - fypos < 0) 997 { 998 m_wb.x = firstY + midY; 999 z = ProjZ(last->trs, midX, 0.0, 1.0); 1000 prevX = ProjX(last->trs, midX, 0.0, z, 1.0); 1001 prevY = ProjY(last->trs, midX, 0.0, z, 1.0); 1002 } 1003 else 1004 { 1005 m_wb.x = firstX - midX; 1006 z = ProjZ(last->trs, midX, last->height - 1.0, 1.0); 1007 prevX = ProjX(last->trs, midX, last->height - 1.0, z, 1.0); 1008 prevY = ProjY(last->trs, midX, last->height - 1.0, z, 1.0); 1009 } 1010 } 1011 dx = prevX - firstX; 1012 if (dx < 0.0) m_wb.direction = 1.0; 1013 else m_wb.direction = -1.0; 1014 m_wb.x = firstX - m_wb.radius * m_wb.direction; 1015 if (dx * m_wb.theta > 0.0) m_wb.width = -m_wb.width; 1016 } 1017 1018 // Calculate the correct correction factor 1019 double deltaX = prevX - m_wb.x; 1020 double deltaY = prevY - m_wb.y; 1021 double length = sqrt(deltaX * deltaX + deltaY * deltaY); 1022 double deltaTheta = (m_wb.horizontal) ? deltaX : deltaY; 1023 deltaTheta = asin(deltaTheta / length); 1024 m_wb.correction = ((m_wb.radius - length) * m_wb.direction) / 1025 (deltaTheta / m_wb.theta); 1026 } 1027} 1028