Trying to do tiling noise with a 2D noise function limits you somewhat. There are a few ways you can do it, but all have their problems or issues.
1) Ensure that each individual octave tiles with itself (which, I suspect, is what your noise2 function does). This can be done by modulo-ing the coordinates at each layer with their repeating period size. The limitation of this, of course, is that it leaves you unable to perform certain common artifact-reduction tricks such as rotating the domain of a noise layer around an arbitrary axis before calling the noise function. This trick is helpful for reducing or eliminating grid line artifacts in the final result. However, since this rotates the noise layer off the grid, the edges will no longer line up and tile the way they need to. Additionally, this imposes the rigid requirement that your lacunarity be integral, so that the entire grid is represented in the output, preventing you from using a non-integral lacunarity to again try to reduce grid-based artifacts. Integral lacunarity means that the grid boundaries of each successive layer line up exactly with the grid boundaries in preceding layers, thus amplifying the results of grid artifacts in certain types of noise, such as ridged noise.
2) Blend 4 regions of noise arranged in a 2x2 pattern, with continuous weighting of the 4 blended values based on distance from the horizontal and vertical edges. This has the problem of skewing the distribution of the function output toward neutral, especially at the center of the image. Consider this image (the one on the left uses this method of seamless tiling)
Because the final value of a pixel in the middle is averaged from 4 different pixels, the result is an averaged blend, with much of the contrast of the function destroyed. Additionally, since the blending is weighted by distance-to-edge, you introduce grid-based artifacts since the averaging becomes less pronounced as you approach the edge. The edge regions have higher contrast (closer to the un-blended base function) than the middle. Applied tiling on a plane, this pattern becomes very obvious.
If you can increase to higher dimensional noise, other methods open up to you:
3) Perform a spherical mapping from a 3D function. This, of course, has the large drawback of introducing a lot of distortion into the final mapping, as singularities are formed at the poles of the sphere, as demonstrated by this example:
If the noise were being mapped back onto the surface of a sphere the distortion would vanish, but since it is being mapped to a 2D rectangle it becomes very pronounced. Set to tile on a plane, this type of mapping would produce severe banding artifacts.
4) In a journal entry
I wrote some time back (which was later re-written as an article for Game Developer Magazine) I discussed an alternative technique: using a 4D noise function to generate seamless 2D noise. The X axis of the final image is mapped to a 2D circle, and the Y axis is mapped to a second 2D circle. The first circle lies upon the X/Y plane in 4D space, and the second circle lies upon the Z/W plane. The result is a 4-tuple that is used to index the noise function at a given (Ximg
) coordinate. In the image of the high-contrast cellular noise posted above, the image on the right is of tiling cellular noise using this method. You can see that this method preserves the overall contrast of the function, unlike the blending method. While some distortion is introduced by the non-linearity of the axis->circle mappings, it is a much more tolerable and even distortion than what is generated by the spherical mapping.
This method works for all sorts of functions, as long as a 4D variant can be constructed. It has the advantage of preserving contrast and working with arbitrarily-rotated domains at any level of the layering. The disadvantages include performance (requiring 4D adds additional computational complexity, making it difficult to justify implementation in some cases, such as inside shader code) and the subtle distortion introduced by the mapping to circle, which makes it difficult if you are working with regular repeating patterns such as checkerboards. (The checkerboard will be warped beyond recognition.)
noise2() implies 2-dimensional noise, yet you're calling it as if it were noise4()? Or do px and py do something other than act as coordinates when passed to your noise2?