git.sesse.net Git - ffmpeg/blob - libavfilter/vf_v360.c

   1 /*
   2  * Copyright (c) 2019 Eugene Lyapustin
   3  *
   4  * This file is part of FFmpeg.
   5  *
   6  * FFmpeg is free software; you can redistribute it and/or
   7  * modify it under the terms of the GNU Lesser General Public
   8  * License as published by the Free Software Foundation; either
   9  * version 2.1 of the License, or (at your option) any later version.
  10  *
  11  * FFmpeg is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14  * Lesser General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU Lesser General Public
  17  * License along with FFmpeg; if not, write to the Free Software
  18  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  19  */
  20
  21 /**
  22  * @file
  23  * 360 video conversion filter.
  24  * Principle of operation:
  25  *
  26  * (for each pixel in output frame)
  27  * 1) Calculate OpenGL-like coordinates (x, y, z) for pixel position (i, j)
  28  * 2) Apply 360 operations (rotation, mirror) to (x, y, z)
  29  * 3) Calculate pixel position (u, v) in input frame
  30  * 4) Calculate interpolation window and weight for each pixel
  31  *
  32  * (for each frame)
  33  * 5) Remap input frame to output frame using precalculated data
  34  */
  35
  36 #include <math.h>
  37
  38 #include "libavutil/avassert.h"
  39 #include "libavutil/imgutils.h"
  40 #include "libavutil/pixdesc.h"
  41 #include "libavutil/opt.h"
  42 #include "avfilter.h"
  43 #include "formats.h"
  44 #include "internal.h"
  45 #include "video.h"
  46 #include "v360.h"
  47
  48 typedef struct ThreadData {
  49     AVFrame *in;
  50     AVFrame *out;
  51 } ThreadData;
  52
  53 #define OFFSET(x) offsetof(V360Context, x)
  54 #define FLAGS AV_OPT_FLAG_FILTERING_PARAM|AV_OPT_FLAG_VIDEO_PARAM
  55
  56 static const AVOption v360_options[] = {
  57     {     "input", "set input projection",              OFFSET(in), AV_OPT_TYPE_INT,    {.i64=EQUIRECTANGULAR}, 0,    NB_PROJECTIONS-1, FLAGS, "in" },
  58     {         "e", "equirectangular",                            0, AV_OPT_TYPE_CONST,  {.i64=EQUIRECTANGULAR}, 0,                   0, FLAGS, "in" },
  59     {  "equirect", "equirectangular",                            0, AV_OPT_TYPE_CONST,  {.i64=EQUIRECTANGULAR}, 0,                   0, FLAGS, "in" },
  60     {      "c3x2", "cubemap 3x2",                                0, AV_OPT_TYPE_CONST,  {.i64=CUBEMAP_3_2},     0,                   0, FLAGS, "in" },
  61     {      "c6x1", "cubemap 6x1",                                0, AV_OPT_TYPE_CONST,  {.i64=CUBEMAP_6_1},     0,                   0, FLAGS, "in" },
  62     {       "eac", "equi-angular cubemap",                       0, AV_OPT_TYPE_CONST,  {.i64=EQUIANGULAR},     0,                   0, FLAGS, "in" },
  63     {  "dfisheye", "dual fisheye",                               0, AV_OPT_TYPE_CONST,  {.i64=DUAL_FISHEYE},    0,                   0, FLAGS, "in" },
  64     {    "barrel", "barrel facebook's 360 format",               0, AV_OPT_TYPE_CONST,  {.i64=BARREL},          0,                   0, FLAGS, "in" },
  65     {        "fb", "barrel facebook's 360 format",               0, AV_OPT_TYPE_CONST,  {.i64=BARREL},          0,                   0, FLAGS, "in" },
  66     {      "c1x6", "cubemap 1x6",                                0, AV_OPT_TYPE_CONST,  {.i64=CUBEMAP_1_6},     0,                   0, FLAGS, "in" },
  67     {    "output", "set output projection",            OFFSET(out), AV_OPT_TYPE_INT,    {.i64=CUBEMAP_3_2},     0,    NB_PROJECTIONS-1, FLAGS, "out" },
  68     {         "e", "equirectangular",                            0, AV_OPT_TYPE_CONST,  {.i64=EQUIRECTANGULAR}, 0,                   0, FLAGS, "out" },
  69     {  "equirect", "equirectangular",                            0, AV_OPT_TYPE_CONST,  {.i64=EQUIRECTANGULAR}, 0,                   0, FLAGS, "out" },
  70     {      "c3x2", "cubemap 3x2",                                0, AV_OPT_TYPE_CONST,  {.i64=CUBEMAP_3_2},     0,                   0, FLAGS, "out" },
  71     {      "c6x1", "cubemap 6x1",                                0, AV_OPT_TYPE_CONST,  {.i64=CUBEMAP_6_1},     0,                   0, FLAGS, "out" },
  72     {       "eac", "equi-angular cubemap",                       0, AV_OPT_TYPE_CONST,  {.i64=EQUIANGULAR},     0,                   0, FLAGS, "out" },
  73     {      "flat", "regular video",                              0, AV_OPT_TYPE_CONST,  {.i64=FLAT},            0,                   0, FLAGS, "out" },
  74     {"rectilinear", "regular video",                             0, AV_OPT_TYPE_CONST,  {.i64=FLAT},            0,                   0, FLAGS, "out" },
  75     {  "gnomonic", "regular video",                              0, AV_OPT_TYPE_CONST,  {.i64=FLAT},            0,                   0, FLAGS, "out" },
  76     {    "barrel", "barrel facebook's 360 format",               0, AV_OPT_TYPE_CONST,  {.i64=BARREL},          0,                   0, FLAGS, "out" },
  77     {        "fb", "barrel facebook's 360 format",               0, AV_OPT_TYPE_CONST,  {.i64=BARREL},          0,                   0, FLAGS, "out" },
  78     {      "c1x6", "cubemap 1x6",                                0, AV_OPT_TYPE_CONST,  {.i64=CUBEMAP_1_6},     0,                   0, FLAGS, "out" },
  79     {        "sg", "stereographic",                              0, AV_OPT_TYPE_CONST,  {.i64=STEREOGRAPHIC},   0,                   0, FLAGS, "out" },
  80     {    "interp", "set interpolation method",      OFFSET(interp), AV_OPT_TYPE_INT,    {.i64=BILINEAR},        0, NB_INTERP_METHODS-1, FLAGS, "interp" },
  81     {      "near", "nearest neighbour",                          0, AV_OPT_TYPE_CONST,  {.i64=NEAREST},         0,                   0, FLAGS, "interp" },
  82     {   "nearest", "nearest neighbour",                          0, AV_OPT_TYPE_CONST,  {.i64=NEAREST},         0,                   0, FLAGS, "interp" },
  83     {      "line", "bilinear interpolation",                     0, AV_OPT_TYPE_CONST,  {.i64=BILINEAR},        0,                   0, FLAGS, "interp" },
  84     {    "linear", "bilinear interpolation",                     0, AV_OPT_TYPE_CONST,  {.i64=BILINEAR},        0,                   0, FLAGS, "interp" },
  85     {      "cube", "bicubic interpolation",                      0, AV_OPT_TYPE_CONST,  {.i64=BICUBIC},         0,                   0, FLAGS, "interp" },
  86     {     "cubic", "bicubic interpolation",                      0, AV_OPT_TYPE_CONST,  {.i64=BICUBIC},         0,                   0, FLAGS, "interp" },
  87     {      "lanc", "lanczos interpolation",                      0, AV_OPT_TYPE_CONST,  {.i64=LANCZOS},         0,                   0, FLAGS, "interp" },
  88     {   "lanczos", "lanczos interpolation",                      0, AV_OPT_TYPE_CONST,  {.i64=LANCZOS},         0,                   0, FLAGS, "interp" },
  89     {         "w", "output width",                   OFFSET(width), AV_OPT_TYPE_INT,    {.i64=0},               0,           INT16_MAX, FLAGS, "w"},
  90     {         "h", "output height",                 OFFSET(height), AV_OPT_TYPE_INT,    {.i64=0},               0,           INT16_MAX, FLAGS, "h"},
  91     { "in_forder", "input cubemap face order",   OFFSET(in_forder), AV_OPT_TYPE_STRING, {.str="rludfb"},        0,     NB_DIRECTIONS-1, FLAGS, "in_forder"},
  92     {"out_forder", "output cubemap face order", OFFSET(out_forder), AV_OPT_TYPE_STRING, {.str="rludfb"},        0,     NB_DIRECTIONS-1, FLAGS, "out_forder"},
  93     {   "in_frot", "input cubemap face rotation",  OFFSET(in_frot), AV_OPT_TYPE_STRING, {.str="000000"},        0,     NB_DIRECTIONS-1, FLAGS, "in_frot"},
  94     {  "out_frot", "output cubemap face rotation",OFFSET(out_frot), AV_OPT_TYPE_STRING, {.str="000000"},        0,     NB_DIRECTIONS-1, FLAGS, "out_frot"},
  95     {    "in_pad", "input cubemap pads",            OFFSET(in_pad), AV_OPT_TYPE_FLOAT,  {.dbl=0.f},           0.f,                 1.f, FLAGS, "in_pad"},
  96     {   "out_pad", "output cubemap pads",          OFFSET(out_pad), AV_OPT_TYPE_FLOAT,  {.dbl=0.f},           0.f,                 1.f, FLAGS, "out_pad"},
  97     {       "yaw", "yaw rotation",                     OFFSET(yaw), AV_OPT_TYPE_FLOAT,  {.dbl=0.f},        -180.f,               180.f, FLAGS, "yaw"},
  98     {     "pitch", "pitch rotation",                 OFFSET(pitch), AV_OPT_TYPE_FLOAT,  {.dbl=0.f},        -180.f,               180.f, FLAGS, "pitch"},
  99     {      "roll", "roll rotation",                   OFFSET(roll), AV_OPT_TYPE_FLOAT,  {.dbl=0.f},        -180.f,               180.f, FLAGS, "roll"},
 100     {    "rorder", "rotation order",                OFFSET(rorder), AV_OPT_TYPE_STRING, {.str="ypr"},           0,                   0, FLAGS, "rorder"},
 101     {     "h_fov", "horizontal field of view",       OFFSET(h_fov), AV_OPT_TYPE_FLOAT,  {.dbl=90.f},     0.00001f,               180.f, FLAGS, "h_fov"},
 102     {     "v_fov", "vertical field of view",         OFFSET(v_fov), AV_OPT_TYPE_FLOAT,  {.dbl=45.f},     0.00001f,                90.f, FLAGS, "v_fov"},
 103     {    "h_flip", "flip out video horizontally",   OFFSET(h_flip), AV_OPT_TYPE_BOOL,   {.i64=0},               0,                   1, FLAGS, "h_flip"},
 104     {    "v_flip", "flip out video vertically",     OFFSET(v_flip), AV_OPT_TYPE_BOOL,   {.i64=0},               0,                   1, FLAGS, "v_flip"},
 105     {    "d_flip", "flip out video indepth",        OFFSET(d_flip), AV_OPT_TYPE_BOOL,   {.i64=0},               0,                   1, FLAGS, "d_flip"},
 106     {   "ih_flip", "flip in video horizontally",   OFFSET(ih_flip), AV_OPT_TYPE_BOOL,   {.i64=0},               0,                   1, FLAGS, "ih_flip"},
 107     {   "iv_flip", "flip in video vertically",     OFFSET(iv_flip), AV_OPT_TYPE_BOOL,   {.i64=0},               0,                   1, FLAGS, "iv_flip"},
 108     {  "in_trans", "transpose video input",   OFFSET(in_transpose), AV_OPT_TYPE_BOOL,   {.i64=0},               0,                   1, FLAGS, "in_transpose"},
 109     { "out_trans", "transpose video output", OFFSET(out_transpose), AV_OPT_TYPE_BOOL,   {.i64=0},               0,                   1, FLAGS, "out_transpose"},
 110     { NULL }
 111 };
 112
 113 AVFILTER_DEFINE_CLASS(v360);
 114
 115 static int query_formats(AVFilterContext *ctx)
 116 {
 117     static const enum AVPixelFormat pix_fmts[] = {
 118         // YUVA444
 119         AV_PIX_FMT_YUVA444P,   AV_PIX_FMT_YUVA444P9,
 120         AV_PIX_FMT_YUVA444P10, AV_PIX_FMT_YUVA444P12,
 121         AV_PIX_FMT_YUVA444P16,
 122
 123         // YUVA422
 124         AV_PIX_FMT_YUVA422P,   AV_PIX_FMT_YUVA422P9,
 125         AV_PIX_FMT_YUVA422P10, AV_PIX_FMT_YUVA422P12,
 126         AV_PIX_FMT_YUVA422P16,
 127
 128         // YUVA420
 129         AV_PIX_FMT_YUVA420P,   AV_PIX_FMT_YUVA420P9,
 130         AV_PIX_FMT_YUVA420P10, AV_PIX_FMT_YUVA420P16,
 131
 132         // YUVJ
 133         AV_PIX_FMT_YUVJ444P, AV_PIX_FMT_YUVJ440P,
 134         AV_PIX_FMT_YUVJ422P, AV_PIX_FMT_YUVJ420P,
 135         AV_PIX_FMT_YUVJ411P,
 136
 137         // YUV444
 138         AV_PIX_FMT_YUV444P,   AV_PIX_FMT_YUV444P9,
 139         AV_PIX_FMT_YUV444P10, AV_PIX_FMT_YUV444P12,
 140         AV_PIX_FMT_YUV444P14, AV_PIX_FMT_YUV444P16,
 141
 142         // YUV440
 143         AV_PIX_FMT_YUV440P, AV_PIX_FMT_YUV440P10,
 144         AV_PIX_FMT_YUV440P12,
 145
 146         // YUV422
 147         AV_PIX_FMT_YUV422P,   AV_PIX_FMT_YUV422P9,
 148         AV_PIX_FMT_YUV422P10, AV_PIX_FMT_YUV422P12,
 149         AV_PIX_FMT_YUV422P14, AV_PIX_FMT_YUV422P16,
 150
 151         // YUV420
 152         AV_PIX_FMT_YUV420P,   AV_PIX_FMT_YUV420P9,
 153         AV_PIX_FMT_YUV420P10, AV_PIX_FMT_YUV420P12,
 154         AV_PIX_FMT_YUV420P14, AV_PIX_FMT_YUV420P16,
 155
 156         // YUV411
 157         AV_PIX_FMT_YUV411P,
 158
 159         // YUV410
 160         AV_PIX_FMT_YUV410P,
 161
 162         // GBR
 163         AV_PIX_FMT_GBRP,   AV_PIX_FMT_GBRP9,
 164         AV_PIX_FMT_GBRP10, AV_PIX_FMT_GBRP12,
 165         AV_PIX_FMT_GBRP14, AV_PIX_FMT_GBRP16,
 166
 167         // GBRA
 168         AV_PIX_FMT_GBRAP,   AV_PIX_FMT_GBRAP10,
 169         AV_PIX_FMT_GBRAP12, AV_PIX_FMT_GBRAP16,
 170
 171         // GRAY
 172         AV_PIX_FMT_GRAY8,  AV_PIX_FMT_GRAY9,
 173         AV_PIX_FMT_GRAY10, AV_PIX_FMT_GRAY12,
 174         AV_PIX_FMT_GRAY14, AV_PIX_FMT_GRAY16,
 175
 176         AV_PIX_FMT_NONE
 177     };
 178
 179     AVFilterFormats *fmts_list = ff_make_format_list(pix_fmts);
 180     if (!fmts_list)
 181         return AVERROR(ENOMEM);
 182     return ff_set_common_formats(ctx, fmts_list);
 183 }
 184
 185 #define DEFINE_REMAP1_LINE(bits, div)                                                           \
 186 static void remap1_##bits##bit_line_c(uint8_t *dst, int width, const uint8_t *src,              \
 187                                       ptrdiff_t in_linesize,                                    \
 188                                       const uint16_t *u, const uint16_t *v, const int16_t *ker) \
 189 {                                                                                               \
 190     const uint##bits##_t *s = (const uint##bits##_t *)src;                                      \
 191     uint##bits##_t *d = (uint##bits##_t *)dst;                                                  \
 192                                                                                                 \
 193     in_linesize /= div;                                                                         \
 194                                                                                                 \
 195     for (int x = 0; x < width; x++)                                                             \
 196         d[x] = s[v[x] * in_linesize + u[x]];                                                    \
 197 }
 198
 199 DEFINE_REMAP1_LINE( 8, 1)
 200 DEFINE_REMAP1_LINE(16, 2)
 201
 202 typedef struct XYRemap {
 203     uint16_t u[4][4];
 204     uint16_t v[4][4];
 205     float ker[4][4];
 206 } XYRemap;
 207
 208 /**
 209  * Generate remapping function with a given window size and pixel depth.
 210  *
 211  * @param ws size of interpolation window
 212  * @param bits number of bits per pixel
 213  */
 214 #define DEFINE_REMAP(ws, bits)                                                                             \
 215 static int remap##ws##_##bits##bit_slice(AVFilterContext *ctx, void *arg, int jobnr, int nb_jobs)          \
 216 {                                                                                                          \
 217     ThreadData *td = (ThreadData*)arg;                                                                     \
 218     const V360Context *s = ctx->priv;                                                                      \
 219     const AVFrame *in = td->in;                                                                            \
 220     AVFrame *out = td->out;                                                                                \
 221                                                                                                            \
 222     for (int plane = 0; plane < s->nb_planes; plane++) {                                                   \
 223         const int in_linesize  = in->linesize[plane];                                                      \
 224         const int out_linesize = out->linesize[plane];                                                     \
 225         const uint8_t *src = in->data[plane];                                                              \
 226         uint8_t *dst = out->data[plane];                                                                   \
 227         const int width = s->planewidth[plane];                                                            \
 228         const int height = s->planeheight[plane];                                                          \
 229                                                                                                            \
 230         const int slice_start = (height *  jobnr     ) / nb_jobs;                                          \
 231         const int slice_end   = (height * (jobnr + 1)) / nb_jobs;                                          \
 232                                                                                                            \
 233         for (int y = slice_start; y < slice_end; y++) {                                                    \
 234             const unsigned map = s->map[plane];                                                            \
 235             const uint16_t *u = s->u[map] + y * width * ws * ws;                                           \
 236             const uint16_t *v = s->v[map] + y * width * ws * ws;                                           \
 237             const int16_t *ker = s->ker[map] + y * width * ws * ws;                                        \
 238                                                                                                            \
 239             s->remap_line(dst + y * out_linesize, width, src, in_linesize, u, v, ker);                     \
 240         }                                                                                                  \
 241     }                                                                                                      \
 242                                                                                                            \
 243     return 0;                                                                                              \
 244 }
 245
 246 DEFINE_REMAP(1,  8)
 247 DEFINE_REMAP(2,  8)
 248 DEFINE_REMAP(4,  8)
 249 DEFINE_REMAP(1, 16)
 250 DEFINE_REMAP(2, 16)
 251 DEFINE_REMAP(4, 16)
 252
 253 #define DEFINE_REMAP_LINE(ws, bits, div)                                                                   \
 254 static void remap##ws##_##bits##bit_line_c(uint8_t *dst, int width, const uint8_t *src,                    \
 255                                            ptrdiff_t in_linesize,                                          \
 256                                            const uint16_t *u, const uint16_t *v, const int16_t *ker)       \
 257 {                                                                                                          \
 258     const uint##bits##_t *s = (const uint##bits##_t *)src;                                                 \
 259     uint##bits##_t *d = (uint##bits##_t *)dst;                                                             \
 260                                                                                                            \
 261     in_linesize /= div;                                                                                    \
 262                                                                                                            \
 263     for (int x = 0; x < width; x++) {                                                                      \
 264         const uint16_t *uu = u + x * ws * ws;                                                              \
 265         const uint16_t *vv = v + x * ws * ws;                                                              \
 266         const int16_t *kker = ker + x * ws * ws;                                                           \
 267         int tmp = 0;                                                                                       \
 268                                                                                                            \
 269         for (int i = 0; i < ws; i++) {                                                                     \
 270             for (int j = 0; j < ws; j++) {                                                                 \
 271                 tmp += kker[i * ws + j] * s[vv[i * ws + j] * in_linesize + uu[i * ws + j]];                \
 272             }                                                                                              \
 273         }                                                                                                  \
 274                                                                                                            \
 275         d[x] = av_clip_uint##bits(tmp >> 14);                                                              \
 276     }                                                                                                      \
 277 }
 278
 279 DEFINE_REMAP_LINE(2,  8, 1)
 280 DEFINE_REMAP_LINE(4,  8, 1)
 281 DEFINE_REMAP_LINE(2, 16, 2)
 282 DEFINE_REMAP_LINE(4, 16, 2)
 283
 284 void ff_v360_init(V360Context *s, int depth)
 285 {
 286     switch (s->interp) {
 287     case NEAREST:
 288         s->remap_line = depth <= 8 ? remap1_8bit_line_c : remap1_16bit_line_c;
 289         break;
 290     case BILINEAR:
 291         s->remap_line = depth <= 8 ? remap2_8bit_line_c : remap2_16bit_line_c;
 292         break;
 293     case BICUBIC:
 294     case LANCZOS:
 295         s->remap_line = depth <= 8 ? remap4_8bit_line_c : remap4_16bit_line_c;
 296         break;
 297     }
 298
 299     if (ARCH_X86)
 300         ff_v360_init_x86(s, depth);
 301 }
 302
 303 /**
 304  * Save nearest pixel coordinates for remapping.
 305  *
 306  * @param du horizontal relative coordinate
 307  * @param dv vertical relative coordinate
 308  * @param r_tmp calculated 4x4 window
 309  * @param u u remap data
 310  * @param v v remap data
 311  * @param ker ker remap data
 312  */
 313 static void nearest_kernel(float du, float dv, const XYRemap *r_tmp,
 314                            uint16_t *u, uint16_t *v, int16_t *ker)
 315 {
 316     const int i = roundf(dv) + 1;
 317     const int j = roundf(du) + 1;
 318
 319     u[0] = r_tmp->u[i][j];
 320     v[0] = r_tmp->v[i][j];
 321 }
 322
 323 /**
 324  * Calculate kernel for bilinear interpolation.
 325  *
 326  * @param du horizontal relative coordinate
 327  * @param dv vertical relative coordinate
 328  * @param r_tmp calculated 4x4 window
 329  * @param u u remap data
 330  * @param v v remap data
 331  * @param ker ker remap data
 332  */
 333 static void bilinear_kernel(float du, float dv, const XYRemap *r_tmp,
 334                             uint16_t *u, uint16_t *v, int16_t *ker)
 335 {
 336     int i, j;
 337
 338     for (i = 0; i < 2; i++) {
 339         for (j = 0; j < 2; j++) {
 340             u[i * 2 + j] = r_tmp->u[i + 1][j + 1];
 341             v[i * 2 + j] = r_tmp->v[i + 1][j + 1];
 342         }
 343     }
 344
 345     ker[0] = (1.f - du) * (1.f - dv) * 16384;
 346     ker[1] =        du  * (1.f - dv) * 16384;
 347     ker[2] = (1.f - du) *        dv  * 16384;
 348     ker[3] =        du  *        dv  * 16384;
 349 }
 350
 351 /**
 352  * Calculate 1-dimensional cubic coefficients.
 353  *
 354  * @param t relative coordinate
 355  * @param coeffs coefficients
 356  */
 357 static inline void calculate_bicubic_coeffs(float t, float *coeffs)
 358 {
 359     const float tt  = t * t;
 360     const float ttt = t * t * t;
 361
 362     coeffs[0] =     - t / 3.f + tt / 2.f - ttt / 6.f;
 363     coeffs[1] = 1.f - t / 2.f - tt       + ttt / 2.f;
 364     coeffs[2] =       t       + tt / 2.f - ttt / 2.f;
 365     coeffs[3] =     - t / 6.f            + ttt / 6.f;
 366 }
 367
 368 /**
 369  * Calculate kernel for bicubic interpolation.
 370  *
 371  * @param du horizontal relative coordinate
 372  * @param dv vertical relative coordinate
 373  * @param r_tmp calculated 4x4 window
 374  * @param u u remap data
 375  * @param v v remap data
 376  * @param ker ker remap data
 377  */
 378 static void bicubic_kernel(float du, float dv, const XYRemap *r_tmp,
 379                            uint16_t *u, uint16_t *v, int16_t *ker)
 380 {
 381     int i, j;
 382     float du_coeffs[4];
 383     float dv_coeffs[4];
 384
 385     calculate_bicubic_coeffs(du, du_coeffs);
 386     calculate_bicubic_coeffs(dv, dv_coeffs);
 387
 388     for (i = 0; i < 4; i++) {
 389         for (j = 0; j < 4; j++) {
 390             u[i * 4 + j] = r_tmp->u[i][j];
 391             v[i * 4 + j] = r_tmp->v[i][j];
 392             ker[i * 4 + j] = du_coeffs[j] * dv_coeffs[i] * 16384;
 393         }
 394     }
 395 }
 396
 397 /**
 398  * Calculate 1-dimensional lanczos coefficients.
 399  *
 400  * @param t relative coordinate
 401  * @param coeffs coefficients
 402  */
 403 static inline void calculate_lanczos_coeffs(float t, float *coeffs)
 404 {
 405     int i;
 406     float sum = 0.f;
 407
 408     for (i = 0; i < 4; i++) {
 409         const float x = M_PI * (t - i + 1);
 410         if (x == 0.f) {
 411             coeffs[i] = 1.f;
 412         } else {
 413             coeffs[i] = sinf(x) * sinf(x / 2.f) / (x * x / 2.f);
 414         }
 415         sum += coeffs[i];
 416     }
 417
 418     for (i = 0; i < 4; i++) {
 419         coeffs[i] /= sum;
 420     }
 421 }
 422
 423 /**
 424  * Calculate kernel for lanczos interpolation.
 425  *
 426  * @param du horizontal relative coordinate
 427  * @param dv vertical relative coordinate
 428  * @param r_tmp calculated 4x4 window
 429  * @param u u remap data
 430  * @param v v remap data
 431  * @param ker ker remap data
 432  */
 433 static void lanczos_kernel(float du, float dv, const XYRemap *r_tmp,
 434                            uint16_t *u, uint16_t *v, int16_t *ker)
 435 {
 436     int i, j;
 437     float du_coeffs[4];
 438     float dv_coeffs[4];
 439
 440     calculate_lanczos_coeffs(du, du_coeffs);
 441     calculate_lanczos_coeffs(dv, dv_coeffs);
 442
 443     for (i = 0; i < 4; i++) {
 444         for (j = 0; j < 4; j++) {
 445             u[i * 4 + j] = r_tmp->u[i][j];
 446             v[i * 4 + j] = r_tmp->v[i][j];
 447             ker[i * 4 + j] = du_coeffs[j] * dv_coeffs[i] * 16384;
 448         }
 449     }
 450 }
 451
 452 /**
 453  * Modulo operation with only positive remainders.
 454  *
 455  * @param a dividend
 456  * @param b divisor
 457  *
 458  * @return positive remainder of (a / b)
 459  */
 460 static inline int mod(int a, int b)
 461 {
 462     const int res = a % b;
 463     if (res < 0) {
 464         return res + b;
 465     } else {
 466         return res;
 467     }
 468 }
 469
 470 /**
 471  * Convert char to corresponding direction.
 472  * Used for cubemap options.
 473  */
 474 static int get_direction(char c)
 475 {
 476     switch (c) {
 477     case 'r':
 478         return RIGHT;
 479     case 'l':
 480         return LEFT;
 481     case 'u':
 482         return UP;
 483     case 'd':
 484         return DOWN;
 485     case 'f':
 486         return FRONT;
 487     case 'b':
 488         return BACK;
 489     default:
 490         return -1;
 491     }
 492 }
 493
 494 /**
 495  * Convert char to corresponding rotation angle.
 496  * Used for cubemap options.
 497  */
 498 static int get_rotation(char c)
 499 {
 500     switch (c) {
 501     case '0':
 502         return ROT_0;
 503     case '1':
 504         return ROT_90;
 505     case '2':
 506         return ROT_180;
 507     case '3':
 508         return ROT_270;
 509     default:
 510         return -1;
 511     }
 512 }
 513
 514 /**
 515  * Convert char to corresponding rotation order.
 516  */
 517 static int get_rorder(char c)
 518 {
 519     switch (c) {
 520     case 'Y':
 521     case 'y':
 522         return YAW;
 523     case 'P':
 524     case 'p':
 525         return PITCH;
 526     case 'R':
 527     case 'r':
 528         return ROLL;
 529     default:
 530         return -1;
 531     }
 532 }
 533
 534 /**
 535  * Prepare data for processing cubemap input format.
 536  *
 537  * @param ctx filter context
 538  *
 539  * @return error code
 540  */
 541 static int prepare_cube_in(AVFilterContext *ctx)
 542 {
 543     V360Context *s = ctx->priv;
 544
 545     for (int face = 0; face < NB_FACES; face++) {
 546         const char c = s->in_forder[face];
 547         int direction;
 548
 549         if (c == '\0') {
 550             av_log(ctx, AV_LOG_ERROR,
 551                    "Incomplete in_forder option. Direction for all 6 faces should be specified.\n");
 552             return AVERROR(EINVAL);
 553         }
 554
 555         direction = get_direction(c);
 556         if (direction == -1) {
 557             av_log(ctx, AV_LOG_ERROR,
 558                    "Incorrect direction symbol '%c' in in_forder option.\n", c);
 559             return AVERROR(EINVAL);
 560         }
 561
 562         s->in_cubemap_face_order[direction] = face;
 563     }
 564
 565     for (int face = 0; face < NB_FACES; face++) {
 566         const char c = s->in_frot[face];
 567         int rotation;
 568
 569         if (c == '\0') {
 570             av_log(ctx, AV_LOG_ERROR,
 571                    "Incomplete in_frot option. Rotation for all 6 faces should be specified.\n");
 572             return AVERROR(EINVAL);
 573         }
 574
 575         rotation = get_rotation(c);
 576         if (rotation == -1) {
 577             av_log(ctx, AV_LOG_ERROR,
 578                    "Incorrect rotation symbol '%c' in in_frot option.\n", c);
 579             return AVERROR(EINVAL);
 580         }
 581
 582         s->in_cubemap_face_rotation[face] = rotation;
 583     }
 584
 585     return 0;
 586 }
 587
 588 /**
 589  * Prepare data for processing cubemap output format.
 590  *
 591  * @param ctx filter context
 592  *
 593  * @return error code
 594  */
 595 static int prepare_cube_out(AVFilterContext *ctx)
 596 {
 597     V360Context *s = ctx->priv;
 598
 599     for (int face = 0; face < NB_FACES; face++) {
 600         const char c = s->out_forder[face];
 601         int direction;
 602
 603         if (c == '\0') {
 604             av_log(ctx, AV_LOG_ERROR,
 605                    "Incomplete out_forder option. Direction for all 6 faces should be specified.\n");
 606             return AVERROR(EINVAL);
 607         }
 608
 609         direction = get_direction(c);
 610         if (direction == -1) {
 611             av_log(ctx, AV_LOG_ERROR,
 612                    "Incorrect direction symbol '%c' in out_forder option.\n", c);
 613             return AVERROR(EINVAL);
 614         }
 615
 616         s->out_cubemap_direction_order[face] = direction;
 617     }
 618
 619     for (int face = 0; face < NB_FACES; face++) {
 620         const char c = s->out_frot[face];
 621         int rotation;
 622
 623         if (c == '\0') {
 624             av_log(ctx, AV_LOG_ERROR,
 625                    "Incomplete out_frot option. Rotation for all 6 faces should be specified.\n");
 626             return AVERROR(EINVAL);
 627         }
 628
 629         rotation = get_rotation(c);
 630         if (rotation == -1) {
 631             av_log(ctx, AV_LOG_ERROR,
 632                    "Incorrect rotation symbol '%c' in out_frot option.\n", c);
 633             return AVERROR(EINVAL);
 634         }
 635
 636         s->out_cubemap_face_rotation[face] = rotation;
 637     }
 638
 639     return 0;
 640 }
 641
 642 static inline void rotate_cube_face(float *uf, float *vf, int rotation)
 643 {
 644     float tmp;
 645
 646     switch (rotation) {
 647     case ROT_0:
 648         break;
 649     case ROT_90:
 650         tmp =  *uf;
 651         *uf = -*vf;
 652         *vf =  tmp;
 653         break;
 654     case ROT_180:
 655         *uf = -*uf;
 656         *vf = -*vf;
 657         break;
 658     case ROT_270:
 659         tmp = -*uf;
 660         *uf =  *vf;
 661         *vf =  tmp;
 662         break;
 663     default:
 664         av_assert0(0);
 665     }
 666 }
 667
 668 static inline void rotate_cube_face_inverse(float *uf, float *vf, int rotation)
 669 {
 670     float tmp;
 671
 672     switch (rotation) {
 673     case ROT_0:
 674         break;
 675     case ROT_90:
 676         tmp = -*uf;
 677         *uf =  *vf;
 678         *vf =  tmp;
 679         break;
 680     case ROT_180:
 681         *uf = -*uf;
 682         *vf = -*vf;
 683         break;
 684     case ROT_270:
 685         tmp =  *uf;
 686         *uf = -*vf;
 687         *vf =  tmp;
 688         break;
 689     default:
 690         av_assert0(0);
 691     }
 692 }
 693
 694 /**
 695  * Normalize vector.
 696  *
 697  * @param vec vector
 698  */
 699 static void normalize_vector(float *vec)
 700 {
 701     const float norm = sqrtf(vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2]);
 702
 703     vec[0] /= norm;
 704     vec[1] /= norm;
 705     vec[2] /= norm;
 706 }
 707
 708 /**
 709  * Calculate 3D coordinates on sphere for corresponding cubemap position.
 710  * Common operation for every cubemap.
 711  *
 712  * @param s filter context
 713  * @param uf horizontal cubemap coordinate [0, 1)
 714  * @param vf vertical cubemap coordinate [0, 1)
 715  * @param face face of cubemap
 716  * @param vec coordinates on sphere
 717  */
 718 static void cube_to_xyz(const V360Context *s,
 719                         float uf, float vf, int face,
 720                         float *vec)
 721 {
 722     const int direction = s->out_cubemap_direction_order[face];
 723     float l_x, l_y, l_z;
 724
 725     uf /= (1.f - s->out_pad);
 726     vf /= (1.f - s->out_pad);
 727
 728     rotate_cube_face_inverse(&uf, &vf, s->out_cubemap_face_rotation[face]);
 729
 730     switch (direction) {
 731     case RIGHT:
 732         l_x =  1.f;
 733         l_y = -vf;
 734         l_z =  uf;
 735         break;
 736     case LEFT:
 737         l_x = -1.f;
 738         l_y = -vf;
 739         l_z = -uf;
 740         break;
 741     case UP:
 742         l_x =  uf;
 743         l_y =  1.f;
 744         l_z = -vf;
 745         break;
 746     case DOWN:
 747         l_x =  uf;
 748         l_y = -1.f;
 749         l_z =  vf;
 750         break;
 751     case FRONT:
 752         l_x =  uf;
 753         l_y = -vf;
 754         l_z = -1.f;
 755         break;
 756     case BACK:
 757         l_x = -uf;
 758         l_y = -vf;
 759         l_z =  1.f;
 760         break;
 761     }
 762
 763     vec[0] = l_x;
 764     vec[1] = l_y;
 765     vec[2] = l_z;
 766
 767     normalize_vector(vec);
 768 }
 769
 770 /**
 771  * Calculate cubemap position for corresponding 3D coordinates on sphere.
 772  * Common operation for every cubemap.
 773  *
 774  * @param s filter context
 775  * @param vec coordinated on sphere
 776  * @param uf horizontal cubemap coordinate [0, 1)
 777  * @param vf vertical cubemap coordinate [0, 1)
 778  * @param direction direction of view
 779  */
 780 static void xyz_to_cube(const V360Context *s,
 781                         const float *vec,
 782                         float *uf, float *vf, int *direction)
 783 {
 784     const float phi   = atan2f(vec[0], -vec[2]);
 785     const float theta = asinf(-vec[1]);
 786     float phi_norm, theta_threshold;
 787     int face;
 788
 789     if (phi >= -M_PI_4 && phi < M_PI_4) {
 790         *direction = FRONT;
 791         phi_norm = phi;
 792     } else if (phi >= -(M_PI_2 + M_PI_4) && phi < -M_PI_4) {
 793         *direction = LEFT;
 794         phi_norm = phi + M_PI_2;
 795     } else if (phi >= M_PI_4 && phi < M_PI_2 + M_PI_4) {
 796         *direction = RIGHT;
 797         phi_norm = phi - M_PI_2;
 798     } else {
 799         *direction = BACK;
 800         phi_norm = phi + ((phi > 0.f) ? -M_PI : M_PI);
 801     }
 802
 803     theta_threshold = atanf(cosf(phi_norm));
 804     if (theta > theta_threshold) {
 805         *direction = DOWN;
 806     } else if (theta < -theta_threshold) {
 807         *direction = UP;
 808     }
 809
 810     switch (*direction) {
 811     case RIGHT:
 812         *uf =  vec[2] / vec[0];
 813         *vf = -vec[1] / vec[0];
 814         break;
 815     case LEFT:
 816         *uf =  vec[2] / vec[0];
 817         *vf =  vec[1] / vec[0];
 818         break;
 819     case UP:
 820         *uf =  vec[0] / vec[1];
 821         *vf = -vec[2] / vec[1];
 822         break;
 823     case DOWN:
 824         *uf = -vec[0] / vec[1];
 825         *vf = -vec[2] / vec[1];
 826         break;
 827     case FRONT:
 828         *uf = -vec[0] / vec[2];
 829         *vf =  vec[1] / vec[2];
 830         break;
 831     case BACK:
 832         *uf = -vec[0] / vec[2];
 833         *vf = -vec[1] / vec[2];
 834         break;
 835     default:
 836         av_assert0(0);
 837     }
 838
 839     face = s->in_cubemap_face_order[*direction];
 840     rotate_cube_face(uf, vf, s->in_cubemap_face_rotation[face]);
 841
 842     (*uf) *= s->input_mirror_modifier[0];
 843     (*vf) *= s->input_mirror_modifier[1];
 844 }
 845
 846 /**
 847  * Find position on another cube face in case of overflow/underflow.
 848  * Used for calculation of interpolation window.
 849  *
 850  * @param s filter context
 851  * @param uf horizontal cubemap coordinate
 852  * @param vf vertical cubemap coordinate
 853  * @param direction direction of view
 854  * @param new_uf new horizontal cubemap coordinate
 855  * @param new_vf new vertical cubemap coordinate
 856  * @param face face position on cubemap
 857  */
 858 static void process_cube_coordinates(const V360Context *s,
 859                                      float uf, float vf, int direction,
 860                                      float *new_uf, float *new_vf, int *face)
 861 {
 862     /*
 863      *  Cubemap orientation
 864      *
 865      *           width
 866      *         <------->
 867      *         +-------+
 868      *         |       |                              U
 869      *         | up    |                   h       ------->
 870      * +-------+-------+-------+-------+ ^ e      |
 871      * |       |       |       |       | | i    V |
 872      * | left  | front | right | back  | | g      |
 873      * +-------+-------+-------+-------+ v h      v
 874      *         |       |                   t
 875      *         | down  |
 876      *         +-------+
 877      */
 878
 879     *face = s->in_cubemap_face_order[direction];
 880     rotate_cube_face_inverse(&uf, &vf, s->in_cubemap_face_rotation[*face]);
 881
 882     if ((uf < -1.f || uf >= 1.f) && (vf < -1.f || vf >= 1.f)) {
 883         // There are no pixels to use in this case
 884         *new_uf = uf;
 885         *new_vf = vf;
 886     } else if (uf < -1.f) {
 887         uf += 2.f;
 888         switch (direction) {
 889         case RIGHT:
 890             direction = FRONT;
 891             *new_uf =  uf;
 892             *new_vf =  vf;
 893             break;
 894         case LEFT:
 895             direction = BACK;
 896             *new_uf =  uf;
 897             *new_vf =  vf;
 898             break;
 899         case UP:
 900             direction = LEFT;
 901             *new_uf =  vf;
 902             *new_vf = -uf;
 903             break;
 904         case DOWN:
 905             direction = LEFT;
 906             *new_uf = -vf;
 907             *new_vf =  uf;
 908             break;
 909         case FRONT:
 910             direction = LEFT;
 911             *new_uf =  uf;
 912             *new_vf =  vf;
 913             break;
 914         case BACK:
 915             direction = RIGHT;
 916             *new_uf =  uf;
 917             *new_vf =  vf;
 918             break;
 919         default:
 920             av_assert0(0);
 921         }
 922     } else if (uf >= 1.f) {
 923         uf -= 2.f;
 924         switch (direction) {
 925         case RIGHT:
 926             direction = BACK;
 927             *new_uf =  uf;
 928             *new_vf =  vf;
 929             break;
 930         case LEFT:
 931             direction = FRONT;
 932             *new_uf =  uf;
 933             *new_vf =  vf;
 934             break;
 935         case UP:
 936             direction = RIGHT;
 937             *new_uf = -vf;
 938             *new_vf =  uf;
 939             break;
 940         case DOWN:
 941             direction = RIGHT;
 942             *new_uf =  vf;
 943             *new_vf = -uf;
 944             break;
 945         case FRONT:
 946             direction = RIGHT;
 947             *new_uf =  uf;
 948             *new_vf =  vf;
 949             break;
 950         case BACK:
 951             direction = LEFT;
 952             *new_uf =  uf;
 953             *new_vf =  vf;
 954             break;
 955         default:
 956             av_assert0(0);
 957         }
 958     } else if (vf < -1.f) {
 959         vf += 2.f;
 960         switch (direction) {
 961         case RIGHT:
 962             direction = UP;
 963             *new_uf =  vf;
 964             *new_vf = -uf;
 965             break;
 966         case LEFT:
 967             direction = UP;
 968             *new_uf = -vf;
 969             *new_vf =  uf;
 970             break;
 971         case UP:
 972             direction = BACK;
 973             *new_uf = -uf;
 974             *new_vf = -vf;
 975             break;
 976         case DOWN:
 977             direction = FRONT;
 978             *new_uf =  uf;
 979             *new_vf =  vf;
 980             break;
 981         case FRONT:
 982             direction = UP;
 983             *new_uf =  uf;
 984             *new_vf =  vf;
 985             break;
 986         case BACK:
 987             direction = UP;
 988             *new_uf = -uf;
 989             *new_vf = -vf;
 990             break;
 991         default:
 992             av_assert0(0);
 993         }
 994     } else if (vf >= 1.f) {
 995         vf -= 2.f;
 996         switch (direction) {
 997         case RIGHT:
 998             direction = DOWN;
 999             *new_uf = -vf;
1000             *new_vf =  uf;
1001             break;
1002         case LEFT:
1003             direction = DOWN;
1004             *new_uf =  vf;
1005             *new_vf = -uf;
1006             break;
1007         case UP:
1008             direction = FRONT;
1009             *new_uf =  uf;
1010             *new_vf =  vf;
1011             break;
1012         case DOWN:
1013             direction = BACK;
1014             *new_uf = -uf;
1015             *new_vf = -vf;
1016             break;
1017         case FRONT:
1018             direction = DOWN;
1019             *new_uf =  uf;
1020             *new_vf =  vf;
1021             break;
1022         case BACK:
1023             direction = DOWN;
1024             *new_uf = -uf;
1025             *new_vf = -vf;
1026             break;
1027         default:
1028             av_assert0(0);
1029         }
1030     } else {
1031         // Inside cube face
1032         *new_uf = uf;
1033         *new_vf = vf;
1034     }
1035
1036     *face = s->in_cubemap_face_order[direction];
1037     rotate_cube_face(new_uf, new_vf, s->in_cubemap_face_rotation[*face]);
1038 }
1039
1040 /**
1041  * Calculate 3D coordinates on sphere for corresponding frame position in cubemap3x2 format.
1042  *
1043  * @param s filter context
1044  * @param i horizontal position on frame [0, height)
1045  * @param j vertical position on frame [0, width)
1046  * @param width frame width
1047  * @param height frame height
1048  * @param vec coordinates on sphere
1049  */
1050 static void cube3x2_to_xyz(const V360Context *s,
1051                            int i, int j, int width, int height,
1052                            float *vec)
1053 {
1054     const float ew = width  / 3.f;
1055     const float eh = height / 2.f;
1056
1057     const int u_face = floorf(i / ew);
1058     const int v_face = floorf(j / eh);
1059     const int face = u_face + 3 * v_face;
1060
1061     const int u_shift = ceilf(ew * u_face);
1062     const int v_shift = ceilf(eh * v_face);
1063     const int ewi = ceilf(ew * (u_face + 1)) - u_shift;
1064     const int ehi = ceilf(eh * (v_face + 1)) - v_shift;
1065
1066     const float uf = 2.f * (i - u_shift) / ewi - 1.f;
1067     const float vf = 2.f * (j - v_shift) / ehi - 1.f;
1068
1069     cube_to_xyz(s, uf, vf, face, vec);
1070 }
1071
1072 /**
1073  * Calculate frame position in cubemap3x2 format for corresponding 3D coordinates on sphere.
1074  *
1075  * @param s filter context
1076  * @param vec coordinates on sphere
1077  * @param width frame width
1078  * @param height frame height
1079  * @param us horizontal coordinates for interpolation window
1080  * @param vs vertical coordinates for interpolation window
1081  * @param du horizontal relative coordinate
1082  * @param dv vertical relative coordinate
1083  */
1084 static void xyz_to_cube3x2(const V360Context *s,
1085                            const float *vec, int width, int height,
1086                            uint16_t us[4][4], uint16_t vs[4][4], float *du, float *dv)
1087 {
1088     const float ew = width  / 3.f;
1089     const float eh = height / 2.f;
1090     float uf, vf;
1091     int ui, vi;
1092     int ewi, ehi;
1093     int i, j;
1094     int direction, face;
1095     int u_face, v_face;
1096
1097     xyz_to_cube(s, vec, &uf, &vf, &direction);
1098
1099     uf *= (1.f - s->in_pad);
1100     vf *= (1.f - s->in_pad);
1101
1102     face = s->in_cubemap_face_order[direction];
1103     u_face = face % 3;
1104     v_face = face / 3;
1105     ewi = ceilf(ew * (u_face + 1)) - ceilf(ew * u_face);
1106     ehi = ceilf(eh * (v_face + 1)) - ceilf(eh * v_face);
1107
1108     uf = 0.5f * ewi * (uf + 1.f);
1109     vf = 0.5f * ehi * (vf + 1.f);
1110
1111     ui = floorf(uf);
1112     vi = floorf(vf);
1113
1114     *du = uf - ui;
1115     *dv = vf - vi;
1116
1117     for (i = -1; i < 3; i++) {
1118         for (j = -1; j < 3; j++) {
1119             int new_ui = ui + j;
1120             int new_vi = vi + i;
1121             int u_shift, v_shift;
1122             int new_ewi, new_ehi;
1123
1124             if (new_ui >= 0 && new_ui < ewi && new_vi >= 0 && new_vi < ehi) {
1125                 face = s->in_cubemap_face_order[direction];
1126
1127                 u_face = face % 3;
1128                 v_face = face / 3;
1129                 u_shift = ceilf(ew * u_face);
1130                 v_shift = ceilf(eh * v_face);
1131             } else {
1132                 uf = 2.f * new_ui / ewi - 1.f;
1133                 vf = 2.f * new_vi / ehi - 1.f;
1134
1135                 uf /= (1.f - s->in_pad);
1136                 vf /= (1.f - s->in_pad);
1137
1138                 process_cube_coordinates(s, uf, vf, direction, &uf, &vf, &face);
1139
1140                 uf *= (1.f - s->in_pad);
1141                 vf *= (1.f - s->in_pad);
1142
1143                 u_face = face % 3;
1144                 v_face = face / 3;
1145                 u_shift = ceilf(ew * u_face);
1146                 v_shift = ceilf(eh * v_face);
1147                 new_ewi = ceilf(ew * (u_face + 1)) - u_shift;
1148                 new_ehi = ceilf(eh * (v_face + 1)) - v_shift;
1149
1150                 new_ui = av_clip(roundf(0.5f * new_ewi * (uf + 1.f)), 0, new_ewi - 1);
1151                 new_vi = av_clip(roundf(0.5f * new_ehi * (vf + 1.f)), 0, new_ehi - 1);
1152             }
1153
1154             us[i + 1][j + 1] = u_shift + new_ui;
1155             vs[i + 1][j + 1] = v_shift + new_vi;
1156         }
1157     }
1158 }
1159
1160 /**
1161  * Calculate 3D coordinates on sphere for corresponding frame position in cubemap1x6 format.
1162  *
1163  * @param s filter context
1164  * @param i horizontal position on frame [0, height)
1165  * @param j vertical position on frame [0, width)
1166  * @param width frame width
1167  * @param height frame height
1168  * @param vec coordinates on sphere
1169  */
1170 static void cube1x6_to_xyz(const V360Context *s,
1171                            int i, int j, int width, int height,
1172                            float *vec)
1173 {
1174     const float ew = width;
1175     const float eh = height / 6.f;
1176
1177     const int face = floorf(j / eh);
1178
1179     const int v_shift = ceilf(eh * face);
1180     const int ehi = ceilf(eh * (face + 1)) - v_shift;
1181
1182     const float uf = 2.f *  i            / ew  - 1.f;
1183     const float vf = 2.f * (j - v_shift) / ehi - 1.f;
1184
1185     cube_to_xyz(s, uf, vf, face, vec);
1186 }
1187
1188 /**
1189  * Calculate 3D coordinates on sphere for corresponding frame position in cubemap6x1 format.
1190  *
1191  * @param s filter context
1192  * @param i horizontal position on frame [0, height)
1193  * @param j vertical position on frame [0, width)
1194  * @param width frame width
1195  * @param height frame height
1196  * @param vec coordinates on sphere
1197  */
1198 static void cube6x1_to_xyz(const V360Context *s,
1199                            int i, int j, int width, int height,
1200                            float *vec)
1201 {
1202     const float ew = width / 6.f;
1203     const float eh = height;
1204
1205     const int face = floorf(i / ew);
1206
1207     const int u_shift = ceilf(ew * face);
1208     const int ewi = ceilf(ew * (face + 1)) - u_shift;
1209
1210     const float uf = 2.f * (i - u_shift) / ewi - 1.f;
1211     const float vf = 2.f *  j            / eh  - 1.f;
1212
1213     cube_to_xyz(s, uf, vf, face, vec);
1214 }
1215
1216 /**
1217  * Calculate frame position in cubemap1x6 format for corresponding 3D coordinates on sphere.
1218  *
1219  * @param s filter context
1220  * @param vec coordinates on sphere
1221  * @param width frame width
1222  * @param height frame height
1223  * @param us horizontal coordinates for interpolation window
1224  * @param vs vertical coordinates for interpolation window
1225  * @param du horizontal relative coordinate
1226  * @param dv vertical relative coordinate
1227  */
1228 static void xyz_to_cube1x6(const V360Context *s,
1229                            const float *vec, int width, int height,
1230                            uint16_t us[4][4], uint16_t vs[4][4], float *du, float *dv)
1231 {
1232     const float eh = height / 6.f;
1233     const int ewi = width;
1234     float uf, vf;
1235     int ui, vi;
1236     int ehi;
1237     int i, j;
1238     int direction, face;
1239
1240     xyz_to_cube(s, vec, &uf, &vf, &direction);
1241
1242     uf *= (1.f - s->in_pad);
1243     vf *= (1.f - s->in_pad);
1244
1245     face = s->in_cubemap_face_order[direction];
1246     ehi = ceilf(eh * (face + 1)) - ceilf(eh * face);
1247
1248     uf = 0.5f * ewi * (uf + 1.f);
1249     vf = 0.5f * ehi * (vf + 1.f);
1250
1251     ui = floorf(uf);
1252     vi = floorf(vf);
1253
1254     *du = uf - ui;
1255     *dv = vf - vi;
1256
1257     for (i = -1; i < 3; i++) {
1258         for (j = -1; j < 3; j++) {
1259             int new_ui = ui + j;
1260             int new_vi = vi + i;
1261             int v_shift;
1262             int new_ehi;
1263
1264             if (new_ui >= 0 && new_ui < ewi && new_vi >= 0 && new_vi < ehi) {
1265                 face = s->in_cubemap_face_order[direction];
1266
1267                 v_shift = ceilf(eh * face);
1268             } else {
1269                 uf = 2.f * new_ui / ewi - 1.f;
1270                 vf = 2.f * new_vi / ehi - 1.f;
1271
1272                 uf /= (1.f - s->in_pad);
1273                 vf /= (1.f - s->in_pad);
1274
1275                 process_cube_coordinates(s, uf, vf, direction, &uf, &vf, &face);
1276
1277                 uf *= (1.f - s->in_pad);
1278                 vf *= (1.f - s->in_pad);
1279
1280                 v_shift = ceilf(eh * face);
1281                 new_ehi = ceilf(eh * (face + 1)) - v_shift;
1282
1283                 new_ui = av_clip(roundf(0.5f *     ewi * (uf + 1.f)), 0,     ewi - 1);
1284                 new_vi = av_clip(roundf(0.5f * new_ehi * (vf + 1.f)), 0, new_ehi - 1);
1285             }
1286
1287             us[i + 1][j + 1] =           new_ui;
1288             vs[i + 1][j + 1] = v_shift + new_vi;
1289         }
1290     }
1291 }
1292
1293 /**
1294  * Calculate frame position in cubemap6x1 format for corresponding 3D coordinates on sphere.
1295  *
1296  * @param s filter context
1297  * @param vec coordinates on sphere
1298  * @param width frame width
1299  * @param height frame height
1300  * @param us horizontal coordinates for interpolation window
1301  * @param vs vertical coordinates for interpolation window
1302  * @param du horizontal relative coordinate
1303  * @param dv vertical relative coordinate
1304  */
1305 static void xyz_to_cube6x1(const V360Context *s,
1306                            const float *vec, int width, int height,
1307                            uint16_t us[4][4], uint16_t vs[4][4], float *du, float *dv)
1308 {
1309     const float ew = width / 6.f;
1310     const int ehi = height;
1311     float uf, vf;
1312     int ui, vi;
1313     int ewi;
1314     int i, j;
1315     int direction, face;
1316
1317     xyz_to_cube(s, vec, &uf, &vf, &direction);
1318
1319     uf *= (1.f - s->in_pad);
1320     vf *= (1.f - s->in_pad);
1321
1322     face = s->in_cubemap_face_order[direction];
1323     ewi = ceilf(ew * (face + 1)) - ceilf(ew * face);
1324
1325     uf = 0.5f * ewi * (uf + 1.f);
1326     vf = 0.5f * ehi * (vf + 1.f);
1327
1328     ui = floorf(uf);
1329     vi = floorf(vf);
1330
1331     *du = uf - ui;
1332     *dv = vf - vi;
1333
1334     for (i = -1; i < 3; i++) {
1335         for (j = -1; j < 3; j++) {
1336             int new_ui = ui + j;
1337             int new_vi = vi + i;
1338             int u_shift;
1339             int new_ewi;
1340
1341             if (new_ui >= 0 && new_ui < ewi && new_vi >= 0 && new_vi < ehi) {
1342                 face = s->in_cubemap_face_order[direction];
1343
1344                 u_shift = ceilf(ew * face);
1345             } else {
1346                 uf = 2.f * new_ui / ewi - 1.f;
1347                 vf = 2.f * new_vi / ehi - 1.f;
1348
1349                 uf /= (1.f - s->in_pad);
1350                 vf /= (1.f - s->in_pad);
1351
1352                 process_cube_coordinates(s, uf, vf, direction, &uf, &vf, &face);
1353
1354                 uf *= (1.f - s->in_pad);
1355                 vf *= (1.f - s->in_pad);
1356
1357                 u_shift = ceilf(ew * face);
1358                 new_ewi = ceilf(ew * (face + 1)) - u_shift;
1359
1360                 new_ui = av_clip(roundf(0.5f * new_ewi * (uf + 1.f)), 0, new_ewi - 1);
1361                 new_vi = av_clip(roundf(0.5f *     ehi * (vf + 1.f)), 0,     ehi - 1);
1362             }
1363
1364             us[i + 1][j + 1] = u_shift + new_ui;
1365             vs[i + 1][j + 1] =           new_vi;
1366         }
1367     }
1368 }
1369
1370 /**
1371  * Calculate 3D coordinates on sphere for corresponding frame position in equirectangular format.
1372  *
1373  * @param s filter context
1374  * @param i horizontal position on frame [0, height)
1375  * @param j vertical position on frame [0, width)
1376  * @param width frame width
1377  * @param height frame height
1378  * @param vec coordinates on sphere
1379  */
1380 static void equirect_to_xyz(const V360Context *s,
1381                             int i, int j, int width, int height,
1382                             float *vec)
1383 {
1384     const float phi   = ((2.f * i) / width  - 1.f) * M_PI;
1385     const float theta = ((2.f * j) / height - 1.f) * M_PI_2;
1386
1387     const float sin_phi   = sinf(phi);
1388     const float cos_phi   = cosf(phi);
1389     const float sin_theta = sinf(theta);
1390     const float cos_theta = cosf(theta);
1391
1392     vec[0] =  cos_theta * sin_phi;
1393     vec[1] = -sin_theta;
1394     vec[2] = -cos_theta * cos_phi;
1395 }
1396
1397 /**
1398  * Calculate 3D coordinates on sphere for corresponding frame position in stereographic format.
1399  *
1400  * @param s filter context
1401  * @param i horizontal position on frame [0, height)
1402  * @param j vertical position on frame [0, width)
1403  * @param width frame width
1404  * @param height frame height
1405  * @param vec coordinates on sphere
1406  */
1407 static void stereographic_to_xyz(const V360Context *s,
1408                                  int i, int j, int width, int height,
1409                                  float *vec)
1410 {
1411     const float x = ((2.f * i) / width  - 1.f) * (s->h_fov / 180.f) * M_PI;
1412     const float y = ((2.f * j) / height - 1.f) * (s->v_fov /  90.f) * M_PI_2;
1413     const float xy = x * x + y * y;
1414
1415     vec[0] = 2.f * x / (1.f + xy);
1416     vec[1] = (-1.f + xy) / (1.f + xy);
1417     vec[2] = 2.f * y / (1.f + xy);
1418
1419     normalize_vector(vec);
1420 }
1421
1422 /**
1423  * Calculate frame position in equirectangular format for corresponding 3D coordinates on sphere.
1424  *
1425  * @param s filter context
1426  * @param vec coordinates on sphere
1427  * @param width frame width
1428  * @param height frame height
1429  * @param us horizontal coordinates for interpolation window
1430  * @param vs vertical coordinates for interpolation window
1431  * @param du horizontal relative coordinate
1432  * @param dv vertical relative coordinate
1433  */
1434 static void xyz_to_equirect(const V360Context *s,
1435                             const float *vec, int width, int height,
1436                             uint16_t us[4][4], uint16_t vs[4][4], float *du, float *dv)
1437 {
1438     const float phi   = atan2f(vec[0], -vec[2]) * s->input_mirror_modifier[0];
1439     const float theta = asinf(-vec[1]) * s->input_mirror_modifier[1];
1440     float uf, vf;
1441     int ui, vi;
1442     int i, j;
1443
1444     uf = (phi   / M_PI   + 1.f) * width  / 2.f;
1445     vf = (theta / M_PI_2 + 1.f) * height / 2.f;
1446     ui = floorf(uf);
1447     vi = floorf(vf);
1448
1449     *du = uf - ui;
1450     *dv = vf - vi;
1451
1452     for (i = -1; i < 3; i++) {
1453         for (j = -1; j < 3; j++) {
1454             us[i + 1][j + 1] = mod(ui + j, width);
1455             vs[i + 1][j + 1] = av_clip(vi + i, 0, height - 1);
1456         }
1457     }
1458 }
1459
1460 /**
1461  * Prepare data for processing equi-angular cubemap input format.
1462  *
1463  * @param ctx filter context
1464
1465  * @return error code
1466  */
1467 static int prepare_eac_in(AVFilterContext *ctx)
1468 {
1469     V360Context *s = ctx->priv;
1470
1471     if (s->ih_flip && s->iv_flip) {
1472         s->in_cubemap_face_order[RIGHT] = BOTTOM_LEFT;
1473         s->in_cubemap_face_order[LEFT]  = BOTTOM_RIGHT;
1474         s->in_cubemap_face_order[UP]    = TOP_LEFT;
1475         s->in_cubemap_face_order[DOWN]  = TOP_RIGHT;
1476         s->in_cubemap_face_order[FRONT] = BOTTOM_MIDDLE;
1477         s->in_cubemap_face_order[BACK]  = TOP_MIDDLE;
1478     } else if (s->ih_flip) {
1479         s->in_cubemap_face_order[RIGHT] = TOP_LEFT;
1480         s->in_cubemap_face_order[LEFT]  = TOP_RIGHT;
1481         s->in_cubemap_face_order[UP]    = BOTTOM_LEFT;
1482         s->in_cubemap_face_order[DOWN]  = BOTTOM_RIGHT;
1483         s->in_cubemap_face_order[FRONT] = TOP_MIDDLE;
1484         s->in_cubemap_face_order[BACK]  = BOTTOM_MIDDLE;
1485     } else if (s->iv_flip) {
1486         s->in_cubemap_face_order[RIGHT] = BOTTOM_RIGHT;
1487         s->in_cubemap_face_order[LEFT]  = BOTTOM_LEFT;
1488         s->in_cubemap_face_order[UP]    = TOP_RIGHT;
1489         s->in_cubemap_face_order[DOWN]  = TOP_LEFT;
1490         s->in_cubemap_face_order[FRONT] = BOTTOM_MIDDLE;
1491         s->in_cubemap_face_order[BACK]  = TOP_MIDDLE;
1492     } else {
1493         s->in_cubemap_face_order[RIGHT] = TOP_RIGHT;
1494         s->in_cubemap_face_order[LEFT]  = TOP_LEFT;
1495         s->in_cubemap_face_order[UP]    = BOTTOM_RIGHT;
1496         s->in_cubemap_face_order[DOWN]  = BOTTOM_LEFT;
1497         s->in_cubemap_face_order[FRONT] = TOP_MIDDLE;
1498         s->in_cubemap_face_order[BACK]  = BOTTOM_MIDDLE;
1499     }
1500
1501     if (s->iv_flip) {
1502         s->in_cubemap_face_rotation[TOP_LEFT]      = ROT_270;
1503         s->in_cubemap_face_rotation[TOP_MIDDLE]    = ROT_90;
1504         s->in_cubemap_face_rotation[TOP_RIGHT]     = ROT_270;
1505         s->in_cubemap_face_rotation[BOTTOM_LEFT]   = ROT_0;
1506         s->in_cubemap_face_rotation[BOTTOM_MIDDLE] = ROT_0;
1507         s->in_cubemap_face_rotation[BOTTOM_RIGHT]  = ROT_0;
1508     } else {
1509         s->in_cubemap_face_rotation[TOP_LEFT]      = ROT_0;
1510         s->in_cubemap_face_rotation[TOP_MIDDLE]    = ROT_0;
1511         s->in_cubemap_face_rotation[TOP_RIGHT]     = ROT_0;
1512         s->in_cubemap_face_rotation[BOTTOM_LEFT]   = ROT_270;
1513         s->in_cubemap_face_rotation[BOTTOM_MIDDLE] = ROT_90;
1514         s->in_cubemap_face_rotation[BOTTOM_RIGHT]  = ROT_270;
1515     }
1516
1517     return 0;
1518 }
1519
1520 /**
1521  * Prepare data for processing equi-angular cubemap output format.
1522  *
1523  * @param ctx filter context
1524  *
1525  * @return error code
1526  */
1527 static int prepare_eac_out(AVFilterContext *ctx)
1528 {
1529     V360Context *s = ctx->priv;
1530
1531     s->out_cubemap_direction_order[TOP_LEFT]      = LEFT;
1532     s->out_cubemap_direction_order[TOP_MIDDLE]    = FRONT;
1533     s->out_cubemap_direction_order[TOP_RIGHT]     = RIGHT;
1534     s->out_cubemap_direction_order[BOTTOM_LEFT]   = DOWN;
1535     s->out_cubemap_direction_order[BOTTOM_MIDDLE] = BACK;
1536     s->out_cubemap_direction_order[BOTTOM_RIGHT]  = UP;
1537
1538     s->out_cubemap_face_rotation[TOP_LEFT]      = ROT_0;
1539     s->out_cubemap_face_rotation[TOP_MIDDLE]    = ROT_0;
1540     s->out_cubemap_face_rotation[TOP_RIGHT]     = ROT_0;
1541     s->out_cubemap_face_rotation[BOTTOM_LEFT]   = ROT_270;
1542     s->out_cubemap_face_rotation[BOTTOM_MIDDLE] = ROT_90;
1543     s->out_cubemap_face_rotation[BOTTOM_RIGHT]  = ROT_270;
1544
1545     return 0;
1546 }
1547
1548 /**
1549  * Calculate 3D coordinates on sphere for corresponding frame position in equi-angular cubemap format.
1550  *
1551  * @param s filter context
1552  * @param i horizontal position on frame [0, height)
1553  * @param j vertical position on frame [0, width)
1554  * @param width frame width
1555  * @param height frame height
1556  * @param vec coordinates on sphere
1557  */
1558 static void eac_to_xyz(const V360Context *s,
1559                        int i, int j, int width, int height,
1560                        float *vec)
1561 {
1562     const float pixel_pad = 2;
1563     const float u_pad = pixel_pad / width;
1564     const float v_pad = pixel_pad / height;
1565
1566     int u_face, v_face, face;
1567
1568     float l_x, l_y, l_z;
1569
1570     float uf = (float)i / width;
1571     float vf = (float)j / height;
1572
1573     // EAC has 2-pixel padding on faces except between faces on the same row
1574     // Padding pixels seems not to be stretched with tangent as regular pixels
1575     // Formulas below approximate original padding as close as I could get experimentally
1576
1577     // Horizontal padding
1578     uf = 3.f * (uf - u_pad) / (1.f - 2.f * u_pad);
1579     if (uf < 0.f) {
1580         u_face = 0;
1581         uf -= 0.5f;
1582     } else if (uf >= 3.f) {
1583         u_face = 2;
1584         uf -= 2.5f;
1585     } else {
1586         u_face = floorf(uf);
1587         uf = fmodf(uf, 1.f) - 0.5f;
1588     }
1589
1590     // Vertical padding
1591     v_face = floorf(vf * 2.f);
1592     vf = (vf - v_pad - 0.5f * v_face) / (0.5f - 2.f * v_pad) - 0.5f;
1593
1594     if (uf >= -0.5f && uf < 0.5f) {
1595         uf = tanf(M_PI_2 * uf);
1596     } else {
1597         uf = 2.f * uf;
1598     }
1599     if (vf >= -0.5f && vf < 0.5f) {
1600         vf = tanf(M_PI_2 * vf);
1601     } else {
1602         vf = 2.f * vf;
1603     }
1604
1605     face = u_face + 3 * v_face;
1606
1607     switch (face) {
1608     case TOP_LEFT:
1609         l_x = -1.f;
1610         l_y = -vf;
1611         l_z = -uf;
1612         break;
1613     case TOP_MIDDLE:
1614         l_x =  uf;
1615         l_y = -vf;
1616         l_z = -1.f;
1617         break;
1618     case TOP_RIGHT:
1619         l_x =  1.f;
1620         l_y = -vf;
1621         l_z =  uf;
1622         break;
1623     case BOTTOM_LEFT:
1624         l_x = -vf;
1625         l_y = -1.f;
1626         l_z =  uf;
1627         break;
1628     case BOTTOM_MIDDLE:
1629         l_x = -vf;
1630         l_y =  uf;
1631         l_z =  1.f;
1632         break;
1633     case BOTTOM_RIGHT:
1634         l_x = -vf;
1635         l_y =  1.f;
1636         l_z = -uf;
1637         break;
1638     default:
1639         av_assert0(0);
1640     }
1641
1642     vec[0] = l_x;
1643     vec[1] = l_y;
1644     vec[2] = l_z;
1645
1646     normalize_vector(vec);
1647 }
1648
1649 /**
1650  * Calculate frame position in equi-angular cubemap format for corresponding 3D coordinates on sphere.
1651  *
1652  * @param s filter context
1653  * @param vec coordinates on sphere
1654  * @param width frame width
1655  * @param height frame height
1656  * @param us horizontal coordinates for interpolation window
1657  * @param vs vertical coordinates for interpolation window
1658  * @param du horizontal relative coordinate
1659  * @param dv vertical relative coordinate
1660  */
1661 static void xyz_to_eac(const V360Context *s,
1662                        const float *vec, int width, int height,
1663                        uint16_t us[4][4], uint16_t vs[4][4], float *du, float *dv)
1664 {
1665     const float pixel_pad = 2;
1666     const float u_pad = pixel_pad / width;
1667     const float v_pad = pixel_pad / height;
1668
1669     float uf, vf;
1670     int ui, vi;
1671     int i, j;
1672     int direction, face;
1673     int u_face, v_face;
1674
1675     xyz_to_cube(s, vec, &uf, &vf, &direction);
1676
1677     face = s->in_cubemap_face_order[direction];
1678     u_face = face % 3;
1679     v_face = face / 3;
1680
1681     uf = M_2_PI * atanf(uf) + 0.5f;
1682     vf = M_2_PI * atanf(vf) + 0.5f;
1683
1684     // These formulas are inversed from eac_to_xyz ones
1685     uf = (uf + u_face) * (1.f - 2.f * u_pad) / 3.f + u_pad;
1686     vf = vf * (0.5f - 2.f * v_pad) + v_pad + 0.5f * v_face;
1687
1688     uf *= width;
1689     vf *= height;
1690
1691     ui = floorf(uf);
1692     vi = floorf(vf);
1693
1694     *du = uf - ui;
1695     *dv = vf - vi;
1696
1697     for (i = -1; i < 3; i++) {
1698         for (j = -1; j < 3; j++) {
1699             us[i + 1][j + 1] = av_clip(ui + j, 0, width  - 1);
1700             vs[i + 1][j + 1] = av_clip(vi + i, 0, height - 1);
1701         }
1702     }
1703 }
1704
1705 /**
1706  * Prepare data for processing flat output format.
1707  *
1708  * @param ctx filter context
1709  *
1710  * @return error code
1711  */
1712 static int prepare_flat_out(AVFilterContext *ctx)
1713 {
1714     V360Context *s = ctx->priv;
1715
1716     const float h_angle = 0.5f * s->h_fov * M_PI / 180.f;
1717     const float v_angle = 0.5f * s->v_fov * M_PI / 180.f;
1718
1719     const float sin_phi   = sinf(h_angle);
1720     const float cos_phi   = cosf(h_angle);
1721     const float sin_theta = sinf(v_angle);
1722     const float cos_theta = cosf(v_angle);
1723
1724     s->flat_range[0] =  cos_theta * sin_phi;
1725     s->flat_range[1] =  sin_theta;
1726     s->flat_range[2] = -cos_theta * cos_phi;
1727
1728     return 0;
1729 }
1730
1731 /**
1732  * Calculate 3D coordinates on sphere for corresponding frame position in flat format.
1733  *
1734  * @param s filter context
1735  * @param i horizontal position on frame [0, height)
1736  * @param j vertical position on frame [0, width)
1737  * @param width frame width
1738  * @param height frame height
1739  * @param vec coordinates on sphere
1740  */
1741 static void flat_to_xyz(const V360Context *s,
1742                         int i, int j, int width, int height,
1743                         float *vec)
1744 {
1745     const float l_x =  s->flat_range[0] * (2.f * i / width  - 1.f);
1746     const float l_y = -s->flat_range[1] * (2.f * j / height - 1.f);
1747     const float l_z =  s->flat_range[2];
1748
1749     vec[0] = l_x;
1750     vec[1] = l_y;
1751     vec[2] = l_z;
1752
1753     normalize_vector(vec);
1754 }
1755
1756 /**
1757  * Calculate frame position in dual fisheye format for corresponding 3D coordinates on sphere.
1758  *
1759  * @param s filter context
1760  * @param vec coordinates on sphere
1761  * @param width frame width
1762  * @param height frame height
1763  * @param us horizontal coordinates for interpolation window
1764  * @param vs vertical coordinates for interpolation window
1765  * @param du horizontal relative coordinate
1766  * @param dv vertical relative coordinate
1767  */
1768 static void xyz_to_dfisheye(const V360Context *s,
1769                             const float *vec, int width, int height,
1770                             uint16_t us[4][4], uint16_t vs[4][4], float *du, float *dv)
1771 {
1772     const float scale = 1.f - s->in_pad;
1773
1774     const float ew = width / 2.f;
1775     const float eh = height;
1776
1777     const float phi   = atan2f(-vec[1], -vec[0]) * s->input_mirror_modifier[0];
1778     const float theta = acosf(fabsf(vec[2])) / M_PI * s->input_mirror_modifier[1];
1779
1780     float uf = (theta * cosf(phi) * scale + 0.5f) * ew;
1781     float vf = (theta * sinf(phi) * scale + 0.5f) * eh;
1782
1783     int ui, vi;
1784     int u_shift;
1785     int i, j;
1786
1787     if (vec[2] >= 0) {
1788         u_shift = 0;
1789     } else {
1790         u_shift = ceilf(ew);
1791         uf = ew - uf;
1792     }
1793
1794     ui = floorf(uf);
1795     vi = floorf(vf);
1796
1797     *du = uf - ui;
1798     *dv = vf - vi;
1799
1800     for (i = -1; i < 3; i++) {
1801         for (j = -1; j < 3; j++) {
1802             us[i + 1][j + 1] = av_clip(u_shift + ui + j, 0, width  - 1);
1803             vs[i + 1][j + 1] = av_clip(          vi + i, 0, height - 1);
1804         }
1805     }
1806 }
1807
1808 /**
1809  * Calculate 3D coordinates on sphere for corresponding frame position in barrel facebook's format.
1810  *
1811  * @param s filter context
1812  * @param i horizontal position on frame [0, height)
1813  * @param j vertical position on frame [0, width)
1814  * @param width frame width
1815  * @param height frame height
1816  * @param vec coordinates on sphere
1817  */
1818 static void barrel_to_xyz(const V360Context *s,
1819                           int i, int j, int width, int height,
1820                           float *vec)
1821 {
1822     const float scale = 0.99f;
1823     float l_x, l_y, l_z;
1824
1825     if (i < 4 * width / 5) {
1826         const float theta_range = M_PI / 4.f;
1827
1828         const int ew = 4 * width / 5;
1829         const int eh = height;
1830
1831         const float phi   = ((2.f * i) / ew - 1.f) * M_PI        / scale;
1832         const float theta = ((2.f * j) / eh - 1.f) * theta_range / scale;
1833
1834         const float sin_phi   = sinf(phi);
1835         const float cos_phi   = cosf(phi);
1836         const float sin_theta = sinf(theta);
1837         const float cos_theta = cosf(theta);
1838
1839         l_x =  cos_theta * sin_phi;
1840         l_y = -sin_theta;
1841         l_z = -cos_theta * cos_phi;
1842     } else {
1843         const int ew = width  / 5;
1844         const int eh = height / 2;
1845
1846         float uf, vf;
1847
1848         if (j < eh) {   // UP
1849             uf = 2.f * (i - 4 * ew) / ew  - 1.f;
1850             vf = 2.f * (j         ) / eh - 1.f;
1851
1852             uf /= scale;
1853             vf /= scale;
1854
1855             l_x =  uf;
1856             l_y =  1.f;
1857             l_z = -vf;
1858         } else {            // DOWN
1859             uf = 2.f * (i - 4 * ew) / ew - 1.f;
1860             vf = 2.f * (j -     eh) / eh - 1.f;
1861
1862             uf /= scale;
1863             vf /= scale;
1864
1865             l_x =  uf;
1866             l_y = -1.f;
1867             l_z =  vf;
1868         }
1869     }
1870
1871     vec[0] = l_x;
1872     vec[1] = l_y;
1873     vec[2] = l_z;
1874
1875     normalize_vector(vec);
1876 }
1877
1878 /**
1879  * Calculate frame position in barrel facebook's format for corresponding 3D coordinates on sphere.
1880  *
1881  * @param s filter context
1882  * @param vec coordinates on sphere
1883  * @param width frame width
1884  * @param height frame height
1885  * @param us horizontal coordinates for interpolation window
1886  * @param vs vertical coordinates for interpolation window
1887  * @param du horizontal relative coordinate
1888  * @param dv vertical relative coordinate
1889  */
1890 static void xyz_to_barrel(const V360Context *s,
1891                           const float *vec, int width, int height,
1892                           uint16_t us[4][4], uint16_t vs[4][4], float *du, float *dv)
1893 {
1894     const float scale = 0.99f;
1895
1896     const float phi   = atan2f(vec[0], -vec[2]) * s->input_mirror_modifier[0];
1897     const float theta = asinf(-vec[1]) * s->input_mirror_modifier[1];
1898     const float theta_range = M_PI / 4.f;
1899
1900     int ew, eh;
1901     int u_shift, v_shift;
1902     float uf, vf;
1903     int ui, vi;
1904     int i, j;
1905
1906     if (theta > -theta_range && theta < theta_range) {
1907         ew = 4 * width / 5;
1908         eh = height;
1909
1910         u_shift = s->ih_flip ? width / 5 : 0;
1911         v_shift = 0;
1912
1913         uf = (phi   / M_PI        * scale + 1.f) * ew / 2.f;
1914         vf = (theta / theta_range * scale + 1.f) * eh / 2.f;
1915     } else {
1916         ew = width  / 5;
1917         eh = height / 2;
1918
1919         u_shift = s->ih_flip ? 0 : 4 * ew;
1920
1921         if (theta < 0.f) {  // UP
1922             uf =  vec[0] / vec[1];
1923             vf = -vec[2] / vec[1];
1924             v_shift = 0;
1925         } else {            // DOWN
1926             uf = -vec[0] / vec[1];
1927             vf = -vec[2] / vec[1];
1928             v_shift = eh;
1929         }
1930
1931         uf *= s->input_mirror_modifier[0] * s->input_mirror_modifier[1];
1932         vf *= s->input_mirror_modifier[1];
1933
1934         uf = 0.5f * ew * (uf * scale + 1.f);
1935         vf = 0.5f * eh * (vf * scale + 1.f);
1936     }
1937
1938     ui = floorf(uf);
1939     vi = floorf(vf);
1940
1941     *du = uf - ui;
1942     *dv = vf - vi;
1943
1944     for (i = -1; i < 3; i++) {
1945         for (j = -1; j < 3; j++) {
1946             us[i + 1][j + 1] = u_shift + av_clip(ui + j, 0, ew - 1);
1947             vs[i + 1][j + 1] = v_shift + av_clip(vi + i, 0, eh - 1);
1948         }
1949     }
1950 }
1951
1952 static void multiply_matrix(float c[3][3], const float a[3][3], const float b[3][3])
1953 {
1954     for (int i = 0; i < 3; i++) {
1955         for (int j = 0; j < 3; j++) {
1956             float sum = 0;
1957
1958             for (int k = 0; k < 3; k++)
1959                 sum += a[i][k] * b[k][j];
1960
1961             c[i][j] = sum;
1962         }
1963     }
1964 }
1965
1966 /**
1967  * Calculate rotation matrix for yaw/pitch/roll angles.
1968  */
1969 static inline void calculate_rotation_matrix(float yaw, float pitch, float roll,
1970                                              float rot_mat[3][3],
1971                                              const int rotation_order[3])
1972 {
1973     const float yaw_rad   = yaw   * M_PI / 180.f;
1974     const float pitch_rad = pitch * M_PI / 180.f;
1975     const float roll_rad  = roll  * M_PI / 180.f;
1976
1977     const float sin_yaw   = sinf(-yaw_rad);
1978     const float cos_yaw   = cosf(-yaw_rad);
1979     const float sin_pitch = sinf(pitch_rad);
1980     const float cos_pitch = cosf(pitch_rad);
1981     const float sin_roll  = sinf(roll_rad);
1982     const float cos_roll  = cosf(roll_rad);
1983
1984     float m[3][3][3];
1985     float temp[3][3];
1986
1987     m[0][0][0] =  cos_yaw;  m[0][0][1] = 0;          m[0][0][2] =  sin_yaw;
1988     m[0][1][0] =  0;        m[0][1][1] = 1;          m[0][1][2] =  0;
1989     m[0][2][0] = -sin_yaw;  m[0][2][1] = 0;          m[0][2][2] =  cos_yaw;
1990
1991     m[1][0][0] = 1;         m[1][0][1] = 0;          m[1][0][2] =  0;
1992     m[1][1][0] = 0;         m[1][1][1] = cos_pitch;  m[1][1][2] = -sin_pitch;
1993     m[1][2][0] = 0;         m[1][2][1] = sin_pitch;  m[1][2][2] =  cos_pitch;
1994
1995     m[2][0][0] = cos_roll;  m[2][0][1] = -sin_roll;  m[2][0][2] =  0;
1996     m[2][1][0] = sin_roll;  m[2][1][1] =  cos_roll;  m[2][1][2] =  0;
1997     m[2][2][0] = 0;         m[2][2][1] =  0;         m[2][2][2] =  1;
1998
1999     multiply_matrix(temp, m[rotation_order[0]], m[rotation_order[1]]);
2000     multiply_matrix(rot_mat, temp, m[rotation_order[2]]);
2001 }
2002
2003 /**
2004  * Rotate vector with given rotation matrix.
2005  *
2006  * @param rot_mat rotation matrix
2007  * @param vec vector
2008  */
2009 static inline void rotate(const float rot_mat[3][3],
2010                           float *vec)
2011 {
2012     const float x_tmp = vec[0] * rot_mat[0][0] + vec[1] * rot_mat[0][1] + vec[2] * rot_mat[0][2];
2013     const float y_tmp = vec[0] * rot_mat[1][0] + vec[1] * rot_mat[1][1] + vec[2] * rot_mat[1][2];
2014     const float z_tmp = vec[0] * rot_mat[2][0] + vec[1] * rot_mat[2][1] + vec[2] * rot_mat[2][2];
2015
2016     vec[0] = x_tmp;
2017     vec[1] = y_tmp;
2018     vec[2] = z_tmp;
2019 }
2020
2021 static inline void set_mirror_modifier(int h_flip, int v_flip, int d_flip,
2022                                        float *modifier)
2023 {
2024     modifier[0] = h_flip ? -1.f : 1.f;
2025     modifier[1] = v_flip ? -1.f : 1.f;
2026     modifier[2] = d_flip ? -1.f : 1.f;
2027 }
2028
2029 static inline void mirror(const float *modifier, float *vec)
2030 {
2031     vec[0] *= modifier[0];
2032     vec[1] *= modifier[1];
2033     vec[2] *= modifier[2];
2034 }
2035
2036 static int allocate_plane(V360Context *s, int sizeof_uv, int sizeof_ker, int p)
2037 {
2038     s->u[p] = av_calloc(s->planewidth[p] * s->planeheight[p], sizeof_uv);
2039     s->v[p] = av_calloc(s->planewidth[p] * s->planeheight[p], sizeof_uv);
2040     if (!s->u[p] || !s->v[p])
2041         return AVERROR(ENOMEM);
2042     if (sizeof_ker) {
2043         s->ker[p] = av_calloc(s->planewidth[p] * s->planeheight[p], sizeof_ker);
2044         if (!s->ker[p])
2045             return AVERROR(ENOMEM);
2046     }
2047
2048     return 0;
2049 }
2050
2051 static int config_output(AVFilterLink *outlink)
2052 {
2053     AVFilterContext *ctx = outlink->src;
2054     AVFilterLink *inlink = ctx->inputs[0];
2055     V360Context *s = ctx->priv;
2056     const AVPixFmtDescriptor *desc = av_pix_fmt_desc_get(inlink->format);
2057     const int depth = desc->comp[0].depth;
2058     int sizeof_uv;
2059     int sizeof_ker;
2060     int elements;
2061     int err;
2062     int p, h, w;
2063     float hf, wf;
2064     float output_mirror_modifier[3];
2065     void (*in_transform)(const V360Context *s,
2066                          const float *vec, int width, int height,
2067                          uint16_t us[4][4], uint16_t vs[4][4], float *du, float *dv);
2068     void (*out_transform)(const V360Context *s,
2069                           int i, int j, int width, int height,
2070                           float *vec);
2071     void (*calculate_kernel)(float du, float dv, const XYRemap *r_tmp,
2072                              uint16_t *u, uint16_t *v, int16_t *ker);
2073     float rot_mat[3][3];
2074
2075     s->input_mirror_modifier[0] = s->ih_flip ? -1.f : 1.f;
2076     s->input_mirror_modifier[1] = s->iv_flip ? -1.f : 1.f;
2077
2078     switch (s->interp) {
2079     case NEAREST:
2080         calculate_kernel = nearest_kernel;
2081         s->remap_slice = depth <= 8 ? remap1_8bit_slice : remap1_16bit_slice;
2082         elements = 1;
2083         sizeof_uv = sizeof(uint16_t) * elements;
2084         sizeof_ker = 0;
2085         break;
2086     case BILINEAR:
2087         calculate_kernel = bilinear_kernel;
2088         s->remap_slice = depth <= 8 ? remap2_8bit_slice : remap2_16bit_slice;
2089         elements = 2 * 2;
2090         sizeof_uv = sizeof(uint16_t) * elements;
2091         sizeof_ker = sizeof(uint16_t) * elements;
2092         break;
2093     case BICUBIC:
2094         calculate_kernel = bicubic_kernel;
2095         s->remap_slice = depth <= 8 ? remap4_8bit_slice : remap4_16bit_slice;
2096         elements = 4 * 4;
2097         sizeof_uv = sizeof(uint16_t) * elements;
2098         sizeof_ker = sizeof(uint16_t) * elements;
2099         break;
2100     case LANCZOS:
2101         calculate_kernel = lanczos_kernel;
2102         s->remap_slice = depth <= 8 ? remap4_8bit_slice : remap4_16bit_slice;
2103         elements = 4 * 4;
2104         sizeof_uv = sizeof(uint16_t) * elements;
2105         sizeof_ker = sizeof(uint16_t) * elements;
2106         break;
2107     default:
2108         av_assert0(0);
2109     }
2110
2111     ff_v360_init(s, depth);
2112
2113     for (int order = 0; order < NB_RORDERS; order++) {
2114         const char c = s->rorder[order];
2115         int rorder;
2116
2117         if (c == '\0') {
2118             av_log(ctx, AV_LOG_ERROR,
2119                    "Incomplete rorder option. Direction for all 3 rotation orders should be specified.\n");
2120             return AVERROR(EINVAL);
2121         }
2122
2123         rorder = get_rorder(c);
2124         if (rorder == -1) {
2125             av_log(ctx, AV_LOG_ERROR,
2126                    "Incorrect rotation order symbol '%c' in rorder option.\n", c);
2127             return AVERROR(EINVAL);
2128         }
2129
2130         s->rotation_order[order] = rorder;
2131     }
2132
2133     switch (s->in) {
2134     case EQUIRECTANGULAR:
2135         in_transform = xyz_to_equirect;
2136         err = 0;
2137         wf = inlink->w;
2138         hf = inlink->h;
2139         break;
2140     case CUBEMAP_3_2:
2141         in_transform = xyz_to_cube3x2;
2142         err = prepare_cube_in(ctx);
2143         wf = inlink->w / 3.f * 4.f;
2144         hf = inlink->h;
2145         break;
2146     case CUBEMAP_1_6:
2147         in_transform = xyz_to_cube1x6;
2148         err = prepare_cube_in(ctx);
2149         wf = inlink->w * 4.f;
2150         hf = inlink->h / 3.f;
2151         break;
2152     case CUBEMAP_6_1:
2153         in_transform = xyz_to_cube6x1;
2154         err = prepare_cube_in(ctx);
2155         wf = inlink->w / 3.f * 2.f;
2156         hf = inlink->h * 2.f;
2157         break;
2158     case EQUIANGULAR:
2159         in_transform = xyz_to_eac;
2160         err = prepare_eac_in(ctx);
2161         wf = inlink->w;
2162         hf = inlink->h / 9.f * 8.f;
2163         break;
2164     case FLAT:
2165         av_log(ctx, AV_LOG_ERROR, "Flat format is not accepted as input.\n");
2166         return AVERROR(EINVAL);
2167     case DUAL_FISHEYE:
2168         in_transform = xyz_to_dfisheye;
2169         err = 0;
2170         wf = inlink->w;
2171         hf = inlink->h;
2172         break;
2173     case BARREL:
2174         in_transform = xyz_to_barrel;
2175         err = 0;
2176         wf = inlink->w / 5.f * 4.f;
2177         hf = inlink->h;
2178         break;
2179     default:
2180         av_log(ctx, AV_LOG_ERROR, "Specified input format is not handled.\n");
2181         return AVERROR_BUG;
2182     }
2183
2184     if (err != 0) {
2185         return err;
2186     }
2187
2188     switch (s->out) {
2189     case EQUIRECTANGULAR:
2190         out_transform = equirect_to_xyz;
2191         err = 0;
2192         w = roundf(wf);
2193         h = roundf(hf);
2194         break;
2195     case CUBEMAP_3_2:
2196         out_transform = cube3x2_to_xyz;
2197         err = prepare_cube_out(ctx);
2198         w = roundf(wf / 4.f * 3.f);
2199         h = roundf(hf);
2200         break;
2201     case CUBEMAP_1_6:
2202         out_transform = cube1x6_to_xyz;
2203         err = prepare_cube_out(ctx);
2204         w = roundf(wf / 4.f);
2205         h = roundf(hf * 3.f);
2206         break;
2207     case CUBEMAP_6_1:
2208         out_transform = cube6x1_to_xyz;
2209         err = prepare_cube_out(ctx);
2210         w = roundf(wf / 2.f * 3.f);
2211         h = roundf(hf / 2.f);
2212         break;
2213     case EQUIANGULAR:
2214         out_transform = eac_to_xyz;
2215         err = prepare_eac_out(ctx);
2216         w = roundf(wf);
2217         h = roundf(hf / 8.f * 9.f);
2218         break;
2219     case FLAT:
2220         out_transform = flat_to_xyz;
2221         err = prepare_flat_out(ctx);
2222         w = roundf(wf * s->flat_range[0] / s->flat_range[1] / 2.f);
2223         h = roundf(hf);
2224         break;
2225     case DUAL_FISHEYE:
2226         av_log(ctx, AV_LOG_ERROR, "Dual fisheye format is not accepted as output.\n");
2227         return AVERROR(EINVAL);
2228     case BARREL:
2229         out_transform = barrel_to_xyz;
2230         err = 0;
2231         w = roundf(wf / 4.f * 5.f);
2232         h = roundf(hf);
2233         break;
2234     case STEREOGRAPHIC:
2235         out_transform = stereographic_to_xyz;
2236         err = 0;
2237         w = FFMAX(roundf(wf), roundf(hf));
2238         h = w;
2239         break;
2240     default:
2241         av_log(ctx, AV_LOG_ERROR, "Specified output format is not handled.\n");
2242         return AVERROR_BUG;
2243     }
2244
2245     if (err != 0) {
2246         return err;
2247     }
2248
2249     // Override resolution with user values if specified
2250     if (s->width > 0 && s->height > 0) {
2251         w = s->width;
2252         h = s->height;
2253     } else if (s->width > 0 || s->height > 0) {
2254         av_log(ctx, AV_LOG_ERROR, "Both width and height values should be specified.\n");
2255         return AVERROR(EINVAL);
2256     } else {
2257         if (s->out_transpose)
2258             FFSWAP(int, w, h);
2259
2260         if (s->in_transpose)
2261             FFSWAP(int, w, h);
2262     }
2263
2264     s->planeheight[1] = s->planeheight[2] = FF_CEIL_RSHIFT(h, desc->log2_chroma_h);
2265     s->planeheight[0] = s->planeheight[3] = h;
2266     s->planewidth[1]  = s->planewidth[2] = FF_CEIL_RSHIFT(w, desc->log2_chroma_w);
2267     s->planewidth[0]  = s->planewidth[3] = w;
2268
2269     outlink->h = h;
2270     outlink->w = w;
2271
2272     s->inplaneheight[1] = s->inplaneheight[2] = FF_CEIL_RSHIFT(inlink->h, desc->log2_chroma_h);
2273     s->inplaneheight[0] = s->inplaneheight[3] = inlink->h;
2274     s->inplanewidth[1]  = s->inplanewidth[2]  = FF_CEIL_RSHIFT(inlink->w, desc->log2_chroma_w);
2275     s->inplanewidth[0]  = s->inplanewidth[3]  = inlink->w;
2276     s->nb_planes = av_pix_fmt_count_planes(inlink->format);
2277
2278     if (desc->log2_chroma_h == desc->log2_chroma_w && desc->log2_chroma_h == 0) {
2279         s->nb_allocated = 1;
2280         s->map[0] = s->map[1] = s->map[2] = s->map[3] = 0;
2281         allocate_plane(s, sizeof_uv, sizeof_ker, 0);
2282     } else if (desc->log2_chroma_h == desc->log2_chroma_w) {
2283         s->nb_allocated = 2;
2284         s->map[0] = 0;
2285         s->map[1] = s->map[2] = 1;
2286         s->map[3] = 0;
2287         allocate_plane(s, sizeof_uv, sizeof_ker, 0);
2288         allocate_plane(s, sizeof_uv, sizeof_ker, 1);
2289     } else {
2290         s->nb_allocated = 3;
2291         s->map[0] = 0;
2292         s->map[1] = 1;
2293         s->map[2] = 2;
2294         s->map[3] = 0;
2295         allocate_plane(s, sizeof_uv, sizeof_ker, 0);
2296         allocate_plane(s, sizeof_uv, sizeof_ker, 1);
2297         allocate_plane(s, sizeof_uv, sizeof_ker, 2);
2298     }
2299
2300     calculate_rotation_matrix(s->yaw, s->pitch, s->roll, rot_mat, s->rotation_order);
2301     set_mirror_modifier(s->h_flip, s->v_flip, s->d_flip, output_mirror_modifier);
2302
2303     // Calculate remap data
2304     for (p = 0; p < s->nb_allocated; p++) {
2305         const int width = s->planewidth[p];
2306         const int height = s->planeheight[p];
2307         const int in_width = s->inplanewidth[p];
2308         const int in_height = s->inplaneheight[p];
2309         float du, dv;
2310         float vec[3];
2311         XYRemap r_tmp;
2312         int i, j;
2313
2314         for (i = 0; i < width; i++) {
2315             for (j = 0; j < height; j++) {
2316                 uint16_t *u = s->u[p] + (j * width + i) * elements;
2317                 uint16_t *v = s->v[p] + (j * width + i) * elements;
2318                 int16_t *ker = s->ker[p] + (j * width + i) * elements;
2319
2320                 if (s->out_transpose)
2321                     out_transform(s, j, i, height, width, vec);
2322                 else
2323                     out_transform(s, i, j, width, height, vec);
2324                 av_assert1(!isnan(vec[0]) && !isnan(vec[1]) && !isnan(vec[2]));
2325                 rotate(rot_mat, vec);
2326                 av_assert1(!isnan(vec[0]) && !isnan(vec[1]) && !isnan(vec[2]));
2327                 normalize_vector(vec);
2328                 mirror(output_mirror_modifier, vec);
2329                 if (s->in_transpose)
2330                     in_transform(s, vec, in_height, in_width, r_tmp.v, r_tmp.u, &du, &dv);
2331                 else
2332                     in_transform(s, vec, in_width, in_height, r_tmp.u, r_tmp.v, &du, &dv);
2333                 av_assert1(!isnan(du) && !isnan(dv));
2334                 calculate_kernel(du, dv, &r_tmp, u, v, ker);
2335             }
2336         }
2337     }
2338
2339     return 0;
2340 }
2341
2342 static int filter_frame(AVFilterLink *inlink, AVFrame *in)
2343 {
2344     AVFilterContext *ctx = inlink->dst;
2345     AVFilterLink *outlink = ctx->outputs[0];
2346     V360Context *s = ctx->priv;
2347     AVFrame *out;
2348     ThreadData td;
2349
2350     out = ff_get_video_buffer(outlink, outlink->w, outlink->h);
2351     if (!out) {
2352         av_frame_free(&in);
2353         return AVERROR(ENOMEM);
2354     }
2355     av_frame_copy_props(out, in);
2356
2357     td.in = in;
2358     td.out = out;
2359
2360     ctx->internal->execute(ctx, s->remap_slice, &td, NULL, FFMIN(outlink->h, ff_filter_get_nb_threads(ctx)));
2361
2362     av_frame_free(&in);
2363     return ff_filter_frame(outlink, out);
2364 }
2365
2366 static av_cold void uninit(AVFilterContext *ctx)
2367 {
2368     V360Context *s = ctx->priv;
2369     int p;
2370
2371     for (p = 0; p < s->nb_allocated; p++) {
2372         av_freep(&s->u[p]);
2373         av_freep(&s->v[p]);
2374         av_freep(&s->ker[p]);
2375     }
2376 }
2377
2378 static const AVFilterPad inputs[] = {
2379     {
2380         .name         = "default",
2381         .type         = AVMEDIA_TYPE_VIDEO,
2382         .filter_frame = filter_frame,
2383     },
2384     { NULL }
2385 };
2386
2387 static const AVFilterPad outputs[] = {
2388     {
2389         .name         = "default",
2390         .type         = AVMEDIA_TYPE_VIDEO,
2391         .config_props = config_output,
2392     },
2393     { NULL }
2394 };
2395
2396 AVFilter ff_vf_v360 = {
2397     .name          = "v360",
2398     .description   = NULL_IF_CONFIG_SMALL("Convert 360 projection of video."),
2399     .priv_size     = sizeof(V360Context),
2400     .uninit        = uninit,
2401     .query_formats = query_formats,
2402     .inputs        = inputs,
2403     .outputs       = outputs,
2404     .priv_class    = &v360_class,
2405     .flags         = AVFILTER_FLAG_SLICE_THREADS,
2406 };