001/**
002 * FaceTracker Licence
003 * -------------------
004 * (Academic, non-commercial, not-for-profit licence)
005 *
006 * Copyright (c) 2010 Jason Mora Saragih
007 * All rights reserved.
008 *
009 * Redistribution and use in source and binary forms, with or without
010 * modification, are permitted provided that the following conditions are met:
011 *
012 *     * The software is provided under the terms of this licence stricly for
013 *       academic, non-commercial, not-for-profit purposes.
014 *     * Redistributions of source code must retain the above copyright notice,
015 *       this list of conditions (licence) and the following disclaimer.
016 *     * Redistributions in binary form must reproduce the above copyright
017 *       notice, this list of conditions (licence) and the following disclaimer
018 *       in the documentation and/or other materials provided with the
019 *       distribution.
020 *     * The name of the author may not be used to endorse or promote products
021 *       derived from this software without specific prior written permission.
022 *     * As this software depends on other libraries, the user must adhere to and
023 *       keep in place any licencing terms of those libraries.
024 *     * Any publications arising from the use of this software, including but
025 *       not limited to academic journal and conference publications, technical
026 *       reports and manuals, must cite the following work:
027 *
028 *       J. M. Saragih, S. Lucey, and J. F. Cohn. Face Alignment through Subspace
029 *       Constrained Mean-Shifts. International Journal of Computer Vision
030 *       (ICCV), September, 2009.
031 *
032 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR "AS IS" AND ANY EXPRESS OR IMPLIED
033 * WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
034 * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO
035 * EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
036 * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
037 * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
038 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY
039 * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
040 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
041 * EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
042 */
043package com.jsaragih;
044
045import java.io.BufferedReader;
046import java.io.BufferedWriter;
047import java.io.FileNotFoundException;
048import java.io.FileReader;
049import java.io.FileWriter;
050import java.io.IOException;
051import java.util.Scanner;
052
053import org.openimaj.image.FImage;
054import org.openimaj.image.processing.transform.RemapProcessor;
055
056import Jama.Matrix;
057
058/**
059 * Piecewise affine warp
060 * 
061 * @author Jason Mora Saragih
062 * @author Jonathon Hare (jsh2@ecs.soton.ac.uk)
063 */
064public class PAW {
065        static {
066                Tracker.init();
067        }
068
069        /** Number of pixels */
070        int _nPix;
071
072        /** Minimum x-coord for src */
073        double _xmin;
074
075        /** Minimum y-coord for src */
076        double _ymin;
077
078        /** Source points */
079        Matrix _src;
080
081        /** destination points */
082        Matrix _dst;
083
084        /** Triangulation */
085        int[][] _tri;
086
087        /** Triangle for each valid pixel */
088        int[][] _tridx;
089
090        /** Valid region mask */
091        FImage _mask;
092
093        /** affine coeffs for all triangles */
094        Matrix _coeff;
095
096        /** matrix of (c,x,y) coeffs for alpha */
097        Matrix _alpha;
098
099        /** matrix of (c,x,y) coeffs for alpha */
100        Matrix _beta;
101
102        /** x-destination of warped points */
103        FImage _mapx;
104
105        /** y-destination of warped points */
106        FImage _mapy;
107
108        boolean sameSide(double x0, double y0, double x1, double y1, double x2,
109                        double y2, double x3, double y3)
110        {
111                final double x = (x3 - x2) * (y0 - y2) - (x0 - x2) * (y3 - y2);
112                final double y = (x3 - x2) * (y1 - y2) - (x1 - x2) * (y3 - y2);
113
114                if (x * y >= 0)
115                        return true;
116                return false;
117        }
118
119        int isWithinTri(double x, double y, int[][] tri, Matrix shape) {
120                final int n = tri.length;
121                final int p = shape.getRowDimension() / 2;
122
123                for (int t = 0; t < n; t++) {
124                        final int i = tri[t][0];
125                        final int j = tri[t][1];
126                        final int k = tri[t][2];
127
128                        final double s11 = shape.get(i, 0);
129                        final double s21 = shape.get(j, 0);
130                        final double s31 = shape.get(k, 0);
131                        final double s12 = shape.get(i + p, 0);
132                        final double s22 = shape.get(j + p, 0);
133                        final double s32 = shape.get(k + p, 0);
134
135                        if (sameSide(x, y, s11, s12, s21, s22, s31, s32)
136                                        && sameSide(x, y, s21, s22, s11, s12, s31, s32)
137                                        && sameSide(x, y, s31, s32, s11, s12, s21, s22))
138                                return t;
139                }
140                return -1;
141        }
142
143        static PAW load(final String fname) throws FileNotFoundException {
144                BufferedReader br = null;
145                try {
146                        br = new BufferedReader(new FileReader(fname));
147                        final Scanner sc = new Scanner(br);
148                        return read(sc, true);
149                } finally {
150                        try {
151                                br.close();
152                        } catch (final IOException e) {
153                        }
154                }
155        }
156
157        void save(final String fname) throws IOException {
158                BufferedWriter bw = null;
159                try {
160                        bw = new BufferedWriter(new FileWriter(fname));
161
162                        write(bw);
163                } finally {
164                        try {
165                                if (bw != null)
166                                        bw.close();
167                        } catch (final IOException e) {
168                        }
169                }
170        }
171
172        void write(BufferedWriter s) throws IOException {
173                s.write(IO.Types.PAW.ordinal() + " " + _nPix + " " + _xmin + " "
174                                + _ymin + " ");
175
176                IO.writeMat(s, _src);
177                IO.writeIntArray(s, _tri);
178                IO.writeIntArray(s, _tridx);
179                IO.writeImg(s, _mask);
180                IO.writeMat(s, _alpha);
181                IO.writeMat(s, _beta);
182        }
183
184        static PAW read(Scanner s, boolean readType) {
185                if (readType) {
186                        final int type = s.nextInt();
187                        assert (type == IO.Types.PAW.ordinal());
188                }
189
190                final PAW paw = new PAW();
191                paw._nPix = s.nextInt();
192                paw._xmin = s.nextDouble();
193                paw._ymin = s.nextDouble();
194
195                paw._src = IO.readMat(s);
196                paw._tri = IO.readIntArray(s);
197                paw._tridx = IO.readIntArray(s);
198                paw._mask = IO.readImgByte(s);
199                paw._alpha = IO.readMat(s);
200                paw._beta = IO.readMat(s);
201
202                paw._mapx = new FImage(paw._mask.width, paw._mask.height);
203                paw._mapy = new FImage(paw._mask.width, paw._mask.height);
204
205                paw._coeff = new Matrix(paw.nTri(), 6);
206                paw._dst = paw._src;
207
208                return paw;
209        }
210
211        int nPoints() {
212                return _src.getRowDimension() / 2;
213        }
214
215        int nTri() {
216                return _tri.length;
217        }
218
219        int width() {
220                return _mask.width;
221        }
222
223        int height() {
224                return _mask.height;
225        }
226
227        PAW(Matrix src, int[][] tri) {
228                assert (src.getColumnDimension() == 1);
229                assert (tri[0].length == 3);
230
231                _src = src.copy();
232                _tri = tri.clone();
233
234                final int n = nPoints();
235
236                _alpha = new Matrix(nTri(), 3);
237                _beta = new Matrix(nTri(), 3);
238
239                for (int i = 0; i < nTri(); i++) {
240                        final int j = _tri[i][0];
241                        final int k = _tri[i][1];
242                        final int l = _tri[i][2];
243
244                        final double c1 = _src.get(l + n, 0) - _src.get(j + n, 0);
245                        final double c2 = _src.get(l, 0) - _src.get(j, 0);
246                        final double c4 = _src.get(k + n, 0) - _src.get(j + n, 0);
247                        final double c3 = _src.get(k, 0) - _src.get(j, 0);
248                        final double c5 = c3 * c1 - c2 * c4;
249
250                        _alpha.set(i, 0, (_src.get(j + n, 0) * c2 - _src.get(j, 0) * c1)
251                                        / c5);
252                        _alpha.set(i, 1, c1 / c5);
253                        _alpha.set(i, 2, -c2 / c5);
254
255                        _beta.set(i, 0, (_src.get(j, 0) * c4 - _src.get(j + n, 0) * c3)
256                                        / c5);
257                        _beta.set(i, 1, -c4 / c5);
258                        _beta.set(i, 2, c3 / c5);
259                }
260
261                double xmax, ymax, xmin, ymin;
262                xmax = xmin = _src.get(0, 0);
263                ymax = ymin = _src.get(n, 0);
264
265                for (int i = 0; i < n; i++) {
266                        final double vx = _src.get(i, 0);
267                        final double vy = _src.get(i + n, 0);
268
269                        xmax = Math.max(xmax, vx);
270                        ymax = Math.max(ymax, vy);
271                        xmin = Math.min(xmin, vx);
272                        ymin = Math.min(ymin, vy);
273                }
274
275                final int w = (int) (xmax - xmin + 1.0);
276                final int h = (int) (ymax - ymin + 1.0);
277                _mask = new FImage(w, h);
278                _tridx = new int[h][w];
279
280                for (int i = 0; i < h; i++) {
281                        for (int j = 0; j < w; j++) {
282                                if ((_tridx[i][j] = isWithinTri(j + xmin, i + ymin, tri, _src)) == -1) {
283                                        _mask.pixels[i][j] = 0;
284                                } else {
285                                        _mask.pixels[i][j] = 0;
286                                }
287                        }
288                }
289
290                _mapx = new FImage(_mask.width, _mask.height);
291                _mapy = new FImage(_mask.width, _mask.height);
292                _coeff = new Matrix(nTri(), 6);
293
294                _dst = _src;
295                _xmin = xmin;
296                _ymin = ymin;
297        }
298
299        PAW() {
300        }
301
302        void crop(FImage src, FImage dst, Matrix s) {
303                assert ((s.getRowDimension() == _src.getRowDimension()) && (s
304                                .getColumnDimension() == 1));
305
306                _dst = s;
307
308                calcCoeff();
309
310                warpRegion(_mapx, _mapy);
311
312                RemapProcessor.remap(src, dst, _mapx, _mapy);
313        }
314
315        void calcCoeff() {
316                final int p = nPoints();
317
318                for (int l = 0; l < nTri(); l++) {
319                        final int i = _tri[l][0];
320                        final int j = _tri[l][1];
321                        final int k = _tri[l][2];
322
323                        final double c1 = _dst.get(i, 0);
324                        final double c2 = _dst.get(j, 0) - c1;
325                        final double c3 = _dst.get(k, 0) - c1;
326                        final double c4 = _dst.get(i + p, 0);
327                        final double c5 = _dst.get(j + p, 0) - c4;
328                        final double c6 = _dst.get(k + p, 0) - c4;
329
330                        final double[] coeff = _coeff.getArray()[l];
331                        final double[] alpha = _alpha.getArray()[l];
332                        final double[] beta = _beta.getArray()[l];
333
334                        coeff[0] = c1 + c2 * alpha[0] + c3 * beta[0];
335                        coeff[1] = c2 * alpha[1] + c3 * beta[1];
336                        coeff[2] = c2 * alpha[2] + c3 * beta[2];
337                        coeff[3] = c4 + c5 * alpha[0] + c6 * beta[0];
338                        coeff[4] = c5 * alpha[1] + c6 * beta[1];
339                        coeff[5] = c5 * alpha[2] + c6 * beta[2];
340                }
341        }
342
343        void warpRegion(FImage mapx, FImage mapy) {
344                if ((mapx.height != _mask.height) || (mapx.width != _mask.width))
345                        _mapx.internalAssign(new FImage(_mask.width, _mask.height));
346
347                if ((mapy.height != _mask.height) || (mapy.width != _mask.width))
348                        _mapy.internalAssign(new FImage(_mask.width, _mask.height));
349
350                int k = -1;
351                double[] a = null, ap;
352
353                final float[][] xp = mapx.pixels;
354                final float[][] yp = mapy.pixels;
355                final float[][] mp = _mask.pixels;
356
357                for (int y = 0; y < _mask.height; y++) {
358                        final double yi = y + _ymin;
359
360                        for (int x = 0; x < _mask.width; x++) {
361                                final double xi = x + _xmin;
362
363                                if (mp[y][x] == 0) {
364                                        xp[y][x] = -1;
365                                        yp[y][x] = -1;
366                                } else {
367                                        final int j = _tridx[y][x];
368
369                                        if (j != k) {
370                                                a = _coeff.getArray()[j];
371                                                k = j;
372                                        }
373                                        ap = a;
374                                        double xo = ap[0];
375                                        xo += ap[1] * xi;
376                                        xp[y][x] = (float) (xo + ap[2] * yi);
377
378                                        double yo = ap[3];
379                                        yo += ap[4] * xi;
380                                        yp[y][x] = (float) (yo + ap[5] * yi);
381                                }
382                        }
383                }
384        }
385}