[FFmpeg-devel,2/2] libavfilter/vf_colorconstancy.c : Adding weighted greyedge

Signed-off-by: Yatendra Singh <yatendra1999luffy@gmail.com>
---
 doc/filters.texi                |  36 ++++++
 libavfilter/Makefile            |   1 +
 libavfilter/allfilters.c        |   1 +
 libavfilter/vf_colorconstancy.c | 215 ++++++++++++++++++++++++++++++++
 4 files changed, 253 insertions(+)

On 6/18/20 3:52 PM, Yatendra Singh wrote:
> Signed-off-by: Yatendra Singh <yatendra1999luffy@gmail.com>
> ---
>   doc/filters.texi                |  36 ++++++
>   libavfilter/Makefile            |   1 +
>   libavfilter/allfilters.c        |   1 +
>   libavfilter/vf_colorconstancy.c | 215 ++++++++++++++++++++++++++++++++
>   4 files changed, 253 insertions(+)
>
> diff --git a/doc/filters.texi b/doc/filters.texi
> index 85a511b205..2946b5b9e6 100644
> --- a/doc/filters.texi
> +++ b/doc/filters.texi
> @@ -20250,6 +20250,42 @@ separatefields,select=eq(mod(n,4),0)+eq(mod(n,4),3),weave
>   @end example
>   @end itemize
>   
> +@section weighted_greyedge
> +Apply the color constancy filter which estimates illumination and updates the
> +image colors accordingly.
> +
> +It accepts the following options:
> +
> +@table @option
> +@item minknorm
> +The Minkowski parameter to be used for calculating the Minkowski distance. Must
> +be chosen in the range [0,20] and default value is 1. Set to 0 for getting
> +max value instead of calculating Minkowski distance.
> +
> +@item sigma
> +The standard deviation of Gaussian blur to be applied on the scene. Must be
> +chosen in the range [0,1024.0] and default value = 1. floor( @var{sigma} * break_off_sigma(3) )
> +can't be equal to 0 if @var{difford} is greater than 0.
> +
> +@item min_err
> +The error angle between the estimated illumination and white light at which the algorithm stops even
> +if it hasn't reached the required number of iterations @code{max_iters}. Must be chosen in the range
> +[0.02,PI] radians with default of 0.1.
> +
> +@item max_iters
> +The number of iterations at which the algorithm stops even if it hasn't reached the required
> +error angle between the estimated illumination and white light @code{min_err}. Must be chosen in the
> +range [1,100] with a default value of 10.
> +
> +@item kappa
> +The power which is applied to the spectral weights to change the impact of weights on illuminant
> +estimation @code{kappa}. Must be chosen in the range [1,25] with a default value of 10.
> +@end table
> +
> +@example
> +ffmpeg -i mondrian.tif -vf "weighted_greyedge=minknorm=0:sigma=1:max_iters=50:min_err=0.02:kappa=10" mondrian_out.tif
> +@end example
> +
>   @section xbr
>   Apply the xBR high-quality magnification filter which is designed for pixel
>   art. It follows a set of edge-detection rules, see
> diff --git a/libavfilter/Makefile b/libavfilter/Makefile
> index 994a4172a3..6973452f8d 100644
> --- a/libavfilter/Makefile
> +++ b/libavfilter/Makefile
> @@ -453,6 +453,7 @@ OBJS-$(CONFIG_VSTACK_FILTER)                 += vf_stack.o framesync.o
>   OBJS-$(CONFIG_W3FDIF_FILTER)                 += vf_w3fdif.o
>   OBJS-$(CONFIG_WAVEFORM_FILTER)               += vf_waveform.o
>   OBJS-$(CONFIG_WEAVE_FILTER)                  += vf_weave.o
> +OBJS-$(CONFIG_WEIGHTED_GREYEDGE_FILTER)      += vf_colorconstancy.o
>   OBJS-$(CONFIG_XBR_FILTER)                    += vf_xbr.o
>   OBJS-$(CONFIG_XFADE_FILTER)                  += vf_xfade.o
>   OBJS-$(CONFIG_XFADE_OPENCL_FILTER)           += vf_xfade_opencl.o opencl.o opencl/xfade.o
> diff --git a/libavfilter/allfilters.c b/libavfilter/allfilters.c
> index f2a44b0090..ad2e07f9c5 100644
> --- a/libavfilter/allfilters.c
> +++ b/libavfilter/allfilters.c
> @@ -432,6 +432,7 @@ extern AVFilter ff_vf_vstack;
>   extern AVFilter ff_vf_w3fdif;
>   extern AVFilter ff_vf_waveform;
>   extern AVFilter ff_vf_weave;
> +extern AVFilter ff_vf_weighted_greyedge;
>   extern AVFilter ff_vf_xbr;
>   extern AVFilter ff_vf_xfade;
>   extern AVFilter ff_vf_xfade_opencl;
> diff --git a/libavfilter/vf_colorconstancy.c b/libavfilter/vf_colorconstancy.c
> index d36400bd35..e2e32b7ca3 100644
> --- a/libavfilter/vf_colorconstancy.c
> +++ b/libavfilter/vf_colorconstancy.c
> @@ -1,5 +1,6 @@
>   /*
>    * Copyright (c) 2018 Mina Sami
> + * Copyright (c) 2020 Yatendra Singh
>    *
>    * This file is part of FFmpeg.
>    *
> @@ -26,6 +27,14 @@
>    *
>    * @cite
>    * J. van de Weijer, Th. Gevers, A. Gijsenij "Edge-Based Color Constancy".
> + *
> + * @cite
> + * J. van de Weijer, Th. Gevers, and J. Geusebroek,
> + * “Edge and corner detection by photometric quasi-invariants”.
> + *
> + * @cite
> + * A. Gijsenij, Th. Gevers, J. van de Weijer,
> + * "Improving Color Constancy by Photometric Edge Weighting".
>    */
>   
>   #include "libavutil/imgutils.h"
> @@ -40,8 +49,10 @@
>   #include <math.h>
>   
>   #define GREY_EDGE "greyedge"
> +#define WEIGHTED_GREY_EDGE "weighted_greyedge"
>   
>   #define SQRT3 1.73205080757
> +#define NORMAL_WHITE 1/SQRT3
>   
>   #define NUM_PLANES    3
>   #define MAX_DIFF_ORD  2
> @@ -77,12 +88,16 @@ typedef struct ColorConstancyContext {
>   
>       int difford;
>       int minknorm; /**< @minknorm = 0 : getMax instead */
> +    int kappa;
>       double sigma;
>   
>       int nb_threads;
>       int planeheight[4];
>       int planewidth[4];
>   
> +    double min_err;
> +    int max_iters;
> +
>       int filtersize;
>       double *gauss[MAX_DIFF_ORD+1];
>   
> @@ -608,6 +623,162 @@ static void chromatic_adaptation(AVFilterContext *ctx, AVFrame *in, AVFrame *out
>       ctx->internal->execute(ctx, diagonal_transformation, &td, NULL, nb_jobs);
>   }
>   
> +/**
> + * Slice function for weighted grey edge algorithm that does partial summing/maximizing
> + * of gaussian derivatives and applies the pixel wise weights.
> + *
> + * @param ctx the filter context.
> + * @param arg data to be passed between threads.
> + * @param jobnr current job nubmer.
> + * @param nb_jobs total number of jobs.
> + *
> + * @return 0.
> + */
> +static int filter_slice_weighted_greyedge(AVFilterContext* ctx, void *arg, int jobnr, int nb_jobs)
> +{
> +    ColorConstancyContext *s    = ctx->priv;
> +    ThreadData *td              = arg;
> +    AVFrame *in                 = td->in;
> +    int minknorm                = s->minknorm;
> +    const uint8_t thresh        = 255;
> +
> +    int height  = s->planeheight[0];
> +    int width   = s->planewidth[0];
> +
> +    int linesize    = in->linesize[0];
> +    int slice_start = (height * jobnr) / nb_jobs;
> +    int slice_end   = (height * (jobnr+1)) / nb_jobs;
> +
> +    for(int h = slice_start; h < slice_end; h++)
> +    {
> +        for(int w = 0; w < width; w++)
> +        {
> +            double spectral_weight, spvar, norm, norm_temp;
> +            double spvar_x, spvar_y;
> +            spvar_x = 0;
> +            spvar_y = 0;
> +            norm    = 0;
> +            for(int plane = 0; plane < NUM_PLANES; plane++)
> +            {
> +                spvar_x     += td->data[INDEX_DX][plane][INDX2D(h,w,width)];
> +                spvar_y     += td->data[INDEX_DY][plane][INDX2D(h,w,width)];
> +                norm_temp   = pow( td->data[INDEX_DX][plane][INDX2D(h,w,width)], 2);
> +                norm_temp   += pow( td->data[INDEX_DY][plane][INDX2D(h,w,width)], 2);
> +                norm        += sqrt( norm_temp);
> +            }
> +            spvar           = sqrt( pow( spvar_x * NORMAL_WHITE, 2) + pow( spvar_y * NORMAL_WHITE, 2) );
> +            norm            = sqrt(norm);
> +            spectral_weight = pow(spvar/norm, s->kappa);
> +
> +            for(int plane = 0; plane < NUM_PLANES; plane++)
> +            {
> +                const uint8_t *img_data = in->data[plane];
> +                const double *src       = td->data[INDEX_NORM][plane];
> +                double *dst             = td->data[INDEX_DST][plane];
> +
> +                if(!minknorm)
> +                {
> +                    spectral_weight = 1.0;
> +                    dst[jobnr] = FFMAX( dst[jobnr], fabs(src[INDX2D(h, w, width)]) * spectral_weight
> +                                        * (img_data[INDX2D(h, w, linesize)] < thresh) );
> +                }
> +                else
> +                {
> +                    dst[jobnr] += ( pow( fabs(src[INDX2D(h, w, width)] / 255.) * spectral_weight, minknorm)
> +                                    * (img_data[INDX2D(h, w, linesize)] < thresh) );
> +                }
> +            }
> +        }
> +    }
> +    return 0;
> +}
> +
> +/**
> + * Main driver function for weighted grey edge algorithm.
> + *
> + * @param ctx the filter context.
> + * @param in holds the input frame.
> + * @param out holds the output frame.
> + *
> + * AVERROR code if any error has occured.
> + */
> +static int filter_weighted_greyedge(AVFilterContext *ctx, AVFrame *in, AVFrame *out)
> +{
> +    ColorConstancyContext *s = ctx->priv;
> +    ThreadData td;
> +    int minknorm  = s->minknorm;
> +    int difford   = s->difford;
> +    double white[NUM_PLANES];
> +    int nb_jobs   = FFMIN3(s->planeheight[1], s->planewidth[1], s->nb_threads);
> +    int num_iters = 0;
> +    int plane, job, ret;
> +    double rep_err, dot;
> +    double W[NUM_PLANES];
> +
> +    td.in = in;
> +    ret = setup_derivative_buffers(ctx, &td);
> +    if (ret) {
> +        return ret;
> +    }
> +
> +    while( num_iters < s->max_iters )
> +    {
> +        get_derivative(ctx, &td);
> +        if (difford > 0) {
> +            ctx->internal->execute(ctx, slice_normalize, &td, NULL, nb_jobs);
> +        }
> +
> +        ctx->internal->execute(ctx, filter_slice_weighted_greyedge, &td, NULL, nb_jobs);
> +        if (!minknorm) {
> +            for (plane = 0; plane < NUM_PLANES; plane++) {
> +                white[plane] = 0; // All values are absolute
> +                for (job = 0; job < nb_jobs; job++) {
> +                    white[plane] = FFMAX(white[plane] , td.data[INDEX_DST][plane][job]);
> +                }
> +            }
> +        } else {
> +            for (plane = 0; plane < NUM_PLANES; ++plane) {
> +                white[plane] = 0;
> +                for (job = 0; job < nb_jobs; job++) {
> +                    white[plane] += td.data[INDEX_DST][plane][job];
> +                }
> +                white[plane] = pow(white[plane], 1./minknorm);
> +            }
> +        }
> +
> +        for( int i = 0; i < NUM_PLANES; i++ ){
> +            s->white[i] = s->white[i] * white[i];
> +        }
> +        normalize_light(s->white);
> +        chromatic_adaptation(ctx, in, out);
> +
> +        num_iters++;
> +
> +        rep_err = 0.0;
> +        dot = 0.0;
> +
> +        for(int i = 0; i < NUM_PLANES; i++)
> +        {
> +            W[i] = 1/s->white[i];
> +        }
> +        normalize_light(W);
> +
> +        for(int i = 0; i < NUM_PLANES; i++ )
> +        {
> +            W[i] *= s->white[i];
> +            dot += W[i];
> +        }
> +        rep_err = acosf(dot);
> +        if (rep_err <= s->min_err)
> +        {
> +            break;
> +        }
> +    }

Somehow cosmetic: maybe move angle error calculation to is own function?


> +
> +    cleanup_derivative_buffers(&td, difford + 1, NUM_PLANES);
> +    return 0;
> +}
> +
>   static int query_formats(AVFilterContext *ctx)
>   {
>       static const enum AVPixelFormat pix_fmts[] = {
> @@ -629,6 +800,8 @@ static int config_props(AVFilterLink *inlink)
>       double sigma = s->sigma;
>       int ret;
>   
> +    s->difford = 1;

Why are you forcing this value?


> +
>       if (!floor(break_off_sigma * sigma + 0.5) && s->difford) {
>           av_log(ctx, AV_LOG_ERROR, "floor(%f * sigma) must be > 0 when difford > 0.\n", break_off_sigma);
>           return AVERROR(EINVAL);
> @@ -681,6 +854,21 @@ static int filter_frame(AVFilterLink *inlink, AVFrame *in)
>           }
>           chromatic_adaptation(ctx, in, out);
>       }
> +    else if (!strcmp(ctx->filter->name, WEIGHTED_GREY_EDGE))
> +    {
> +        ColorConstancyContext *s = ctx->priv;
> +        for(int i = 0; i < NUM_PLANES; i++)
> +        {
> +            s->white[i] = NORMAL_WHITE;
> +        }
> +        chromatic_adaptation(ctx, in, out);

If you are setting initial estimation to normal white, does it make any 
difference to do a chromatic adaptation?


> +        ret = filter_weighted_greyedge(ctx, in, out);
> +        if (ret)
> +        {
> +            av_frame_free(&in);
> +            return ret;
> +        }
> +    }
>   
>       if (!direct)
>           av_frame_free(&in);
> @@ -741,3 +929,30 @@ AVFilter ff_vf_greyedge = {
>   };
>   
>   #endif /* CONFIG_GREY_EDGE_FILTER */
> +
> +#if CONFIG_WEIGHTED_GREYEDGE_FILTER
> +
> +static const AVOption weighted_greyedge_options[] = {
> +    { "minknorm",  "set Minkowski norm",         OFFSET(minknorm),  AV_OPT_TYPE_INT,    {.i64=1},   0,    20,     FLAGS },
> +    { "sigma",     "set sigma",                  OFFSET(sigma),     AV_OPT_TYPE_DOUBLE, {.dbl=1},   0.0,  1024.0, FLAGS },
> +    { "min_err",   "set minimum angular error",  OFFSET(min_err),   AV_OPT_TYPE_DOUBLE, {.dbl=0.1}, 0.02, M_PI,   FLAGS },
> +    { "max_iters", "set the maximum iterations", OFFSET(max_iters), AV_OPT_TYPE_INT,    {.i64=10},  1,    100,    FLAGS },
> +    { "kappa",     "set the kappa for weights",  OFFSET(kappa),     AV_OPT_TYPE_INT,    {.i64=10},  1,    25,    FLAGS },

Why 25?

Cosmetic: align FLAGS


> +    { NULL }
> +};
> +
> +AVFILTER_DEFINE_CLASS(weighted_greyedge);
> +
> +AVFilter ff_vf_weighted_greyedge = {
> +    .name          = WEIGHTED_GREY_EDGE,
> +    .description   = NULL_IF_CONFIG_SMALL("Estimates scene illumination by grey edge assumption."),
> +    .priv_size     = sizeof(ColorConstancyContext),
> +    .priv_class    = &weighted_greyedge_class,
> +    .query_formats = query_formats,
> +    .uninit        = uninit,
> +    .inputs        = colorconstancy_inputs,
> +    .outputs       = colorconstancy_outputs,
> +    .flags         = AVFILTER_FLAG_SUPPORT_TIMELINE_GENERIC | AVFILTER_FLAG_SLICE_THREADS,
> +};
> +
> +#endif /* CONFIG_WEIGHTED_GREY_EDGE_FILTER */

On Thu, Jun 18, 2020 at 19:22:56 +0530, Yatendra Singh wrote:
> +    for(int h = slice_start; h < slice_end; h++)
> +    {

Space between "for" and "(", and the opening bracket '{' goes on the
same line.

> +        for(int w = 0; w < width; w++)
> +        {

Same here, and elsewhere (also for "if" and "while").

> +                norm_temp   = pow( td->data[INDEX_DX][plane][INDX2D(h,w,width)], 2);
> +                norm_temp   += pow( td->data[INDEX_DY][plane][INDX2D(h,w,width)], 2);
> +                norm        += sqrt( norm_temp);

No space after the function calls' opening brackets.

> +    .description   = NULL_IF_CONFIG_SMALL("Estimates scene illumination by grey edge assumption."),

Imperative: "Estimate".

Moritz

>
> > +    s->difford = 1;
>
> Why are you forcing this value?
>
The algorithm only uses first order differentials for the calculation of
the specular variant and it does not seem to be an option to use other
orders as far as I have understood from the paper and the official code.

>
> If you are setting initial estimation to normal white, does it make any
> difference to do a chromatic adaptation?
>
 I was trying to reproduce the same values as the official code in matlab
and hence this was done.
 It might provide a very minor improvement in speed but having an estimate
beforehand might help converge faster.
 I would update the next patch with chromatic adaptation.

>
> > +
> > +static const AVOption weighted_greyedge_options[] = {
> > +    { "minknorm",  "set Minkowski norm",         OFFSET(minknorm),
> AV_OPT_TYPE_INT,    {.i64=1},   0,    20,     FLAGS },
> > +    { "sigma",     "set sigma",                  OFFSET(sigma),
>  AV_OPT_TYPE_DOUBLE, {.dbl=1},   0.0,  1024.0, FLAGS },
> > +    { "min_err",   "set minimum angular error",  OFFSET(min_err),
>  AV_OPT_TYPE_DOUBLE, {.dbl=0.1}, 0.02, M_PI,   FLAGS },
> > +    { "max_iters", "set the maximum iterations", OFFSET(max_iters),
> AV_OPT_TYPE_INT,    {.i64=10},  1,    100,    FLAGS },
> > +    { "kappa",     "set the kappa for weights",  OFFSET(kappa),
>  AV_OPT_TYPE_INT,    {.i64=10},  1,    25,    FLAGS },
>
> Why 25?
>
 Honestly, I have no idea as to what the actual upper limit on the power of
weights should be, so I used a random upper limit while testing.
 What would you suggest I should use?


I would also update the new patch with all the rest of the changes.

>
> If you are setting initial estimation to normal white, does it make any
>> difference to do a chromatic adaptation?
>>
>  I was trying to reproduce the same values as the official code in matlab
> and hence this was done.
>  It might provide a very minor improvement in speed but having an estimate
> beforehand might help converge faster.
>  I would update the next patch with chromatic adaptation.
>
 Sorry, I kind of misread this part while replying earlier.
 I was following the algorithm which was laid out in the paper which does
chromatic adaptation as the first step in the loop.
 So according to the algorithm it first performs correction based on the
initial estimate and then calculates the weights and new value of estimated
light.

>