SDM For Face Alignment 流程介紹及Matlab代碼實現之訓練篇

SDM 訓練階段的任務如下：

載入標準化的數據（包括400*400的正臉及特征點）
對每一張標準化的圖片，模擬人臉檢測儀，產生10個擾動的人臉框及相應的初始特征點 $x_0$ 。
求解 $\Delta x$ , $\Phi$ ,其中 $\Delta x=x_*-x_0$ , $x_*$ 表示true shape, $\Phi$ 表示每個特征點的特征向量
求解最小二乘問題，得到一系列 $\{R_k\}$

下面分別來說明：

載入數據

載入811個訓練數據，按照上一章預備篇關于第一幅圖片的裁剪方法裁剪這811張圖片。
matlab代碼如下：

function [Data] = load_single_data2 ( dbpath_img, dbpath_pts,image_index, options )%% output format
%{
DATA.
- width_orig: the width of the original image.
- height_orig: the height of the original image.
- img_gray: the crop image.
- height: the height of crop image.
- wdith: the width of crop image.
- shape_gt: ground-truth landmark.
- bbox_gt: bounding box of ground-truth.
%}
slash = options.slash;
dbname = options.datasetName;imlist = dir([dbpath_img slash '*.*g']);%% load imagesimg = im2uint8(imread([dbpath_img slash imlist(image_index).name]));Data.width_orig  = size(img,2);Data.height_orig = size(img,1);%% load shapeData.shape_gt = double(annotation_load(...[dbpath_pts slash imlist(image_index).name(1:end-3) 'pts'] , dbname));if 0figure(1); imshow(img); hold on;draw_shape(Data.shape_gt(:,1),...Data.shape_gt(:,2),'y');hold off;pause;end    %% get bounding boxData.bbox_gt = getbbox(Data.shape_gt);%% enlarge region of faceregion     = enlargingbbox(Data.bbox_gt, 2.0);region(2)  = double(max(region(2), 1));%這里主要是為了防止求出的包圍盒超過圖像，因此一旦超過，則region(2)必然小于0，因此此時取1即可。region(1)  = double(max(region(1), 1));bottom_y   = double(min(region(2) + region(4) - 1, ...Data.height_orig));right_x    = double(min(region(1) + region(3) - 1, ...Data.width_orig));%防止長和寬超過圖片大小，因此取二者最小值img_region = img(region(2):bottom_y, region(1):right_x, :);%取人臉區域%% recalculate(重新計算) the location of groundtruth shape and bounding boxData.shape_gt = bsxfun(@minus, Data.shape_gt,...double([region(1) region(2)]));%等價于Data{iimgs}.shape_gt-repeat( double([region(1) region(2)]),size(Data{iimgs}.shape_gt,1),1)%將圖像的坐標原點移到人臉包圍盒的左上角,并因此得以重新計算新的特征點Data.bbox_gt = getbbox(Data.shape_gt);%新的特征點的包圍盒的左上角坐標發生了改變，但是寬和高沒有變化if size(img_region, 3) == 1Data.img_gray = img_region;elseData.img_gray = rgb2gray(img_region);endData.width    = size(img_region, 2);Data.height   = size(img_region, 1);if 0figure(2); imshow(Data.img_gray); hold on;draw_shape(Data.shape_gt(:,1),...Data.shape_gt(:,2),'y');hold off;pause;end%% normalized the image to the mean-shapesr = options.canvasSize(1)/Data.width;sc = options.canvasSize(2)/Data.height;Data.img_gray = imresize(Data.img_gray,options.canvasSize);Data.width    = options.canvasSize(1);Data.height   = options.canvasSize(2);Data.shape_gt = bsxfun(@times, Data.shape_gt, [sr sc]); Data.bbox_gt(1:2) = bsxfun(@times, Data.bbox_gt(1:2), [sr sc]);%補充Data.bbox_gt(3:4) = bsxfun(@times, Data.bbox_gt(3:4), [sr sc]);%補充if 0figure(3); imshow(Data.img_gray); hold on;draw_shape(Data.shape_gt(:,1),...Data.shape_gt(:,2),'r');hold on; rectangle('Position',  Data.bbox_gt, 'EdgeColor', 'k');pause;end       endfunction region = enlargingbbox(bbox, scale)
%同前面一樣，初始時刻這里得到僅僅是特征點盒子，而我們如果想要包住整個人臉，就必須先將原始盒子的左上角平移一半的寬高，然后再放大兩倍。這個在前面求解
%rect = get_correct_region( boxes, shape,Dataa(i).img, 1 );中也用到過
%因此這里得到的盒子是包住全部人臉的盒子。    
region(1) = floor(bbox(1) - (scale - 1)/2*bbox(3));
region(2) = floor(bbox(2) - (scale - 1)/2*bbox(4));region(3) = floor(scale*bbox(3));
region(4) = floor(scale*bbox(4));end

模擬人臉檢測，產生10個初始值

事實上，每張圖片都有一個ground-truth poins,因此可以求出它的包圍盒，也可以通過opencv或其他的檢測器可以檢測出這樣的框來。但兩者有點不一樣。如下：

，我們可以對opencv的檢測盒做一些變換就可以得到近似的box gt了。
我們需要對包圍盒擾動，以產生更多的盒子。怎么擾動呢？

對于一個盒子，有四個屬性：x,y,width,height.因此我們只要產生10種屬性即可。或者，也可以從另外一種角度來考慮這個問題。假設新的盒子已產生，那么它與原來的盒子之間就會產生4個方向的偏差，因此我們只需要對這些偏差做估計即可。

事實上，我們通過對811張圖片的init shape 與ground truth shape求解偏差的均值與方差，以此可以產生兩個（分別是(x,y),(width,height)）二維正太分布，因此就可以產生正太分布的隨機數，于是10種屬性的偏差就產生了，然后加上原來盒子的屬性，就產生了10個擾動的盒子。再將mean shape對齊到10個盒子上產生了10個初始值。
do_learn_variation.m：用來產生偏差的均值和方差

function do_learn_variation( options )%% loading learned shape model
load([options.modelPath options.slash options.datasetName '_ShapeModel.mat']);imgDir = options.trainingImageDataPath;
ptsDir = options.trainingTruthDataPath;%% loading data
Data = load_data( imgDir, ptsDir, options );n = length(Data);transVec   = zeros(n,2);
scaleVec   = zeros(n,2);debug = 0;%% computing the translation and scale vectors %%%%%%%%%%%%%%%%%%%%%%%%%%%%for i = 1 : n%% the information of i-th imagedisp(Data(i).img);img   = imread(Data(i).img);shape = Data(i).shape;%% if detect face using viola opencv% boxes = detect_face( img , options );%% if using ground-truthboxes = [];%% predict the face boxrect = get_correct_region( boxes, shape,img, 1 );%% predict initial location[initX,initY,width,height] = init_face_location( rect );%注意：上面算出的人臉框比較大，一般是特征點包圍盒的4倍，因此上面算出的width和height分別是rect寬和高的一半，實際上從bounding_box的計算中可以看出，%特征點的包圍盒分別向左上和右下延伸了一半的寬和高，導致人臉的包圍盒的面積是特征點包圍盒的4倍.init_shape = align_init_shape(ShapeModel.MeanShape, ...initX, initY, width, height);if debugfigure(1); imshow(img); hold on;rectangle('Position',  rect, 'EdgeColor', 'g');draw_shape(init_shape.XY(1:2:end), init_shape.XY(2:2:end), 'y');%繪制每幅人臉圖上的平均人臉點hold on;plot(initX, initY, 'b*');%中心點draw_shape(shape(:,1), shape(:,2), 'r');hold off;pause;end[aligned_shape, cropIm] = align_to_mean_shape( ShapeModel, img , ...vec_2_shape(init_shape.XY) , options );%vec_2_shape將一維向量轉化為二維向量,獲取400*400下的圖像和在此標準下的真實人臉點和初始化人臉點[aligned_true_shape] = align_shape(aligned_shape.TransM,shape_2_vec(shape));if debugfigure(1); imshow(cropIm); hold on;draw_shape(aligned_shape.XY(1:2:end), ...aligned_shape.XY(2:2:end), 'y');draw_shape(aligned_true_shape(1:2:end), ...aligned_true_shape(2:2:end), 'r');%hold off;pause;end    initVector = vec_2_shape(aligned_shape.XY);trueVector = vec_2_shape(aligned_true_shape);%compute mean and covariance matrices of translation.%計算平移的平均值和協方差矩陣meanInitVector  = mean(initVector);meanTrueVector  = mean(trueVector);%compute bounding box sizeinitLeftTop     = min(initVector);initRightBottom = max(initVector);initFaceSize = abs(initLeftTop - initRightBottom);trueLeftTop     = min(trueVector);trueRightBottom = max(trueVector);trueFaceSize = abs(trueLeftTop - trueRightBottom);transVec(i,:) = (meanInitVector - meanTrueVector)./initFaceSize;%平移要除以一個標準的人臉大小是為了消除人臉大小帶來的不一致scaleVec(i,:) = initFaceSize./trueFaceSize;clear img;clear xy;%    endend%compute mean and covariance matrices of scale.%計算縮放的平均值和協方差矩陣
[mu_trans,cov_trans] = mean_covariance_of_data ( transVec );
[mu_scale,cov_scale] = mean_covariance_of_data ( scaleVec );DataVariation.mu_trans  = mu_trans;
DataVariation.cov_trans = cov_trans;
DataVariation.mu_scale  = mu_scale;
DataVariation.cov_scale = cov_scale;save([options.modelPath options.slash options.datasetName ...'_DataVariation.mat'], 'DataVariation');clear Data;end

random_init_position.m:產生10個盒子

function [rbbox] = random_init_position( bbox, ...DataVariation, nRandInit,options)rbbox(1,:) = bbox;    if nRandInit > 1center = bbox(1:2) + bbox(3:4)/2;                                            mu_trans  = DataVariation.mu_trans;
cov_trans = DataVariation.cov_trans;
mu_scale  = DataVariation.mu_scale;
cov_scale = DataVariation.cov_scale;rInit_trans = mvnrnd(mu_trans,cov_trans,nRandInit-1);
%rInit_trans = zeros(nRandInit-1,2);rCenter = repmat(center,nRandInit-1,1) + ...rInit_trans.*repmat([bbox(3) bbox(4)],nRandInit-1,1);rInit_scale = mvnrnd(mu_scale,cov_scale,nRandInit-1);%r = mvnrnd(MU,SIGMA,cases)——從均值為MU(1*d)，協方差矩陣為SIGMA(d*d)的正態分布中隨機抽取cases個樣本，返回cases*d的矩陣r。
%rInit_scale = ones(nRandInit-1,2);rWidth  = zeros(nRandInit-1,1);
rHeight = zeros(nRandInit-1,1);for i = 1 : nRandInit - 1rWidth(i)  = bbox(3)*rInit_scale(i,1);%原始是除rHeight(i) = bbox(4)*rInit_scale(i,2);
endrbbox(2:nRandInit,1:2) = rCenter - [rWidth(:,1) rHeight(:,1)]/2;
rbbox(2:nRandInit,3:4) = [rWidth(:,1) rHeight(:,1)];
%補充項，防止擾動超過圖片的邊界
rbbox(1:nRandInit,1:2)=max(rbbox(1:nRandInit,1:2),1);
rbbox(1:nRandInit,1:2)=min(rbbox(1:nRandInit,1:2)+rbbox(1:nRandInit,3:4),options.canvasSize(1) )-rbbox(1:nRandInit,3:4);
endend

resetshape.m：將shape_union對齊到bbox

function [shape_initial] = resetshape(bbox, shape_union)
%RESETSHAPE Summary of this function goes here
%   Function: reset the initial shape according to the groundtruth shape and union shape for all faces
%   Detailed explanation goes here
%   Input: 
%       bbox: bbounding box of groundtruth shape
%       shape_union: uniionshape
%   Output:
%       shape_initial: reset initial shape
%       bbox: bounding box of face image% get the bounding box according to the ground truth shape
width_union = (max(shape_union(:, 1)) - min(shape_union(:, 1)));
height_union = (max(shape_union(:, 2)) - min(shape_union(:, 2)));shape_union = bsxfun(@minus, (shape_union), (min(shape_union)));shape_initial = bsxfun(@times, shape_union, [(bbox(3)/width_union) (bbox(4)/height_union)]);
shape_initial = bsxfun(@plus, shape_initial, double([bbox(1) bbox(2)]));end

求解特征點之差和特征向量

上面我們對每幅圖片求得了10個初始特征點，這樣我們就很容易求解 $\Delta x$ 了。同樣對于特征向量 $\Phi$ ，我們也可以很容易地求出來。關于特征向量，又名描述子。我們可以選擇Sift特征或者Hog特征。
local_descriptors:求解特征向量

function [desc] = local_descriptors( img, xy, dsize, dbins, options )%計算描述子featType = options.descType;stage = options.current_cascade;dsize = options.descScale(stage) * size(img,1);if strcmp(featType,'raw')if size(img,3) == 3im = im2double(rgb2gray(uint8(img)));elseim = im2double(uint8(img));endfor ipts = 1 : nptsdesc(ipts,:) = raw(im,xy(ipts,:),desc_scale,desc_size);endelseif strcmp(featType,'xx_sift')%     i = randi([1 68],1,1);
%     rect = [xy(18,:) - [dsize/2 dsize/2] dsize dsize];
%     
%     if 1
%         figure(2); imshow(img); hold on;
%         rectangle('Position',  rect, 'EdgeColor', 'g');
%         hold off;
%         pause;
%     endif size(img,3) == 3im = im2double(rgb2gray(uint8(img)));elseim = im2double(uint8(img));endxy = xy - repmat(dsize/2,size(xy,1),2);desc = xx_sift(im,xy,'nsb',dbins,'winsize',dsize);elseif strcmp(featType,'hog')if size(img,3) == 3im = im2double(rgb2gray(uint8(img)));elseim = im2double(uint8(img));endnpts = size(xy,1);for ipts = 1 : npts%disp(ipts);if isempty(im)disp('empty im');endif isempty(dsize)disp('empty dsize');enddesc(ipts,:) = hog(im,xy(ipts,:),dsize);endendend

求解最小二乘問題

問題：

M i n | | Δ X ? R Φ | | 22

$Min ||\Delta X-R\Phi||_2^2$
其中

ΔX∈R(68?2)×n,Φ∈R(128?68)×n $\Delta X\in R^{(68*2)\times n},\Phi\in R^{(128*68)\times n}$
這里68為特征點的個數，128為每個特征點的特征向量的維數，n為樣本量，這里為811.
顯然這是個最小二乘問題，可以直接求解。

Δ X = R Φ Δ X Φ' = R Φ Φ' Δ X Φ' (Φ Φ' + λ I) ? 1 = R Φ Φ' (Φ Φ' + λ I) ? 1 Δ X Φ' (Φ Φ' + λ I) ? 1 = R

$\Delta X=R\Phi\\ \Delta X\Phi'=R\Phi\Phi'\\ \Delta X\Phi'(\Phi\Phi'+\lambda I)^{-1}=R\Phi\Phi'(\Phi\Phi'+\lambda I)^{-1}\\ \Delta X\Phi'(\Phi\Phi'+\lambda I)^{-1}=R$
也可以通過SVM方法求解，這里我們調用了 liblinear的SVR方法來求解。
linreg.m:求解最小二乘問題

function [R,lambda] = linreg( X , Y , lambda )%X = [ones(size(X,1),1) X];%% method 1: soving linear regression using close-form solution %%%%%%%%%%%% R = (X'*X+eye(size(X,2))*lambda)\X'*Y;%先是X'*Y,再是除法%% method 2: using SVR in liblinear %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
featdim   = size(X,2);
shapedim  = size(Y,2);param = sprintf('-s 12 -p 0 -c %f -q', lambda);
%param = sprintf('-s 12 -p 0 -c 0.3 -q');
R_tmp = zeros( featdim, shapedim );
tic;
for o = 1 : shapedimdisp(['Training landmarks ' num2str(o)]);model = train(Y(:,o),sparse(X),param);R_tmp(:,o) = model.w';
end
toc;R = R_tmp;end

后續的話，我們還需要根據求解的R來更新 $x_0$ ,進而更新 $\Delta x,\Phi$ ,
最后求解新的最小二乘問題，得到新的R，以此下去，迭代5步即可。
這時產生的 $\{R_k\}$ 就可以用來進行下一步的test了。如下為5次的迭代的特征點效果圖：
1st step
2nd step
3rd step
4th step
5th step
我們可以看到越往后迭代，產生的新的特征點就越接近true shape.