Perlin noise is an accumulation of many different wave functions, so the obvious thing to do is to separate each of these out and work on them one at a time. Each one needs a set of unique keypoints at regular intervals over the domain of the function. Because the noise generator is two dimensional, each keypoint needs to have an x coordinate, a y coordinate, and a magnitude. The amount of the keypoints should not be fixed, as the different functions that will need to be generated will have to have different frequencies and therefore different amounts of keypoints. Luckily, I wrote the one-dimensional keypoint generator with this in mind - now all I had to do to run this in two dimensions was to simply use tail recursion to call the one-dimensional keypoint generator many times (making sure to update the random generator each time, or all the random magnitudes will be the same) and append all the results to each other.
The next interesting part - figuring out which keypoints are the ones bounding a given (x, y) coordinate. The output of this function serves as input to an interpolation function, which uses the bounding box magnitudes to figure out what the magnitude of a given point would be. When I was thinking about how to do this, I was wondering what would happen if an (x, y) coordinate lied directly on a keypoint, or what if it doesn't lie in between a square of 4 but actually lies on a border of the square, in between exactly 2 keypoints. Based on this reasoning, I had this function return either 1, 2, or 4 keypoints; if the function returns 1 coordinate, the given point lies exactly on a keypoint, and the interpolation function can just use the magnitude of that point. If two points are returned, the function can use the linear interpolation that I already wrote in the last version of the noise generator. If four points are returned, then the full bicubic interpolation must be run.
That being said, my implementation of this function is not the best. In order to find if the given point lies on a keypoint, I filter the list of keypoints for one that lies on the given point. If this filtered list is empty, I proceed to filter the list of keypoints for all the ones that have only one coordinate equal to that coordinate in the given point (for the linear case). If this list isn't empty, I filter based on the greatest value that is less than the given point and the smalles value that is less than the given point, and return those. If none of those cases were true, I filter the list again twice to get all the points that are in the four quadrants around the given point. The top-left point is simply the last point in the top-left quadrant, and the bottom-right point is simply the first coordinate in the bottom-right quadrant (assuming the points come sorted, and filter doesn't rearrange the list). In order to get the top-right point and the bottom-left point, I had to do folds across the list to find one with minimum in one coordinate and maximum in the other. That makes for a ton of filters, each of which is O(n), and even 2 folds thrown into the mix. After I was done, running the program took a very long time.
Looking back, I feel that I should change the output of the function so that it always returns 4 points. If the given point lies on a point, I should return a specified square (say, to the bottom-right of the given point), and if a point lies on a line, I should return the square either below or to the right of the given point. This way, I can unify all the multiple cases into one case, limit the number of filters, and the interpolation function can deal with when a coordinate value = 0. Cubic and bicubic interpolation do this for us, so no additional code would have to be added.
The interpolation function itself is also a fairly interesting function. After reading the Wikipedia article on bicubic interpolation, it kindly provides a matrix that can be used. The use of this matrix is fairly simple - make a vector of all the values and partial derivatives as specified, multiply them by the given matrix, and use the output vector as coefficients into the simple double-sigma formula at the top of the page (Under "Bicubic spline interpolation"). Because Haskell is a very math-y language, implementing the double-sigma formula was extremely simple. After thinking for a second or two about how to do the matrix multiplication, I realized that Haskell's standard formulas were enough. Matrix multiplication is simply multiplying each entry i in row j with entry i in the vector. So, I am going to be doing some computation over ever row of the matrix - this looks like a map. The computation that I'm going to be doing is a dot product - which looks like a sum over a zip. All in all - the multiplication looks like this:
map (sum . (zipWith (*) vector)) matrix
Of course, this assumes that the matrix is in row-major form.
All in all, everything is mostly the same as the one-dimensional case, except with another dimension artificially added in. My code can be found here.
I did have a neat idea about making an image. All the two-dimensional noise images that I had seen have been grayscale, so that it looks kind of like fog. My idea was to simply run the algorithm 3 times, mapping the output into each of the red, green, and bue channels of the image, respectively. This would create a bright, more appealing image. Here is one such image, with 3 noise functions added together: