On Wed, 2006-10-11 at 21:33 +0200, Geert Jordaens wrote:

In order not to rewrite interpolation fnctionality for every operation I would like to propose following API.

GeglInterpolator +- GeglInterpolatorLinear
|
+ GeglInterpolatorCubic
|
+ GeglInterpolatorLanczos

usage :

/* first create the interpolator */ GeglInterpolatorLinear *interpolator = g_object_new (GEGL_TYPE_INTERPLATOR_LINEAR,
"source", buffer, "format", babl_format ("RaGaBaA float"),
NULL);

gegl_interpolator_prepare(interpolator); /* prepares a linear buffer in the format suitable for the interpolation routine. */

gegl_interpolator_get(interpolator, x, y, dst); /* returns a interpolated pixel in dst */

sllooooooooooow. I've encountered this oh-so-wonderful abstraction before (JAI uses it, among others). Making a method call in the inner loop of your pixel calculation routine is a big-nono. The first optimization you will make is to remove that. Typically, the common case will not use the abstraction, and an operator will be provided that can accept the general interpolation object, but the operator will be so slow no one will use it. It is wonderful when developing your foundations, to be sure, but it will eventually be ripped out as soon as someone tries to speed it up.

If you can provide versions that accept arrays as input, then the general case can be used everywhere. Actually that can be made a general rule. Anything that manipulates pixels much accept arrays as input (possibly masks as well) or they will be ripped out at the earliest opportunity.

On 10/12/06, Daniel Rogers wrote:

sllooooooooooow. I've encountered this oh-so-wonderful abstraction before (JAI uses it, among others). Making a method call in the inner loop of your pixel calculation routine is a big-nono. The first optimization you will make is to remove that. Typically, the common case will not use the abstraction, and an operator will be provided that can accept the general interpolation object, but the operator will be so slow no one will use it. It is wonderful when developing your foundations, to be sure, but it will eventually be ripped out as soon as someone tries to speed it up.

Any such optimizations are premature optimizations at the moment. At the current stage making the source code of the included operations general as well as easy to read, understand and debug is more important than efficient code.

I guess the overhead of a function call is smaller in C than in Java, it might also be possible to inline the actual interpolating sampling functions.

The operation that geert is currently working on is a displacement map, what you are suggesting would lead to writing a separate displacement map operation for each interpolation method. This would lead to very unwieldy code.

If you can provide versions that accept arrays as input, then the general case can be used everywhere.

The proposed solution is the closest you will get to operate directly on a linear buffer.

Actually that can be made a
general rule. Anything that manipulates pixels much accept arrays as input (possibly masks as well) or they will be ripped out at the earliest opportunity.

My guess is that the proposed solution is a good compromise, at a later stage more optimized version can be added. And then I actually mean adding code, not replacing the current code. The reference version of the base operation set int floating point should remain readable code, and optimized version of 8bit/16bit/floating point that are regression tested against the reference version.

On Oct 12, 2006, at 2:03 AM, Øyvind Kolås wrote:

On 10/12/06, Daniel Rogers wrote:

sllooooooooooow. I've encountered this oh-so-wonderful abstraction before (JAI uses it, among others). Making a method call in the inner loop of your pixel calculation routine is a big-nono. The first optimization you will make is to remove that. Typically, the common case will not use the abstraction, and an operator will be provided that
can accept the general interpolation object, but the operator will be so
slow no one will use it. It is wonderful when developing your foundations, to be sure, but it will eventually be ripped out as soon as
someone tries to speed it up.

Any such optimizations are premature optimizations at the moment. At the
current stage making the source code of the included operations general as well
as easy to read, understand and debug is more important than efficient code.

I guess the overhead of a function call is smaller in C than in Java, it might
also be possible to inline the actual interpolating sampling functions.

In this particular case it will not be possible to inline. Since the relevant GObject function must be virtual it won't be able to be inlined. You can actually manage to inline it anyway in a more managable language, since you can declare a class final and then specify the final class as the concrete type.

The operation that geert is currently working on is a displacement map, what
you are suggesting would lead to writing a separate displacement map operation
for each interpolation method. This would lead to very unwieldy code.

Yeah, the IP libraries I've worked with all do something similar, unfortunately.

If you can provide versions that accept arrays as input, then the general case can be used everywhere.

The proposed solution is the closest you will get to operate directly on a linear buffer.

I'm not going to be pushy about this. Just declare my experience.

On 10/12/06, Daniel Rogers wrote:

On Oct 12, 2006, at 2:03 AM, Øyvind Kolås wrote:

On 10/12/06, Daniel Rogers wrote:

sllooooooooooow. I've encountered this oh-so-wonderful abstraction before (JAI uses it, among others). Making a method call in the inner loop of your pixel calculation routine is a big-nono. The first optimization you will make is to remove that. Typically, the common

Any such optimizations are premature optimizations at the moment. At the
current stage making the source code of the included operations general as well
as easy to read, understand and debug is more important than efficient code.
The operation that geert is currently working on is a displacement map, what
you are suggesting would lead to writing a separate displacement map operation
for each interpolation method. This would lead to very unwieldy code.

Yeah, the IP libraries I've worked with all do something similar, unfortunately.

Something similar to what? Having a proxy object for doing interpolation or implement separate displacement etc. operation per interpolator?

If you mean the first, I do not find readable code unfortunate, a problem with the current GIMP code base IMHO is that some algorithms are optimized beyond recognition.

/Øyvind K.

On Fri, 2006-10-13 at 16:29 +0200, Øyvind Kolås wrote:

On 10/12/06, Daniel Rogers wrote:

On Oct 12, 2006, at 2:03 AM, Øyvind Kolås wrote:

On 10/12/06, Daniel Rogers wrote:

sllooooooooooow. I've encountered this oh-so-wonderful

abstraction

before (JAI uses it, among others). Making a method call in the

inner

loop of your pixel calculation routine is a big-nono. The first optimization you will make is to remove that. Typically, the

common

Any such optimizations are premature optimizations at the moment. At the
current stage making the source code of the included operations general as well
as easy to read, understand and debug is more important than efficient code.
The operation that geert is currently working on is a displacement map, what
you are suggesting would lead to writing a separate displacement map operation
for each interpolation method. This would lead to very unwieldy

code.

Yeah, the IP libraries I've worked with all do something similar, unfortunately.

Something similar to what? Having a proxy object for doing

interpolation or

implement separate displacement etc. operation per interpolator?

If you mean the first, I do not find readable code unfortunate, a problem with the current GIMP code base IMHO is that some algorithms are optimized beyond recognition.

No, I meant the later. I am all in favor of readable code, actually. I like the particular interface we are discussing. Too slow is still too slow, however. My experience has shown that when performance matters, that abstraction is the first to go. I'm not saying it's a bad idea, I am just bringing up what I've found to be a thorny issue.

On 10/15/06, Geert Jordaens wrote:

I've added the implementation to bug report : Bug 360888 – Interpolation

snip

usage :

snip

Maybe something along the following lines would be a better API?

/Øyvind K.

void displacementmap (GeglBuffer *src,
GeglBuffer *aux,
GeglBuffer *dst,
gdouble scale,
gint cx,
gint cy)
{
gint i, j;
gint dst_pos=0;
gfloat *buf = g_malloc0 (4 * 4 * 4); gfloat *aux_buf = g_malloc0 (aux->width * aux->height * 4 * 4); gfloat *dst_buf = g_malloc0 (dst->width * dst->height * 4 * 4); gegl_buffer_get_fmt (aux, aux_buf, babl_format ("RGBA float")); gegl_buffer_get_fmt (dst, dst_buf, babl_format ("RaGaBaA float"));

*GeglSampler *sampler = g_object_new (GEGL_TYPE_SAMPLER_CUBIC, "input", src,
"format", babl_format ("RaGaBaA float"), "x", 50,
"y", 50,
"width", 200, /* or a rect could be passed, or a subbuffer of "height", 200, buffer expected to be created before making a sampler */
NULL);

for (i = 0; i < src->height; i++) for (j = 0; j < src->width; j++) {
gdouble dx, dy, xc, yc;

/* Calculate x/y - offset relative to origin */ dx = (gdouble) j / (src->width-1.0); dy = (gdouble) i / (src->height-1.0); xc = j + (gint) (displacement (aux, aux_buf, dx, dy, cx) * scale); yc = i + (gint) (displacement (aux, aux_buf, dx, dy, cy) * scale);

/* returns a freshly sampled pixel in dst */ gegl_sampler_get (sampler, xc, yc, buf);

#if 0 /* using this form, might make more sense, since it would allow accurate sampling as needed by affine/perspective transforms that are scaling down as well as when scaling up/interpolating */
gegl_sampler_get_quad (sampler, xc-0.5, yc-0.5, xc+0.5, yc-0.5, xc+0.5, yc+0.5, xc-0.5, yc+0.5, buf); #endif

{
int component;
for (component = 0 ; component < 4 ; component++) dst_buf[dst_pos+component]= buf[component]; }
dst_pos+=4;
}

gegl_buffer_set_fmt (dst, dst_buf, babl_format ("RaGaBaA float")); g_object_unref (sampler);
g_free (aux_buf);
g_free (dst_buf);
}

On 10/16/06, Øyvind Kolås wrote:

#if 0
/* using this form, might make more sense, since it would allow accurate sampling as needed by affine/perspective transforms that are scaling down as well as when scaling up/interpolating */
gegl_sampler_get_quad (sampler, xc-0.5, yc-0.5, xc+0.5, yc-0.5, xc+0.5, yc+0.5, xc-0.5, yc+0.5, buf); #endif

I think this is wrong, each of the seperate corners of the output pixel should probably have their coordinates computed separatly, before being passed in.

/Øyvind K.

The change below i think would complicate the usage. Since every time accessing the interpolator the coordinates for accessing should be ajusted by the value of x, y.

*GeglSampler *sampler = g_object_new (GEGL_TYPE_SAMPLER_CUBIC, "input", src, "format", babl_format ("RaGaBaA float"), "x", 50, "y", 50, "width", 200, /* or a rect could be passed, or a subbuffer of "height", 200, buffer expected to be created before making a sampler */ NULL);

The api below would introduce more complexity because each implementation would need other arguments.

gegl_sampler_get_quad (sampler, xc-0.5, yc-0.5, xc+0.5, yc-0.5, xc+0.5, yc+0.5, xc-0.5, yc+0.5, buf);
The curren API proposal works on the complete passed buffer and should alos work for scaling and transforming. Since the fractional part of x,y determines the offset.

void affine_cubic (GeglBuffer *dest,
GeglBuffer *src,
Matrix3 matrix)
{
gdouble arecip;

gint x, y, i, j,
pos,
src_w = src->width,
src_h = src->height,
u,
u1 = src->x,
u2 = src_w - 1,
pu,
v,
v1 = src->y,
v2 = src_h - 1,
pv;
gfloat *src_buf,
*dest_buf,
*src_ptr,
*dest_ptr,
abyss = 0.;
gdouble newval[4];
gdouble data[64];
gdouble du, dv, fu, fv;
Matrix3 inverse;

if (gegl_buffer_pixels (src) == 0 || gegl_buffer_pixels (dest) == 0) return;
dest_buf = g_new (gfloat, gegl_buffer_pixels (dest) << 2); g_assert (src_buf && dest_buf);
GeglInterpolatorCubic sampler = g_object_new (GEGL_TYPE_INTERPOLATOR_CUBIC, "input", src, "format", babl_format ("RaGaBaA float"), NULL);

matrix3_copy (inverse, matrix); matrix3_invert (inverse);

fu = du = inverse[0][0] * dest->x + inverse[0][1] * dest->y + inverse[0][2] - src->x; fv = dv = inverse[1][0] * dest->x + inverse[1][1] * dest->y + inverse[1][2] - src->y;

for (dest_ptr = dest_buf, y = 0; y < dest->height; y++) {
for (x = 0; x < dest->width; x++) {
gegl_interpolator_get (interpolator, fu, fv,dest_buf); fu += inverse [0][0];
fv += inverse [1][0];
}
du += inverse [0][1];
dv += inverse [1][1];
fu = du;
fv = dv;
}
gegl_buffer_set_fmt (dest, dest_buf, babl_format ("RaGaBaA float")); g_free (dest_buf);
}

On 10/16/06, geert.jordaens-CNXmb7IdZIWZIoH1IeqzKA@public.gmane.org wrote:

The change below i think would complicate the usage. Since every time accessing the interpolator the coordinates for accessing should be ajusted by the value of x, y.

*GeglSampler *sampler = g_object_new (GEGL_TYPE_SAMPLER_CUBIC, "input", src, "format", babl_format ("RaGaBaA float"), "x", 50, "y", 50, "width", 200, /* or a rect could be passed, or a subbuffer of "height", 200, buffer expected to be created before making a sampler */ NULL);

The purpose of this change would be to indicate which window within the buffer needs to be prepared for access (the specific sampler would then need to know exactly which data must be prepared). But as I wrote in my comment to those API additions relying on the buffer to be pre-clipped is probably better.

The api below would introduce more complexity because each implementation would need other arguments.

gegl_sampler_get_quad (sampler, xc-0.5, yc-0.5, xc+0.5, yc-0.5, xc+0.5, yc+0.5, xc-0.5, yc+0.5, buf);

The curren API proposal works on the complete passed buffer and should alos work for scaling and transforming. Since the fractional part of x,y determines the offset.

Nope, the point is that interpolation is wrong when scaling down you as will happen when scaling down using an affine transform or a perspective transform. By transforming the corners of a pixel, one would get an idea about the size needed for the resampling kernel.

If you scale a image to 10% of the original size using cubic, you have a situation where the data for each destination pixel is taken from a region of 4x4pixel, whilst it should at least be taken from a region of 10x10pixels, 84% of the image data is thrown away.

/Øyvind K.

Nope, the point is that interpolation is wrong when scaling down you as will happen when scaling down using an affine transform or a perspective transform. By transforming the corners of a pixel, one would get an idea about the size needed for the resampling kernel.

If you scale a image to 10% of the original size using cubic, you have a situation where the data for each destination pixel is taken from a region of 4x4pixel, whilst it should at least be taken from a region of 10x10pixels, 84% of the image data is thrown away.

OK, the handling of scaling down is not yet in the proposition

Could we not just add the scale factor to the API ?

On 10/17/06, geert.jordaens-CNXmb7IdZIWZIoH1IeqzKA@public.gmane.org wrote:

If you scale a image to 10% of the original size using cubic, you have a situation where the data for each destination pixel is taken from a region of 4x4pixel, whilst it should at least be taken from a region of 10x10pixels, 84% of the image data is thrown away.

OK, the handling of scaling down is not yet in the proposition

Could we not just add the scale factor to the API ?

The scale factor is not enough, if we scale it to 10% horizontally and 70% vertically (or add some kind of rotation as well). A fixed scale factor would no longer be correct. Doing a reverse transform of the "corners" of the destination pixel would give us all the information we need, and work for perspective transforms as well, hence the method I suggested.

/Øyvind K.

On Tue, 2006-10-17 at 09:27 +0200, Øyvind Kolås wrote:

On 10/17/06, geert.jordaens-CNXmb7IdZIWZIoH1IeqzKA@public.gmane.org wrote:

If you scale a image to 10% of the original size using cubic, you have a situation where the data for each destination pixel is taken from a region of 4x4pixel, whilst it should at least be taken from a region of 10x10pixels, 84% of the image data is thrown away.

OK, the handling of scaling down is not yet in the proposition

Could we not just add the scale factor to the API ?

The scale factor is not enough, if we scale it to 10% horizontally and 70% vertically (or add some kind of rotation as well). A fixed scale factor would no longer be correct. Doing a reverse transform of the "corners" of the destination pixel would give us all the information we need, and work for perspective transforms as well, hence the method I suggested.

I think you are mixing up interpolation and resampling. They are similar (abstractly anyway), but separate operations.

Interpolation is designed to add points between two or more nearby points using linear or cubic interpolating curves. Translating by fractional coordinates uses interpolation. Rotation (without scaling) also uses interpolation.

Resampling is designed to remove high frequency components while down sampling an image. Without resampling, when you downscale you get "jaggies" (aliasing) throughout the image. Examples of resampling are taking the average of the input area (box filter), gaussian weighted average, or applying a sync filter. Resampling is a convolution.

These operations are best implemented seperately. First interpolate to get a set of pixels that can be input into a resampling function to complete your transformation.

For example, suppose your x-scale is .1 (matrix X), your y scale is 2 (matrix Y) and your rotation angle is 45 degrees (matrix R). They are combined like so: YXR. First the rotation, do the rotation by 45 degrees by interpolating from nearby points. Then do the X scale by resampling 10 or more of the interpolated pixels (say, with a sync function, which is the best filter). Then complete the transformation by interpolating one extra Y point between every existing Y point using an interpolation function to complete your Y scale.

Geert's interpolation API is sufficient. Interpolation rarely takes or needs more than 16 pixels. A rescaling API needs to be able to declare the size of the neighborhood around the pixel (an interpolation API could probably use this). They probably both need another object for extending the region at the edges (copy, reflect, zeros, wrap, etc).

Both these operations are similar in that they need a neighborhood and produce a new value based on that neighborhood, but mathematicaly they are acomplishing two very different goals. Interpolation doesn't affect resolution and is used for shifting and scaling up. Resampling is for downsampling (aka decimating) a signal and is a type of anti-aliasing. They are used in very different cases, and should probably have different interfaces to reflect that.

On 10/17/06, Daniel Rogers wrote:

I think you are mixing up interpolation and resampling. They are similar (abstractly anyway), but separate operations.

My claim is that interpolation and decimation both are resampling operations, and thus shuld be grouped together.

Geert's interpolation API is sufficient. Interpolation rarely takes or needs more than 16 pixels. A rescaling API needs to be able to declare the size of the neighborhood around the pixel (an interpolation API could probably use this). They probably both need another object for extending the region at the edges (copy, reflect, zeros, wrap, etc).

As I see it this needs to be handled by the buffer implementation and not the resamplers, each operation will only operate on a smaller region of the image any abyss behavior should already be applied to the data outside when data is mapped from storage to a linear buffer to be operated on. (currently the behavior is always zeroed out data)

Both these operations are similar in that they need a neighborhood and produce a new value based on that neighborhood, but mathematicaly they are acomplishing two very different goals. Interpolation doesn't affect resolution and is used for shifting and scaling up. Resampling is for downsampling (aka decimating) a signal and is a type of anti-aliasing. They are used in very different cases, and should probably have different interfaces to reflect that.

Both scaling up and scaling down is resampling, changing the spatial sampling rate of the image raster. Increasing the samling rate is interpolation, decreasing it is decimation, both of them are resampling.

One place in the GIMP code where these are used is the transformation tools, the transformation code in GIMP uses a matrix to represent the entire transformation. In the case of a perspective transform you need to do both interpolation and decimation with the same resampler. What I was suggesting is to keep the same cache of pixels available both for the interpolation and decimation and use different functions for retrieving the data.

/OEyvind K.