Nonrigid image deformation using moving regularized least quares

MLS算法是基于最小二乘法的移动网格变形算法,该算法依赖于网格大小的选择,而本文介绍的这篇论文,就是对网格变形的改进,使他不依赖于网格的选择。

为了方便更多的人看懂,本文主要讲算法实现的流程,不侧重于数学公式的推导之类。所以这里直接给出算法实现过程:

其中:

上面的算法流程经作者实践,比较耗时,于是作者做了优化,对公式8做了改进,如下:

论文效果图如下:

本人代码实现这个算法如下:

// 对矩阵求逆,结果保存在a中
int f_GaussInverseMatrix(double A[], double B[], int nColumns)
{  memcpy(B, A, sizeof(double) * nColumns * nColumns);int *is,*js,i,j,k,l,u,v;  double d,p;  is = new int[nColumns];  js = new int[nColumns];  for (k=0; k<=nColumns-1; k++)  {  d=0.0;  for (i=k; i<=nColumns-1; i++)  for (j=k; j<=nColumns-1; j++)  {  l=i*nColumns+j; p=fabs(B[l]);  if (p>d) { d=p; is[k]=i; js[k]=j;}  }  if (d+1.0==1.0)  {  free(is); free(js); printf("err**not inv\n");  return(1);  }  if (is[k]!=k)  for (j=0; j<=nColumns-1; j++)  {  u=k*nColumns+j; v=is[k]*nColumns+j;  p=B[u]; B[u]=B[v]; B[v]=p;  }  if (js[k]!=k)  for (i=0; i<=nColumns-1; i++)  {  u=i*nColumns+k; v=i*nColumns+js[k];  p=B[u]; B[u]=B[v]; B[v]=p;  }  l=k*nColumns+k;  B[l]=1.0f/B[l];  for (j=0; j<=nColumns-1; j++)  if (j!=k)  { u=k*nColumns+j; B[u]=B[u]*B[l];}  for (i=0; i<=nColumns-1; i++)  if (i!=k)  for (j=0; j<=nColumns-1; j++)  if (j!=k)  {  u=i*nColumns+j;  B[u] -= B[i*nColumns+k]*B[k*nColumns+j];  }  for (i=0; i<=nColumns-1; i++)  if (i!=k)  {  u=i*nColumns+k;  B[u] = -B[u]*B[l];  }  }  for (k=nColumns-1; k>=0; k--)  {  if (js[k]!=k)  for (j=0; j<=nColumns-1; j++)  {  u=k*nColumns+j; v=js[k]*nColumns+j;  p=B[u]; B[u]=B[v]; B[v]=p;  }  if (is[k]!=k)  for (i=0; i<=nColumns-1; i++)  {   u=i*nColumns+k; v=i*nColumns+is[k];  p=B[u]; B[u]=B[v]; B[v]=p;  }  }  free(is); free(js);  return 0;
}  // 求矩阵的乘积 C = A.*B
void f_matrixMul(double a[], double b[], int m, int n, int k, double c[])
{  int i, j, l, u;  for (i = 0; i < m; i++)  {for (j = 0; j < k; j++)  {  u = i * k + j; c[u] = 0.0;  for (l = 0; l < n; l++)  {c[u] += a[i * n + l] * b[l * k + j];}}  }
}
#define MIN2(a, b) ((a) < (b) ? (a) : (b))
#define MAX2(a, b) ((a) > (b) ? (a) : (b))
#define CLIP3(x, a, b) MIN2(MAX2(a,x), b)
//------------------------------------------------------------------
//Function: MRLS
//Params:
//return: 0 or false
//Reference: Nonrigid image deformation using moving regularized least squares. Jiayi Ma,2013
//-------------------------------------------------------------------
int f_MLSR(unsigned char *srcData, int width, int height, int stride, int srcPoint[], int dstPoint[], int pointNum)
{int ret = 0;unsigned char* tempData = (unsigned char*)malloc(sizeof(unsigned char)* height * stride);memcpy(tempData, srcData, sizeof(unsigned char) * height * stride);//Processdouble *srcPx = (double*)malloc(sizeof(double) * pointNum);double *srcPy = (double*)malloc(sizeof(double) * pointNum);double *dstPx = (double*)malloc(sizeof(double) * pointNum);double *dstPy = (double*)malloc(sizeof(double) * pointNum);double* GammaIJ= (double*)malloc(sizeof(double) * pointNum * pointNum); double* GammaIJT= (double*)malloc(sizeof(double) * pointNum * pointNum); double* GammaP = (double*)malloc(sizeof(double) * pointNum); double* S      = (double*)malloc(sizeof(double) * pointNum); double* W      = (double*)malloc(sizeof(double) * pointNum * pointNum); double* WInv   = (double*)malloc(sizeof(double) * pointNum * pointNum); double* tempIJ = (double*)malloc(sizeof(double) * pointNum * pointNum); double belta = 2.0;double alpha = 2.0;double lambda = 1.0;double fbelta = 1.0 / (belta * belta);for(int i = 0; i < pointNum; i++){srcPx[i] = (double)srcPoint[2 * i] / width;srcPy[i] = (double)srcPoint[2 * i + 1] / height;dstPx[i] = (double)dstPoint[2 * i] / width;dstPy[i] = (double)dstPoint[2 * i + 1] / height;}int w = pointNum;for(int j = 0; j < pointNum; j++){for(int i = 0; i < pointNum; i++){GammaIJ[i + j * w] = exp(-((double)(srcPx[i] - srcPx[j]) * (srcPx[i] - srcPx[j]) + (srcPy[i] - srcPy[j]) * (srcPy[i] - srcPy[j])) * fbelta);}}unsigned char* pSrc = srcData;int aa, bb, cc, dd, pos, pos1;double r1, r2, r3;unsigned char *pSrcL1;unsigned char *pSrcL2;unsigned char* p = tempData;for(int j = 0; j < height; j++){for(int i = 0; i < width; i++){double px = (double)i / width, py = (double)j / height;for(int n = 0; n < pointNum; n++){for(int m = 0; m < pointNum; m++){if(n == m)W[m + n * w] = pow(((double)(px - srcPx[m]) * (px - srcPx[m]) + (py - srcPy[m]) * (py - srcPy[m])), -alpha);elseW[m + n * w] = 0;}}//compute inverse matrix of Wf_GaussInverseMatrix(W, WInv, pointNum);//compute Gamma + lambda * WInvfor(int n = 0; n < pointNum; n++){for(int m = 0; m < pointNum; m++){//W is a diagonal matrix with the i-th entry Wi(p),other's 0GammaIJT[m + n * w] = GammaIJ[m + n * w] + lambda * WInv[m + n * w];}}//compute inverse matrix of (Gamma + lambda * WInv)f_GaussInverseMatrix(GammaIJT, tempIJ, pointNum);//compute GammaP   exp(-(p-xi) * (p - xi) / (belta * belta))for(int m = 0; m < pointNum; m++)GammaP[m] = exp((-((double)(px - srcPx[m]) * (px - srcPx[m]) + (py - srcPy[m]) * (py - srcPy[m])) * fbelta));//compute Sf_matrixMul(GammaP, tempIJ, 1, pointNum, pointNum, S);//compute fp   T + SYdouble sumx = 0, sumy = 0;for(int m = 0; m < pointNum; m++){sumx += S[m] * srcPx[m];sumy += S[m] * srcPy[m];}px = px - sumx;py = py - sumy;sumx = 0, sumy = 0;for(int m = 0; m < pointNum; m++){sumx += S[m] * dstPx[m];sumy += S[m] * dstPy[m];}px = px + sumx;py = py + sumy;double x_in = CLIP3(px * width, 0, width - 1);double y_in = CLIP3(py * height, 0, height - 1);//please change interpolation to get better effects.int xx = (int)x_in;int yy = (int)y_in;     pSrcL1 = p + yy * stride;pSrcL2 = pSrcL1 + stride;pos = (xx << 2);aa = pSrcL1[pos];bb = pSrcL1[pos + 4];cc = pSrcL2[pos];dd = pSrcL2[pos + 4];r1 = aa + (bb - aa) * (x_in - xx);r2 = cc + (dd - cc) * (x_in - xx);r3 = r1 + (r2 - r1) * (y_in - yy);pSrc[0]=(unsigned char)(CLIP3(r3, 0, 255));//B分量aa = pSrcL1[pos + 1];bb = pSrcL1[pos + 4 +1];cc = pSrcL2[pos + 1];dd = pSrcL2[pos + 4 + 1];r1 = aa + (bb - aa) * (x_in - xx);r2 = cc + (dd - cc) * (x_in - xx);r3 = r1 + (r2 - r1) * (y_in - yy);pSrc[1]=(unsigned char)(CLIP3(r3, 0, 255));//G分量aa = pSrcL1[pos + 2];bb = pSrcL1[pos + 4 + 2];cc = pSrcL2[pos + 2];dd = pSrcL2[pos + 4 + 2];r1 = aa + (bb - aa) * (x_in - xx);r2 = cc + (dd - cc) * (x_in - xx);r3 = r1 + (r2 - r1) * (y_in - yy);pSrc[2]=(unsigned char)(CLIP3(r3, 0, 255));//R分量aa = pSrcL1[pos + 3];bb = pSrcL1[pos + 4 + 3];cc = pSrcL2[pos + 3];dd = pSrcL2[pos + 4 + 3];r1=aa + (bb - aa) * (x_in - xx);r2=cc + (dd - cc) * (x_in - xx);r3=r1 + (r2 - r1) * (y_in - yy);pSrc[3]=(unsigned char)(CLIP3(r3, 0, 255));//A分量pSrc += 4;}}free(tempData);free(srcPx);free(srcPy);free(dstPx);free(dstPy);free(GammaIJ);free(GammaP);free(S     );free(W     );free(WInv  );free(tempIJ);free(GammaIJT);return ret;
};

具体效果,本人与上一篇博文保持一致仍然以瘦脸效果为例,原图和效果图如下:

原图(蓝色点为原始点,红色点为拖拽点)

瘦脸效果图

图片来自网络,如有侵权敬请告知

该算法优缺点:

算法变形效果非常好,但是在10个点位个数以上,速度会随点位增加而变慢,对于全脸100+个特征点变形而言,速度太慢;算法适用于10个点位左右的人脸五官变形,比如瘦脸,胖脸,小脸,嘴巴变形等等,当然也可以做其他变形功能需求;

本人QQ1358009172

Imagewarping变形算法研究---MLSR(Nonrigid image deformation using moving regularized least quares)相关推荐

  1. 反距离加权matlab算法,ImageWarping变形算法研究---反距离加权插值(IDW)

    参考论文:Image Warping with Scattered Data Interpolation Inverse distance weighted interpolation算法(IDW)实 ...

  2. ImageWarping变形算法研究---反距离加权插值(IDW)

    参考论文:Image Warping with Scattered Data Interpolation Inverse distance weighted interpolation算法(IDW)实 ...

  3. 基于深度学习的场景分割算法研究综述

    基于深度学习的场景分割算法研究综述 人工智能技术与咨询 来自<计算机研究与发展> ,作者张 蕊等 摘 要 场景分割的目标是判断场景图像中每个像素的类别.场景分割是计算机视觉领域重要的基本问 ...

  4. 基于语音的疲劳度检测算法研究

    基于语音的疲劳度检测算法研究 摘 要 疲劳是一种自然现象,是人体的一种自我调节和保护功能.检测疲劳状态对于当今社会从事各行各业都有积极意义.本课题提出了一种基于语音特征参数和概率神经网络的语音疲劳度识 ...

  5. (转)卡马克卷轴算法研究

    卡马克卷轴算法研究 中文摘要 对于J2ME框架下的手机游戏程序的开发,其地图滚动的重绘有多种算法,由于手机性能的限制和开发周期等其他非技术条件,需要根据情况灵活选择所需的技术.但在及其苛刻条件下,如系 ...

  6. matlab人脸识别开题报告,基于人脸识别的出勤点名系统中特征提取算法研究开题报告...

    基于人脸识别的出勤点名系统中特征提取算法研究 一.本课题研究的目的,意义 人脸识别是一项既有科学研究价值,又有广泛应用前景的研究课题.国际上大量研究人员几十年的研究取得了丰硕的研究成果,自动人脸识别技 ...

  7. 人像美妆---妆容迁移算法研究(Makeup transfer)

    对于人像美妆算法,现在的美妆相机.玩美彩妆之类的app已经做的比较成熟了,但是具体算法,基本网络上是杳无可查,今天本人介绍一种自动的人像美妆算法----(Makeup Transfer)妆容迁移 妆容 ...

  8. 照片美妆--人像变老算法研究

    人像变老技术可以把一张小孩子的照片或者年轻人的照片转换为变老以后的样子,目前市面上已有相应的应用,这里本人先放两张效果,然后分析算法: 这个效果是本人算法的效果,现在我们来分析一下人像变老的技术情况. ...

  9. mimo 雷达成像 matlab,MIMO雷达成像算法研究

    摘要: 作为一种新兴雷达技术,多输入多输出(MIMO)雷达受到国内外科研人员的广泛关注.该雷达综合利用阵列与分集技术,能够引入远多于实际物理阵元数目的观测通道和自由度,并在信号层联合进行多通道回波数据 ...

最新文章

  1. Road-SLAM:基于道路标线车道级精度SLAM
  2. java静态代理与动态代理
  3. c语言编程员工管理的代码,员工信息管理完整(含附源代码).doc
  4. OpenCV均值漂移meanshift algorithm算法的实例(附完整代码)
  5. SQL语句 SELECT LIKE用法详解
  6. Java回顾之多线程同步
  7. 使用 ABAP 读取每个月的月份名称和编号
  8. ​北京大学 2022 年博士研究生招生简章
  9. python-ImageDraw
  10. 转:能和LoadRunner匹敌的VS2010/2012Web负载测试
  11. 2023王道计算机考研数据结构第一章-绪论
  12. hp linux 禁用u盘启动项,惠普台式机UEFI BIOS设置U盘启动
  13. 2022年调味品行业研究报告
  14. codeforces 407C Curious Array 数学
  15. python一行输入多个值用空格隔开_Python 实现一行输入多个数字(用空格隔开)
  16. 第六届“强网杯”青少年专项赛 writewp by 楠辞姐姐后援团
  17. 数仓4.0(三)------数据仓库系统
  18. 【Axure10基础教程】第七章 设置文本
  19. 艾森哲面试 Accenture
  20. 计算机专业论文结束语,毕业设计论文的结束语

热门文章

  1. 基于思维导图的研究生创新能力培养
  2. 以全局产业观领航智慧城市建设
  3. VLC for Android源码下载和编译
  4. ASP.NET上传一个文件夹
  5. springdoc swagger3 文件上传API正确写法
  6. Python+Django+sqlite3实现基于内容的音乐推荐系统
  7. 基于 Ng-zorro-antd 的企业后台模板 ng-alain
  8. 计算机使用了休眠 怎么唤醒,电脑睡眠模式怎么唤醒?
  9. matlab半小提琴图,【ggplot2】不同方法画half -小提琴图
  10. 用户输入一个整数,求出它的各个位数,并求各位数之和