Downsampling in gavl is broken

[Nobody/Anonymous]

Downsampling in gavl is broken - wrong number of factors per pixel is taken in gavl_video_scale_table_init().

When downscaling, the proper number of factors should be

  filter_support_size * (source_size / destination_size)

, i.e. for downsampling 1000*1000 -> 100*50 (width * height) and GAVL_SCALE_CUBIC_MITCHELL

  factors_per_pixel_horiz = 4 * (1000 / 100) = 40
  factors_per_pixel_vert  = 4 * (1000 / 50) = 80

In gavl_video_scale_table_init() the number of factors per pixel taken is always equal to filter support size - 4 for GAVL_SCALE_CUBIC_* - but this is right for upsampling only.

Target image, 1000x1000
http://ipicture.ru/uploads/080801/P5Sdz10UR4.png

Downsampling to 200x200 with gavl-1.0.1/src/scaletest, GAVL_SCALE_CUBIC_MITCHELL:
http://ipicture.ru/uploads/080801/lmPXewy7Uh.png

Downsampling to 200x200 with ImageMagick, (convert Target_1000x1000.png -filter Mitchell -resize 200x200 Scaled_200x200_IM_Mitchell.png)

http://ipicture.ru/uploads/080801/yUSs1VX9Lr.png

[Burkhard Plaum]

Thanks for the notice.

Does the same apply for sinc scaling?

I'll try to fix this.

[Nobody/Anonymous]

Yes, the same applies to all the scaling - from GAVL_SCALE_NEAREST to GAVL_SCALE_SINC_LANCZOS.

My quick-and-dirty fix :

scale_table.c:

void gavl_video_scale_table_init(gavl_video_scale_table_t * tab,
                                 gavl_video_options_t * opt,
                                 double src_off, double src_size,
                                 int dst_size, int src_width)
  {
  double t;
  int i, j, src_index_min, src_index_nearest;
  double src_index_f;
   
  gavl_video_scale_get_weight weight_func;
 
  double scale_factor;
 
  double filter_support;
  int src_taps;
  double scale_k;

  //  src_off = -0.25;
 
  scale_factor = (double)(dst_size) / src_size;
  scale_k = (scale_factor >= 1.0) ? 1.0 : scale_factor;
 
  /* Get the kernel generator */

  weight_func =
    gavl_video_scale_get_weight_func(opt, &(tab->factors_per_pixel));
 
  // +++
  filter_support = tab->factors_per_pixel;
  
  if (dst_size < src_size)  {
      src_taps = filter_support * src_size / dst_size + 0.5;
      if (!(src_taps & 1)) src_taps++;
      tab->factors_per_pixel = src_taps;
  }

  //  fprintf(stderr, "tab->factors_per_pixel: %d, src_width: %d\n",
  //          tab->factors_per_pixel, src_width);
 
  if(tab->factors_per_pixel > src_width)
    {
    switch(src_width)
      {
      case 1:
        opt->scale_mode = GAVL_SCALE_NEAREST;
        //        fprintf(stderr, "gavl: Changing scale mode to nearest (image too small)\n");
        break;
      case 2:
      case 3:
        opt->scale_mode = GAVL_SCALE_BILINEAR;
        //        fprintf(stderr, "gavl: Changing scale mode to bilinear (image too small)\n");
        break;
      default:
        opt->scale_mode = GAVL_SCALE_CUBIC_BSPLINE;
        //        fprintf(stderr, "gavl: Changing scale mode to bspline (image too small)\n");
        break;
      }
    weight_func =
      gavl_video_scale_get_weight_func(opt, &(tab->factors_per_pixel));
    }
 
  //  fprintf(stderr, "gavl_video_scale_table_init: %f %f %d %d\n",
  //          src_off, src_size, dst_size, tab->factors_per_pixel);
 
  /* (Re)allocate memory */

  alloc_table(tab, dst_size);

  for(i = 0; i < dst_size; i++)
    {

    /* src_index_f:       Fractional source index */
    /* src_index_nearest: Nearest integer source index */
   
    src_index_f = DST_TO_SRC((double)i);

    //    fprintf(stderr, "dst: %d -> src: %f (offset: s: %.2f)\n", i, src_index_f,
    //            src_off);
   
    //    if(src_index_f > src_size - 1.0)
    //      src_index_f = src_size - 1.0;
   
    src_index_nearest = ROUND(src_index_f);

    src_index_min = src_index_nearest - tab->factors_per_pixel/2;
   
    if(((double)src_index_nearest < src_index_f) && !(tab->factors_per_pixel % 2))
      {
      src_index_min++;
      }

    tab->pixels[i].index = src_index_min;

    //    fprintf(stderr, "src_index_f: %f, src_index_nearest: %d, src_index_min: %d, dst_index: %d\n",
    //            src_index_f, src_index_nearest, src_index_min, i);
   
    /* For nearest neighbour, we don't need any factors */
    if(tab->factors_per_pixel == 1)
      {
      if(tab->pixels[i].index < 0)
        tab->pixels[i].index = 0;
      if(tab->pixels[i].index > src_width - 1)
        tab->pixels[i].index = src_width - 1;
      continue;
      }
   
    /* Normalized distance of the destination pixel to the first source pixel
       in src coordinates */   
    t = src_index_f - src_index_min;
   
    for(j = 0; j < tab->factors_per_pixel; j++)
      {
      tab->pixels[i].factor_f[j] = weight_func(opt, t*scale_k);
      //      fprintf(stderr, "t: %f, w: %f\n", t, weight_func(opt, t));
      t -= 1.0;
      }
    }

  //  fprintf(stderr, "Before shift\n");
  //  gavl_video_scale_table_dump(tab);
 
  shift_borders(tab, src_width);
 
  // +++ // if(opt->scale_mode == GAVL_SCALE_SINC_LANCZOS)
    normalize_table(tab);
 
  //  fprintf(stderr, "After shift %d\n", src_width);
  //if(deinterlace || (total_fields == 2))
  //  gavl_video_scale_table_dump(tab);
 
  }

scale.c:

// can't use optimized routines for downsampling - number of factors is unknown
void gavl_init_scale_funcs_ex(gavl_scale_funcs_t * tab, gavl_video_options_t * opt,
                           int src_advance, int dst_advance, int shrinking)
  {
   
    if (shrinking) {
    gavl_init_scale_funcs_generic_c(tab);
    } else {
    gavl_init_scale_funcs_ex(tab, opt, src_advance, dst_advance)
        }
  }
   

scale_context.c:

int gavl_video_scale_context_init(gavl_video_scale_context_t*ctx,
                                  gavl_video_options_t * opt, int plane,
                                  const gavl_video_format_t * src_format,
                                  const gavl_video_format_t * dst_format,
                                  int src_field, int dst_field,
                                  int src_fields, int dst_fields)
  {
  int bits_h = 0, bits_v = 0;
  int sub_h_in = 1, sub_v_in = 1, sub_h_out = 1, sub_v_out = 1;
  int scale_x, scale_y;
  gavl_video_options_t tmp_opt, tmp_opt_y;
  double scale_factor_x, scale_factor_y;
  float src_chroma_offset_x, src_chroma_offset_y, dst_chroma_offset_x, dst_chroma_offset_y;
  gavl_rectangle_i_t src_rect_i;

 
  int src_width, src_height; /* Needed for generating the scale table */
  float offset_x, offset_y;
 
  int shrink_x, shrink_y;
 
  gavl_scale_funcs_t funcs;

  ctx->first_scanline = 0;

#if 0 
  fprintf(stderr, "scale_context_init: src_field: %d, dst_field: %d plane: %d\n",
          src_field, dst_field, plane);
#endif 
  gavl_rectangle_f_copy(&(ctx->src_rect), &(opt->src_rect));
  gavl_rectangle_i_copy(&(ctx->dst_rect), &(opt->dst_rect));
 
  scale_factor_x = (float)ctx->dst_rect.w / ctx->src_rect.w;
  scale_factor_y = (float)ctx->dst_rect.h / ctx->src_rect.h;
 
  shrink_x = (scale_factor_x < 1.0); shrink_y = (scale_factor_y < 1.0);
   
  ctx->plane = plane;

 
 
  if(plane)
    {
    /* Get chroma subsampling factors for source and destination */
    gavl_pixelformat_chroma_sub(src_format->pixelformat, &sub_h_in, &sub_v_in);
    gavl_pixelformat_chroma_sub(dst_format->pixelformat, &sub_h_out, &sub_v_out);
   
    ctx->src_rect.w /= sub_h_in;
    ctx->src_rect.x /= sub_h_in;

    ctx->src_rect.h /= sub_v_in;
    ctx->src_rect.y /= sub_v_in;
   
    ctx->dst_rect.w /= sub_h_out;
    ctx->dst_rect.x /= sub_h_out;

    ctx->dst_rect.h /= sub_v_out;
    ctx->dst_rect.y /= sub_v_out;

    src_width = src_format->image_width / sub_h_in;
    src_height = src_format->image_height / sub_v_in;

    }
  else
    {
    src_width = src_format->image_width;
    src_height = src_format->image_height;
    }
 
  if(src_fields == 2)
    {
    ctx->src_rect.h /= 2.0;
    ctx->src_rect.y /= 2.0;
    src_height /= 2;
    }
  if(dst_fields == 2)
    {
    ctx->dst_rect.h /= 2;
    ctx->dst_rect.y /= 2;
    }
 
#if 0
  fprintf(stderr, "gavl_video_scale_context_init\n");
  gavl_rectangle_f_dump(&(ctx->src_rect));
  fprintf(stderr, "\n");
  gavl_rectangle_i_dump(&(ctx->dst_rect));
  fprintf(stderr, "\n");
#endif
 

  /* Calculate chroma offsets  */
   
  gavl_video_format_get_chroma_offset(src_format, src_field, plane,
                                      &src_chroma_offset_x,
                                      &src_chroma_offset_y);
 
  gavl_video_format_get_chroma_offset(dst_format, dst_field, plane,
                                      &dst_chroma_offset_x,
                                      &dst_chroma_offset_y);
 
  offset_x = get_scale_offset(0, 0, 1, 1,
                              scale_factor_x, sub_h_in, sub_h_out,
                              src_chroma_offset_x, dst_chroma_offset_x);
  offset_y = get_scale_offset(src_field, dst_field, src_fields, dst_fields,
                              scale_factor_y, sub_v_in, sub_v_out,
                              src_chroma_offset_y, dst_chroma_offset_y);
 
  if((fabs(ctx->src_rect.w - ctx->dst_rect.w) > EPS) ||
     (fabs(ctx->src_rect.x) > EPS) || (fabs(offset_x) > EPS))
    scale_x = 1;
  else
    scale_x = 0;
 
  if((fabs(ctx->src_rect.h - ctx->dst_rect.h) > EPS) ||
     (fabs(ctx->src_rect.y) > EPS) || (fabs(offset_y) > EPS))
    scale_y = 1;
  else
    scale_y = 0;
#if 0
  fprintf(stderr, "Offsets: %f %f, scale_factors: %f %f\n",
          offset_x, offset_y, ctx->dst_rect.w / ctx->src_rect.w,
          ctx->dst_rect.h / ctx->src_rect.h);
#endif 
  offset_x += ctx->src_rect.x;
  offset_y += ctx->src_rect.y;

  ctx->func1 = NULL;
  ctx->func2 = NULL;
 
  ctx->num_directions = 0;
  if(scale_x)
    ctx->num_directions++;
  if(scale_y)
    ctx->num_directions++;

  /* Set source and destination frame planes */

  if(gavl_pixelformat_is_planar(src_format->pixelformat))
    ctx->src_frame_plane = plane;
  else
    ctx->src_frame_plane = 0;

  if(gavl_pixelformat_is_planar(dst_format->pixelformat))
    ctx->dst_frame_plane = plane;
  else
    ctx->dst_frame_plane = 0;

  /* Set bytes per line. We need this for copying */
  ctx->bytes_per_line = gavl_pixelformat_is_planar(src_format->pixelformat) ?
    ctx->dst_rect.w * gavl_pixelformat_bytes_per_component(src_format->pixelformat) :
    ctx->dst_rect.w * gavl_pixelformat_bytes_per_pixel(src_format->pixelformat);
 
  /* Set source and destination offsets */
 
  if(ctx->num_directions == 1)
    {
    get_offset_internal(src_format->pixelformat,
                plane, &ctx->offset1.src_advance, &ctx->offset1.src_offset);
    get_offset_internal(dst_format->pixelformat,
                plane, &ctx->offset1.dst_advance, &ctx->offset1.dst_offset);

    /* We set this once here */

    ctx->offset = &(ctx->offset1);
    ctx->dst_size = ctx->dst_rect.w;
    }
  else if(ctx->num_directions == 2)
    {
    get_offset_internal(src_format->pixelformat,
                plane, &ctx->offset1.src_advance, &ctx->offset1.src_offset);
    get_offset_internal(dst_format->pixelformat,
                plane, &ctx->offset2.dst_advance, &ctx->offset2.dst_offset);

    ctx->offset1.dst_offset = 0;
    ctx->offset1.dst_advance = ctx->offset1.src_advance;

    if((src_format->pixelformat == GAVL_YUY2) ||
       (src_format->pixelformat == GAVL_UYVY))
      {
      ctx->offset1.dst_advance = 1;
      }
    ctx->offset2.src_advance = ctx->offset1.dst_advance;
    ctx->offset2.src_offset  = ctx->offset1.dst_offset;
    }
  else if(!ctx->num_directions)
    {
    if((src_format->pixelformat == GAVL_YUY2) ||
       (src_format->pixelformat == GAVL_UYVY) ||
       (dst_format->pixelformat == GAVL_YUY2) ||
       (dst_format->pixelformat == GAVL_UYVY))
      ctx->func1 = copy_scanline_advance;
    else
      ctx->func1 = copy_scanline_noadvance;

    /* Set source and destination offsets */
    get_offset_internal(src_format->pixelformat,
                plane, &ctx->offset1.src_advance, &ctx->offset1.src_offset);
    get_offset_internal(dst_format->pixelformat,
                plane, &ctx->offset1.dst_advance, &ctx->offset1.dst_offset);

    /* We set this once here */

    ctx->offset = &(ctx->offset1);
    ctx->dst_size = ctx->dst_rect.w;

    ctx->num_directions = 1;
    return 1;
    }

  if(scale_x && scale_y)
    {
    // fprintf(stderr, "Initializing x table %f\n", ctx->src_rect.x + offset_x);
    gavl_video_options_copy(&tmp_opt, opt);
    gavl_video_scale_table_init(&(ctx->table_h), &tmp_opt,
                                offset_x,
                                ctx->src_rect.w, ctx->dst_rect.w, src_width);
    //    fprintf(stderr, "Initializing x table done\n");

    //    fprintf(stderr, "Initializing y table %f\n",
    //            ctx->src_rect.y + offset_y);
    gavl_video_options_copy(&tmp_opt_y, opt);
    gavl_video_scale_table_init(&(ctx->table_v), &tmp_opt_y,
                                offset_y,
                                ctx->src_rect.h, ctx->dst_rect.h, src_height);
   
    //    fprintf(stderr, "Initializing y table done\n");
   
    /* Check if we can scale in x and y-directions at once */

    if((tmp_opt.scale_mode == tmp_opt_y.scale_mode) &&
       (tmp_opt.scale_order == tmp_opt_y.scale_order) &&
       (shrink_x == shrink_y)
       )
      {
      memset(&funcs, 0, sizeof(funcs));
      gavl_init_scale_funcs_ex(&funcs, &tmp_opt, ctx->offset1.src_advance,
                            ctx->offset2.dst_advance, shrink_x);
      ctx->func1 = get_func(&(funcs.funcs_xy), src_format->pixelformat, &bits_h);
      //      fprintf(stderr, "X AND Y %d\n");
      }
   
    if(ctx->func1) /* Scaling routines for x-y are there, good */
      {
      ctx->num_directions = 1;
           
      gavl_video_scale_table_init_int(&(ctx->table_h), bits_h);
      /* Must be bits_h since we have only one function (and thus one accuracy) */
      gavl_video_scale_table_init_int(&(ctx->table_v), bits_h);
      ctx->offset = &(ctx->offset1);
      ctx->dst_size = ctx->dst_rect.w;
     
      /* Switch offsets back */
      ctx->offset1.dst_advance = ctx->offset2.dst_advance;
      ctx->offset1.dst_offset  = ctx->offset2.dst_offset;
      }
    else
      {
      gavl_video_scale_table_get_src_indices(&(ctx->table_h),
                                             &src_rect_i.x, &src_rect_i.w);
      gavl_video_scale_table_get_src_indices(&(ctx->table_v),
                                             &src_rect_i.y, &src_rect_i.h);

      /* Decide the scale order depending on whats computationally less expensive */
      /* We calculate the sizes (in pixels) of the temporary frame for both options and
         take the smaller one */
#if 0
      fprintf(stderr, "Src rect: ");
      gavl_rectangle_i_dump(&src_rect_i);
      fprintf(stderr, "\nDst rect: ");
      gavl_rectangle_i_dump(&ctx->dst_rect);
      fprintf(stderr, "\n");
#endif
      if(src_rect_i.h * ctx->dst_rect.w < ctx->dst_rect.h * src_rect_i.w)
        {
        //        fprintf(stderr, "X then Y\n");
        /* X then Y */

        ctx->buffer_width  = ctx->dst_rect.w;
        ctx->buffer_height = src_rect_i.h;
       
        gavl_video_scale_table_shift_indices(&(ctx->table_v), -src_rect_i.y);
        ctx->first_scanline = src_rect_i.y;
       
        memset(&funcs, 0, sizeof(funcs));
        gavl_init_scale_funcs_ex(&funcs, &tmp_opt,
                              ctx->offset1.src_advance,
                              ctx->offset1.dst_advance, shrink_x);
        ctx->func1 = get_func(&funcs.funcs_x, src_format->pixelformat, &bits_h);
       
        gavl_video_scale_table_init_int(&(ctx->table_h), bits_h);

        memset(&funcs, 0, sizeof(funcs));
        gavl_init_scale_funcs_ex(&funcs, &tmp_opt_y,
                              ctx->offset2.src_advance,
                              ctx->offset2.dst_advance, shrink_y);
        ctx->func2 = get_func(&funcs.funcs_y,
                              src_format->pixelformat, &bits_v);

        gavl_video_scale_table_init_int(&(ctx->table_v), bits_v);
        }
      else
        {
        //        fprintf(stderr, "Y then X\n");
        /* Y then X */

        ctx->buffer_width = src_rect_i.w;
        ctx->buffer_height = ctx->dst_rect.h;
       
        ctx->offset1.src_offset += src_rect_i.x * ctx->offset1.src_advance;
       
        gavl_video_scale_table_shift_indices(&(ctx->table_h), -src_rect_i.x);

        memset(&funcs, 0, sizeof(funcs));
        gavl_init_scale_funcs_ex(&funcs, &tmp_opt_y,
                              ctx->offset1.src_advance,
                              ctx->offset1.dst_advance, shrink_y);
        ctx->func1 = get_func(&funcs.funcs_y, src_format->pixelformat, &bits_v);

        gavl_video_scale_table_init_int(&(ctx->table_v), bits_v);

        memset(&funcs, 0, sizeof(funcs));
        gavl_init_scale_funcs_ex(&funcs, &tmp_opt,
                              ctx->offset2.src_advance,
                              ctx->offset2.dst_advance, shrink_x);
        ctx->func2 = get_func(&(funcs.funcs_x),
                              src_format->pixelformat, &bits_h);
       
        gavl_video_scale_table_init_int(&(ctx->table_h), bits_h);
        }
     
      /* Allocate temporary buffer */
      alloc_temp(ctx, src_format->pixelformat);
      }
    }
  else if(scale_x)
    {
    gavl_video_options_copy(&tmp_opt, opt);
    gavl_video_scale_table_init(&(ctx->table_h), &tmp_opt,
                                offset_x,
                                ctx->src_rect.w, ctx->dst_rect.w, src_width);
    //    fprintf(stderr, "Initializing x table done\n");

    memset(&funcs, 0, sizeof(funcs));
    gavl_init_scale_funcs_ex(&funcs, &tmp_opt,
                          ctx->offset1.src_advance,
                          ctx->offset1.dst_advance, shrink_x);
    ctx->func1 = get_func(&(funcs.funcs_x), src_format->pixelformat, &bits_h);

    gavl_video_scale_table_init_int(&(ctx->table_h), bits_h);
    }
  else if(scale_y)
    {
    //    fprintf(stderr, "Initializing y table\n");
    gavl_video_options_copy(&tmp_opt, opt);

    gavl_video_scale_table_init(&(ctx->table_v), &tmp_opt, offset_y,
                                ctx->src_rect.h, ctx->dst_rect.h, src_height);
    //    fprintf(stderr, "Initializing y table done\n");
    memset(&funcs, 0, sizeof(funcs));
    gavl_init_scale_funcs_ex(&funcs, &tmp_opt,
                          ctx->offset1.src_advance,
                          ctx->offset1.dst_advance, shrink_y);
    ctx->func1 = get_func(&(funcs.funcs_y), src_format->pixelformat, &bits_v);
   
    gavl_video_scale_table_init_int(&(ctx->table_v), bits_v);
    }

  if(!ctx->func1 || ((ctx->num_directions == 2) && !ctx->func2))
    {
    //    fprintf(stderr, "Initializing scale context failed: %d %p %p\n",
    //            ctx->num_directions,
    //            ctx->func1, ctx->func2);
    return 0;
    }
#if 0 
  /* Dump final scale tables */
  fprintf(stderr, "Horizontal table:\n");
  gavl_video_scale_table_dump(&(ctx->table_h));
  fprintf(stderr, "Vertical table:\n");
  gavl_video_scale_table_dump(&(ctx->table_v));
#endif

 
  get_minmax(src_format->pixelformat, ctx->min_values_h, ctx->max_values_h, ctx->min_values_f, ctx->max_values_f);
  get_minmax(src_format->pixelformat, ctx->min_values_v, ctx->max_values_v, ctx->min_values_f, ctx->max_values_f);

#if 0
  fprintf(stderr, "Min: %d %d %d, max: %d %d %d\n",
          ctx->min_values_h[0],
          ctx->min_values_h[1],
          ctx->min_values_h[2],
          ctx->max_values_h[0],
          ctx->max_values_h[1],
          ctx->max_values_h[2]);
#endif
#if 0
  for(i = 0; i < 4; i++)
    {
    ctx->min_values_h[i] <<= bits_h;
    ctx->max_values_h[i] <<= bits_h;

    ctx->min_values_v[i] <<= bits_v;
    ctx->max_values_v[i] <<= bits_v;
    }
#endif
#if 0
  fprintf(stderr, "Min: %d %d %d, max: %d %d %d\n",
          ctx->min_values_h[0],
          ctx->min_values_h[1],
          ctx->min_values_h[2],
          ctx->max_values_h[0],
          ctx->max_values_h[1],
          ctx->max_values_h[2]);
#endif
   
#if 0
  for(i = 0; i < 4; i++)
    {
    fprintf(stderr, "Channel %d: min: %08x, max: %08x\n", i,
            ctx->min_values[i], ctx->max_values[i]);
    }
#endif
  return 1;
  }

And gavl_video_scale_context_init_convolve() needs to be changed like gavl_video_scale_context_init()

[Burkhard Plaum]

I located the respective code imagemagick, and I'm thinking how
this can be cleanly done.

2 questions:

- For subsampled chroma planes, do the filter taps for the chroma planes
  have to be downsampled as well (currently they aren't)?

- Could you send me your changes as a unified diff?

[Nobody/Anonymous]

1. I think the filter taps don't have to be downsampled :

subsampled_chroma_taps = filter_support * (source_chroma_samples / dest_chroma_samples) = filter_support * ( (source_size / subsampling_factor) / (dest_size / subsampling_factor) ) = filter_support * (source_size / dest_size)

2.

diff -u M:\Soft\Linux\gavl-1.0.0\gavl\/scale.c E:\temp\gavl-1.0.0.fixed\/scale.c
--- M:\Soft\Linux\gavl-1.0.0\gavl\/scale.c    2008-03-11 16:08:46.000000000 +0300
+++ E:\temp\gavl-1.0.0.fixed\/scale.c    2008-08-06 15:47:25.794854400 +0400
@@ -256,6 +256,16 @@
     }
   }

+void gavl_init_scale_funcs_ex(gavl_scale_funcs_t * tab, gavl_video_options_t * opt,
+                           int src_advance, int dst_advance, int shrinking)
+{
+    if (shrinking) {
+    gavl_init_scale_funcs_generic_c(tab);
+    } else {
+    gavl_init_scale_funcs(tab, opt, src_advance, dst_advance);
+    }   
+}
+

int gavl_video_scaler_init(gavl_video_scaler_t * scaler,
                            const gavl_video_format_t * src_format,

diff -u M:\Soft\Linux\gavl-1.0.0\gavl\/scale_context.c E:\temp\gavl-1.0.0.fixed\/scale_context.c
--- M:\Soft\Linux\gavl-1.0.0\gavl\/scale_context.c    2008-03-11 14:50:22.000000000 +0300
+++ E:\temp\gavl-1.0.0.fixed\/scale_context.c    2008-07-29 12:47:01.000000000 +0400
@@ -571,6 +571,8 @@
   int src_width, src_height; /* Needed for generating the scale table */
   float offset_x, offset_y;
  
+  int shrink_x, shrink_y;
+ 
   gavl_scale_funcs_t funcs;

   ctx->first_scanline = 0;
@@ -584,6 +586,8 @@
  
   scale_factor_x = (float)ctx->dst_rect.w / ctx->src_rect.w;
   scale_factor_y = (float)ctx->dst_rect.h / ctx->src_rect.h;
+ 
+  shrink_x = (scale_factor_x < 1.0); shrink_y = (scale_factor_y < 1.0);
    
   ctx->plane = plane;

@@ -778,11 +782,13 @@
     /* Check if we can scale in x and y-directions at once */

     if((tmp_opt.scale_mode == tmp_opt_y.scale_mode) &&
-       (tmp_opt.scale_order == tmp_opt_y.scale_order))
+       (tmp_opt.scale_order == tmp_opt_y.scale_order) &&
+       (shrink_x == shrink_y)
+       )
       {
       memset(&funcs, 0, sizeof(funcs));
-      gavl_init_scale_funcs(&funcs, &tmp_opt, ctx->offset1.src_advance,
-                            ctx->offset2.dst_advance);
+      gavl_init_scale_funcs_ex(&funcs, &tmp_opt, ctx->offset1.src_advance,
+                            ctx->offset2.dst_advance, shrink_x);
       ctx->func1 = get_func(&(funcs.funcs_xy), src_format->pixelformat, &bits_h);
       //      fprintf(stderr, "X AND Y %d\n");
       }
@@ -830,17 +836,17 @@
         ctx->first_scanline = src_rect_i.y;
        
         memset(&funcs, 0, sizeof(funcs));
-        gavl_init_scale_funcs(&funcs, &tmp_opt,
+        gavl_init_scale_funcs_ex(&funcs, &tmp_opt,
                               ctx->offset1.src_advance,
-                              ctx->offset1.dst_advance);
+                              ctx->offset1.dst_advance, shrink_x);
         ctx->func1 = get_func(&funcs.funcs_x, src_format->pixelformat, &bits_h);
        
         gavl_video_scale_table_init_int(&(ctx->table_h), bits_h);

         memset(&funcs, 0, sizeof(funcs));
-        gavl_init_scale_funcs(&funcs, &tmp_opt_y,
+        gavl_init_scale_funcs_ex(&funcs, &tmp_opt_y,
                               ctx->offset2.src_advance,
-                              ctx->offset2.dst_advance);
+                              ctx->offset2.dst_advance, shrink_y);
         ctx->func2 = get_func(&funcs.funcs_y,
                               src_format->pixelformat, &bits_v);

@@ -859,17 +865,17 @@
         gavl_video_scale_table_shift_indices(&(ctx->table_h), -src_rect_i.x);

         memset(&funcs, 0, sizeof(funcs));
-        gavl_init_scale_funcs(&funcs, &tmp_opt_y,
+        gavl_init_scale_funcs_ex(&funcs, &tmp_opt_y,
                               ctx->offset1.src_advance,
-                              ctx->offset1.dst_advance);
+                              ctx->offset1.dst_advance, shrink_y);
         ctx->func1 = get_func(&funcs.funcs_y, src_format->pixelformat, &bits_v);

         gavl_video_scale_table_init_int(&(ctx->table_v), bits_v);

         memset(&funcs, 0, sizeof(funcs));
-        gavl_init_scale_funcs(&funcs, &tmp_opt,
+        gavl_init_scale_funcs_ex(&funcs, &tmp_opt,
                               ctx->offset2.src_advance,
-                              ctx->offset2.dst_advance);
+                              ctx->offset2.dst_advance, shrink_x);
         ctx->func2 = get_func(&(funcs.funcs_x),
                               src_format->pixelformat, &bits_h);
        
@@ -889,9 +895,9 @@
     //    fprintf(stderr, "Initializing x table done\n");

     memset(&funcs, 0, sizeof(funcs));
-    gavl_init_scale_funcs(&funcs, &tmp_opt,
+    gavl_init_scale_funcs_ex(&funcs, &tmp_opt,
                           ctx->offset1.src_advance,
-                          ctx->offset1.dst_advance);
+                          ctx->offset1.dst_advance, shrink_x);
     ctx->func1 = get_func(&(funcs.funcs_x), src_format->pixelformat, &bits_h);

@@ -906,9 +912,9 @@
                                 ctx->src_rect.h, ctx->dst_rect.h, src_height);
     //    fprintf(stderr, "Initializing y table done\n");
     memset(&funcs, 0, sizeof(funcs));
-    gavl_init_scale_funcs(&funcs, &tmp_opt,
+    gavl_init_scale_funcs_ex(&funcs, &tmp_opt,
                           ctx->offset1.src_advance,
-                          ctx->offset1.dst_advance);
+                          ctx->offset1.dst_advance, shrink_y);
     ctx->func1 = get_func(&(funcs.funcs_y), src_format->pixelformat, &bits_v);
    
     gavl_video_scale_table_init_int(&(ctx->table_v), bits_v);

diff -u M:\Soft\Linux\gavl-1.0.0\gavl\/scale_table.c E:\temp\gavl-1.0.0.fixed\/scale_table.c
--- M:\Soft\Linux\gavl-1.0.0\gavl\/scale_table.c    2008-01-23 14:24:16.000000000 +0300
+++ E:\temp\gavl-1.0.0.fixed\/scale_table.c    2008-08-06 15:42:13.956452800 +0400
@@ -88,13 +88,29 @@
   gavl_video_scale_get_weight weight_func;
  
   double scale_factor;
+ 
+  double filter_support;
+  int src_taps;
+  double scale_k;

   //  src_off = -0.25;
  
+  scale_factor = (double)(dst_size) / src_size;
+  scale_k = (scale_factor >= 1.0) ? 1.0 : scale_factor;
+ 
   /* Get the kernel generator */

   weight_func =
     gavl_video_scale_get_weight_func(opt, &(tab->factors_per_pixel));
+ 
+  // +++
+  filter_support = tab->factors_per_pixel;
+  
+  if (dst_size < src_size)  {
+      src_taps = filter_support * src_size / dst_size + 0.5;
+      if (!(src_taps & 1)) src_taps++;
+      tab->factors_per_pixel = src_taps;
+  }

   //  fprintf(stderr, "tab->factors_per_pixel: %d, src_width: %d\n",
   //          tab->factors_per_pixel, src_width);
@@ -128,8 +144,6 @@

   alloc_table(tab, dst_size);
 
-  scale_factor = (double)(dst_size) / src_size;
- 
   for(i = 0; i < dst_size; i++)
     {

@@ -174,7 +188,7 @@
    
     for(j = 0; j < tab->factors_per_pixel; j++)
       {
-      tab->pixels[i].factor_f[j] = weight_func(opt, t);
+      tab->pixels[i].factor_f[j] = weight_func(opt, t*scale_k);
       //      fprintf(stderr, "t: %f, w: %f\n", t, weight_func(opt, t));
       t -= 1.0;
       }
@@ -185,7 +199,7 @@
  
   shift_borders(tab, src_width);
  
-  if(opt->scale_mode == GAVL_SCALE_SINC_LANCZOS)
+  // +++ // if(opt->scale_mode == GAVL_SCALE_SINC_LANCZOS)
     normalize_table(tab);
  
   //  fprintf(stderr, "After shift %d\n", src_width);

[Burkhard Plaum]

Thanks a lot.

While researching a bit, I found that the method you suggest
(increasing filter taps according to downscaling factor) is
certainly better than it's done now, but mathematically not
100% correct.

What should be done, is to:

1. Apply a lowpass filter, whose cutoff frequency corresponds
   to the downsampling factor (that's missing right now, hence the
   aliasing effects).

2. Do the actual downsampling by some interpolation method.

If e.g. lowpass filtering is done by a Gaussian filter, interpolation
is done with a sinc filter, the overall filter coefficinents can be
obtained from convolution of the sinc coefficiens with the Gauss
coefficients.

This will always lead to more filter taps:

taps_total = taps_sinc + 2 * (taps_gauss - 1)

but this does not necessarily correspond to a factor of 1/shrink_factor.

What do you think?

[Burkhard Plaum]

taps_total = taps_sinc + 2 * (taps_gauss - 1)

must of course be

taps_total = taps_sinc + taps_gauss - 1

I implemented both methods, and they look almost the same. Mine is a bit better
for nearest neighbor scaling. There are 2 new functions:

gavl_video_options_set_downscale_filter();

and

gavl_video_options_set_downscale_blur();

which control the process.

Thanks for pointing this issue out. I always knew the lowpass before downsampling
was missing, but normal video images never have such high-frequency components, so
there were no visible aliasing effects.