I used GPUImage framework (GPUImageToneCurveFilter) to read Adobe acv file and generate the LUT for shader texture recently. After successfully using this utility to render customized color tone in my image, I'm curious what kind of spline interpolation algorithm was used to achieve that. In the sample code, I could find out there're some second derivative calculation. But I can't understand where the math parameters came from. Is it possible to tell me some theoretical reference for further interpolation studying? Thanks.
[Edit]
As Spektre mentioned below, I paste the codes of spline calculation in GPUImageToneCurveFilter.m. It was Objective-C version.
After getting control points from ACV file and convert them from (0, 1) to (0, 255), send the points to splineCurve function for interpolation as below:
NSMutableArray *splinePoints = [self splineCurve:convertedPoints];
Two interpolating functions as:
- (NSMutableArray *)splineCurve:(NSArray *)points
{
NSMutableArray *sdA = [self secondDerivative:points];
// [points count] is equal to [sdA count]
NSInteger n = [sdA count];
if (n < 1)
{
return nil;
}
double sd[n];
// From NSMutableArray to sd[n];
for (int i=0; i<n; i++)
{
sd[i] = [[sdA objectAtIndex:i] doubleValue];
}
NSMutableArray *output = [NSMutableArray arrayWithCapacity:(n+1)];
for(int i=0; i<n-1 ; i++)
{
#if TARGET_IPHONE_SIMULATOR || TARGET_OS_IPHONE
CGPoint cur = [[points objectAtIndex:i] CGPointValue];
CGPoint next = [[points objectAtIndex:(i+1)] CGPointValue];
#else
NSPoint cur = [[points objectAtIndex:i] pointValue];
NSPoint next = [[points objectAtIndex:(i+1)] pointValue];
#endif
for(int x=cur.x;x<(int)next.x;x++)
{
double t = (double)(x-cur.x)/(next.x-cur.x);
double a = 1-t;
double b = t;
double h = next.x-cur.x;
double y= a*cur.y + b*next.y + (h*h/6)*( (a*a*a-a)*sd[i]+ (b*b*b-b)*sd[i+1] );
if (y > 255.0)
{
y = 255.0;
}
else if (y < 0.0)
{
y = 0.0;
}
#if TARGET_IPHONE_SIMULATOR || TARGET_OS_IPHONE
[output addObject:[NSValue valueWithCGPoint:CGPointMake(x, y)]];
#else
[output addObject:[NSValue valueWithPoint:NSMakePoint(x, y)]];
#endif
}
}
// The above always misses the last point because the last point is the last next, so we approach but don't equal it.
[output addObject:[points lastObject]];
return output;
}
- (NSMutableArray *)secondDerivative:(NSArray *)points
{
const NSInteger n = [points count];
if ((n <= 0) || (n == 1))
{
return nil;
}
double matrix[n][3];
double result[n];
matrix[0][1]=1;
// What about matrix[0][1] and matrix[0][0]? Assuming 0 for now (Brad L.)
matrix[0][0]=0;
matrix[0][2]=0;
for(int i=1;i<n-1;i++)
{
#if TARGET_IPHONE_SIMULATOR || TARGET_OS_IPHONE
CGPoint P1 = [[points objectAtIndex:(i-1)] CGPointValue];
CGPoint P2 = [[points objectAtIndex:i] CGPointValue];
CGPoint P3 = [[points objectAtIndex:(i+1)] CGPointValue];
#else
NSPoint P1 = [[points objectAtIndex:(i-1)] pointValue];
NSPoint P2 = [[points objectAtIndex:i] pointValue];
NSPoint P3 = [[points objectAtIndex:(i+1)] pointValue];
#endif
matrix[i][0]=(double)(P2.x-P1.x)/6;
matrix[i][1]=(double)(P3.x-P1.x)/3;
matrix[i][2]=(double)(P3.x-P2.x)/6;
result[i]=(double)(P3.y-P2.y)/(P3.x-P2.x) - (double)(P2.y-P1.y)/(P2.x-P1.x);
}
// What about result[0] and result[n-1]? Assuming 0 for now (Brad L.)
result[0] = 0;
result[n-1] = 0;
matrix[n-1][1]=1;
// What about matrix[n-1][0] and matrix[n-1][2]? For now, assuming they are 0 (Brad L.)
matrix[n-1][0]=0;
matrix[n-1][2]=0;
// solving pass1 (up->down)
for(int i=1;i<n;i++)
{
double k = matrix[i][0]/matrix[i-1][1];
matrix[i][1] -= k*matrix[i-1][2];
matrix[i][0] = 0;
result[i] -= k*result[i-1];
}
// solving pass2 (down->up)
for(NSInteger i=n-2;i>=0;i--)
{
double k = matrix[i][2]/matrix[i+1][1];
matrix[i][1] -= k*matrix[i+1][0];
matrix[i][2] = 0;
result[i] -= k*result[i+1];
}
double y2[n];
for(int i=0;i<n;i++) y2[i]=result[i]/matrix[i][1];
NSMutableArray *output = [NSMutableArray arrayWithCapacity:n];
for (int i=0;i<n;i++)
{
[output addObject:[NSNumber numberWithDouble:y2[i]]];
}
return output;
}
Hope these explanations could help to figure out what kind of algorithm that Brad L. used to interpolate the curve that makes the result as Photoshop did.
I do nor recognize the language (apart the C++ like stuff) so it is hard to read it (for me) but I think this
double y= a*cur.y + b*next.y + (h*h/6)*( (a*a*a-a)*sd[i]+ (b*b*b-b)*sd[i+1] );
is what you are asking about. The term:
a*cur.y + b*next.y
is basic linear interpolation so when t=0.0 the result is current point and if t=1.0 the result is next point. The rest is some kind of superponated curve with specific parameters probably producing some kind of bump pattern for the gradient scaled by the distance of points (that is why deltas are used for it). To reverse it would be hard (and I am too lazy for that) but the therms are reminding me of Bernstein polynomials (Used for BEZIER and SPLINE). The a^3 and b^3 suggest cubic curve. You should start with plotting graph which it produce and than try to create a polynomial interpolation that will have similar properties. If the result match the term you found you answer. For more info see this:
And all the links there (especially the How to construct own interpolation ones).