diff --git a/doc/filters.texi b/doc/filters.texi
index b4f888c14d..d066041e21 100644
--- a/doc/filters.texi
+++ b/doc/filters.texi
@@ -7712,6 +7712,61 @@ tblend=all_mode=grainextract
 @subsection Commands
 This filter supports same @ref{commands} as options.
 
+@anchor{blurriness}
+@section blurriness
+
+Determines blurriness of frames without altering the input frames.
+
+Based on Marziliano, Pina, et al. "A no-reference perceptual blur metric."
+Allows for a block-based abbreviation.
+
+The filter accepts the following options:
+
+@table @option
+@item low
+@item high
+Set low and high threshold values used by the Canny thresholding
+algorithm.
+
+The high threshold selects the "strong" edge pixels, which are then
+connected through 8-connectivity with the "weak" edge pixels selected
+by the low threshold.
+
+@var{low} and @var{high} threshold values must be chosen in the range
+[0,1], and @var{low} should be lesser or equal to @var{high}.
+
+Default value for @var{low} is @code{20/255}, and default value for @var{high}
+is @code{50/255}.
+
+@item radius
+Define the radius to search around an edge pixel for local maxima.
+
+@item block_pct
+Determine blurriness only for the most significant blocks, given in percentage.
+
+@item block_width
+Determine blurriness for blocks of width @var{block_width}. If set to any value smaller 1, no blocks are used and the whole image is processed as one no matter of @var{block_height}.
+
+@item block_height
+Determine blurriness for blocks of height @var{block_height}. If set to any value smaller 1, no blocks are used and the whole image is processed as one no matter of @var{block_width}.
+
+@item stats_file, f
+If specified the filter will use the named file to save the PSNR of
+each individual frame. When filename equals "-" the data is sent to
+standard output.
+@end table
+
+@subsection Examples
+
+@itemize
+@item
+Determine blurriness for 80% of most significant 32x32 blocks and write to stdout:
+@example
+blurriness=block_width=32:block_height=32:block_pct=80:f=-
+@end example
+@end itemize
+
+
 @section bm3d
 
 Denoise frames using Block-Matching 3D algorithm.
diff --git a/libavfilter/Makefile b/libavfilter/Makefile
index 08edc92d8c..be70635c82 100644
--- a/libavfilter/Makefile
+++ b/libavfilter/Makefile
@@ -188,6 +188,7 @@ OBJS-$(CONFIG_BITPLANENOISE_FILTER)          += vf_bitplanenoise.o
 OBJS-$(CONFIG_BLACKDETECT_FILTER)            += vf_blackdetect.o
 OBJS-$(CONFIG_BLACKFRAME_FILTER)             += vf_blackframe.o
 OBJS-$(CONFIG_BLEND_FILTER)                  += vf_blend.o framesync.o
+OBJS-$(CONFIG_BLURRINESS_FILTER)             += vf_blurriness.o
 OBJS-$(CONFIG_BM3D_FILTER)                   += vf_bm3d.o framesync.o
 OBJS-$(CONFIG_BOXBLUR_FILTER)                += vf_boxblur.o boxblur.o
 OBJS-$(CONFIG_BOXBLUR_OPENCL_FILTER)         += vf_avgblur_opencl.o opencl.o \
diff --git a/libavfilter/allfilters.c b/libavfilter/allfilters.c
index f250020159..121e37391d 100644
--- a/libavfilter/allfilters.c
+++ b/libavfilter/allfilters.c
@@ -180,6 +180,7 @@ extern const AVFilter ff_vf_bitplanenoise;
 extern const AVFilter ff_vf_blackdetect;
 extern const AVFilter ff_vf_blackframe;
 extern const AVFilter ff_vf_blend;
+extern const AVFilter ff_vf_blurriness;
 extern const AVFilter ff_vf_bm3d;
 extern const AVFilter ff_vf_boxblur;
 extern const AVFilter ff_vf_boxblur_opencl;
diff --git a/libavfilter/vf_blurriness.c b/libavfilter/vf_blurriness.c
new file mode 100644
index 0000000000..19cd2675f6
--- /dev/null
+++ b/libavfilter/vf_blurriness.c
@@ -0,0 +1,435 @@
+/*
+ * Copyright (c) 2021 Thilo Borgmann <thilo.borgmann _at_ mail.de>
+ *
+ * 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
+ * No-reference blurriness metric filter
+ *
+ * Implementing:
+ * Marziliano, Pina, et al. "A no-reference perceptual blur metric." Proceedings.
+ * International conference on image processing. Vol. 3. IEEE, 2002.
+ * https://infoscience.epfl.ch/record/111802/files/14%20A%20no-reference%20perceptual%20blur%20metric.pdf
+ *
+ * @author Thilo Borgmann <thilo.borgmann _at_ mail.de>
+ */
+
+#include "libavutil/avassert.h"
+#include "libavutil/imgutils.h"
+#include "libavutil/opt.h"
+#include "libavutil/pixelutils.h"
+#include "libavutil/motion_vector.h"
+#include "avfilter.h"
+#include "formats.h"
+#include "internal.h"
+#include "video.h"
+
+static int comp(const void *a,const void *b)
+{
+    float x = *(float*)a;
+    float y = *(float*)b;
+    if (x > y) return  1;
+    if (x < y) return -1;
+    return 0;
+}
+
+typedef struct BLRContext {
+    const AVClass *class;
+
+    FILE *stats_file;
+    char *stats_file_str;
+
+    float   low, high;
+    uint8_t low_u8, high_u8;
+    int     radius;        // radius during local maxima detection
+    int     block_pct;     // percentage of "sharpest" blocks in the image to use for bluriness calculation
+    int     block_width;   // width for block abbreviation
+    int     block_height;  // height for block abbreviation
+
+    float    blur_total;
+    uint64_t nb_frames;
+
+    float    *blks;
+    uint8_t  *filterbuf;
+    uint8_t  *tmpbuf;
+    uint16_t *gradients;
+    char     *directions;
+} BLRContext;
+
+#define OFFSET(x) offsetof(BLRContext, x)
+#define FLAGS AV_OPT_FLAG_FILTERING_PARAM|AV_OPT_FLAG_VIDEO_PARAM
+static const AVOption blr_options[] = {
+    { "high",          "set high threshold", OFFSET(high), AV_OPT_TYPE_FLOAT, {.dbl=30/255.}, 0, 1, FLAGS },
+    { "low",           "set low threshold",  OFFSET(low),  AV_OPT_TYPE_FLOAT, {.dbl=15/255.}, 0, 1, FLAGS },
+    { "radius",        "search radius for maxima detection", OFFSET(radius), AV_OPT_TYPE_INT, {.i64=50}, 1, 100, FLAGS },
+    { "block_pct",     "block pooling threshold when calculating blurriness", OFFSET(block_pct), AV_OPT_TYPE_INT, {.i64=80}, 1, 100, FLAGS },
+    { "block_width",   "block size for block-based abbreviation of blurriness", OFFSET(block_width), AV_OPT_TYPE_INT, {.i64=-1}, -1, INT_MAX, FLAGS },
+    { "block_height",  "block size for block-based abbreviation of blurriness", OFFSET(block_height), AV_OPT_TYPE_INT, {.i64=-1}, -1, INT_MAX, FLAGS },
+    { "stats_file",    "Set file where to store per-frame blur information", OFFSET(stats_file_str), AV_OPT_TYPE_STRING, {.str=NULL}, 0, 0, FLAGS },
+    { "f",             "Set file where to store per-frame blur information", OFFSET(stats_file_str), AV_OPT_TYPE_STRING, {.str=NULL}, 0, 0, FLAGS },
+    { NULL }
+};
+
+AVFILTER_DEFINE_CLASS(blr);
+
+static av_cold int blr_init(AVFilterContext *ctx)
+{
+    BLRContext *blr = ctx->priv;
+
+    blr->low_u8  = blr->low  * 255. + .5;
+    blr->high_u8 = blr->high * 255. + .5;
+
+    if (blr->stats_file_str) {
+        if (!strcmp(blr->stats_file_str, "-")) {
+            blr->stats_file = stdout;
+        } else {
+            blr->stats_file = fopen(blr->stats_file_str, "w");
+            if (!blr->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",
+                       blr->stats_file_str, buf);
+                return err;
+            }
+        }
+    }
+
+    return 0;
+}
+
+static int blr_config_input(AVFilterLink *inlink)
+{
+    AVFilterContext *ctx    = inlink->dst;
+    BLRContext      *blr    = ctx->priv;
+    const int       bufsize = inlink->w * inlink->h;
+
+    if (blr->block_width  < 1 || blr->block_height < 1) {
+        blr->block_width  = inlink->w;
+        blr->block_height = inlink->h;
+    }
+
+    blr->tmpbuf     = av_malloc(bufsize);
+    blr->filterbuf  = av_malloc(bufsize);
+    blr->gradients  = av_calloc(bufsize, sizeof(*blr->gradients));
+    blr->directions = av_malloc(bufsize);
+    if (!blr->tmpbuf || !blr->filterbuf || !blr->gradients || !blr->directions)
+        return AVERROR(ENOMEM);
+
+    blr->blks = av_calloc((inlink->w / blr->block_width) * (inlink->h / blr->block_height),
+                          sizeof(*blr->blks));
+
+    return 0;
+}
+
+static void gaussian_blur(int w, int h,
+                          uint8_t *dst, int dst_linesize,
+                          const uint8_t *src, int src_linesize)
+{
+    int i, j;
+
+    memcpy(dst, src, w); dst += dst_linesize; src += src_linesize;
+    memcpy(dst, src, w); dst += dst_linesize; src += src_linesize;
+    for (j = 2; j < h - 2; j++) {
+        dst[0] = src[0];
+        dst[1] = src[1];
+        for (i = 2; i < w - 2; i++) {
+            /* Gaussian mask of size 5x5 with sigma = 1.4 */
+            dst[i] = ((src[-2*src_linesize + i-2] + src[2*src_linesize + i-2]) * 2
+                      + (src[-2*src_linesize + i-1] + src[2*src_linesize + i-1]) * 4
+                      + (src[-2*src_linesize + i  ] + src[2*src_linesize + i  ]) * 5
+                      + (src[-2*src_linesize + i+1] + src[2*src_linesize + i+1]) * 4
+                      + (src[-2*src_linesize + i+2] + src[2*src_linesize + i+2]) * 2
+
+                      + (src[  -src_linesize + i-2] + src[  src_linesize + i-2]) *  4
+                      + (src[  -src_linesize + i-1] + src[  src_linesize + i-1]) *  9
+                      + (src[  -src_linesize + i  ] + src[  src_linesize + i  ]) * 12
+                      + (src[  -src_linesize + i+1] + src[  src_linesize + i+1]) *  9
+                      + (src[  -src_linesize + i+2] + src[  src_linesize + i+2]) *  4
+
+                      + src[i-2] *  5
+                      + src[i-1] * 12
+                      + src[i  ] * 15
+                      + src[i+1] * 12
+                      + src[i+2] *  5) / 159;
+        }
+        dst[i    ] = src[i    ];
+        dst[i + 1] = src[i + 1];
+
+        dst += dst_linesize;
+        src += src_linesize;
+    }
+    memcpy(dst, src, w); dst += dst_linesize; src += src_linesize;
+    memcpy(dst, src, w);
+}
+
+// edge width is defined as the distance between surrounding maxima of the edge pixel
+static float edge_width(BLRContext *blr, int i, int j, int8_t dir, int w, int h,
+                        int edge, const uint8_t *src, int src_linesize)
+{
+    float width = 0;
+    int dX, dY;
+    int sign;
+    int tmp;
+    int p1;
+    int p2;
+    int k, x, y;
+    int edge1;
+    int edge2;
+    float luma1 = 0.0; // average luma difference per edge pixel
+    float luma2 = 0.0;
+    int radius = blr->radius;
+
+    switch(dir)
+    {
+    case DIRECTION_HORIZONTAL: dX = 1; dY =  0; break;
+    case DIRECTION_VERTICAL:   dX = 0; dY =  1; break;
+    case DIRECTION_45UP:       dX = 1; dY = -1; break;
+    case DIRECTION_45DOWN:     dX = 1; dY =  1; break;
+    }
+    if (dir == DIRECTION_HORIZONTAL) return 0;
+
+    // determines if search in direction dX/dY is looking for a maximum or minimum
+    sign = src[j * src_linesize + i] > src[(j - dY) * src_linesize + i - dX] ? 1 : -1;
+
+    // search in -(dX/dY) direction
+    for (k = 0; k < radius; k++)
+    {
+        x = i - k*dX;
+        y = j - k*dY;
+        p1 = y * src_linesize + x;
+        x -= dX;
+        y -= dY;
+        p2 = y * src_linesize + x;
+        if (x < 0 || x >= w || y < 0 || y >= h)
+            return 0;
+
+        tmp = (src[p1] - src[p2]) * sign;
+
+        if (tmp <= 0) // local maximum found
+            break;
+
+        luma1 += tmp;
+    }
+    if (k > 0) luma1 /= k;
+    edge1 = k;
+    width += k;
+
+    // search in +(dX/dY) direction
+    for (k = 0; k < radius; k++)
+    {
+        x = i + k * dX;
+        y = j + k * dY;
+        p1 = y * src_linesize + x;
+        x += dX;
+        y += dY;
+        p2 = y * src_linesize + x;
+        if (x < 0 || x >= w || y < 0 || y >= h)
+            return 0;
+
+        tmp = (src[p1] - src[p2]) * sign;
+
+        if (tmp >= 0) // local maximum found
+            break;
+
+        luma2 -= tmp;
+    }
+    if (k > 0) luma2 /= k;
+    edge2 = k;
+    width += k;
+
+    // for 45 degree directions approximate edge width in pixel units: 0.7 ~= sqrt(2)/2
+    if (dir == DIRECTION_45UP || dir == DIRECTION_45DOWN)
+        width *= 0.7;
+
+    return width;
+}
+
+static float calculate_blur(BLRContext *blr, int w, int h,
+                            uint8_t* dir, int dir_linesize,
+                            uint8_t* dst, int dst_linesize,
+                            uint8_t* src, int src_linesize)
+{
+    float total_width = 0.0;
+    int block_count;
+    float block_total_width;
+
+    int i, j;
+    int blkcnt = 0;
+
+    float *blks = blr->blks;
+    int block_pct = blr->block_pct;
+    float block_pool_threshold = block_pct / 100.0;
+
+    int block_width  = blr->block_width;
+    int block_height = blr->block_height;
+    int brows = h / block_height;
+    int bcols = w / block_width;
+
+    for (int blkj = 0; blkj < brows; blkj++) {
+        for (int blki = 0; blki < bcols; blki++) {
+            block_total_width = 0;
+            block_count = 0;
+            for (int inj = 0; inj < block_height; inj++) {
+                for (int ini = 0; ini < block_width; ini++) {
+                    i = blki * block_width + ini;
+                    j = blkj * block_height + inj;
+
+                    if (dst[j * dst_linesize + i] > 0) {
+                        float width = edge_width(blr, i, j, dir[j*dir_linesize+i],
+                                                 w, h, dst[j*dst_linesize+i],
+                                                 src, src_linesize);
+                        if (width > 0.001) { // throw away zeros
+                            block_count++;
+                            block_total_width += width;
+                        }
+                    }
+                }
+            }
+            // if not enough edge pixels in a block, consider it smooth
+            if (block_total_width >= 2) {
+                blks[blkcnt] = block_total_width / block_count;
+                blkcnt++;
+            }
+        }
+    }
+
+    // simple block pooling by sorting and keeping the sharper blocks
+    qsort(blks, blkcnt, sizeof(*blks), comp);
+    blkcnt = ceil(blkcnt * block_pool_threshold);
+    for (int i = 0; i < blkcnt; i++) {
+        total_width += blks[i];
+    }
+
+    return  total_width / blkcnt;
+}
+
+static void set_meta(AVDictionary **metadata, const char *key, float d)
+{
+    char value[128];
+    snprintf(value, sizeof(value), "%f", d);
+    av_dict_set(metadata, key, value, 0);
+}
+
+static int blr_filter_frame(AVFilterLink *inlink, AVFrame *in)
+{
+    AVFilterContext *ctx = inlink->dst;
+    BLRContext *blr      = ctx->priv;
+
+    const int w = inlink->w;
+    const int h = inlink->h;
+
+    uint8_t *tmpbuf     = blr->tmpbuf;
+    uint8_t *filterbuf  = blr->filterbuf;
+    uint16_t *gradients = blr->gradients;
+    int8_t *directions  = blr->directions;
+
+    float blur;
+    AVDictionary **metadata;
+    metadata = &in->metadata;
+
+    // gaussian filter to reduce noise
+    gaussian_blur(w, h,
+                  filterbuf,  w,
+                  in->data[0], in->linesize[0]);
+
+    // compute the 16-bits gradients and directions for the next step
+    av_image_sobel(w, h, gradients, w, directions, w, filterbuf, w);
+
+    // non_maximum_suppression() will actually keep & clip what's necessary and
+    // ignore the rest, so we need a clean output buffer
+    memset(tmpbuf, 0, w * h);
+    av_image_non_maximum_suppression(w, h, tmpbuf, w, directions, w, gradients, w);
+
+
+    // keep high values, or low values surrounded by high values
+    av_image_double_threshold(blr->low_u8, blr->high_u8, w, h,
+                              tmpbuf, w, tmpbuf, w);
+
+
+    blur = calculate_blur(blr, w, h, directions, w,
+                          tmpbuf, w, filterbuf, w);
+
+    blr->blur_total += blur;
+
+    // write stats
+    av_log(ctx, AV_LOG_VERBOSE, "blur: %.7f\n", blur);
+
+    set_meta(metadata, "lavfi.blur", blur);
+
+    if (blr->stats_file) {
+        fprintf(blr->stats_file, "n:%"PRId64" ", blr->nb_frames);
+        fprintf(blr->stats_file, "blur:%.7f\n", blur);
+    }
+
+    blr->nb_frames++;
+
+    return ff_filter_frame(ctx->outputs[0], in);
+}
+
+static av_cold void blr_uninit(AVFilterContext *ctx)
+{
+    BLRContext *blr = ctx->priv;
+
+    if (blr->stats_file && blr->stats_file != stdout)
+        fclose(blr->stats_file);
+
+    if (blr->nb_frames > 0) {
+        int nb_frames = blr->nb_frames;
+        av_log(ctx, AV_LOG_INFO, "blur mean: %.7f\n",
+               blr->blur_total / nb_frames);
+    }
+
+    av_freep(&blr->tmpbuf);
+    av_freep(&blr->filterbuf);
+    av_freep(&blr->gradients);
+    av_freep(&blr->directions);
+    av_freep(&blr->blks);
+}
+
+static const enum AVPixelFormat pix_fmts[] = { AV_PIX_FMT_YUV420P, AV_PIX_FMT_GRAY8, AV_PIX_FMT_NONE };
+
+static const AVFilterPad blr_inputs[] = {
+    {
+        .name         = "default",
+        .type         = AVMEDIA_TYPE_VIDEO,
+        .config_props = blr_config_input,
+        .filter_frame = blr_filter_frame,
+    },
+};
+
+static const AVFilterPad blr_outputs[] = {
+    {
+        .name = "default",
+        .type = AVMEDIA_TYPE_VIDEO,
+    },
+};
+
+AVFilter ff_vf_blurriness = {
+    .name          = "blurriness",
+    .description   = NULL_IF_CONFIG_SMALL("Blurriness metric filter."),
+    .priv_size     = sizeof(BLRContext),
+    .init          = blr_init,
+    .uninit        = blr_uninit,
+    FILTER_PIXFMTS_ARRAY(pix_fmts),
+    FILTER_INPUTS(blr_inputs),
+    FILTER_OUTPUTS(blr_outputs),
+    .priv_class    = &blr_class,
+    .flags         = AVFILTER_FLAG_SUPPORT_TIMELINE_GENERIC,
+};
