# Why my cosine interpolation of a cube's face doesn't work?

NB please : executable use cases are available at the end of this question. I begin this question by showing you the problems of my program, then I explain how the latter works, and finally I end up with these both executable use cases.

# What I want to do

I want to draw the bottom face of a cube. To do that, I interpolate between the 4 vertices of this face. I use two kinds of interpolation :

1. Linear

2. Cosine

I know it's typically the aim of GPUs, but I wanted to implement it myself to learn and try the "deep" mechanisms of 3D rendering.

# Results

Each point of the following both pictures are interpolated between the 4 vertices of the bottom face (except the 8 vertices of my cube).

## With linear interpolation

It's OK, no problem.

## With cosine interpolation

Error. As you can see, there is some weird problem : face's points are inclined towards left.

# How do I interpolate between 4 points ?

I test all the combinations of interpolation's weights and, if the sum of the weights equals 1, I draw the resulting point. I know this method is not very good (a lot of consumption of CPU and RAM), and there are P duplicates. But this problem is (1) secondary and (2) out of the range of this StackOverflow question.

In other words : the definition of the interpolation between 4 points is : P = aA + bB + cC + dD. a + b + c + d must equals 1 to draw P, otherwise it's not drawn. (a;b;c;d) \in [0;1]^4 by the way.

A concrete example is (step = 1) :

P = 0a + 0b + 0c + 0d is NOT drawn

P = 0a + 0b + 0c + 1d is drawn

P = 0a + 0b + 1c + 0d is drawn

P = 0a + 0b + 1c + 1d is NOT drawn

etc. ("I test all the combinations"). (Last is : P = 1a + 1b + 1c + 1d, which is not drawn).

To test all these combinations, I iterates on the weights, with a step, using recursion.

## Case of the linear interpolation

The above definition doesn't change : P = aA + bB + cC + dD.

## Case of the cosine interpolation

The above definition does change : P = aA + bB + cC + dD becomes P = a'A + b'B + c'C + d'D, with : a' = 0.5(1 - cos(a * \pi)) and the idea is the same for the three other weights.

This new definition is not random, I could explain how I found it, but it's out of the range of this topic (in résumé : a simple remap of the cosine function).

# Question

Why does it works well with the linear interpolation, and not the cosine one ? Indeed, the cosine face's points are inclined towards left.

# By the way : interpolations between two points work well

Linear and cosine interpolations both work very well when used between two points, as you can see below.

## Linear interpolation

It's OK, no problem.

## Cosine interpolation

It's OK, no problem.

# Implementation and Executable

I wrote the program in Scala. First I show you the functions, then two use cases you can execute along with the given functions.

## Tests all the combinations of interpolations between n points of k coordinates

def computeAllPossibleInterpolatedPoints(step : Double, points : Seq[Seq[Double]], transform: (Double) => Double) : Seq[Seq[Double]] = {
var returned_object : Seq[Seq[Double]] = Seq.empty[Seq[Double]]
recursiveInterpolation(0, Seq.empty[Double])

def recursiveInterpolation(current_weight_id : Int, building_weights : Seq[Double]) : Unit = {

(.0 to 1.0 by step).foreach(current_step => {
if (current_weight_id < points.size - 1) {
recursiveInterpolation(current_weight_id + 1, building_weights :+ current_step)
} else {
val found_solution = (building_weights :+ current_step).map(transform)
if(BigDecimal(found_solution.sum).setScale(5, BigDecimal.RoundingMode.HALF_UP).toDouble == 1.0) {
returned_object = returned_object :+ interpolation(found_solution, points)
}
}
})
}

returned_object
}


## Interpolates between n points of k coordinates

def interpolation(weights: Seq[Double], points: Seq[Seq[Double]], transform: (Double) => Double = null) : Seq[Double] = {
var transformed_weights = weights
if(transform != null) {
transformed_weights = weights.map(transform)
}
if(BigDecimal(transformed_weights.sum).setScale(5, BigDecimal.RoundingMode.HALF_UP).toDouble != 1.0) {
println("ERROR : SUM(weights) != 1. Returning null.")
return null
}
if(transformed_weights.exists(weight => BigDecimal(weight).setScale(5, BigDecimal.RoundingMode.HALF_UP).toDouble < 0)
||
transformed_weights.exists(weight => BigDecimal(weight).setScale(5, BigDecimal.RoundingMode.HALF_UP).toDouble > 1)) {
println("ERROR : EXISTS(weight) / weight < -1 OR weight > 1. Returning null.")
return null
}

transformed_weights.zip(points).map(
weight_point => weight_point._2.map(coordinate => weight_point._1 * coordinate)
).reduce((point_a : Seq[Double], point_b : Seq[Double]) => point_a.zip(point_b).map(coordinate_points => coordinate_points._1 + coordinate_points._2))
}


## Interpolation functions

def linear(weight : Double) : Double = {
weight
}

def cosine(weight : Double) : Double = {
(1 - Math.cos(weight * Math.PI)) * 0.5
}


## Two use cases : between 2 points, each made of 2 coordinates (first use case : linear interpolation, second use case : cosine interpolation)

Note that you can add 2 new points to each of the following functions, to interpolate a face instead of a segment (for now, the below code indeed interpolates points of a segment).

val interpolated_points : Seq[Seq[Int]] = interpolator.computeAllPossibleInterpolatedPoints(0.05, Seq(Seq(22, 22), Seq(33, 33)), interpolator.linear).map(_.map(_.intValue()))

val interpolated_points : Seq[Seq[Int]] = interpolator.computeAllPossibleInterpolatedPoints(0.05, Seq(Seq(22, 22), Seq(33, 33)), interpolator.cosine).map(_.map(_.intValue()))

• Not sure if this suits this cite (and not stackoverflow), but I definitely love the pictures. – lisyarus Jun 8 '18 at 14:30
• Ahah thank you:)! – JarsOfJam-Scheduler Jun 8 '18 at 14:39
• I think it can be an aliasing problem. It happens when sampling periodic functions. – Cesareo Jun 8 '18 at 14:43
• Maybe this aliasing would occurs for linear interpolation between 2 and 4 points, and for cosine interpolation for 2 points. However, these 3 cases work well :/ In other words : why aliasing would only concern cosine interpolation between 4 points ? – JarsOfJam-Scheduler Jun 8 '18 at 15:17
• I haven't read your implementation yet, but why do you draw $P = 0a + 0b + 0c + 0d$? According to your explanation, you should only draw points if the sum of weights is 1, but this point has weight sum 0. – N.Bach Jun 8 '18 at 15:42

I only remember after struggling to read through code, but I'd recommend posting pseudo-code rather than code here. Regardless, I tried to understand your implementation, and as far as I can tell, you test weights $a,b,c,d$ such that they sum to 1, and for these $a+b+c+d=1$ you directly plug them into your cosine interpolation to compute $a',b',c',d'$. In general, this will not give you a point in the face because you actually want the constraint $a'+b'+c'+d'=1$ to perform an interpolation.
Now, why does this work with a segment interpolation and not a face interpolation? The tl;dr is that $a+b=1$ will ensure $a'+b'=1$. Indeed $$\cos(a\pi)=\cos(\pi-b\pi)=-\cos(-b\pi)=-\cos(b\pi)$$ and this yields $$a'+b' = \frac 12(1-\cos(a\pi)) +\frac 12(1-\cos(b\pi))=1$$ When you have more coefficients however, this is unlikely to hold, so your interpolation is not really an interpolation.
• ...able to figure out which is which quickly. Now regarding how to correct your implementation, I suggested to use transform earlier in your code, but that implies you have to change more than just putting transform before testing your weights. If I was doing the implementation, I'd store the transformed weights, check they are valid, and then perform a "plain" interpolation using these valid weights. I don't understand the last couple lines of interpolation, but I'd get rid of map(transform) or change it to map(linear). In computeAllPossibleInterpolatedPoints I'd put ... – N.Bach Jun 8 '18 at 16:57
• ... building_weights :+ transform(current_step) each time you sample a new weight. – N.Bach Jun 8 '18 at 16:57
• @JarsOfJam-Scheduler I haven't fully read your update, but I tested an implementation of my own in another language. I obtained mostly the same result, and this is where I wish I had more programming experience. From what I observed, the main issue stems from testing the sum of indices being equal to 1.0. As a rule of thumb, you don't want to test equality with floating point arithmetic, especially when sampling something. Due to rounding errors in floating point arithmetic, two floats/double are rarely really equal to one another. It depends on the implementation, but ... – N.Bach Jun 9 '18 at 1:11