diff --git a/doc/filters.texi b/doc/filters.texi
index 6a66b0ed11..09bf5f8d32 100644
--- a/doc/filters.texi
+++ b/doc/filters.texi
@@ -8941,6 +8941,30 @@ colorlevels=romin=0.5:gomin=0.5:bomin=0.5
 
 This filter supports the all above options as @ref{commands}.
 
+@section colormap
+
+Apply color maps with custom color mappings.
+
+The filter accepts the following options:
+
+@table @option
+@item source
+Set the list of source colors that will be mapped to target colors.
+
+@item target
+Set the list of target colors that will be mapped from source colors.
+
+@item kernel
+Set the kernel used to measure color differences between mapped colors.
+
+The accepted values are:
+@table @samp
+@item absproduct
+@item euclidean
+@item logarithmic
+@end table
+@end table
+
 @section colormatrix
 
 Convert color matrix.
diff --git a/libavfilter/Makefile b/libavfilter/Makefile
index 3f905bced4..d9865ef25e 100644
--- a/libavfilter/Makefile
+++ b/libavfilter/Makefile
@@ -219,6 +219,7 @@ OBJS-$(CONFIG_COLORKEY_OPENCL_FILTER)        += vf_colorkey_opencl.o opencl.o \
                                                 opencl/colorkey.o
 OBJS-$(CONFIG_COLORHOLD_FILTER)              += vf_colorkey.o
 OBJS-$(CONFIG_COLORLEVELS_FILTER)            += vf_colorlevels.o
+OBJS-$(CONFIG_COLORMAP_FILTER)               += vf_colormap.o
 OBJS-$(CONFIG_COLORMATRIX_FILTER)            += vf_colormatrix.o
 OBJS-$(CONFIG_COLORSPACE_FILTER)             += vf_colorspace.o colorspace.o colorspacedsp.o
 OBJS-$(CONFIG_COLORTEMPERATURE_FILTER)       += vf_colortemperature.o
diff --git a/libavfilter/allfilters.c b/libavfilter/allfilters.c
index b278423884..8574d62e85 100644
--- a/libavfilter/allfilters.c
+++ b/libavfilter/allfilters.c
@@ -205,6 +205,7 @@ extern const AVFilter ff_vf_colorkey;
 extern const AVFilter ff_vf_colorkey_opencl;
 extern const AVFilter ff_vf_colorhold;
 extern const AVFilter ff_vf_colorlevels;
+extern const AVFilter ff_vf_colormap;
 extern const AVFilter ff_vf_colormatrix;
 extern const AVFilter ff_vf_colorspace;
 extern const AVFilter ff_vf_colortemperature;
diff --git a/libavfilter/vf_colormap.c b/libavfilter/vf_colormap.c
new file mode 100644
index 0000000000..6fcd1e34b5
--- /dev/null
+++ b/libavfilter/vf_colormap.c
@@ -0,0 +1,510 @@
+/*
+ * Copyright (c) 2022 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
+ * Compute a look-up table from map of colors.
+ */
+
+#include "config_components.h"
+
+#include "libavutil/attributes.h"
+#include "libavutil/avstring.h"
+#include "libavutil/common.h"
+#include "libavutil/opt.h"
+#include "libavutil/pixdesc.h"
+#include "avfilter.h"
+#include "formats.h"
+#include "internal.h"
+#include "video.h"
+
+#define MAX_SIZE 64
+
+enum KernelType {
+    ABSPRODUCT,
+    EUCLIDEAN,
+    LOGARITHMIC,
+    NB_KERNELS,
+};
+
+typedef struct ColorMapContext {
+    const AVClass *class;
+
+    int nb_maps;
+    char *source_str;
+    char *target_str;
+
+    float source[MAX_SIZE][4];
+    float target[MAX_SIZE][4];
+    float icoeff[4][4];
+    float coeff[MAX_SIZE][4];
+
+    int kernel_type;
+    float (*kernel)(const float *x, const float *y);
+} ColorMapContext;
+
+#define OFFSET(x) offsetof(ColorMapContext, x)
+#define FLAGS AV_OPT_FLAG_FILTERING_PARAM|AV_OPT_FLAG_VIDEO_PARAM|AV_OPT_FLAG_RUNTIME_PARAM
+
+static const AVOption colormap_options[] = {
+    { "source", "set source map color values",   OFFSET(source_str),  AV_OPT_TYPE_STRING, { .str = "0" }, .flags = FLAGS },
+    { "target", "set target map color values",   OFFSET(target_str),  AV_OPT_TYPE_STRING, { .str = "0" }, .flags = FLAGS },
+    { "kernel", "set the kernel used for measuring color difference",  OFFSET(kernel_type), AV_OPT_TYPE_INT, {.i64=1}, 0, NB_KERNELS-1, FLAGS, "kernel" },
+    { "absproduct",  "absolute product of differences",            0, AV_OPT_TYPE_CONST, {.i64=ABSPRODUCT},  0, 0, FLAGS, "kernel" },
+    { "euclidean",   "square root of sum of squared differences",  0, AV_OPT_TYPE_CONST, {.i64=EUCLIDEAN},   0, 0, FLAGS, "kernel" },
+    { "logarithmic", "modified euclidean with logarithmic factor", 0, AV_OPT_TYPE_CONST, {.i64=LOGARITHMIC}, 0, 0, FLAGS, "kernel" },
+    { NULL }
+};
+
+static int parse_color_str(AVFilterContext *ctx, float colors[64][4], char *str, int *N)
+{
+    char *p, *arg, *saveptr;
+    int n = MAX_SIZE;
+
+    p = str;
+    for (int i = 0; i < MAX_SIZE; i++) {
+        float r, g, b;
+
+        if (!(arg = av_strtok(p, "|", &saveptr))) {
+            n = i;
+            break;
+        }
+
+        p = NULL;
+
+        if (av_sscanf(arg, "%f %f %f", &r, &g, &b) != 3) {
+            av_log(ctx, AV_LOG_ERROR, "Invalid syntax for color.\n");
+            return AVERROR(EINVAL);
+        }
+
+        colors[i][0] = r;
+        colors[i][1] = g;
+        colors[i][2] = b;
+    }
+
+    *N = n;
+
+    return 0;
+}
+
+static int gauss_make_triangular(double *A, int *p, int n)
+{
+    p[n - 1] = n - 1;
+    for (int k = 0; k < n; k++) {
+        double t1;
+        int m = k;
+
+        for (int i = k + 1; i < n; i++)
+            if (fabs(A[k + n * i]) > fabs(A[k + n * m]))
+                m = i;
+        p[k] = m;
+        t1 = A[k + n * m];
+        A[k + n * m] = A[k + n * k];
+        A[k + n * k] = t1;
+        if (t1 != 0) {
+            for (int i = k + 1; i < n; i++)
+                A[k + n * i] /= -t1;
+            if (k != m)
+                for (int i = k + 1; i < n; i++) {
+                    double t2 = A[i + n * m];
+                    A[i + n * m] = A[i + n * k];
+                    A[i + n * k] = t2;
+                }
+            for (int j = k + 1; j < n; j++)
+                for (int i = k + 1; i < n; i++)
+                    A[i + n * j] += A[k + j * n] * A[i + k * n];
+        } else {
+            return 0;
+        }
+    }
+
+    return 1;
+}
+
+static void gauss_solve_triangular(const double *A, const int *p, double *b, int n)
+{
+    for(int k = 0; k < n - 1; k++) {
+        int m = p[k];
+        double t = b[m];
+        b[m] = b[k];
+        b[k] = t;
+        for (int i = k + 1; i < n; i++)
+            b[i] += A[k + n * i] * t;
+    }
+
+    for(int k = n - 1; k > 0; k--) {
+        b[k] /= A[k + n * k];
+        double t = b[k];
+        for (int i = 0; i < k; i++)
+            b[i] -= A[k + n * i] * t;
+    }
+
+    b[0] /= A[0 + 0 * n];
+}
+
+static int gauss_solve(double *A, double *b, int n)
+{
+    int *p = av_calloc(n, sizeof(*p));
+
+    if (!p)
+        return 1;
+
+    if (!gauss_make_triangular(A, p, n)) {
+        av_freep(&p);
+        return 1;
+    }
+
+    gauss_solve_triangular(A, p, b, n);
+
+    av_freep(&p);
+
+    return 0;
+}
+
+#define P2(x) ((x)*(x))
+
+static float absproduct(const float *x, const float *y)
+{
+    const float d = (x[0]-y[0]) *
+                    (x[1]-y[1]) *
+                    (x[2]-y[2]);
+    return fabsf(d);
+}
+
+static float euclidean_kernel(const float *x, const float *y)
+{
+    const float d2 = P2(x[0]-y[0]) +
+                     P2(x[1]-y[1]) +
+                     P2(x[2]-y[2]);
+    return sqrtf(d2);
+}
+
+static float logarithmic_kernel(const float *x, const float *y)
+{
+    const float d2 = P2(x[0]-y[0]) +
+                     P2(x[1]-y[1]) +
+                     P2(x[2]-y[2]);
+    return sqrtf(d2) * log(fmaxf(1e-8f, d2));
+}
+
+static int config_props(AVFilterLink *inlink)
+{
+    AVFilterContext *ctx = inlink->dst;
+    ColorMapContext *s = ctx->priv;
+    int ret, n = 0;
+
+    switch (s->kernel_type) {
+    case ABSPRODUCT:
+        s->kernel = absproduct;
+        break;
+    case EUCLIDEAN:
+        s->kernel = euclidean_kernel;
+        break;
+    case LOGARITHMIC:
+        s->kernel = logarithmic_kernel;
+        break;
+    default:
+        return AVERROR_BUG;
+    }
+
+    ret = parse_color_str(ctx, s->source, s->source_str, &n);
+    if (ret < 0)
+        return ret;
+
+    ret = parse_color_str(ctx, s->target, s->target_str, &n);
+    if (ret < 0)
+        return ret;
+
+    s->nb_maps = n;
+
+    for (int c = 0; c < 3; c++) {
+        for (int j = 0; j < s->nb_maps; j++)
+            s->coeff[j][c] = 0.f;
+
+        for (int j = 0; j < 4; j++) {
+            s->icoeff[j][c] = 0;
+            s->icoeff[j][c] = 0;
+            s->icoeff[j][c] = 0;
+        }
+
+        s->icoeff[c+1][c] = 1.f;
+
+        switch (n) {
+        case 1:
+            {
+                float div = fabsf(s->source[0][c]) < 1e-6f ? 1e-6f : s->source[0][c];
+                s->icoeff[c][1+c] = s->target[0][c] / div;
+            }
+            break;
+        case 2:
+            {
+                double A[2 * 2] = { 1, s->source[0][c],
+                                    1, s->source[1][c] };
+                double b[2] = { s->target[0][c], s->target[1][c] };
+
+                if (gauss_solve(A, b, 2))
+                    continue;
+
+                s->icoeff[0  ][c] = b[0];
+                s->icoeff[1+c][c] = b[1];
+            }
+            break;
+        case 3:
+            {
+                const uint8_t idx[3][3] = {{ 0, 1, 2 },
+                                           { 1, 0, 2 },
+                                           { 2, 0, 1 }};
+                const uint8_t didx[3][4] = {{ 0, 1, 2, 2 },
+                                            { 0, 2, 1, 2 },
+                                            { 0, 2, 2, 1 }};
+                const int C0 = idx[c][0];
+                const int C1 = idx[c][1];
+                const int C2 = idx[c][2];
+                double A[3 * 3] = { 1, s->source[0][C0], s->source[0][C1] + s->source[0][C2],
+                                    1, s->source[1][C0], s->source[1][C1] + s->source[1][C2],
+                                    1, s->source[2][C0], s->source[2][C1] + s->source[2][C2] };
+                double b[3] = { s->target[0][c], s->target[1][c], s->target[2][c] };
+
+                if (gauss_solve(A, b, 3))
+                    continue;
+
+                s->icoeff[0][c] = b[didx[c][0]];
+                s->icoeff[1][c] = b[didx[c][1]];
+                s->icoeff[2][c] = b[didx[c][2]];
+                s->icoeff[3][c] = b[didx[c][3]];
+            }
+            break;
+        case 4:
+            {
+                double A[4 * 4] = { 1, s->source[0][0], s->source[0][1], s->source[0][2],
+                                    1, s->source[1][0], s->source[1][1], s->source[1][2],
+                                    1, s->source[2][0], s->source[2][1], s->source[2][2],
+                                    1, s->source[3][0], s->source[3][1], s->source[3][2] };
+                double b[4] = { s->target[0][c], s->target[1][c], s->target[2][c], s->target[3][c] };
+                int pivot[4];
+
+                if (!gauss_make_triangular(A, pivot, 4))
+                    continue;
+                gauss_solve_triangular(A, pivot, b, 4);
+
+                s->icoeff[0][c] = b[0];
+                s->icoeff[1][c] = b[1];
+                s->icoeff[2][c] = b[2];
+                s->icoeff[3][c] = b[3];
+            }
+            break;
+        default:
+            {
+                const int N = s->nb_maps;
+                const int N4 = N + 4;
+                double *A = av_calloc(sizeof(*A), N4 * N4);
+                double *b = av_calloc(sizeof(*b), N4);
+                int *pivot = NULL;
+
+                if (!A || !b)
+                    goto error;
+
+                for (int j = 0; j < N; j++)
+                    for (int i = j; i < N; i++)
+                        A[j*N4+i] = A[i*N4+j] = s->kernel(s->source[i], s->source[j]);
+
+                for (int i = 0; i < N; i++)
+                    A[i*N4+N+0] = A[(N+0)*N4+i] = 1;
+                for (int i = 0; i < N; i++)
+                    A[i*N4+N+1] = A[(N+1)*N4+i] = s->source[i][0];
+                for (int i = 0; i < N; i++)
+                    A[i*N4+N+2] = A[(N+2)*N4+i] = s->source[i][1];
+                for (int i = 0; i < N; i++)
+                    A[i*N4+N+3] = A[(N+3)*N4+i] = s->source[i][2];
+
+                for (int j = N; j < N4; j++)
+                    for (int i = N;i < N4; i++)
+                        A[j * N4 + i] = 0.;
+
+                pivot = av_calloc(N4, sizeof(*pivot));
+                if (!pivot)
+                    goto error;
+
+                if (gauss_make_triangular(A, pivot, N4)) {
+                    for (int i = 0;i < N; i++)
+                        b[i] = s->target[i][c];
+                    for (int i = N; i < N + 4; i++)
+                        b[i] = 0;
+
+                    gauss_solve_triangular(A, pivot, b, N4);
+
+                    for (int i = 0; i < N; i++)
+                        s->coeff[i][c] = b[i];
+
+                    for (int i = 0; i < 4; i++)
+                        s->icoeff[i][c] = b[N + i];
+                }
+error:
+                av_free(pivot);
+                av_free(b);
+                av_free(A);
+            }
+        }
+    }
+
+    return 0;
+}
+
+typedef struct ThreadData {
+    AVFrame *in, *out;
+} ThreadData;
+
+static int colormap_slice(AVFilterContext *ctx, void *arg, int jobnr, int nb_jobs)
+{
+    ColorMapContext *s = ctx->priv;
+    ThreadData *td = arg;
+    AVFrame *in = td->in;
+    AVFrame *out = td->out;
+    const int maps = s->nb_maps;
+    const int width = out->width;
+    const int height = out->height;
+    const int slice_start = (height * jobnr) / nb_jobs;
+    const int slice_end = (height * (jobnr + 1)) / nb_jobs;
+    const int sr_linesize = in->linesize[2] / 4;
+    const int dr_linesize = out->linesize[2] / 4;
+    const int sg_linesize = in->linesize[0] / 4;
+    const int dg_linesize = out->linesize[0] / 4;
+    const int sb_linesize = in->linesize[1] / 4;
+    const int db_linesize = out->linesize[1] / 4;
+    const float *sr = (float *)in->data[2] + slice_start * sr_linesize;
+    const float *sg = (float *)in->data[0] + slice_start * sg_linesize;
+    const float *sb = (float *)in->data[1] + slice_start * sb_linesize;
+    float *r = (float *)out->data[2] + slice_start * dr_linesize;
+    float *g = (float *)out->data[0] + slice_start * dg_linesize;
+    float *b = (float *)out->data[1] + slice_start * db_linesize;
+
+    for (int y = slice_start; y < slice_end; y++) {
+        for (int x = 0; x < width; x++) {
+            const float input[3] = { sr[x], sg[x], sb[x] };
+            float srv, sgv, sbv;
+            float rv, gv, bv;
+
+            srv = sr[x];
+            sgv = sg[x];
+            sbv = sb[x];
+
+            rv = s->icoeff[0][0];
+            gv = s->icoeff[0][1];
+            bv = s->icoeff[0][2];
+
+            rv += s->icoeff[1][0] * srv + s->icoeff[2][0] * sgv + s->icoeff[3][0] * sbv;
+            gv += s->icoeff[1][1] * srv + s->icoeff[2][1] * sgv + s->icoeff[3][1] * sbv;
+            bv += s->icoeff[1][2] * srv + s->icoeff[2][2] * sgv + s->icoeff[3][2] * sbv;
+
+            for (int z = 0; z < maps && maps > 4; z++) {
+                const float cr = s->coeff[z][0];
+                const float cg = s->coeff[z][1];
+                const float cb = s->coeff[z][2];
+                const float f = s->kernel(input, s->source[z]);
+
+                rv += f * cr;
+                gv += f * cg;
+                bv += f * cb;
+            }
+
+            r[x] = rv;
+            g[x] = gv;
+            b[x] = bv;
+        }
+
+        sg += sg_linesize;
+        g += dg_linesize;
+        sb += sb_linesize;
+        b += db_linesize;
+        sr += sr_linesize;
+        r += dr_linesize;
+    }
+
+    return 0;
+}
+
+static int filter_frame(AVFilterLink *inlink, AVFrame *in)
+{
+    AVFilterContext *ctx = inlink->dst;
+    AVFilterLink *outlink = ctx->outputs[0];
+    ThreadData td;
+    AVFrame *out;
+
+    if (av_frame_is_writable(in)) {
+        out = in;
+    } else {
+        out = ff_get_video_buffer(outlink, outlink->w, outlink->h);
+        if (!out) {
+            av_frame_free(&in);
+            return AVERROR(ENOMEM);
+        }
+        av_frame_copy_props(out, in);
+    }
+
+    td.in = in;
+    td.out = out;
+    ff_filter_execute(ctx, colormap_slice, &td, NULL,
+                      FFMIN(in->height, ff_filter_get_nb_threads(ctx)));
+
+    if (out != in)
+        av_frame_free(&in);
+
+    return ff_filter_frame(outlink, out);
+}
+
+static int process_command(AVFilterContext *ctx, const char *cmd, const char *args,
+                           char *res, int res_len, int flags)
+{
+    int ret = ff_filter_process_command(ctx, cmd, args, res, res_len, flags);
+
+    if (ret < 0)
+        return ret;
+
+    return config_props(ctx->inputs[0]);
+}
+
+static const AVFilterPad inputs[] = {
+    { .name         = "default",
+      .type         = AVMEDIA_TYPE_VIDEO,
+      .filter_frame = filter_frame,
+      .config_props = config_props,
+    },
+};
+
+static const AVFilterPad outputs[] = {
+    { .name = "default",
+      .type = AVMEDIA_TYPE_VIDEO,
+    },
+};
+
+AVFILTER_DEFINE_CLASS(colormap);
+
+const AVFilter ff_vf_colormap = {
+    .name          = "colormap",
+    .description   = NULL_IF_CONFIG_SMALL("Apply custom Color Maps to video stream."),
+    .priv_class    = &colormap_class,
+    .priv_size     = sizeof(ColorMapContext),
+    FILTER_INPUTS(inputs),
+    FILTER_OUTPUTS(outputs),
+    FILTER_PIXFMTS(AV_PIX_FMT_GBRPF32, AV_PIX_FMT_GBRAPF32),
+    .flags         = AVFILTER_FLAG_SUPPORT_TIMELINE_GENERIC |
+                     AVFILTER_FLAG_SLICE_THREADS,
+    .process_command = process_command,
+};
