Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse files

initial commit

  • Loading branch information...
commit da6446e1c6650bbd3a2eddf96972238d0f573264 0 parents
@pkmital authored
Showing with 58,362 additions and 0 deletions.
  1. 0  README
  2. +64 −0 model/face.con
  3. +1 −0  model/face.tracker
  4. +94 −0 model/face.tri
  5. +1 −0  model/face2.tracker
  6. +23,550 −0 model/haarcascade_frontalface_alt2.xml
  7. BIN  src/.DS_Store
  8. +385 −0 src/aam-library/AAM_Basic.cpp
  9. +99 −0 src/aam-library/AAM_Basic.h
  10. +454 −0 src/aam-library/AAM_CAM.cpp
  11. +118 −0 src/aam-library/AAM_CAM.h
  12. +73 −0 src/aam-library/AAM_Config.h
  13. +613 −0 src/aam-library/AAM_IC.cpp
  14. +118 −0 src/aam-library/AAM_IC.h
  15. +57 −0 src/aam-library/AAM_MovieAVI.cpp
  16. +39 −0 src/aam-library/AAM_MovieAVI.h
  17. +644 −0 src/aam-library/AAM_PAW.cpp
  18. +121 −0 src/aam-library/AAM_PAW.h
  19. +291 −0 src/aam-library/AAM_PDM.cpp
  20. +85 −0 src/aam-library/AAM_PDM.h
  21. +611 −0 src/aam-library/AAM_Shape.cpp
  22. +97 −0 src/aam-library/AAM_Shape.h
  23. +316 −0 src/aam-library/AAM_TDM.cpp
  24. +89 −0 src/aam-library/AAM_TDM.h
  25. +530 −0 src/aam-library/AAM_Util.cpp
  26. +114 −0 src/aam-library/AAM_Util.h
  27. +69 −0 src/aam-library/VJfacedetect.cpp
  28. +40 −0 src/aam-library/VJfacedetect.h
  29. +11 −0 src/app.cpp
  30. +121 −0 src/app.h
  31. BIN  src/exe/.DS_Store
  32. +23 −0 src/exe/main.cpp
  33. +312 −0 src/facetracker/CLM.cc
  34. +86 −0 src/facetracker/CLM.h
  35. +40 −0 src/facetracker/COPYRIGHT
  36. +136 −0 src/facetracker/FCheck.cc
  37. +93 −0 src/facetracker/FCheck.h
  38. +208 −0 src/facetracker/FDet.cc
  39. +84 −0 src/facetracker/FDet.h
  40. +134 −0 src/facetracker/IO.cc
  41. +60 −0 src/facetracker/IO.h
  42. +214 −0 src/facetracker/PAW.cc
  43. +86 −0 src/facetracker/PAW.h
  44. +362 −0 src/facetracker/PDM.cc
  45. +82 −0 src/facetracker/PDM.h
  46. +244 −0 src/facetracker/Patch.cc
  47. +98 −0 src/facetracker/Patch.h
  48. +52 −0 src/facetracker/README
  49. +185 −0 src/facetracker/Tracker.cc
  50. +116 −0 src/facetracker/Tracker.h
  51. +11 −0 src/pkmFaceModeler.cpp
  52. +142 −0 src/pkmFaceModeler.h
  53. +11 −0 src/pkmFaceTracker.cpp
  54. +181 −0 src/pkmFaceTracker.h
  55. BIN  xcode/.DS_Store
  56. BIN  xcode/build/.DS_Store
  57. +64 −0 xcode/build/model/face.con
  58. +1 −0  xcode/build/model/face.tracker
  59. +94 −0 xcode/build/model/face.tri
  60. +1 −0  xcode/build/model/face2.tracker
  61. +23,550 −0 xcode/build/model/haarcascade_frontalface_alt2.xml
  62. +79 −0 xcode/xcode.1
  63. +764 −0 xcode/xcode.xcodeproj/pkmital.pbxuser
  64. +1,611 −0 xcode/xcode.xcodeproj/pkmital.perspectivev3
  65. +433 −0 xcode/xcode.xcodeproj/project.pbxproj
0  README
No changes.
64 model/face.con
@@ -0,0 +1,64 @@
+n_connections: 61
+{
+0 1
+1 2
+2 3
+3 4
+4 5
+5 6
+6 7
+7 8
+8 9
+9 10
+10 11
+11 12
+12 13
+13 14
+14 15
+15 16
+17 18
+18 19
+19 20
+20 21
+22 23
+23 24
+24 25
+25 26
+27 28
+28 29
+29 30
+31 32
+32 33
+33 34
+34 35
+36 37
+37 38
+38 39
+39 40
+40 41
+41 36
+42 43
+43 44
+44 45
+45 46
+46 47
+47 42
+48 49
+49 50
+50 51
+51 52
+52 53
+53 54
+54 55
+55 56
+56 57
+57 58
+58 59
+59 48
+60 65
+60 61
+61 62
+62 63
+63 64
+64 65
+}
1  model/face.tracker
1 addition, 0 deletions not shown
94 model/face.tri
@@ -0,0 +1,94 @@
+n_tri: 91
+{
+20 21 23
+21 22 23
+0 1 36
+15 16 45
+0 17 36
+16 26 45
+17 18 37
+25 26 44
+17 36 37
+26 44 45
+18 19 38
+24 25 43
+18 37 38
+25 43 44
+19 20 38
+23 24 43
+20 21 39
+22 23 42
+20 38 39
+23 42 43
+21 22 27
+21 27 39
+22 27 42
+27 28 42
+27 28 39
+28 42 47
+28 39 40
+1 36 41
+15 45 46
+1 2 41
+14 15 46
+28 29 40
+28 29 47
+2 40 41
+14 46 47
+2 29 40
+14 29 47
+2 3 29
+13 14 29
+29 30 31
+29 30 35
+3 29 31
+13 29 35
+30 32 33
+30 33 34
+30 31 32
+30 34 35
+3 4 31
+12 13 35
+4 5 48
+11 12 54
+5 6 48
+10 11 54
+6 48 59
+10 54 55
+6 7 59
+9 10 55
+7 58 59
+9 55 56
+8 57 58
+8 56 57
+7 8 58
+8 9 56
+4 31 48
+12 35 54
+31 48 49
+35 53 54
+31 49 50
+35 52 53
+31 32 50
+34 35 52
+32 33 50
+33 34 52
+33 50 51
+33 51 52
+48 49 60
+49 60 50
+50 60 61
+50 51 61
+51 52 61
+61 62 52
+52 53 62
+53 54 62
+54 55 63
+55 56 63
+56 63 64
+56 57 64
+64 65 57
+57 58 65
+58 59 65
+48 59 65
+}
1  model/face2.tracker
1 addition, 0 deletions not shown
23,550 model/haarcascade_frontalface_alt2.xml
23,550 additions, 0 deletions not shown
BIN  src/.DS_Store
Binary file not shown
385 src/aam-library/AAM_Basic.cpp
@@ -0,0 +1,385 @@
+/****************************************************************************
+* AAMLibrary
+* http://code.google.com/p/aam-library
+* Copyright (c) 2008-2009 by GreatYao, all rights reserved.
+****************************************************************************/
+
+#include "AAM_Basic.h"
+
+using namespace std;
+
+//============================================================================
+AAM_Basic::AAM_Basic()
+{
+ __G = 0;
+ __current_c_q = 0;
+ __update_c_q = 0;
+ __delta_c_q = 0;
+ __c = 0;
+ __p = 0;
+ __q = 0;
+ __lamda = 0;
+ __s = 0;
+ __t_m = 0;
+ __t_s = 0;
+ __delta_t = 0;
+}
+
+//============================================================================
+AAM_Basic::~AAM_Basic()
+{
+ cvReleaseMat(&__G); cvReleaseMat(&__current_c_q);
+ cvReleaseMat(&__update_c_q); cvReleaseMat(&__delta_c_q);
+ cvReleaseMat(&__p); cvReleaseMat(&__q);
+ cvReleaseMat(&__c); cvReleaseMat(&__lamda);
+ cvReleaseMat(&__s); cvReleaseMat(&__t_s);
+ cvReleaseMat(&__t_m); cvReleaseMat(&__delta_t);
+}
+
+void AAM_Basic::Build(const file_lists& pts_files, const file_lists& img_files,
+ double scale /* = 1.0 */, int nPlane /* = 3 */)
+{
+ Train(pts_files, img_files, scale, nPlane);
+}
+
+//============================================================================
+void AAM_Basic::Train(const file_lists& pts_files, const file_lists& img_files,
+ double scale /* = 1.0 */, int nPlane /* = 3 */,
+ double shape_percentage /* = 0.975 */,
+ double texture_percentage /* = 0.975 */,
+ double appearance_percentage /* = 0.975 */)
+{
+ if(pts_files.size() != img_files.size())
+ {
+ AAM_FormatMSG("#Shapes != #Images");
+ AAM_ERROR(errmsg);
+ }
+
+ printf("################################################\n");
+ printf("Build Fixed Jocobian Active Appearance Model ...\n");
+
+ __cam.Train(pts_files, img_files, scale, nPlane, shape_percentage,
+ texture_percentage, appearance_percentage);
+
+ printf("Build Jacobian Matrix...\n");
+ __G = cvCreateMat(__cam.nModes()+4, __cam.__texture.nPixels()*__cam.nPlanes(), CV_64FC1);
+ CalcJacobianMatrix(pts_files, img_files);
+
+ //allocate memory for on-line fitting
+ __current_c_q = cvCreateMat(1, __cam.nModes()+4, CV_64FC1);
+ __update_c_q = cvCreateMat(1, __cam.nModes()+4, CV_64FC1);
+ __delta_c_q = cvCreateMat(1, __cam.nModes()+4, CV_64FC1);
+ __c = cvCreateMat(1, __cam.nModes(), CV_64FC1);
+ __p = cvCreateMat(1, __cam.__shape.nModes(), CV_64FC1);
+ __q = cvCreateMat(1, 4, CV_64FC1);
+ __lamda = cvCreateMat(1, __cam.__texture.nModes(), CV_64FC1);
+ __s = cvCreateMat(1, __cam.__shape.nPoints()*2, CV_64FC1);
+ __t_s = cvCreateMat(1, __cam.__texture.nPixels()*__cam.nPlanes(), CV_64FC1);
+ __t_m = cvCreateMat(1, __cam.__texture.nPixels()*__cam.nPlanes(), CV_64FC1);
+ __delta_t = cvCreateMat(1, __cam.__texture.nPixels()*__cam.nPlanes(), CV_64FC1);
+
+ printf("################################################\n\n");
+}
+
+//============================================================================
+static double rand_in_between(double a, double b)
+{
+ int A = rand() % 50;
+ return a + (b-a)*A/49;
+}
+
+//============================================================================
+void AAM_Basic::CalcJacobianMatrix(const file_lists& pts_files,
+ const file_lists& img_files,
+ double disp_scale /* = 0.2 */,
+ double disp_angle /* = 20 */,
+ double disp_trans /* = 5.0 */,
+ double disp_std /* = 1.0 */,
+ int nExp /* = 30 */)
+{
+ CvMat* J = cvCreateMat(__cam.nModes()+4, __cam.__texture.nPixels()*__cam.nPlanes(), CV_64FC1);
+ CvMat* d = cvCreateMat(1, __cam.nModes()+4, CV_64FC1);
+ CvMat* o = cvCreateMat(1, __cam.nModes()+4, CV_64FC1);
+ CvMat* oo = cvCreateMat(1, __cam.nModes()+4, CV_64FC1);
+ CvMat* t = cvCreateMat(1, __cam.__texture.nPixels()*__cam.nPlanes(), CV_64FC1);
+ CvMat* t_m = cvCreateMat(1, __cam.__texture.nPixels()*__cam.nPlanes(), CV_64FC1);
+ CvMat* t_s = cvCreateMat(1, __cam.__texture.nPixels()*__cam.nPlanes(), CV_64FC1);
+ CvMat* t1 = cvCreateMat(1, __cam.__texture.nPixels()*__cam.nPlanes(), CV_64FC1);
+ CvMat* t2 = cvCreateMat(1, __cam.__texture.nPixels()*__cam.nPlanes(), CV_64FC1);
+ CvMat* u = cvCreateMat(1, __cam.nModes()+4, CV_64FC1);
+ CvMat* c = cvCreateMat(1, __cam.nModes(), CV_64FC1);
+ CvMat* s = cvCreateMat(1, __cam.__shape.nPoints()*2, CV_64FC1);
+ CvMat* q = cvCreateMat(1, 4, CV_64FC1);
+ CvMat* p = cvCreateMat(1, __cam.__shape.nModes(),CV_64FC1);
+ CvMat* lamda = cvCreateMat(1, __cam.__texture.nModes(), CV_64FC1);
+
+ double theta = disp_angle * CV_PI / 180;
+ double aa = MAX(fabs(disp_scale*cos(theta)), fabs(disp_scale*sin(theta)));
+ cvmSet(d,0,0,aa); cvmSet(d,0,1,aa); cvmSet(d,0,2,disp_trans); cvmSet(d,0,3,disp_trans);
+ for(int nmode = 0; nmode < __cam.nModes(); nmode++)
+ cvmSet(d,0,4+nmode,disp_std*sqrt(__cam.Var(nmode)));
+
+ srand(unsigned(time(0)));
+ cvSetZero(u);cvSetZero(J);
+ for(int i = 0; i < (int)pts_files.size(); i++)
+ {
+ IplImage* image = cvLoadImage(img_files[i].c_str(), -1);
+ AAM_Shape Shape;
+ if(!Shape.ReadAnnotations(pts_files[i]))
+ Shape.ScaleXY(image->width, image->height);
+ Shape.Point2Mat(s);
+ AAM_Common::CheckShape(s, image->width, image->height);
+
+ //calculate current texture vector
+ __cam.__paw.CalcWarpTexture(s, image, t, __cam.__texture.nPlanes());
+ __cam.__texture.NormalizeTexture(__cam.__MeanG, t);
+
+ //calculate appearance parameters
+ __cam.__shape.CalcParams(s, p, q);
+ __cam.__texture.CalcParams(t, lamda);
+ __cam.CalcParams(c, p, lamda);
+
+ //update appearance and pose parameters
+ CvMat subo;
+ cvGetCols(o, &subo, 0, 4); cvCopy(q, &subo);
+ cvGetCols(o, &subo, 4, 4+__cam.nModes()); cvCopy(c, &subo);
+
+ //get optimal EstResidual
+ EstResidual(image, o, s, t_m, t_s, t1);
+
+ for(int j = 0; j < nExp; j++)
+ {
+ printf("Pertubing (%d/%d) for image (%d/%d)...\r", j, nExp, i, pts_files.size());
+ fflush(stdout);
+
+ for(int l = 0; l < 4+__cam.nModes(); l++)
+ {
+ double D = cvmGet(d,0,l);
+ double v = rand_in_between(-D, D);
+ cvCopy(o, oo); CV_MAT_ELEM(*oo,double,0,l) += v;
+ EstResidual(image, oo, s, t_m, t_s, t2);
+
+ cvSub(t1, t2, t2);
+ cvConvertScale(t2, t2, 1.0/v);
+
+ //accumulate into l-th row
+ CvMat Jl; cvGetRow(J, &Jl, l);
+ cvAdd(&Jl, t2, &Jl);
+
+ CV_MAT_ELEM(*u, double, 0, l) += 1.0;
+ }
+ }
+ cvReleaseImage(&image);
+ }
+
+ //normalize
+ for(int l = 0; l < __cam.nModes()+4; l++)
+ {
+ CvMat Jl; cvGetRow(J, &Jl, l);
+ cvConvertScale(&Jl, &Jl, 1.0/cvmGet(u,0,l));
+ }
+
+ CvMat* JtJ = cvCreateMat(__cam.nModes()+4, __cam.nModes()+4, CV_64FC1);
+ CvMat* InvJtJ = cvCreateMat(__cam.nModes()+4, __cam.nModes()+4, CV_64FC1);
+
+ cvGEMM(J, J, 1, NULL, 0, JtJ, CV_GEMM_B_T);
+ cvInvert(JtJ, InvJtJ, CV_SVD);
+ cvMatMul(InvJtJ, J, __G);
+
+ cvReleaseMat(&J); cvReleaseMat(&d); cvReleaseMat(&o);
+ cvReleaseMat(&oo); cvReleaseMat(&t); cvReleaseMat(&t_s);
+ cvReleaseMat(&t_m); cvReleaseMat(&t1); cvReleaseMat(&t2);
+ cvReleaseMat(&u); cvReleaseMat(&c); cvReleaseMat(&s);
+ cvReleaseMat(&q); cvReleaseMat(&p); cvReleaseMat(&lamda);
+ cvReleaseMat(&JtJ); cvReleaseMat(&InvJtJ);
+}
+
+
+//============================================================================
+double AAM_Basic::EstResidual(const IplImage* image, const CvMat* c_q,
+ CvMat* s, CvMat* t_m,
+ CvMat* t_s, CvMat* deltat)
+{
+ CvMat c, q;
+ cvGetCols(c_q, &q, 0, 4);
+ cvGetCols(c_q, &c, 4, 4+__cam.nModes());
+
+ // generate model texture
+ __cam.CalcTexture(t_m, &c);
+
+ // generate model shape
+ __cam.CalcShape(s, c_q);
+
+ // generate warped texture
+ AAM_Common::CheckShape(s, image->width, image->height);
+ __cam.__paw.CalcWarpTexture(s, image, t_s, __cam.__texture.nPlanes());
+ __cam.__texture.NormalizeTexture(__cam.__MeanG, t_s);
+
+ // calculate pixel difference
+ cvSub(t_m, t_s, deltat);
+
+ return cvNorm(deltat);
+}
+
+//============================================================================
+void AAM_Basic::SetAllParamsZero()
+{
+ cvSetZero(__q);
+ cvSetZero(__c);
+}
+
+//============================================================================
+void AAM_Basic::InitParams(const IplImage* image)
+{
+ //shape parameter
+ __cam.__shape.CalcParams(__s, __p, __q);
+
+ //texture parameter
+ __cam.__paw.CalcWarpTexture(__s, image, __t_s, __cam.__texture.nPlanes());
+ __cam.__texture.NormalizeTexture(__cam.__MeanG, __t_s);
+ __cam.__texture.CalcParams(__t_s, __lamda);
+
+ //combined appearance parameter
+ __cam.CalcParams(__c, __p, __lamda);
+}
+
+//============================================================================
+void AAM_Basic::Fit(const IplImage* image, AAM_Shape& Shape,
+ int max_iter /* = 30 */,bool showprocess /* = false */)
+{
+ //intial some stuff
+ double t = AAM_GetTime();
+ double e1, e2;
+ int k;
+ const int np = 5;
+ double k_values[np] = {1, 0.5, 0.25, 0.125, 0.0625};
+ IplImage* Drawimg = 0;
+
+ Shape.Point2Mat(__s);
+ InitParams(image);
+ CvMat subcq;
+ cvGetCols(__current_c_q, &subcq, 0, 4); cvCopy(__q, &subcq);
+ cvGetCols(__current_c_q, &subcq, 4, 4+__cam.nModes()); cvCopy(__c, &subcq);
+
+ //calculate error
+ e1 = EstResidual(image, __current_c_q, __s, __t_m, __t_s, __delta_t);
+
+ //do a number of iteration until convergence
+ for(int iter = 0; iter <max_iter; iter++)
+ {
+ if(showprocess)
+ {
+ if(Drawimg == 0) Drawimg = cvCloneImage(image);
+ else cvCopy(image, Drawimg);
+ __cam.CalcShape(__s, __current_c_q);
+ Shape.Mat2Point(__s);
+ Draw(Drawimg, Shape, 2);
+ AAM_MKDIR("result");
+ char filename[100];
+ sprintf(filename, "result/Iter-%02d.jpg", iter);
+ cvSaveImage(filename, Drawimg);
+ }
+
+ // predict parameter update
+ cvGEMM(__delta_t, __G, 1, NULL, 0, __delta_c_q, CV_GEMM_B_T);
+
+ //force first iteration
+ if(iter == 0)
+ {
+ cvAdd(__current_c_q, __delta_c_q, __current_c_q);
+ CvMat c; cvGetCols(__current_c_q, &c, 4, 4+__cam.nModes());
+ //constrain parameters
+ __cam.Clamp(&c);
+ e1 = EstResidual(image, __current_c_q, __s, __t_m, __t_s, __delta_t);
+ }
+
+ //find largest step size which reduces texture EstResidual
+ else
+ {
+ for(k = 0; k < np; k++)
+ {
+ cvScaleAdd(__delta_c_q, cvScalar(k_values[k]), __current_c_q, __update_c_q);
+ //constrain parameters
+ CvMat c; cvGetCols(__update_c_q, &c, 4, 4+__cam.nModes());
+ __cam.Clamp(&c);
+
+ e2 = EstResidual(image, __update_c_q, __s, __t_m, __t_s, __delta_t);
+ if(e2 <= e1) break;
+ }
+
+ //check for convergence
+ if(iter > 0)
+ {
+ if(k == np)
+ {
+ e1 = e2;
+ cvCopy(__update_c_q, __current_c_q);
+ }
+
+ else if(fabs(e2-e1)<0.001*e1) break;
+ else if (cvNorm(__delta_c_q)<0.001) break;
+ else
+ {
+ cvCopy(__update_c_q, __current_c_q);
+ e1 = e2;
+ }
+ }
+ }
+ }
+
+ cvReleaseImage(&Drawimg);
+ __cam.CalcShape(__s, __current_c_q);
+ Shape.Mat2Point(__s);
+ t = AAM_GetTime() - t;
+ printf("AAM-Basic Fitting time cost: %.3f millisec\n", t);
+}
+
+
+//===========================================================================
+void AAM_Basic::Draw(IplImage* image, const AAM_Shape& Shape, int type)
+{
+ if(type == 0) AAM_Common::DrawPoints(image, Shape);
+ else if(type == 1) AAM_Common::DrawTriangles(image, Shape, __cam.__paw.__tri);
+ else if(type == 2)
+ {
+ double minV, maxV;
+ cvMinMaxLoc(__t_m, &minV, &maxV);
+ cvConvertScale(__t_m, __t_m, 255/(maxV-minV), -minV*255/(maxV-minV));
+ AAM_PAW paw;
+ paw.Train(Shape, __cam.__Points, __cam.__Storage, __cam.__paw.GetTri(), false);
+ AAM_Common::DrawAppearance(image, Shape, __t_m, paw, __cam.__paw, __cam.nPlanes());
+ }
+ else
+ {
+ AAM_FormatMSG("Unsupported drawing type\n");
+ AAM_ERROR(errmsg);
+ }
+}
+
+//===========================================================================
+void AAM_Basic::Write(std::ofstream& os)
+{
+ __cam.Write(os);
+ os << __G;
+}
+
+//===========================================================================
+void AAM_Basic::Read(std::ifstream& is)
+{
+ __cam.Read(is);
+ __G = cvCreateMat(__cam.nModes()+4, __cam.__texture.nPixels()*__cam.nPlanes(), CV_64FC1);
+ is >> __G;
+
+ //allocate memory for on-line fitting
+ __current_c_q = cvCreateMat(1, __cam.nModes()+4, CV_64FC1);
+ __update_c_q = cvCreateMat(1, __cam.nModes()+4, CV_64FC1);
+ __delta_c_q = cvCreateMat(1, __cam.nModes()+4, CV_64FC1);
+ __c = cvCreateMat(1, __cam.nModes(), CV_64FC1);
+ __p = cvCreateMat(1, __cam.__shape.nModes(), CV_64FC1);
+ __q = cvCreateMat(1, 4, CV_64FC1);
+ __lamda = cvCreateMat(1, __cam.__texture.nModes(), CV_64FC1);
+ __s = cvCreateMat(1, __cam.__shape.nPoints()*2, CV_64FC1);
+ __t_s = cvCreateMat(1, __cam.__texture.nPixels()*__cam.nPlanes(), CV_64FC1);
+ __t_m = cvCreateMat(1, __cam.__texture.nPixels()*__cam.nPlanes(), CV_64FC1);
+ __delta_t = cvCreateMat(1, __cam.__texture.nPixels()*__cam.nPlanes(), CV_64FC1);
+}
99 src/aam-library/AAM_Basic.h
@@ -0,0 +1,99 @@
+/****************************************************************************
+* AAMLibrary
+* http://code.google.com/p/aam-library
+* Copyright (c) 2008-2009 by GreatYao, all rights reserved.
+****************************************************************************/
+
+#ifndef AAM_BASIC_H
+#define AAM_BASIC_H
+
+#include "AAM_Config.h"
+#include "AAM_Util.h"
+#include "AAM_Shape.h"
+#include "AAM_PDM.h"
+#include "AAM_TDM.h"
+#include "AAM_PAW.h"
+#include "AAM_CAM.h"
+
+class AAM_Pyramid;
+
+/**
+ The basic Active appearace model building and fitting method.
+ Refer to: T.F. Cootes, G.J. Edwards and C.J. Taylor. "Active Appearance Models". ECCV 1998
+*/
+class AAM_EXPORTS AAM_Basic:public AAM
+{
+public:
+ AAM_Basic();
+ ~AAM_Basic();
+
+ virtual const int GetType()const { return TYPE_AAM_BASIC; }
+
+ // Build cootes's basis aam
+ void Train(const file_lists& pts_files, const file_lists& img_files,
+ double scale = 1.0, int nPlane = 3,
+ double shape_percentage = 0.975,
+ double texture_percentage = 0.975,
+ double appearance_percentage = 0.975);
+
+ virtual void Build(const file_lists& pts_files, const file_lists& img_files,
+ double scale = 1.0, int nPlane = 3);
+
+ // Fit the image using aam basic.
+ virtual void Fit(const IplImage* image, AAM_Shape& Shape,
+ int max_iter = 30, bool showprocess = false);
+
+ // Virtual void SetAllParamsZero();
+
+ virtual const AAM_Shape GetMeanShape()const{ return __cam.__shape.GetMeanShape(); }
+ const AAM_Shape GetReferenceShape()const{ return __cam.__paw.__referenceshape; }
+
+ // Draw the image according the searching result(0:point, 1:triangle, 2:appearance)
+ virtual void Draw(IplImage* image, const AAM_Shape& Shape, int type);
+
+ // Read data from stream
+ virtual void Read(std::ifstream& is);
+
+ // Write data to stream
+ virtual void Write(std::ofstream& os);
+
+ // Set search parameters zero
+ virtual void SetAllParamsZero();
+
+ // Init search parameters
+ virtual void InitParams(const IplImage* image);
+
+private:
+ // Calculates the pixel difference from a model instance and an image
+ double EstResidual(const IplImage* image, const CvMat* c_q,
+ CvMat* s, CvMat* t_m, CvMat* t_s, CvMat* deltat);
+
+ // Estimate jacobian matrix during fitting
+ void CalcJacobianMatrix(const file_lists& pts_files,
+ const file_lists& img_files,
+ double disp_scale = 0.2, double disp_angle = 20,
+ double disp_trans = 5.0, double disp_std = 1.0,
+ int nExp = 30);
+
+private:
+ AAM_CAM __cam;
+ CvMat* __G;
+
+private:
+ //speed up for on-line alignment
+ CvMat* __current_c_q; //current appearance+pose parameters
+ CvMat* __update_c_q; //update appearance+pose parameters
+ CvMat* __delta_c_q; //displacement appearance+pose parameters
+ CvMat* __c; //appearance parameters
+ CvMat* __p; //shape parameters
+ CvMat* __q; //pose parameters
+ CvMat* __lamda; //texture parameters
+ CvMat* __s; //current shape
+ CvMat* __t_m; //model texture instance
+ CvMat* __t_s; //warped texture at current shape instance
+ CvMat* __delta_t; //differnce between __ts and __tm
+
+ friend class AAM_Pyramid;
+};
+
+#endif //
454 src/aam-library/AAM_CAM.cpp
@@ -0,0 +1,454 @@
+/****************************************************************************
+* AAMLibrary
+* http://code.google.com/p/aam-library
+* Copyright (c) 2008-2009 by GreatYao, all rights reserved.
+****************************************************************************/
+
+#include "AAM_CAM.h"
+
+using namespace std;
+
+//============================================================================
+AAM_CAM::AAM_CAM()
+{
+ __MeanAppearance = 0;
+ __AppearanceEigenValues = 0;
+ __AppearanceEigenVectors = 0;
+ __Qs = 0;
+ __Qg = 0;
+ __MeanS = 0;
+ __MeanG = 0;
+
+ __Points = 0;
+ __Storage = 0;
+ __pq = 0;
+
+ __a = 0;
+}
+
+//============================================================================
+AAM_CAM::~AAM_CAM()
+{
+ cvReleaseMat(&__MeanAppearance);
+ cvReleaseMat(&__AppearanceEigenValues);
+ cvReleaseMat(&__AppearanceEigenVectors);
+ cvReleaseMat(&__Qg);
+ cvReleaseMat(&__Qs);
+ cvReleaseMat(&__MeanS);
+ cvReleaseMat(&__MeanG);
+
+ cvReleaseMat(&__Points);
+ cvReleaseMemStorage(&__Storage);
+ cvReleaseMat(&__pq);
+
+ cvReleaseMat(&__a);
+}
+
+//============================================================================
+void AAM_CAM::Train(const file_lists& pts_files,
+ const file_lists& img_files,
+ double scale /* = 1.0 */,
+ int nPlane /* = 3 */,
+ double shape_percentage /* = 0.975 */,
+ double texture_percentage /* = 0.975 */,
+ double appearance_percentage /* = 0.975 */)
+{
+ //building shape and texture distribution model
+ std::vector<AAM_Shape> AllShapes;
+ for(int ii = 0; ii < (int)pts_files.size(); ii++)
+ {
+ AAM_Shape Shape;
+ bool flag = Shape.ReadAnnotations(pts_files[ii]);
+ if(!flag)
+ {
+ IplImage* image = cvLoadImage(img_files[ii].c_str(), -1);
+ Shape.ScaleXY(image->width, image->height);
+ cvReleaseImage(&image);
+ }
+ //Shape.Write(std::cout);
+ AllShapes.push_back(Shape);
+ }
+
+ printf("Build point distribution model...\n");
+ __shape.Train(AllShapes, scale, shape_percentage);
+
+ printf("Build warp information of mean shape mesh...");
+ __Points = cvCreateMat (1, __shape.nPoints(), CV_32FC2);
+ __Storage = cvCreateMemStorage(0);
+ AAM_Shape refShape = __shape.__AAMRefShape/* * scale */;
+ //if(refShape.GetWidth() > 120)
+ // refShape.Scale(120/refShape.GetWidth());
+
+ __paw.Train(refShape, __Points, __Storage);
+ printf("[%d by %d, %d triangles, %d*%d pixels]\n",
+ __paw.Width(), __paw.Height(), __paw.nTri(), __paw.nPix(), nPlane);
+
+ printf("Build texture distribution model...\n");
+ __texture.Train(pts_files, img_files, __paw, nPlane, texture_percentage, true);
+ __pq = cvCreateMat(1, __shape.nModes()+4, CV_64FC1);
+
+ printf("Build combined appearance model...\n");
+ int nsamples = pts_files.size();
+ int npointsby2 = __shape.nPoints()*2;
+ int npixels = __texture.nPixels();
+ int nfeatures = __shape.nModes() + __texture.nModes();
+ CvMat* AllAppearances = cvCreateMat(nsamples, nfeatures, CV_64FC1);
+ CvMat* s = cvCreateMat(1, npointsby2, CV_64FC1);
+ CvMat* t = cvCreateMat(1, npixels*nPlane, CV_64FC1);
+ __MeanS = cvCreateMat(1, npointsby2, CV_64FC1);
+ __MeanG = cvCreateMat(1, npixels*nPlane, CV_64FC1);
+ cvCopy(__shape.GetMean(), __MeanS);
+ cvCopy(__texture.GetMean(), __MeanG);
+
+ //calculate ratio of shape to appearance
+ CvScalar Sum1 = cvSum(__shape.__ShapesEigenValues);
+ CvScalar Sum2 = cvSum(__texture.__TextureEigenValues);
+ __WeightsS2T = sqrt(Sum2.val[0] / Sum1.val[0]);
+
+ printf("Combine shape and texture parameters...\n");
+ for(int i = 0; i < nsamples; i++)
+ {
+ //Get Shape and Texture respectively
+ IplImage* image = cvLoadImage(img_files[i].c_str(), -1);
+
+ AAM_Shape Shape;
+ if(!Shape.ReadAnnotations(pts_files[i]))
+ Shape.ScaleXY(image->width, image->height);
+ Shape.Point2Mat(s);
+ AAM_Common::CheckShape(s, image->width, image->height);
+
+ __paw.CalcWarpTexture(s, image, t, __texture.nPlanes());
+ __texture.NormalizeTexture(__MeanG, t);
+
+ //combine shape and texture parameters
+ CvMat OneAppearance;
+ cvGetRow(AllAppearances, &OneAppearance, i);
+ ShapeTexture2Combined(s, t, &OneAppearance);
+
+ cvReleaseImage(&image);
+ }
+
+ //Do PCA of appearances
+ DoPCA(AllAppearances, appearance_percentage);
+
+ int np = __AppearanceEigenVectors->rows;
+
+ printf("Extracting the shape and texture part of the combined eigen vectors..\n");
+
+ // extract the shape part of the combined eigen vectors
+ CvMat Ps;
+ cvGetCols(__AppearanceEigenVectors, &Ps, 0, __shape.nModes());
+ __Qs = cvCreateMat(np, npointsby2, CV_64FC1);
+ cvMatMul(&Ps, __shape.GetBases(), __Qs);
+ cvConvertScale(__Qs, __Qs, 1.0/__WeightsS2T);
+
+ // extract the texture part of the combined eigen vectors
+ CvMat Pg;
+ cvGetCols(__AppearanceEigenVectors, &Pg, __shape.nModes(), nfeatures);
+ __Qg = cvCreateMat(np, npixels*nPlane, CV_64FC1);
+ cvMatMul(&Pg, __texture.GetBases(), __Qg);
+
+ __a = cvCreateMat(1, __AppearanceEigenVectors->cols, CV_64FC1);
+}
+
+//============================================================================
+void AAM_CAM::ShapeTexture2Combined(const CvMat* Shape, const CvMat* Texture,
+ CvMat* Appearance)
+{
+ __shape.CalcParams(Shape, __pq);
+ CvMat mat1, mat2;
+ cvGetCols(__pq, &mat1, 4, 4+__shape.nModes());
+ cvGetCols(Appearance, &mat2, 0, __shape.nModes());
+ cvCopy(&mat1, &mat2);
+ cvConvertScale(&mat2, &mat2, __WeightsS2T);
+
+ cvGetCols(Appearance, &mat2, __shape.nModes(), __shape.nModes()+__texture.nModes());
+ __texture.CalcParams(Texture, &mat2);
+}
+
+//============================================================================
+void AAM_CAM::DoPCA(const CvMat* AllAppearances, double percentage)
+{
+ printf("Doing PCA of appearance datas...");
+
+ int nSamples = AllAppearances->rows;
+ int nfeatures = AllAppearances->cols;
+ int nEigenAtMost = MIN(nSamples, nfeatures);
+
+ CvMat* tmpEigenValues = cvCreateMat(1, nEigenAtMost, CV_64FC1);
+ CvMat* tmpEigenVectors = cvCreateMat(nEigenAtMost, nfeatures, CV_64FC1);
+ __MeanAppearance = cvCreateMat(1, nfeatures, CV_64FC1 );
+
+ cvCalcPCA(AllAppearances, __MeanAppearance,
+ tmpEigenValues, tmpEigenVectors, CV_PCA_DATA_AS_ROW);
+
+ double allSum = cvSum(tmpEigenValues).val[0];
+ double partSum = 0.0;
+ int nTruncated = 0;
+ double largesteigval = cvmGet(tmpEigenValues, 0, 0);
+ for(int i = 0; i < nEigenAtMost; i++)
+ {
+ double thiseigval = cvmGet(tmpEigenValues, 0, i);
+ if(thiseigval / largesteigval < 0.0001) break; // firstly check
+ partSum += thiseigval;
+ ++ nTruncated;
+ if(partSum/allSum >= percentage) break; //secondly check
+ }
+
+ __AppearanceEigenValues = cvCreateMat(1, nTruncated, CV_64FC1);
+ __AppearanceEigenVectors = cvCreateMat(nTruncated, nfeatures, CV_64FC1);
+
+ CvMat G;
+ cvGetCols(tmpEigenValues, &G, 0, nTruncated);
+ cvCopy(&G, __AppearanceEigenValues);
+
+ cvGetRows(tmpEigenVectors, &G, 0, nTruncated);
+ cvCopy(&G, __AppearanceEigenVectors);
+
+ cvReleaseMat(&tmpEigenVectors);
+ cvReleaseMat(&tmpEigenValues);
+ printf("Done (%d/%d)\n", nTruncated, nEigenAtMost);
+}
+
+//============================================================================
+void AAM_CAM::CalcLocalShape(CvMat* s, const CvMat* c)
+{
+ cvMatMul(c, __Qs, s);
+ cvAdd(s, __MeanS, s);
+}
+
+//============================================================================
+void AAM_CAM::CalcGlobalShape(CvMat* s, const CvMat* pose)
+{
+ int npoints = s->cols/2;
+ double* fasts = s->data.db;
+ double a=cvmGet(pose,0,0)+1, b=cvmGet(pose,0,1),
+ tx=cvmGet(pose,0,2), ty=cvmGet(pose,0,3);
+ double x, y;
+ for(int i = 0; i < npoints; i++)
+ {
+ x = fasts[2*i ];
+ y = fasts[2*i+1];
+
+ fasts[2*i ] = a*x-b*y+tx;
+ fasts[2*i+1] = b*x+a*y+ty;
+ }
+}
+
+//============================================================================
+void AAM_CAM::CalcTexture(CvMat* t, const CvMat* c)
+{
+ cvMatMul(c, __Qg, t);
+ cvAdd(t, __MeanG, t);
+}
+
+//============================================================================
+void AAM_CAM::CalcParams(CvMat* c, const CvMat* bs, const CvMat* bg)
+{
+ double* fasta = __a->data.db;
+ double* fastbs = bs->data.db;
+ double* fastbg = bg->data.db;
+
+ int i;
+ for(i = 0; i < bs->cols; i++) fasta[i] = __WeightsS2T * fastbs[i];
+ for(i = 0; i < bg->cols; i++) fasta[i+bs->cols] = fastbg[i];
+
+ cvProjectPCA(__a, __MeanAppearance, __AppearanceEigenVectors, c);
+}
+
+//============================================================================
+void AAM_CAM::Clamp(CvMat* c, double s_d /* = 3.0 */)
+{
+ double* fastc = c->data.db;
+ double* fastv = __AppearanceEigenValues->data.db;
+ int nmodes = nModes();
+ double limit;
+
+ for(int i = 0; i < nmodes; i++)
+ {
+ limit = s_d*sqrt(fastv[i]);
+ if(fastc[i] > limit) fastc[i] = limit;
+ else if(fastc[i] < -limit) fastc[i] = -limit;
+ }
+}
+
+//============================================================================
+void AAM_CAM::DrawAppearance(IplImage* image, const AAM_Shape& Shape, CvMat* Texture)
+{
+ AAM_PAW paw;
+ int nPlane = __texture.nPlanes();
+ int x1, x2, y1, y2, idx1 = 0, idx2 = 0;
+ int tri_idx, v1, v2, v3;
+ int minx, miny, maxx, maxy;
+ paw.Train(Shape, __Points, __Storage, __paw.GetTri(), false);
+ AAM_Shape refShape = __paw.__referenceshape;
+ double minV, maxV;
+ cvMinMaxLoc(Texture, &minV, &maxV);
+ cvConvertScale(Texture, Texture, 1/(maxV-minV)*255, -minV*255/(maxV-minV));
+
+ minx = (int)Shape.MinX(); miny = (int)Shape.MinY();
+ maxx = (int)Shape.MaxX(); maxy = (int)Shape.MaxY();
+ for(int y = miny; y < maxy; y++)
+ {
+ y1 = y-miny;
+ for(int x = minx; x < maxx; x++)
+ {
+ x1 = x-minx;
+ idx1 = paw.Rect(y1, x1);
+ if(idx1 >= 0)
+ {
+ tri_idx = paw.PixTri(idx1);
+ v1 = paw.Tri(tri_idx, 0);
+ v2 = paw.Tri(tri_idx, 1);
+ v3 = paw.Tri(tri_idx, 2);
+
+ x2 = (int)(paw.Alpha(idx1)*refShape[v1].x + paw.Belta(idx1)*refShape[v2].x +
+ paw.Gamma(idx1)*refShape[v3].x);
+ y2 = (int)(paw.Alpha(idx1)*refShape[v1].y + paw.Belta(idx1)*refShape[v2].y +
+ paw.Gamma(idx1)*refShape[v3].y);
+
+ idx2 = __paw.Rect(y2, x2);
+ if(idx2 < 0) continue;
+
+ if(nPlane == 3)
+ {
+ CV_IMAGE_ELEM(image, byte, y, 3*x) = (byte)cvmGet(Texture, 0, 3*idx2);
+ CV_IMAGE_ELEM(image, byte, y, 3*x+1) = (byte)cvmGet(Texture, 0, 3*idx2+1);
+ CV_IMAGE_ELEM(image, byte, y, 3*x+2) = (byte)cvmGet(Texture, 0, 3*idx2+2);
+ }
+ else
+ {
+ CV_IMAGE_ELEM(image, byte, y, 3*x) =
+ CV_IMAGE_ELEM(image, byte, y, 3*x+1) =
+ CV_IMAGE_ELEM(image, byte, y, 3*x+2) = (byte)cvmGet(Texture, 0, idx2);
+ }
+ }
+ }
+ }
+}
+
+//============================================================================
+void AAM_CAM::Write(std::ofstream& os)
+{
+ __shape.Write(os);
+ __texture.Write(os);
+ __paw.Write(os);
+
+ os << __AppearanceEigenVectors->rows << " " << __AppearanceEigenVectors->cols
+ << " " << __WeightsS2T << std::endl;
+ os << __MeanAppearance;
+ os << __AppearanceEigenValues;
+ os << __AppearanceEigenVectors;
+ os << __Qs;
+ os << __Qg;
+ os << __MeanS;
+ os << __MeanG;
+}
+
+//============================================================================
+void AAM_CAM::Read(std::ifstream& is)
+{
+ __shape.Read(is);
+ __texture.Read(is);
+ __paw.Read(is);
+
+ int np, nfeatures;
+ is >> np >> nfeatures >> __WeightsS2T;
+
+ __MeanAppearance = cvCreateMat(1, nfeatures, CV_64FC1);
+ __AppearanceEigenValues = cvCreateMat(1, np, CV_64FC1);
+ __AppearanceEigenVectors = cvCreateMat(np, nfeatures, CV_64FC1);
+ __Qs = cvCreateMat(np, __shape.nPoints()*2, CV_64FC1);
+ __Qg = cvCreateMat(np, __texture.nPixels()*__texture.nPlanes(), CV_64FC1);
+ __MeanS = cvCreateMat(1, __shape.nPoints()*2, CV_64FC1);
+ __MeanG = cvCreateMat(1, __texture.nPixels()*__texture.nPlanes(), CV_64FC1);
+
+ is >> __MeanAppearance;
+ is >> __AppearanceEigenValues;
+ is >> __AppearanceEigenVectors;
+ is >> __Qs;
+ is >> __Qg;
+ is >> __MeanS;
+ is >> __MeanG;
+
+ __Points = cvCreateMat (1, __shape.nPoints(), CV_32FC2);
+ __Storage = cvCreateMemStorage(0);
+ __pq = cvCreateMat(1, __shape.nModes()+4, CV_64FC1);
+ __a = cvCreateMat(1, __AppearanceEigenVectors->cols, CV_64FC1);
+}
+
+//============================================================================
+static AAM_CAM *g_cam;
+static const int n = 6;//appearance modes
+static int b_c[n];
+static const int offset = 40;
+static CvMat* c = 0;//conbined appearance parameters
+static CvMat* s = 0;//shape instance
+static CvMat* t = 0;//texture instance
+static IplImage* image = 0;
+static AAM_Shape aam_s;
+//============================================================================
+
+//============================================================================
+void ontrackcam(int pos)
+{
+ if(c == 0)
+ {
+ c = cvCreateMat(1, g_cam->nModes(), CV_64FC1);cvZero(c);
+ s = cvCreateMat(1, g_cam->__shape.nPoints()*2, CV_64FC1);
+ t = cvCreateMat(1, g_cam->__texture.nPixels()*g_cam->nPlanes(), CV_64FC1);
+ }
+
+ double var;
+ //registrate appearance parameters
+ for(int i = 0; i < n; i++)
+ {
+ var = 3*sqrt(g_cam->Var(i))*(double(b_c[i])/offset-1.0);
+ cvmSet(c, 0, i, var);
+ }
+
+ //generate shape and texture instance
+ g_cam->CalcLocalShape(s, c);
+ g_cam->CalcTexture(t, c);
+
+ //warp texture instance from base mesh to current shape instance
+ aam_s.Mat2Point(s);
+ int w = (int)aam_s.GetWidth(), h = (int)(aam_s.MaxY()-aam_s.MinY());
+ aam_s.Translate(w, h);
+ if(image == 0)image = cvCreateImage(cvSize(w*2,h*2), 8, 3);
+ cvSet(image, cvScalar(128, 128, 128));
+ g_cam->DrawAppearance(image, aam_s, t);
+
+ cvNamedWindow("Combined Appearance Model",1);
+ cvShowImage("Combined Appearance Model", image);
+
+ if(cvWaitKey(10) == 27)
+ {
+ cvReleaseImage(&image);
+ cvReleaseMat(&s);
+ cvReleaseMat(&t);
+ cvReleaseMat(&c);
+ cvDestroyWindow("Parameters");
+ cvDestroyWindow("Combined Appearance Model");
+ }
+}
+
+//============================================================================
+void AAM_CAM::ShowVariation()
+{
+ printf("Show modes of appearance variations...\n");
+ cvNamedWindow("Parameters",1);
+
+ //create trackbars for appearance
+ for(int i = 0; i < n; i++)
+ {
+ char barname[100];
+ sprintf(barname, "a %d", i);
+ b_c[i] = offset;
+ cvCreateTrackbar(barname, "Parameters", &b_c[i], 2*offset+1, ontrackcam);
+ }
+
+ g_cam = this;
+ ontrackcam(1);
+ cvWaitKey(0);
+}
118 src/aam-library/AAM_CAM.h
@@ -0,0 +1,118 @@
+/****************************************************************************
+* AAMLibrary
+* http://code.google.com/p/aam-library
+* Copyright (c) 2008-2009 by GreatYao, all rights reserved.
+****************************************************************************/
+
+#ifndef AAM_CAM_H
+#define AAM_CAM_H
+
+#include "AAM_Config.h"
+#include "AAM_TDM.h"
+#include "AAM_PDM.h"
+#include "AAM_Shape.h"
+#include "AAM_Util.h"
+
+class AAM_Basic;
+class AAM_Pyramid;
+
+//combined appearance model
+class AAM_EXPORTS AAM_CAM
+{
+ friend class AAM_Basic;
+ friend class AAM_Pyramid;
+
+public:
+ AAM_CAM();
+ ~AAM_CAM();
+
+ // Build combined appearance model
+ void Train(const file_lists& pts_files, const file_lists& img_files,
+ double scale = 1.0, int nPlane = 3,
+ double shape_percentage = 0.975,
+ double texture_percentage = 0.975, double appearance_percentage = 0.975);
+
+ // Get dimension of combined appearance vector
+ inline const int nParameters()const { return __AppearanceEigenVectors->cols;}
+
+ // Get number of modes of combined appearance variation
+ inline const int nModes()const { return __AppearanceEigenVectors->rows;}
+
+ // Get image planes of models
+ inline const int nPlanes()const { return __texture.nPlanes(); }
+
+ // Get variance of i-th mode of combined appearance variation
+ inline double Var(int i)const { return cvmGet(__AppearanceEigenValues, 0, i); }
+
+ // Get mean combined appearance
+ inline const CvMat* GetMean()const{ return __MeanAppearance; }
+
+ // Get combined appearance eigen-vectors of PCA (appearance modes)
+ inline const CvMat* GetBases()const{ return __AppearanceEigenVectors; }
+
+ // Show Model Variation according to various of parameters
+ void ShowVariation();
+ // function used in ShowVariation
+ friend void ontrackcam(int pos);
+
+ // Draw the image according the searching result
+ void DrawAppearance(IplImage* image, const AAM_Shape& Shape, CvMat* Texture);
+
+ // Calculate shape according to appearance parameters
+ void CalcLocalShape(CvMat* s, const CvMat* c);
+ void CalcGlobalShape(CvMat* s, const CvMat* pose);
+ inline void CalcShape(CvMat* s, const CvMat* c_q)
+ { CvMat c; cvGetCols(c_q, &c, 4, 4+nModes()); CalcLocalShape(s, &c);
+ CvMat q; cvGetCols(c_q, &q, 0, 4); CalcGlobalShape(s, &q); }
+ inline void CalcShape(CvMat* s, const CvMat* c, const CvMat* pose)
+ { CalcLocalShape(s, c); CalcGlobalShape(s, pose); }
+
+
+ // Calculate texture according to appearance parameters
+ void CalcTexture(CvMat* t, const CvMat* c);
+
+ //Calculate combined appearance parameters from shape and texture params.
+ void CalcParams(CvMat* c, const CvMat* p, const CvMat* lamda);
+
+ // Limit appearance parameters.
+ void Clamp(CvMat* c, double s_d = 3.0);
+
+ // Read data from stream
+ void Read(std::ifstream& is);
+
+ // Write data to stream
+ void Write(std::ofstream& os);
+
+private:
+ // Do PCA of appearance datas
+ void DoPCA(const CvMat* AllAppearances, double percentage);
+
+ // Convert shape and texture instance to appearance parameters
+ void ShapeTexture2Combined(const CvMat* Shape, const CvMat* Texture,
+ CvMat* Appearance);
+
+private:
+ AAM_PDM __shape; /*shape distribution model*/
+ AAM_TDM __texture; /*texture distribution model*/
+ AAM_PAW __paw; /*piecewise affine warp*/
+ double __WeightsS2T; /*ratio between shape and texture model*/
+
+ CvMat* __MeanAppearance;
+ CvMat* __AppearanceEigenValues;
+ CvMat* __AppearanceEigenVectors;
+
+ CvMat* __Qs;
+ CvMat* __Qg;
+ CvMat* __MeanS;
+ CvMat* __MeanG;
+
+ CvMat* __a;
+
+private:
+ //these cached variables are used for speed up
+ CvMat* __Points;
+ CvMemStorage* __Storage;
+ CvMat* __pq;
+};
+
+#endif // !AAM_CAM_H
73 src/aam-library/AAM_Config.h
@@ -0,0 +1,73 @@
+/****************************************************************************
+* AAMLibrary
+* http://code.google.com/p/aam-library
+* Copyright (c) 2008-2009 by GreatYao, all rights reserved.
+****************************************************************************/
+
+#ifndef AAM_CONFIG_H
+#define AAM_CONFIG_H
+
+#include <iostream>
+#include <vector>
+#include <string>
+#include <fstream>
+#include <cstdio>
+#include <cstdarg>
+#include <cstring>
+#include <cstdlib>
+#include <ctime>
+#include <cmath>
+
+#include <opencv2/opencv.hpp>
+
+#ifdef WIN32
+#include <direct.h>
+#include <io.h>
+#define AAM_MKDIR(dir) _mkdir(dir)
+#else
+#include <sys/stat.h>
+#define AAM_MKDIR(dir) mkdir(dir, 0644)
+#endif
+
+#ifdef WIN32
+#define AAM_EXTERN_C extern "C"
+#define AAM_EXPORTS __declspec(dllexport)
+#else
+#define AAM_EXTERN_C
+#define AAM_EXPORTS
+#endif
+#define AAM_API AAM_EXTERN_C AAM_EXPORTS
+
+#ifndef byte
+#define byte unsigned char
+#endif
+
+#define TYPE_AAM_BASIC 0
+#define TYPE_AAM_IC 1
+
+typedef std::vector<std::string> file_lists;
+
+extern char errmsg[];
+
+#define AAM_ERROR(msg) AAM_Error(msg, __FILE__, __LINE__)
+
+#ifdef __cplusplus
+extern "C"{
+#endif
+
+AAM_API double AAM_GetTime();
+
+AAM_API void AAM_FormatMSG(const char* format, ...);
+
+AAM_API void AAM_Error(const char *msg, const char* file, int line);
+
+AAM_API std::ostream& operator<<(std::ostream &os, const CvMat* mat);
+
+AAM_API std::istream& operator>>(std::istream &is, CvMat* mat);
+
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif //AAM_CONFIG_H
613 src/aam-library/AAM_IC.cpp
@@ -0,0 +1,613 @@
+/****************************************************************************
+* AAMLibrary
+* http://code.google.com/p/aam-library
+* Copyright (c) 2008-2009 by GreatYao, all rights reserved.
+****************************************************************************/
+
+#include "AAM_IC.h"
+
+using namespace std;
+
+//============================================================================
+AAM_IC::AAM_IC()
+{
+ __Points = 0;
+ __Storage = 0;
+
+ __update_s0 = 0;
+ __warp_t = 0;
+ __error_t = 0;
+ __search_pq = 0;
+ __delta_pq = 0;
+ __current_s = 0;
+ __update_s = 0;
+}
+
+//============================================================================
+AAM_IC::~AAM_IC()
+{
+ cvReleaseMat(&__Points);
+ cvReleaseMemStorage(&__Storage);
+
+ cvReleaseMat(&__update_s0);
+ cvReleaseMat(&__warp_t);
+ cvReleaseMat(&__error_t);
+ cvReleaseMat(&__search_pq);
+ cvReleaseMat(&__delta_pq);
+ cvReleaseMat(&__current_s);
+ cvReleaseMat(&__update_s);
+}
+
+//============================================================================
+CvMat* AAM_IC::CalcGradIdx()
+{
+ CvMat* pos= cvCreateMat(__paw.nPix(), 4, CV_32SC1);
+
+ int i = 0;
+ int width = __paw.Width(), height = __paw.Height();
+ for(int y = 0; y < height; y++)
+ {
+ for(int x = 0; x < width; x++)
+ {
+ if(__paw.Rect(y, x) >= 0)
+ {
+ int *ppos = (int*)(pos->data.ptr + pos->step*i);
+ ppos[0] = (x-1<0) ?-1:__paw.Rect(y, x-1); // left
+ ppos[1] = (x+1>=width) ?-1:__paw.Rect(y, x+1); // right
+ ppos[2] = (y-1<0) ?-1:__paw.Rect(y-1, x); // top
+ ppos[3] = (y+1>=height) ?-1:__paw.Rect(y+1, x); // bottom
+ i++;
+ }
+ }
+ }
+
+ return pos;
+}
+
+//============================================================================
+void AAM_IC::CalcTexGrad(const CvMat* texture, CvMat* dTx, CvMat* dTy)
+{
+ double* _x = dTx->data.db;
+ double* _y = dTy->data.db;
+ double* t = texture->data.db;
+ CvMat *idx = CalcGradIdx();
+ int nPlane = __texture.nPlanes();
+
+ if(nPlane == 3)
+ {
+ for(int i = 0; i < __paw.nPix(); i++)
+ {
+ int *fastp = (int*)(idx->data.ptr + idx->step*i);
+
+ // x direction
+ if(fastp[0] >= 0 && fastp[1] >= 0)
+ {
+ _x[3*i+0] = (t[3*fastp[1]+0] - t[3*fastp[0]+0])/2;
+ _x[3*i+1] = (t[3*fastp[1]+1] - t[3*fastp[0]+1])/2;
+ _x[3*i+2] = (t[3*fastp[1]+2] - t[3*fastp[0]+2])/2;
+ }
+
+ else if(fastp[0] >= 0 && fastp[1] < 0)
+ {
+ _x[3*i+0] = t[3*i+0] - t[3*fastp[0]+0];
+ _x[3*i+1] = t[3*i+1] - t[3*fastp[0]+1];
+ _x[3*i+2] = t[3*i+2] - t[3*fastp[0]+2];
+ }
+
+ else if(fastp[0] < 0 && fastp[1] >= 0)
+ {
+ _x[3*i+0] = t[3*fastp[1]+0] - t[3*i+0];
+ _x[3*i+1] = t[3*fastp[1]+1] - t[3*i+1];
+ _x[3*i+2] = t[3*fastp[1]+2] - t[3*i+2];
+ }
+ else
+ {
+ _x[3*i+0] = 0;
+ _x[3*i+1] = 0;
+ _x[3*i+2] = 0;
+ }
+
+ // y direction
+ if(fastp[2] >= 0 && fastp[3] >= 0)
+ {
+ _y[3*i+0] = (t[3*fastp[3]+0] - t[3*fastp[2]+0])/2;
+ _y[3*i+1] = (t[3*fastp[3]+1] - t[3*fastp[2]+1])/2;
+ _y[3*i+2] = (t[3*fastp[3]+2] - t[3*fastp[2]+2])/2;
+ }
+
+ else if(fastp[2] >= 0 && fastp[3] < 0)
+ {
+ _y[3*i+0] = t[3*i+0] - t[3*fastp[2]+0];
+ _y[3*i+1] = t[3*i+1] - t[3*fastp[2]+1];
+ _y[3*i+2] = t[3*i+2] - t[3*fastp[2]+2];
+ }
+
+ else if(fastp[2] < 0 && fastp[3] >= 0)
+ {
+ _y[3*i+0] = t[3*fastp[3]+0] - t[3*i+0];
+ _y[3*i+1] = t[3*fastp[3]+1] - t[3*i+1];
+ _y[3*i+2] = t[3*fastp[3]+2] - t[3*i+2];
+ }
+
+ else
+ {
+ _y[3*i+0] = 0;
+ _y[3*i+1] = 0;
+ _y[3*i+2] = 0;
+ }
+ }
+ }
+ else
+ {
+ for(int i = 0; i < __paw.nPix(); i++)
+ {
+ int *fastp = (int*)(idx->data.ptr + idx->step*i);
+
+ // x direction
+ if(fastp[0] >= 0 && fastp[1] >= 0)
+ {
+ _x[i] = (t[fastp[1]] - t[fastp[0]])/2;
+ }
+
+ else if(fastp[0] >= 0 && fastp[1] < 0)
+ {
+ _x[i] = t[i] - t[fastp[0]];
+ }
+
+ else if(fastp[0] < 0 && fastp[1] >= 0)
+ {
+ _x[i] = t[fastp[1]] - t[i];
+ }
+ else
+ {
+ _x[i] = 0;
+ }
+
+ // y direction
+ if(fastp[2] >= 0 && fastp[3] >= 0)
+ {
+ _y[i] = (t[fastp[3]] - t[fastp[2]])/2;
+ }
+
+ else if(fastp[2] >= 0 && fastp[3] < 0)
+ {
+ _y[i] = t[i] - t[fastp[2]];
+ }
+
+ else if(fastp[2] < 0 && fastp[3] >= 0)
+ {
+ _y[i] = t[fastp[3]] - t[i];
+ }
+
+ else
+ {
+ _y[i] = 0;
+ }
+ }
+ }
+ cvReleaseMat(&idx);
+}
+
+//============================================================================
+void AAM_IC::CalcWarpJacobian(CvMat* Jx, CvMat* Jy)
+{
+ int nPoints = __shape.nPoints();
+ __sMean.Mat2Point(__shape.GetMean());
+ __sStar1.resize(nPoints); __sStar2.resize(nPoints);
+ __sStar3.resize(nPoints); __sStar4.resize(nPoints);
+ for(int n = 0; n < nPoints; n++) // Equation (43)
+ {
+ __sStar1[n].x = __sMean[n].x; __sStar1[n].y = __sMean[n].y;
+ __sStar2[n].x = -__sMean[n].y; __sStar2[n].y = __sMean[n].x;
+ __sStar3[n].x = 1; __sStar3[n].y = 0;
+ __sStar4[n].x = 0; __sStar4[n].y = 1;
+ }
+
+ const CvMat* B = __shape.GetBases();
+ //const CvMat* mean = __shape.GetMean();
+ cvZero(Jx); cvZero(Jy);
+ for(int i = 0; i < __paw.nPix(); i++)
+ {
+ int tri_idx = __paw.PixTri(i);
+ int v1 = __paw.Tri(tri_idx, 0);
+ int v2 = __paw.Tri(tri_idx, 1);
+ int v3 = __paw.Tri(tri_idx, 2);
+ double *fastJx = (double*)(Jx->data.ptr + Jx->step*i);
+ double *fastJy = (double*)(Jy->data.ptr + Jy->step*i);
+
+ // Equation (50) dN_dq
+ fastJx[0] = __paw.Alpha(i)*__sStar1[v1].x +
+ __paw.Belta(i)*__sStar1[v2].x + __paw.Gamma(i)*__sStar1[v3].x;
+ fastJy[0] = __paw.Alpha(i)*__sStar1[v1].y +
+ __paw.Belta(i)*__sStar1[v2].y + __paw.Gamma(i)*__sStar1[v3].y;
+
+ fastJx[1] = __paw.Alpha(i)*__sStar2[v1].x +
+ __paw.Belta(i)*__sStar2[v2].x + __paw.Gamma(i)*__sStar2[v3].x;
+ fastJy[1] = __paw.Alpha(i)*__sStar2[v1].y +
+ __paw.Belta(i)*__sStar2[v2].y + __paw.Gamma(i)*__sStar2[v3].y;
+
+ fastJx[2] = __paw.Alpha(i)*__sStar3[v1].x +
+ __paw.Belta(i)*__sStar3[v2].x + __paw.Gamma(i)*__sStar3[v3].x;
+ fastJy[2] = __paw.Alpha(i)*__sStar3[v1].y +
+ __paw.Belta(i)*__sStar3[v2].y + __paw.Gamma(i)*__sStar3[v3].y;
+
+ fastJx[3] = __paw.Alpha(i)*__sStar4[v1].x +
+ __paw.Belta(i)*__sStar4[v2].x + __paw.Gamma(i)*__sStar4[v3].x;
+ fastJy[3] = __paw.Alpha(i)*__sStar4[v1].y +
+ __paw.Belta(i)*__sStar4[v2].y + __paw.Gamma(i)*__sStar4[v3].y;
+
+ // Equation (51) dW_dp
+ for(int j = 0; j < __shape.nModes(); j++)
+ {
+ fastJx[j+4] = __paw.Alpha(i)*cvmGet(B,j,2*v1) +
+ __paw.Belta(i)*cvmGet(B,j,2*v2) + __paw.Gamma(i)*cvmGet(B,j,2*v3);
+
+ fastJy[j+4] = __paw.Alpha(i)*cvmGet(B,j,2*v1+1) +
+ __paw.Belta(i)*cvmGet(B,j,2*v2+1) + __paw.Gamma(i)*cvmGet(B,j,2*v3+1);
+ }
+ }
+}
+
+//============================================================================
+void AAM_IC::CalcModifiedSD(CvMat* SD, const CvMat* dTx, const CvMat* dTy,
+ const CvMat* Jx, const CvMat* Jy)
+{
+ //create steepest descent images
+ int i, j;
+ int nPlane = __texture.nPlanes();
+ double* _x = dTx->data.db;
+ double* _y = dTy->data.db;
+ double temp;
+
+ if(nPlane == 3)
+ {
+ for(i = 0; i < __shape.nModes()+4; i++)
+ {
+ for(j = 0; j < __paw.nPix(); j++)
+ {
+ temp = _x[3*j ]*cvmGet(Jx,j,i) +_y[3*j ]*cvmGet(Jy,j,i);
+ cvmSet(SD,i,3*j,temp);
+
+ temp = _x[3*j+1]*cvmGet(Jx,j,i) +_y[3*j+1]*cvmGet(Jy,j,i);
+ cvmSet(SD,i,3*j+1,temp);
+
+ temp = _x[3*j+2]*cvmGet(Jx,j,i) +_y[3*j+2]*cvmGet(Jy,j,i);
+ cvmSet(SD,i,3*j+2,temp);
+ }
+ }
+ }
+ else
+ {
+ for(i = 0; i < __shape.nModes()+4; i++)
+ {
+ for(j = 0; j < __paw.nPix(); j++)
+ {
+ temp = _x[j]*cvmGet(Jx,j,i) +_y[j]*cvmGet(Jy,j,i);
+ cvmSet(SD,i,j,temp);
+ }
+ }
+ }
+
+ //project out appearance variation (and linear lighting parameters)
+ const CvMat* B = __texture.GetBases();
+ CvMat* V = cvCreateMat(4+__shape.nModes(), __texture.nModes(), CV_64FC1);
+ CvMat SDMat, BMat;
+
+ cvGEMM(SD, B, 1., NULL, 1., V, CV_GEMM_B_T);
+ // Equation (63),(64)
+ for(i = 0; i < __shape.nModes()+4; i++)
+ {
+ for(j = 0; j < __texture.nModes(); j++)
+ {
+ cvGetRow(SD, &SDMat, i);
+ cvGetRow(B, &BMat, j);
+ cvScaleAdd(&BMat, cvScalar(-cvmGet(V,i,j)), &SDMat, &SDMat);
+ }
+ }
+
+ cvReleaseMat(&V);
+}
+
+//============================================================================
+void AAM_IC::CalcHessian(CvMat* H, const CvMat* SD)
+{
+ CvMat* HH = cvCreateMat(H->rows, H->cols, CV_64FC1);
+ cvMulTransposed(SD, HH, 0);// Equation (65)
+ cvInvert(HH, H, CV_SVD);
+ cvReleaseMat(&HH);
+}
+
+//============================================================================
+void AAM_IC::Build(const file_lists& pts_files, const file_lists& img_files,
+ double scale /* = 1.0 */, int nPlane /* = 3 */)
+{
+ Train(pts_files, img_files, scale, nPlane);
+}
+
+//============================================================================
+void AAM_IC::Train(const file_lists& pts_files, const file_lists& img_files,
+ double scale /* = 1.0 */, int nPlane /* = 3 */,
+ double shape_percentage /* = 0.975 */,
+ double texture_percentage /* = 0.975 */)
+{
+ if(pts_files.size() != img_files.size())
+ {
+ AAM_FormatMSG("#Shapes != #Images");
+ AAM_ERROR(errmsg);
+ }
+
+ printf("################################################\n");
+ printf("Build Inverse Compositional Image Alignmennt Model...\n");
+
+ std::vector<AAM_Shape> AllShapes;
+ for(int ii = 0; ii < (int)pts_files.size(); ii++)
+ {
+ AAM_Shape Shape;
+ bool flag = Shape.ReadAnnotations(pts_files[ii]);
+ if(!flag)
+ {
+ IplImage* image = cvLoadImage(img_files[ii].c_str(), -1);
+ Shape.ScaleXY(image->width, image->height);
+ cvReleaseImage(&image);
+ }
+ AllShapes.push_back(Shape);
+ }
+
+ //building shape and texture distribution model
+ printf("Build point distribution model...\n");
+ __shape.Train(AllShapes, scale, shape_percentage);
+
+ printf("Build warp information of mean shape mesh...");
+ __Points = cvCreateMat (1, __shape.nPoints(), CV_32FC2);
+ __Storage = cvCreateMemStorage(0);
+
+ double sp = 1.0;
+ //if(__shape.GetMeanShape().GetWidth() > 120)
+ // sp = 120/__shape.GetMeanShape().GetWidth();
+
+ __paw.Train(__shape.GetMeanShape()*sp, __Points, __Storage);
+ printf("[%d by %d, triangles #%d, pixels #%d*%d]\n",
+ __paw.Width(), __paw.Height(), __paw.nTri(), __paw.nPix(), nPlane);
+
+ printf("Build texture distribution model...\n");
+ __texture.Train(pts_files, img_files, __paw, nPlane, texture_percentage, true);
+
+ //calculate gradient of texture
+ printf("Calculating texture gradient...\n");
+ CvMat* dTx = cvCreateMat(1, __texture.nPixels()*__texture.nPlanes(), CV_64FC1);
+ CvMat* dTy = cvCreateMat(1, __texture.nPixels()*__texture.nPlanes(), CV_64FC1);
+ CalcTexGrad(__texture.GetMean(), dTx, dTy);
+
+ // save gradient image
+ AAM_MKDIR("Modes");
+ __paw.SaveWarpTextureToImage("Modes/dTx.jpg", dTx, nPlane);
+ __paw.SaveWarpTextureToImage("Modes/dTy.jpg", dTy, nPlane);
+
+ //calculate warp Jacobian at base shape
+ printf("Calculating warp Jacobian...\n");
+ CvMat* Jx = cvCreateMat(__paw.nPix(), __shape.nModes()+4, CV_64FC1);
+ CvMat* Jy = cvCreateMat(__paw.nPix(), __shape.nModes()+4, CV_64FC1);
+ CalcWarpJacobian(Jx,Jy);
+
+ //calculate modified steepest descent image
+ printf("Calculating steepest descent images...\n");
+ CvMat* SD = cvCreateMat(__shape.nModes()+4, __texture.nPixels()*__texture.nPlanes(), CV_64FC1);
+ CalcModifiedSD(SD, dTx, dTy, Jx, Jy);
+
+ //calculate inverse Hessian matrix
+ printf("Calculating Hessian inverse matrix...\n");
+ CvMat* H = cvCreateMat(__shape.nModes()+4, __shape.nModes()+4, CV_64FC1);
+ CalcHessian(H, SD);
+
+ //calculate update matrix (multiply inverse Hessian by modified steepest descent image)
+ __G = cvCreateMat(__shape.nModes()+4, __texture.nPixels()*__texture.nPlanes(), CV_64FC1);
+ cvMatMul(H, SD, __G);
+
+ //release
+ cvReleaseMat(&Jx);
+ cvReleaseMat(&Jy);
+ cvReleaseMat(&dTx);
+ cvReleaseMat(&dTy);
+ cvReleaseMat(&SD);
+ cvReleaseMat(&H);
+
+ //alocate memory for on-line fitting stuff
+ __update_s0 = cvCreateMat(1, __shape.nPoints()*2, CV_64FC1);
+ __inv_pq = cvCreateMat(1, __shape.nModes()+4, CV_64FC1);
+ __warp_t = cvCreateMat(1, __texture.nPixels()*__texture.nPlanes(), CV_64FC1);
+ __error_t = cvCreateMat(1, __texture.nPixels()*__texture.nPlanes(), CV_64FC1);
+ __search_pq = cvCreateMat(1, __shape.nModes()+4, CV_64FC1);
+ __delta_pq = cvCreateMat(1, __shape.nModes()+4, CV_64FC1);
+ __current_s = cvCreateMat(1, __shape.nPoints()*2, CV_64FC1);
+ __update_s = cvCreateMat(1, __shape.nPoints()*2, CV_64FC1);
+ __lamda = cvCreateMat(1, __texture.nModes(), CV_64FC1);
+
+ printf("################################################\n\n");
+}
+
+//============================================================================
+void AAM_IC::Fit(const IplImage* image, AAM_Shape& Shape,
+ int max_iter /* = 30 */, bool showprocess /* = false */)
+{
+ //initialize some stuff
+ double t = AAM_GetTime();
+ const CvMat* A0 = __texture.GetMean();
+ CvMat p; cvGetCols(__search_pq, &p, 4, 4+__shape.nModes());
+ Shape.Point2Mat(__current_s);
+ SetAllParamsZero();
+ __shape.CalcParams(__current_s, __search_pq);
+ IplImage* Drawimg = 0;
+
+ for(int iter = 0; iter < max_iter; iter++)
+ {
+ if(showprocess)
+ {
+ if(Drawimg == 0) Drawimg = cvCloneImage(image);
+ else cvCopy(image, Drawimg);
+ Shape.Mat2Point(__current_s);
+ Draw(Drawimg, Shape, 2);
+ AAM_MKDIR("result");
+ char filename[100];
+ sprintf(filename, "result/Iter-%02d.jpg", iter);
+ cvSaveImage(filename, Drawimg);
+
+ }
+
+ //check the current shape
+ AAM_Common::CheckShape(__current_s, image->width, image->height);
+
+ //warp image to mesh shape mesh
+ __paw.CalcWarpTexture(__current_s, image, __warp_t, __texture.nPlanes());
+ AAM_TDM::NormalizeTexture(A0, __warp_t);
+ cvSub(__warp_t, A0, __error_t);
+
+ //calculate updates (and scale to account for linear lighting gain)
+ cvGEMM(__error_t, __G, 1, NULL, 1, __delta_pq, CV_GEMM_B_T);
+
+ //check for parameter convergence
+ if(cvNorm(__delta_pq) < 1e-6) break;
+
+ //apply inverse compositional algorithm to update parameters
+ InverseCompose(__delta_pq, __current_s, __update_s);
+
+ //smooth shape
+ cvAddWeighted(__current_s, 0.4, __update_s, 0.6, 0, __update_s);
+ //update parameters
+ __shape.CalcParams(__update_s, __search_pq);
+ //calculate constrained new shape
+ __shape.CalcShape(__search_pq, __update_s);
+
+ //check for shape convergence
+ if(cvNorm(__current_s, __update_s, CV_L2) < 0.001) break;
+ else cvCopy(__update_s, __current_s);
+ }
+
+ Shape.Mat2Point(__current_s);
+
+ t = AAM_GetTime() - t;
+ printf("AAM IC Fitting time cost %.3f millisec\n", t);
+
+ cvReleaseImage(&Drawimg);
+}
+
+//============================================================================
+void AAM_IC::SetAllParamsZero()
+{
+ cvZero(__warp_t);
+ cvZero(__error_t);
+ cvZero(__search_pq);
+ cvZero(__delta_pq);
+ cvZero(__lamda);
+}
+
+//============================================================================
+void AAM_IC::InverseCompose(const CvMat* dpq, const CvMat* s, CvMat* NewS)
+{
+ // Firstly: Estimate the corresponding changes to the base mesh
+ cvConvertScale(dpq, __inv_pq, -1);
+ __shape.CalcShape(__inv_pq, __update_s0); // __update_s0 = N.W(s0, -delta_p, -delta_q)
+
+ //Secondly: Composing the Incremental Warp with the Current Warp Estimate.
+ double *S0 = __update_s0->data.db;
+ double *S = s->data.db;
+ double *SEst = NewS->data.db;
+ double x, y, xw, yw;
+ int k, tri_idx;
+ int v1, v2, v3;
+ const std::vector<std::vector<int> >& tri = __paw.__tri;
+ const std::vector<std::vector<int> >& vtri = __paw.__vtri;
+
+ for(int i = 0; i < __shape.nPoints(); i++)
+ {
+ x = 0.0; y = 0.0;
+ k = 0;
+ //The only problem with this approach is which triangle do we use?
+ //In general there will be several triangles that share the i-th vertex.
+ for(k = 0; k < (int)vtri[i].size(); k++)// see Figure (11)
+ {
+ tri_idx = vtri[i][k];
+ v1 = tri[tri_idx][0];
+ v2 = tri[tri_idx][1];
+ v3 = tri[tri_idx][2];
+
+ AAM_PAW::Warp(S0[2*i],S0[2*i+1],
+ __sMean[v1].x, __sMean[v1].y,__sMean[v2].x, __sMean[v2].y,__sMean[v3].x, __sMean[v3].y,
+ xw, yw, S[2*v1], S[2*v1+1], S[2*v2], S[2*v2+1], S[2*v3], S[2*v3+1]);
+ x += xw; y += yw;
+ }
+ // average the result so as to smooth the warp at each vertex
+ SEst[2*i] = x/k; SEst[2*i+1] = y/k;
+ }
+}
+
+
+//============================================================================
+void AAM_IC::Draw(IplImage* image, const AAM_Shape& Shape, int type)
+{
+ if(type == 0) AAM_Common::DrawPoints(image, Shape);
+ else if(type == 1) AAM_Common::DrawTriangles(image, Shape, __paw.__tri);
+ else if(type == 2)
+ {
+ cvGEMM(__error_t, __texture.GetBases(), 1, NULL, 1, __lamda, CV_GEMM_B_T);
+ __texture.CalcTexture(__lamda, __warp_t);
+ AAM_PAW paw;
+ double minV, maxV;
+ cvMinMaxLoc(__warp_t, &minV, &maxV);
+ cvConvertScale(__warp_t, __warp_t, 255/(maxV-minV), -minV*255/(maxV-minV));
+ paw.Train(Shape, __Points, __Storage, __paw.GetTri(), false);
+ AAM_Common::DrawAppearance(image, Shape, __warp_t, paw, __paw, __texture.nPlanes());
+ }
+ else
+ {
+ AAM_FormatMSG("Unsupported drawing type\n");
+ AAM_ERROR(errmsg);
+ }
+}
+
+
+//============================================================================
+void AAM_IC::Write(std::ofstream& os)
+{
+ __shape.Write(os);
+ __texture.Write(os);
+ __paw.Write(os);
+
+ __sMean.Write(os);
+ __sStar1.Write(os); __sStar2.Write(os);
+ __sStar3.Write(os); __sStar4.Write(os);
+
+ os << __G << std::endl;
+}
+
+//============================================================================
+void AAM_IC::Read(std::ifstream& is)
+{
+ __shape.Read(is);
+ __texture.Read(is);
+ __paw.Read(is);
+
+ int nPoints = __shape.nPoints();
+ __sMean.resize(nPoints);
+ __sStar1.resize(nPoints); __sStar2.resize(nPoints);
+ __sStar3.resize(nPoints); __sStar4.resize(nPoints);
+ __sMean.Read(is);
+ __sStar1.Read(is); __sStar2.Read(is);
+ __sStar3.Read(is); __sStar4.Read(is);
+
+ __G = cvCreateMat(__shape.nModes()+4, __texture.nPixels()*__texture.nPlanes(), CV_64FC1);
+ is >> __G;
+
+ //alocate memory for on-line fitting stuff
+ __Points = cvCreateMat (1, __shape.nPoints(), CV_32FC2);
+ __Storage = cvCreateMemStorage(0);
+
+ __update_s0 = cvCreateMat(1, __shape.nPoints()*2, CV_64FC1);
+ __inv_pq = cvCreateMat(1, __shape.nModes()+4, CV_64FC1);
+ __warp_t = cvCreateMat(1, __texture.nPixels()*__texture.nPlanes(), CV_64FC1);
+ __error_t = cvCreateMat(1, __texture.nPixels()*__texture.nPlanes(), CV_64FC1);
+ __search_pq = cvCreateMat(1, __shape.nModes()+4, CV_64FC1);
+ __delta_pq = cvCreateMat(1, __shape.nModes()+4, CV_64FC1);
+ __current_s = cvCreateMat(1, __shape.nPoints()*2, CV_64FC1);
+ __update_s = cvCreateMat(1, __shape.nPoints()*2, CV_64FC1);
+ __lamda = cvCreateMat(1, __texture.nModes(), CV_64FC1);
+}
118 src/aam-library/AAM_IC.h
@@ -0,0 +1,118 @@
+/****************************************************************************
+* AAMLibrary
+* http://code.google.com/p/aam-library
+* Copyright (c) 2008-2009 by GreatYao, all rights reserved.
+****************************************************************************/
+
+
+#ifndef AAM_IC_H
+#define AAM_IC_H
+
+#include "AAM_Config.h"
+#include "AAM_Util.h"
+#include "AAM_Shape.h"
+#include "AAM_PDM.h"
+#include "AAM_TDM.h"
+#include "AAM_PAW.h"
+
+/**
+ Active appearace model using project-out inverse compositional fitting method.
+ Refer to: I. Matthews and S. Baker. "Active Appearance Models Revisited". 2004
+*/
+class AAM_EXPORTS AAM_IC : public AAM
+{
+public:
+ AAM_IC();
+ ~AAM_IC();
+
+ virtual const int GetType()const { return TYPE_AAM_IC; }
+
+ // Build aam inverse compositional model
+ virtual void Build(const file_lists& pts_files, const file_lists& img_files,
+ double scale = 1.0, int nPlane = 3);
+
+ void Train(const file_lists& pts_files, const file_lists& img_files,
+ double scale = 1.0, int nPlane = 3,
+ double shape_percentage = 0.975, double texture_percentage = 0.975);
+
+ // Fit the image using Inverse Compositional.
+ virtual void Fit(const IplImage* image, AAM_Shape& Shape,
+ int max_iter = 30, bool showprocess = false);
+
+ // Set all search parameters zero
+ virtual void SetAllParamsZero();
+
+ // Init search parameters
+ virtual void InitParams(const IplImage* image){}
+
+ // Draw the image according the searching result(0:point, 1:triangle, 2:appearance)
+ virtual void Draw(IplImage* image, const AAM_Shape& Shape, int type);
+
+ // Read data from stream
+ virtual void Read(std::ifstream& is);
+
+ // Write data to stream
+ virtual void Write(std::ofstream& os);
+
+ // Get Mean Shape of IC model
+ inline const AAM_Shape GetMeanShape()const{ return __sMean; }
+ const AAM_Shape GetReferenceShape()const{ return __paw.__referenceshape; }
+
+
+private:
+
+ // Calculate the texture parameters project to linear subspace span(A)
+ void CalcAppearanceVariation(const CvMat* error_t, CvMat* lamda);
+
+ // Evaluate the Jacobians dN_dq and dW_dp of piecewise affine warp at(x;0)
+ void CalcWarpJacobian(CvMat* Jx, CvMat* Jy);
+
+ // Calculate index of gradients for every point in texture.
+ // If point is outside texture, set to -1.
+ CvMat* CalcGradIdx();
+
+ // Calculate the gradient of texture template A0.
+ void CalcTexGrad(const CvMat* texture, CvMat* dTx, CvMat* dTy);
+
+ // Calculating the modified steepest descent image.
+ void CalcModifiedSD(CvMat* SD, const CvMat* dTx, const CvMat* dTy,
+ const CvMat* Jx, const CvMat* Jy);
+
+ // Inverse compose current warp with shape parameter update.
+ // Update warp N.W(x;p,q)<-- N.W(x;p,q) . N.W(x;delta_p,delta_q)^-1.
+ void InverseCompose(const CvMat* dpq, const CvMat* s, CvMat* NewS);
+
+ // Compute the Hessian matrix using modified steepest descent image.
+ void CalcHessian(CvMat* H, const CvMat* SD);
+
+private:
+
+ //these variables are used for train PAW
+ CvMat* __Points;
+ CvMemStorage* __Storage;
+
+private:
+
+ AAM_PDM __shape; /*shape distribution model*/
+ AAM_TDM __texture; /*shape distribution model*/
+ AAM_PAW __paw; /*piecewise affine warp*/
+ AAM_Shape __sMean; /*mean shape of model*/
+ AAM_Shape __sStar1, __sStar2, __sStar3, __sStar4;/*global shape transform bases*/
+ CvMat* __G; /*Update matrix*/
+ /*product of inverse Hessian with steepest descent image*/
+
+private:
+ //pre-allocated stuff for online alignment
+ CvMat* __update_s0; /*shape change at the base mesh */
+ CvMat* __inv_pq; /*inverse parameters at the base mesh*/
+
+ CvMat* __warp_t; /*warp image to base mesh*/
+ CvMat* __error_t; /*error between warp image and template image A0*/
+ CvMat* __search_pq; /*search parameters */
+ CvMat* __delta_pq; /*parameters change to be updated*/
+ CvMat* __current_s; /*current search shape*/
+ CvMat* __update_s; /*shape after composing the warp*/
+ CvMat* __lamda; /*appearance parameters*/
+};
+
+#endif
57 src/aam-library/AAM_MovieAVI.cpp
@@ -0,0 +1,57 @@
+/****************************************************************************
+* AAMLibrary
+* http://code.google.com/p/aam-library
+* Copyright (c) 2008-2009 by GreatYao, all rights reserved.
+****************************************************************************/
+
+#include "AAM_MovieAVI.h"
+
+using namespace std;
+
+//============================================================================
+AAM_MovieAVI::AAM_MovieAVI(): capture(0), capimg(0), image(0)
+{
+}
+
+AAM_MovieAVI::~AAM_MovieAVI()
+{
+ Close();
+}
+
+
+void AAM_MovieAVI::Open(const char* videofile)
+{
+ capture = cvCaptureFromAVI(videofile);
+ if(!capture)
+ {
+ AAM_FormatMSG("CANNOT open video file \"%s\" !\n", videofile);
+ AAM_ERROR(errmsg);
+ }
+
+ cvSetCaptureProperty(capture, CV_CAP_PROP_POS_FRAMES, 0);
+ capimg = cvQueryFrame( capture );
+ image = cvCreateImage(cvGetSize(capimg), capimg->depth, capimg->nChannels);
+}
+
+//============================================================================
+void AAM_MovieAVI::Close()
+{
+ cvReleaseCapture(&capture);
+ capture = 0;
+ cvReleaseImage(&image);
+ image = 0;
+}
+
+//============================================================================
+IplImage* AAM_MovieAVI:: ReadFrame(int frame_no )
+{
+ cvSetCaptureProperty(capture, CV_CAP_PROP_POS_FRAMES, frame_no);
+ capimg = cvQueryFrame( capture );
+
+ if(capimg->origin == 0)
+ cvCopy(capimg, image);
+ else
+ cvFlip(capimg, image);
+
+ return image;
+}
39 src/aam-library/AAM_MovieAVI.h
@@ -0,0 +1,39 @@
+/****************************************************************************
+* AAMLibrary
+* http://code.google.com/p/aam-library
+* Copyright (c) 2008-2009 by GreatYao, all rights reserved.
+****************************************************************************/
+
+#ifndef AAM_MOVIEAVI_H
+#define AAM_MOVIEAVI_H
+
+#include "AAM_Config.h"
+
+class AAM_EXPORTS AAM_MovieAVI
+{
+public:
+ AAM_MovieAVI();
+ ~AAM_MovieAVI();
+
+ // Open a AVI file
+ void Open(const char* videofile);
+
+ // Close it
+ void Close();
+
+ // Get concrete frame of the video
+ // Notice: for speed up you have no need to release the returned image
+ IplImage* ReadFrame(int frame_no = -1);
+
+ // frame count of this video
+ const int FrameCount()const
+ {return (int)cvGetCaptureProperty(capture, CV_CAP_PROP_FRAME_COUNT);}
+
+private:
+ CvCapture* capture;
+ IplImage* capimg;//captured from video
+ IplImage *image;
+
+};
+
+#endif // !
644 src/aam-library/AAM_PAW.cpp
@@ -0,0 +1,644 @@
+/****************************************************************************
+* AAMLibrary
+* http://code.google.com/p/aam-library
+* Copyright (c) 2008-2009 by GreatYao, all rights reserved.
+****************************************************************************/
+
+#include <set>
+
+#include "AAM_PAW.h"
+#include "AAM_TDM.h"
+
+#define BINLINEAR 1
+
+using namespace std;
+
+#define free2dvector(vec) \
+{ \
+ for(int i = 0; i < (int)vec.size(); i++) vec[i].clear(); \
+ vec.clear(); \
+}
+
+//============================================================================
+AAM_PAW::AAM_PAW()
+{
+ __nTriangles = 0;
+ __nPixels = 0;
+ __n = 0;
+ __width = 0;
+ __height = 0;
+}
+
+//============================================================================
+AAM_PAW::~AAM_PAW()
+{
+ __pixTri.clear();
+ __alpha.clear();
+ __belta.clear();
+ __gamma.clear();
+
+ free2dvector(__rect);
+ free2dvector(__vtri);
+ free2dvector(__tri);
+}
+
+//============================================================================
+void AAM_PAW::Train(const AAM_Shape& ReferenceShape,
+ CvMat* Points,
+ CvMemStorage* Storage,
+ const std::vector<std::vector<int> >* tri,
+ bool buildVtri)
+{
+ __referenceshape = ReferenceShape;
+
+ __n = __referenceshape.NPoints();// get the number of vertex point
+
+ for(int i = 0; i < __n; i++)
+ CV_MAT_ELEM(*Points, CvPoint2D32f, 0, i) = __referenceshape[i];
+
+ CvMat* ConvexHull = cvCreateMat (1, __n, CV_32FC2);
+ cvConvexHull2(Points, ConvexHull, CV_CLOCKWISE, 0);
+
+ CvRect rect = cvBoundingRect(ConvexHull, 0);
+ CvSubdiv2D* Subdiv = cvCreateSubdivDelaunay2D(rect, Storage);
+ for(int ii = 0; ii < __n; ii++)
+ cvSubdivDelaunay2DInsert(Subdiv, __referenceshape[ii]);
+
+ //firstly: build triangle
+ if(tri == 0) Delaunay(Subdiv, ConvexHull);
+ else __tri = *tri;
+ __nTriangles = __tri.size();// get the number of triangles
+
+ //secondly: build correspondence of Vertex-Triangle
+ if(buildVtri) FindVTri();
+
+ //Thirdly: build pixel point in all triangles
+ if(tri == 0) CalcPixelPoint(rect, ConvexHull);
+ else FastCalcPixelPoint(rect);
+ __nPixels = __pixTri.size();// get the number of pixels
+
+ cvReleaseMat(&ConvexHull);
+}
+
+//============================================================================
+void AAM_PAW::Delaunay(const CvSubdiv2D* Subdiv, const CvMat* ConvexHull)
+{
+ // firstly we build edges
+ int i;
+ CvSeqReader reader;
+ CvQuadEdge2D* edge;
+ CvPoint2D32f org, dst;
+ CvSubdiv2DPoint* org_pt, * dst_pt;
+ std::vector<std::vector<int> > edges;
+ std::vector<int> one_edge; one_edge.resize(2);
+ std::vector<int> one_tri; one_tri.resize(3);
+ int ind1, ind2;
+
+ cvStartReadSeq( (CvSeq*)(Subdiv->edges), &reader, 0 );
+ for(i = 0; i < Subdiv->edges->total; i++)
+ {
+ edge = (CvQuadEdge2D*)(reader.ptr);
+ if( CV_IS_SET_ELEM(edge)){
+ org_pt = cvSubdiv2DEdgeOrg((CvSubdiv2DEdge)edge);
+ dst_pt = cvSubdiv2DEdgeDst((CvSubdiv2DEdge)edge);
+
+ if( org_pt && dst_pt ){
+ org = org_pt->pt;
+ dst = dst_pt->pt;
+ if (cvPointPolygonTest(ConvexHull, org, 0) >= 0 &&
+ cvPointPolygonTest( ConvexHull, dst, 0) >= 0){
+ for (int j = 0; j < __n; j++){
+ if (fabs(org.x-__referenceshape[j].x)<1e-6 &&
+ fabs(org.y-__referenceshape[j].y)<1e-6)
+ {
+ for (int k = 0; k < __n; k++)
+ {
+ if (fabs(dst.x-__referenceshape[k].x)<1e-6
+ &&fabs(dst.y-__referenceshape[k].y)<1e-6)
+ {
+ one_edge[0] = j;
+ one_edge[1] = k;
+ edges.push_back (one_edge);
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ CV_NEXT_SEQ_ELEM( Subdiv->edges->elem_size, reader );
+ }
+
+ // secondly we start to build triangles
+ for (i = 0; i < (int)edges.size(); i++)
+ {
+ ind1 = edges[i][0];
+ ind2 = edges[i][1];
+
+ for (int j = 0; j < __n; j++)
+ {
+ // At most, there are only 2 triangles that can be added
+ if(AAM_PAW::IsEdgeIn(ind1, j, edges) && AAM_PAW::IsEdgeIn(ind2, j, edges) )
+ {
+ one_tri[0] = ind1;
+ one_tri[1] = ind2;
+ one_tri[2] = j;
+ if (AAM_PAW::IsTriangleNotIn(one_tri, __tri) )
+ {
+ __tri.push_back(one_tri);
+ }
+ }
+ }
+ }
+
+ //OK, up to now, we have already builded the triangles!
+}
+
+//============================================================================
+bool AAM_PAW::IsEdgeIn(int ind1, int ind2,
+ const std::vector<std::vector<int> > &edges)
+{
+ for (int i = 0; i < (int)edges.size (); i++)
+ {
+ if ((edges[i][0] == ind1 && edges[i][1] == ind2) ||
+ (edges[i][0] == ind2 && edges[i][1] == ind1) )
+ return true;
+ }
+
+ return false;
+}
+
+//============================================================================
+bool AAM_PAW::IsTriangleNotIn(const std::vector<int>& one_tri,
+ const std::vector<std::vector<int> > &tris)
+{
+ std::set<int> tTriangle;
+ std::set<int> sTriangle;
+
+ for (int i = 0; i < (int)tris.size (); i ++)
+ {
+ tTriangle.clear();
+ sTriangle.clear();
+ for (int j = 0; j < 3; j++ )
+ {
+ tTriangle.insert(tris[i][j]);
+ sTriangle.insert(one_tri[j]);
+ }
+ if (tTriangle == sTriangle) return false;
+ }
+
+ return true;
+}
+
+//============================================================================
+void AAM_PAW::CalcPixelPoint(const CvRect rect, CvMat* ConvexHull)//cost too much time
+{
+ CvPoint2D32f point[3];
+ CvMat tempVert = cvMat(1, 3, CV_32FC2, point);
+ int ll = 0;
+ double alpha, belta, gamma;
+ CvPoint2D32f pt;
+ int ind1, ind2, ind3;
+ int ii, jj;
+ double x, y, x1, y1, x2, y2, x3, y3, c;
+
+ __xmin = rect.x;
+ __ymin = rect.y;
+ __width = rect.width;
+ __height = rect.height;
+ int left = rect.x, right = left + __width;
+ int top = rect.y, bottom = top + __height;
+
+ __rect.resize(__height);
+ for (int i = top; i < bottom; i++)
+ {
+ ii = i - top;
+ __rect[ii].resize(__width);
+ for (int j = left; j < right; j++)
+ {
+ jj = j - left;
+ pt = cvPoint2D32f(j, i);
+ __rect[ii][jj] = -1;
+
+ // firstly: the point (j, i) is located in the ConvexHull
+ if(cvPointPolygonTest(ConvexHull, pt, 0) >= 0 )
+ {
+ // then we find the triangle that the point lies in
+ for (int k = 0; k < __nTriangles; k++)
+ {
+ ind1 = __tri[k][0];
+ ind2 = __tri[k][1];
+ ind3 = __tri[k][2];
+ point[0] = __referenceshape[ind1];
+ point[1] = __referenceshape[ind2];
+ point[2] = __referenceshape[ind3];
+
+ // secondly: the point(j,i) is located in the k-th triangle
+ if(cvPointPolygonTest(&tempVert, pt, 0) >= 0)
+ {
+ __rect[ii][jj] = ll++;
+ __pixTri.push_back(k);
+
+ // calculate alpha and belta for warp
+ x = j; y = i;
+ x1 = point[0].x; y1 = point[0].y;
+ x2 = point[1].x; y2 = point[1].y;
+ x3 = point[2].x; y3 = point[2].y,
+
+ c = 1.0/(+x2*y3-x2*y1-x1*y3-x3*y2+x3*y1+x1*y2);
+ alpha = (y*x3-y3*x+x*y2-x2*y+x2*y3-x3*y2)*c;
+ belta = (-y*x3+x1*y+x3*y1+y3*x-x1*y3-x*y1)*c;
+ gamma = 1 - alpha - belta;
+
+ __alpha.push_back(alpha);
+ __belta.push_back(belta);
+ __gamma.push_back(gamma);
+
+ // make sure each point only located in only one triangle
+ break;
+ }
+
+ }
+ }
+ }
+ }
+}
+
+//============================================================================
+void AAM_PAW::FastCalcPixelPoint(const CvRect rect)
+{
+ CvPoint2D32f point[3];
+ CvMat oneTri = cvMat(1, 3, CV_32FC2, point);
+ double alpha, belta, gamma;
+ CvPoint2D32f pt;
+ int ind1, ind2, ind3;
+ int ll = 0;
+ double x, y, x1, y1, x2, y2, x3, y3, c;
+
+ __xmin = rect.x; __ymin = rect.y;
+ __width = rect.width; __height = rect.height;
+ int left = rect.x, top = rect.y;
+ double aa, bb, cc, dd;
+
+ __rect.resize(__height);
+ for(int i = 0; i < __height; i++)
+ {
+ __rect[i].resize(__width);
+ for(int j = 0; j < __width; j++)
+ __rect[i][j] = -1;
+ }
+
+ for(int k = 0; k < __nTriangles; k++)
+ {
+ ind1 = __tri[k][0];
+ ind2 = __tri[k][1];
+ ind3 = __tri[k][2];
+
+ point[0] = __referenceshape[ind1];
+ point[1] = __referenceshape[ind2];
+ point[2] = __referenceshape[ind3];
+
+ x1 = point[0].x; y1 = point[0].y;
+ x2 = point[1].x; y2 = point[1].y;
+ x3 = point[2].x; y3 = point[2].y;
+ c = 1.0/(+x2*y3-x2*y1-x1*y3-x3*y2+x3*y1+x1*y2);
+
+ aa = MIN(x1, MIN(x2, x3)); //left x
+ bb = MAX(x1, MAX(x2, x3)); //right x
+ cc = MIN(y1, MIN(y2, y3)); //top y
+ dd = MAX(y1, MAX(y2, y3)); //bot y
+
+ for(y=ceil(cc); y<=dd;y +=1)
+ {
+ for(x=ceil(aa); x<=bb;x+=1)
+ {
+ pt = cvPoint2D32f(x, y);
+ //the point is located in the k-th triangle
+ if(cvPointPolygonTest(&oneTri, pt, 0) >= 0)
+ {
+ __rect[y-top][x-left] = ll++;
+ __pixTri.push_back(k);
+
+ alpha = (y*x3-y3*x+x*y2-x2*y+x2*y3-x3*y2)*c;
+ belta = (-y*x3+x1*y+x3*y1+y3*x-x1*y3-x*y1)*c;
+ gamma = 1 - alpha - belta;
+
+ __alpha.push_back(alpha);
+ __belta.push_back(belta);
+ __gamma.push_back(gamma);
+ }
+
+ }
+ }
+ }
+}
+
+//============================================================================
+void AAM_PAW::FindVTri()
+{
+ __vtri.resize(__n);
+ for(int i = 0; i < __n; i++)
+ {
+ for(int j = 0; j < __nTriangles; j++)
+ {
+ if(__tri[j][0] == i || __tri[j][1] == i || __tri[j][2] == i)
+ __vtri[i].push_back(j);
+ }
+ }
+}
+
+//============================================================================
+void AAM_PAW::CalcWarpParameters(double x, double y, double x1, double y1,
+ double x2, double y2, double x3, double y3,
+ double &alpha, double &belta, double &gamma)
+{
+ double c = (+x2*y3-x2*y1-x1*y3-x3*y2+x3*y1+x1*y2);
+ alpha = (y*x3-y3*x+x*y2-x2*y+x2*y3-x3*y2) / c;
+ belta = (-y*x3+x1*y+x3*y1+y3*x-x1*y3-x*y1) / c;
+ gamma = 1 - alpha - belta;
+}
+
+//============================================================================
+void AAM_PAW::Warp(double x, double y,
+ double x1, double y1, double x2, double y2, double x3, double y3,
+ double& X, double& Y,
+ double X1, double Y1, double X2, double Y2, double X3, double Y3)
+{
+ double c = 1.0/(+x2*y3-x2*y1-x1*y3-x3*y2+x3*y1+x1*y2);
+ double alpha = (y*x3-y3*x+x*y2-x2*y+x2*y3-x3*y2)*c;
+ double belta = (-y*x3+x1*y+x3*y1+y3*x-x1*y3-x*y1)*c;
+ double gamma = 1.0 - alpha - belta;
+
+ X = alpha*X1 + belta*X2 + gamma*X3;
+ Y = alpha*Y1 + belta*Y2 + gamma*Y3;
+}
+
+//==========================================================================
+void AAM_PAW::TextureToImage(IplImage* image, const CvMat* t, int nPlane)const
+{
+ CvMat* tt = cvCloneMat(t);
+ double minV, maxV;
+ cvMinMaxLoc(tt, &minV, &maxV);
+ cvConvertScale(tt, tt, 255/(maxV-minV), -minV*255/(maxV-minV));
+
+ int k, x3;
+ register double *T = tt->data.db;
+ byte* p;
+
+ if(nPlane == 3)
+ {
+ for(int y = 0; y < __height; y++)
+ {
+ p = (byte*)(image->imageData + image->widthStep*y);
+ for(int x = 0; x < __width; x++)
+ {
+ k = __rect[y][x];
+ if(k >= 0)
+ {
+ x3 = x+(x<<1); k = k+(k<<1);
+ p[x3 ] = (byte)T[k];
+ p[x3+1] = (byte)T[k+1];
+ p[x3+2] = (byte)T[k+2];
+ }
+ }
+ }
+ }
+ else
+ {
+ for(int y = 0; y < __height; y++)
+ {
+ p = (byte*)(image->imageData + image->widthStep*y);
+ for(int x = 0; x < __width; x++)
+ {
+ k = __rect[y][x];
+ if(k >= 0)
+ {
+ p[x] = (byte)T[k];
+ }
+ }
+ }
+ }
+
+ cvReleaseMat(&tt);
+}
+
+//==========================================================================
+void AAM_PAW::SaveWarpTextureToImage(const char* filename,
+ const CvMat* t, int nPlane)const
+{
+ IplImage* image = cvCreateImage(cvSize(__width, __height), IPL_DEPTH_8U, nPlane);
+ cvSetZero(image);
+ TextureToImage(image, t, nPlane);
+ cvSaveImage(filename, image);
+ cvReleaseImage(&image);
+}
+
+void AAM_PAW::CalcWarpTexture(const CvMat* s, const IplImage* image,
+ CvMat* t, int nPlane)const
+{
+ if(nPlane != 3 && nPlane != 1)
+ {
+ AAM_FormatMSG("Unsuported image plane parameter for CalcWarpTexture, only color/gray is supported!");
+ AAM_ERROR(errmsg);
+ }
+
+ int v1, v2, v3, tri_idx;
+ double x, y;
+ int X, Y, X1, Y1;
+ double s0 , t0, s1, t1;
+ int ixB1, ixG1, ixR1, ixB2, ixG2, ixR2;
+ register double *fastt = t->data.db;
+ register double *ss = s->data.db;
+ register byte* p1, * p2;
+ register char* imgdata = image->imageData;
+ register int step = image->widthStep;
+ register int nchannel = image->nChannels;
+ register int off_g = (nchannel == 3) ? 1 : 0;
+ register int off_r = (nchannel == 3) ? 2 : 0;
+ int i, k;
+ int w = image->width-1, h = image->height-1;
+
+ //std::cout << s << w << " " << h << std::endl;
+
+ switch(nPlane)
+ {
+ case 1:
+ for(i = 0, k = 0; i < __nPixels; i++, k++)
+ {
+ tri_idx = __pixTri[i];
+ v1 = __tri[tri_idx][0];
+ v2 = __tri[tri_idx][1];
+ v3 = __tri[tri_idx][2];
+
+ x = __alpha[i]*ss[v1<<1] + __belta[i]*ss[v2<<1] +
+ __gamma[i]*ss[v3<<1];
+ y = __alpha[i]*ss[1+(v1<<1)] + __belta[i]*ss[1+(v2<<1)] +
+ __gamma[i]*ss[1+(v3<<1)];
+
+ X = cvFloor(x); Y = cvFloor(y); X1 = cvCeil(x); Y1 = cvCeil(y);
+ s0 = x-X; t0 = y-Y; s1 = 1-s0; t1 = 1-t0;
+
+ p1 = (byte*)(imgdata + step*Y);
+ p2 = (byte*)(imgdata + step*Y1);
+
+ if(nchannel == 3)
+ {
+ ixB1 = nchannel*X; ixG1= ixB1+off_g; ixR1 = ixB1+off_r;
+ ixB2 = nchannel*X1; ixG2= ixB2+off_g; ixR2 = ixB2+off_r;
+
+ fastt[k] =
+ (s1*(t1*p1[ixB1]+t0*p2[ixB1])+s0*(t1*p1[ixB2]+t0*p2[ixB2]) +
+ s1*(t1*p1[ixG1]+t0*p2[ixG1])+s0*(t1*p1[ixG2]+t0*p2[ixG2]) +