diff --git a/doc/filters.texi b/doc/filters.texi
index 7400e7dd31..dcb98f772b 100644
--- a/doc/filters.texi
+++ b/doc/filters.texi
@@ -11152,6 +11152,64 @@ geq=lum=255*gauss((X/W-0.5)*3)*gauss((Y/H-0.5)*3)/gauss(0)/gauss(0),format=gray
 @end example
 @end itemize
 
+@section gmsd
+
+Obtain the GMSD (Gradient Magnitude Similarity Deviation) between two input videos.
+
+This filter takes in input two input videos, the first input is
+considered the "main" source and is passed unchanged to the
+output. The second input is used as a "reference" video for computing
+the GMSD.
+
+Both video inputs must have the same resolution and pixel format for
+this filter to work correctly. Also it assumes that both inputs
+have the same number of frames, which are compared one by one.
+
+The filter stores the calculated GMSD of each frame.
+
+The description of the accepted parameters follows.
+
+@table @option
+@item stats_file, f
+If specified the filter will use the named file to save the SSIM of
+each individual frame. When filename equals "-" the data is sent to
+standard output.
+@end table
+
+The file printed if @var{stats_file} is selected, contains a sequence of
+key/value pairs of the form @var{key}:@var{value} for each compared
+couple of frames.
+
+A description of each shown parameter follows:
+
+@table @option
+@item n
+sequential number of the input frame, starting from 1
+
+@item Y, U, V, R, G, B
+GMSD of the compared frames for the component specified by the suffix.
+
+@item All
+GMSD of the compared frames for the whole frame.
+@end table
+
+This filter also supports the @ref{framesync} options.
+
+For example:
+@example
+movie=ref_movie.mpg, setpts=PTS-STARTPTS [main];
+[main][ref] gmsd="stats_file=stats.log" [out]
+@end example
+
+On this example the input file being processed is compared with the
+reference file @file{ref_movie.mpg}. The GMSD of each individual frame
+is stored in @file{stats.log}.
+
+Another example with both gmsd and ssim at same time:
+@example
+ffmpeg -i main.mpg -i ref.mpg -lavfi  "ssim;[0:v][1:v]gmsd" -f null -
+@end example
+
 @section gradfun
 
 Fix the banding artifacts that are sometimes introduced into nearly flat
diff --git a/libavfilter/Makefile b/libavfilter/Makefile
index 63d2fba861..d966b9532b 100644
--- a/libavfilter/Makefile
+++ b/libavfilter/Makefile
@@ -256,6 +256,7 @@ OBJS-$(CONFIG_FREI0R_FILTER)                 += vf_frei0r.o
 OBJS-$(CONFIG_FSPP_FILTER)                   += vf_fspp.o
 OBJS-$(CONFIG_GBLUR_FILTER)                  += vf_gblur.o
 OBJS-$(CONFIG_GEQ_FILTER)                    += vf_geq.o
+OBJS-$(CONFIG_GMSD_FILTER)                   += vf_gmsd.o framesync.o
 OBJS-$(CONFIG_GRADFUN_FILTER)                += vf_gradfun.o
 OBJS-$(CONFIG_GRAPHMONITOR_FILTER)           += f_graphmonitor.o
 OBJS-$(CONFIG_GREYEDGE_FILTER)               += vf_colorconstancy.o
diff --git a/libavfilter/allfilters.c b/libavfilter/allfilters.c
index e4186f93db..d98c5971ac 100644
--- a/libavfilter/allfilters.c
+++ b/libavfilter/allfilters.c
@@ -241,6 +241,7 @@ extern AVFilter ff_vf_frei0r;
 extern AVFilter ff_vf_fspp;
 extern AVFilter ff_vf_gblur;
 extern AVFilter ff_vf_geq;
+extern AVFilter ff_vf_gmsd;
 extern AVFilter ff_vf_gradfun;
 extern AVFilter ff_vf_graphmonitor;
 extern AVFilter ff_vf_greyedge;
diff --git a/libavfilter/vf_gmsd.c b/libavfilter/vf_gmsd.c
new file mode 100644
index 0000000000..894e079fa7
--- /dev/null
+++ b/libavfilter/vf_gmsd.c
@@ -0,0 +1,376 @@
+/*
+ * Copyright (c) 2019 Paul B Mahol
+ *
+ * This file is part of FFmpeg.
+ *
+ * FFmpeg is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * FFmpeg is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with FFmpeg; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+/*
+ * @file
+ * Caculate the GMSD between two input videos.
+ */
+
+#include "libavutil/avstring.h"
+#include "libavutil/opt.h"
+#include "libavutil/pixdesc.h"
+#include "avfilter.h"
+#include "drawutils.h"
+#include "formats.h"
+#include "framesync.h"
+#include "internal.h"
+#include "video.h"
+
+typedef struct GMSDContext {
+    const AVClass *class;
+    FFFrameSync fs;
+    FILE *stats_file;
+    char *stats_file_str;
+    int nb_components;
+    uint64_t nb_frames;
+    double gmsd[4], gmsd_total;
+    char comps[4];
+    float coefs[4];
+    uint8_t rgba_map[4];
+    int planewidth[4];
+    int planeheight[4];
+    float *gms;
+    int is_rgb;
+    float (*gmsd_plane)(const uint8_t *master,
+                        int master_linesize,
+                        const uint8_t *ref, int ref_linesize,
+                        int w, int h, float *gms);
+} GMSDContext;
+
+#define OFFSET(x) offsetof(GMSDContext, x)
+#define FLAGS AV_OPT_FLAG_FILTERING_PARAM|AV_OPT_FLAG_VIDEO_PARAM
+
+static const AVOption gmsd_options[] = {
+    {"stats_file", "Set file where to store per-frame difference information", OFFSET(stats_file_str), AV_OPT_TYPE_STRING, {.str=NULL}, 0, 0, FLAGS },
+    {"f",          "Set file where to store per-frame difference information", OFFSET(stats_file_str), AV_OPT_TYPE_STRING, {.str=NULL}, 0, 0, FLAGS },
+    { NULL }
+};
+
+FRAMESYNC_DEFINE_CLASS(gmsd, GMSDContext, fs);
+
+static void set_meta(AVDictionary **metadata, const char *key, char comp, float d)
+{
+    char value[128];
+    snprintf(value, sizeof(value), "%0.2f", d);
+    if (comp) {
+        char key2[128];
+        snprintf(key2, sizeof(key2), "%s%c", key, comp);
+        av_dict_set(metadata, key2, value, 0);
+    } else {
+        av_dict_set(metadata, key, value, 0);
+    }
+}
+
+#define SQR(x) ((x) * (x))
+
+#define GMSD_PLANE(name, type, div)                                    \
+static float gmsd_plane##name(const uint8_t *mmaster,                  \
+                              int master_linesize,                     \
+                              const uint8_t *rref, int ref_linesize,   \
+                              int w, int h, float *gms)                \
+{                                                                      \
+    const type *master = (const type *)mmaster;                        \
+    const type *ref = (const type *)rref;                              \
+    float gmsm = 0.f;                                                  \
+    float gmsd = 0.f;                                                  \
+    int i = 0;                                                         \
+                                                                       \
+    master_linesize /= div;                                            \
+    ref_linesize /= div;                                               \
+                                                                       \
+    master += master_linesize;                                         \
+    ref    += ref_linesize;                                            \
+                                                                       \
+    for (int y = 1; y < h - 1; y++) {                                  \
+        for (int x = 1; x < w - 1; x++) {                              \
+            float masx = master[x - master_linesize - 1] * 1.f / 3.f + \
+                         master[x                   - 1] * 1.f / 3.f + \
+                         master[x + master_linesize - 1] * 1.f / 3.f + \
+                         master[x - master_linesize + 1] *-1.f / 3.f + \
+                         master[x                   + 1] *-1.f / 3.f + \
+                         master[x + master_linesize + 1] *-1.f / 3.f;  \
+            float masy = master[x - master_linesize - 1] * 1.f / 3.f + \
+                         master[x - master_linesize + 0] * 1.f / 3.f + \
+                         master[x - master_linesize + 1] * 1.f / 3.f + \
+                         master[x + master_linesize - 1] *-1.f / 3.f + \
+                         master[x + master_linesize + 0] *-1.f / 3.f + \
+                         master[x + master_linesize + 1] *-1.f / 3.f;  \
+                                                                       \
+            float refx = ref[x - ref_linesize - 1] * 1.f / 3.f +       \
+                         ref[x                - 1] * 1.f / 3.f +       \
+                         ref[x + ref_linesize - 1] * 1.f / 3.f +       \
+                         ref[x - ref_linesize + 1] *-1.f / 3.f +       \
+                         ref[x                + 1] *-1.f / 3.f +       \
+                         ref[x + ref_linesize + 1] *-1.f / 3.f;        \
+            float refy = ref[x - ref_linesize - 1] * 1.f / 3.f +       \
+                         ref[x - ref_linesize + 0] * 1.f / 3.f +       \
+                         ref[x - ref_linesize + 1] * 1.f / 3.f +       \
+                         ref[x + ref_linesize - 1] *-1.f / 3.f +       \
+                         ref[x + ref_linesize + 0] *-1.f / 3.f +       \
+                         ref[x + ref_linesize + 1] *-1.f / 3.f;        \
+                                                                       \
+            float gm = masx * masx + masy * masy;                      \
+            float gr = refx * refx + refy * refy;                      \
+                                                                       \
+            gms[i] = (2.f * sqrtf(gm * gr) + 0.0026) /                 \
+                     (gm + gr + 0.0026);                               \
+            gmsm += gms[i];                                            \
+            i++;                                                       \
+        }                                                              \
+                                                                       \
+        master += master_linesize;                                     \
+        ref    += ref_linesize;                                        \
+    }                                                                  \
+                                                                       \
+    gmsm /= i;                                                         \
+                                                                       \
+    for (int j = 0; j < i; j++) {                                      \
+        gmsd += SQR(gms[j] - gmsm);                                    \
+    }                                                                  \
+                                                                       \
+    gmsd /= i;                                                         \
+                                                                       \
+    return sqrtf(gmsd);                                                \
+}
+
+GMSD_PLANE(8, uint8_t, 1)
+GMSD_PLANE(16, uint16_t, 2)
+
+static int do_gmsd(FFFrameSync *fs)
+{
+    AVFilterContext *ctx = fs->parent;
+    GMSDContext *s = ctx->priv;
+    AVFrame *master, *ref;
+    AVDictionary **metadata;
+    float c[4], gmsdv = 0.0;
+    int ret, i;
+
+    ret = ff_framesync_dualinput_get(fs, &master, &ref);
+    if (ret < 0)
+        return ret;
+    if (!ref)
+        return ff_filter_frame(ctx->outputs[0], master);
+    metadata = &master->metadata;
+
+    s->nb_frames++;
+
+    for (i = 0; i < s->nb_components; i++) {
+        c[i] = s->gmsd_plane(master->data[i], master->linesize[i],
+                             ref->data[i], ref->linesize[i],
+                             s->planewidth[i], s->planeheight[i], s->gms);
+        gmsdv += s->coefs[i] * c[i];
+        s->gmsd[i] += c[i];
+    }
+    for (i = 0; i < s->nb_components; i++) {
+        int cidx = s->is_rgb ? s->rgba_map[i] : i;
+        set_meta(metadata, "lavfi.gmsd.", s->comps[i], c[cidx]);
+    }
+    s->gmsd_total += gmsdv;
+
+    set_meta(metadata, "lavfi.gmsd.All", 0, gmsdv);
+
+    if (s->stats_file) {
+        fprintf(s->stats_file, "n:%"PRId64" ", s->nb_frames);
+
+        for (i = 0; i < s->nb_components; i++) {
+            int cidx = s->is_rgb ? s->rgba_map[i] : i;
+            fprintf(s->stats_file, "%c:%f ", s->comps[i], c[cidx]);
+        }
+
+        fprintf(s->stats_file, "All:%f\n", gmsdv);
+    }
+
+    return ff_filter_frame(ctx->outputs[0], master);
+}
+
+static av_cold int init(AVFilterContext *ctx)
+{
+    GMSDContext *s = ctx->priv;
+
+    if (s->stats_file_str) {
+        if (!strcmp(s->stats_file_str, "-")) {
+            s->stats_file = stdout;
+        } else {
+            s->stats_file = fopen(s->stats_file_str, "w");
+            if (!s->stats_file) {
+                int err = AVERROR(errno);
+                char buf[128];
+                av_strerror(err, buf, sizeof(buf));
+                av_log(ctx, AV_LOG_ERROR, "Could not open stats file %s: %s\n",
+                       s->stats_file_str, buf);
+                return err;
+            }
+        }
+    }
+
+    s->fs.on_event = do_gmsd;
+
+    return 0;
+}
+
+static int query_formats(AVFilterContext *ctx)
+{
+    static const enum AVPixelFormat pix_fmts[] = {
+        AV_PIX_FMT_GRAY8, AV_PIX_FMT_GRAY9, AV_PIX_FMT_GRAY10,
+        AV_PIX_FMT_GRAY12, AV_PIX_FMT_GRAY14, AV_PIX_FMT_GRAY16,
+        AV_PIX_FMT_YUV420P, AV_PIX_FMT_YUV422P, AV_PIX_FMT_YUV444P,
+        AV_PIX_FMT_YUV440P, AV_PIX_FMT_YUV411P, AV_PIX_FMT_YUV410P,
+        AV_PIX_FMT_YUVJ411P, AV_PIX_FMT_YUVJ420P, AV_PIX_FMT_YUVJ422P,
+        AV_PIX_FMT_YUVJ440P, AV_PIX_FMT_YUVJ444P,
+        AV_PIX_FMT_GBRP,
+#define PF(suf) AV_PIX_FMT_YUV420##suf,  AV_PIX_FMT_YUV422##suf,  AV_PIX_FMT_YUV444##suf, AV_PIX_FMT_GBR##suf
+        PF(P9), PF(P10), PF(P12), PF(P14), PF(P16),
+        AV_PIX_FMT_NONE
+    };
+
+    AVFilterFormats *fmts_list = ff_make_format_list(pix_fmts);
+    if (!fmts_list)
+        return AVERROR(ENOMEM);
+    return ff_set_common_formats(ctx, fmts_list);
+}
+
+static int config_input_ref(AVFilterLink *inlink)
+{
+    const AVPixFmtDescriptor *desc = av_pix_fmt_desc_get(inlink->format);
+    AVFilterContext *ctx  = inlink->dst;
+    GMSDContext *s = ctx->priv;
+    int sum = 0, i;
+
+    s->nb_components = desc->nb_components;
+
+    if (ctx->inputs[0]->w != ctx->inputs[1]->w ||
+        ctx->inputs[0]->h != ctx->inputs[1]->h) {
+        av_log(ctx, AV_LOG_ERROR, "Width and height of input videos must be same.\n");
+        return AVERROR(EINVAL);
+    }
+    if (ctx->inputs[0]->format != ctx->inputs[1]->format) {
+        av_log(ctx, AV_LOG_ERROR, "Inputs must be of same pixel format.\n");
+        return AVERROR(EINVAL);
+    }
+
+    s->is_rgb = ff_fill_rgba_map(s->rgba_map, inlink->format) >= 0;
+    s->comps[0] = s->is_rgb ? 'R' : 'Y';
+    s->comps[1] = s->is_rgb ? 'G' : 'U';
+    s->comps[2] = s->is_rgb ? 'B' : 'V';
+    s->comps[3] = 'A';
+
+    s->planeheight[1] = s->planeheight[2] = AV_CEIL_RSHIFT(inlink->h, desc->log2_chroma_h);
+    s->planeheight[0] = s->planeheight[3] = inlink->h;
+    s->planewidth[1]  = s->planewidth[2]  = AV_CEIL_RSHIFT(inlink->w, desc->log2_chroma_w);
+    s->planewidth[0]  = s->planewidth[3]  = inlink->w;
+    for (i = 0; i < s->nb_components; i++)
+        sum += s->planeheight[i] * s->planewidth[i];
+    for (i = 0; i < s->nb_components; i++)
+        s->coefs[i] = (double) s->planeheight[i] * s->planewidth[i] / sum;
+
+    s->gms = av_calloc(s->planewidth[0] * s->planeheight[0], sizeof(*s->gms));
+    if (!s->gms)
+        return AVERROR(ENOMEM);
+    s->gmsd_plane = (1 << desc->comp[0].depth) <= 8 ? gmsd_plane8 : gmsd_plane16;
+
+    return 0;
+}
+
+static int config_output(AVFilterLink *outlink)
+{
+    AVFilterContext *ctx = outlink->src;
+    GMSDContext *s = ctx->priv;
+    AVFilterLink *mainlink = ctx->inputs[0];
+    int ret;
+
+    ret = ff_framesync_init_dualinput(&s->fs, ctx);
+    if (ret < 0)
+        return ret;
+    outlink->w = mainlink->w;
+    outlink->h = mainlink->h;
+    outlink->time_base = mainlink->time_base;
+    outlink->sample_aspect_ratio = mainlink->sample_aspect_ratio;
+    outlink->frame_rate = mainlink->frame_rate;
+
+    if ((ret = ff_framesync_configure(&s->fs)) < 0)
+        return ret;
+
+    return 0;
+}
+
+static int activate(AVFilterContext *ctx)
+{
+    GMSDContext *s = ctx->priv;
+    return ff_framesync_activate(&s->fs);
+}
+
+static av_cold void uninit(AVFilterContext *ctx)
+{
+    GMSDContext *s = ctx->priv;
+
+    if (s->nb_frames > 0) {
+        char buf[256];
+        int i;
+        buf[0] = 0;
+        for (i = 0; i < s->nb_components; i++) {
+            int c = s->is_rgb ? s->rgba_map[i] : i;
+            av_strlcatf(buf, sizeof(buf), " %c:%f ", s->comps[i], s->gmsd[c] / s->nb_frames);
+        }
+        av_log(ctx, AV_LOG_INFO, "GMSD%s All:%f\n", buf,
+               s->gmsd_total / s->nb_frames);
+    }
+
+    ff_framesync_uninit(&s->fs);
+
+    if (s->stats_file && s->stats_file != stdout)
+        fclose(s->stats_file);
+
+    av_freep(&s->gms);
+}
+
+static const AVFilterPad gmsd_inputs[] = {
+    {
+        .name         = "main",
+        .type         = AVMEDIA_TYPE_VIDEO,
+    },{
+        .name         = "reference",
+        .type         = AVMEDIA_TYPE_VIDEO,
+        .config_props = config_input_ref,
+    },
+    { NULL }
+};
+
+static const AVFilterPad gmsd_outputs[] = {
+    {
+        .name          = "default",
+        .type          = AVMEDIA_TYPE_VIDEO,
+        .config_props  = config_output,
+    },
+    { NULL }
+};
+
+AVFilter ff_vf_gmsd = {
+    .name          = "gmsd",
+    .description   = NULL_IF_CONFIG_SMALL("Calculate the GMSD between two video streams."),
+    .preinit       = gmsd_framesync_preinit,
+    .init          = init,
+    .uninit        = uninit,
+    .query_formats = query_formats,
+    .activate      = activate,
+    .priv_size     = sizeof(GMSDContext),
+    .priv_class    = &gmsd_class,
+    .inputs        = gmsd_inputs,
+    .outputs       = gmsd_outputs,
+};
