Nice solution.  You could also use GDAL utilities to handle this-- which are
written in C++ (and also have APIs for Python, etc).

Best,
Steve


"Adrian" <nomail@nomail> wrote:
> clipka <ano### [at] anonymousorg> wrote:
> > Am 07.11.2011 11:49, schrieb Adrian:
> > > I have an image [1] of a spherical surface that I want to map onto a sphere, but
> > > it is a Mollweide projection. I tried to convert it to a rectangular Mercator
> > > projection (is that even the right one?) with Mathematica [2], but it doesn't
> > > really work. I do get a rectangular image which kind of looks right, but when
> > You want to convert the image to an "equirectangular projection" (see
> > http://en.wikipedia.org/wiki/Equirectangular_projection); that'll allow
> > you to use map type 1.
>
> Thanks, I guess I was over-complicating things! Also thanks to Christian for the
> hint at Matthew's "project". Unfortunately it yields a wrong result for my case.
> I am considering letting him know about it, because that would have been a great
> tool otherwise. I appreciate the hint at the WMAP data as well, but I was
> curious on how to handle this for future projects as well.
>
> For the record, since Mathematica didn't yield a satisfactory result for
> whatever reason, I wrote a small Python script using SciPy that takes a cropped
> *ppm Mollweide projection and converts it to plate carree, which can be mapped
> onto a sphere by Povray just fine. Of course, it is much slower than a C/C++
> solution, but it's fast enough for my purposes (around 2m on my PC for a
> 2048x1024 image). If anyone wants to use it, see below for the source.
>
> Adrian
>
>
>
>
> #!/usr/bin/python
>
> import numpy as np
> from scipy import misc
> import sys
>
> infile = "wmap.ppm"
> outfile = "wmap2.ppm"
>
> arr = misc.imread(infile)
> result = np.empty(arr.shape)
> height, width, depth = arr.shape # Backwards but correct
> print "Dimensions: %ix%i" % (width, height)
>
> class memoized(object):
>    def __init__(self, func):
>       self.func = func
>       self.cache = {}
>    def __call__(self, *args):
>       try:
>          return self.cache[args]
>       except KeyError:
>          value = self.func(*args)
>          self.cache[args] = value
>          return value
>
> @memoized #Improves runtime by factor of 2 easily
> def theta(phi):
>     # First guess for newton method
>     xn = 0
>     old = xn + 1.
>     while abs(xn - old)>1e-6:
>         old = xn
>         xn = xn - (xn + np.sin(xn) - np.pi * np.sin(phi)) / (1 + np.cos(xn))
>     return xn / 2.
>
> def mollweide(x, y):
>     th = theta(y)
>     return np.sqrt(8) * x * np.cos(th) / np.pi, np.sqrt(2) * np.sin(th)
>
> def img2coord(a, b):
>     return (a * 1. / width - .5) * 2 * np.pi, (b * 1. / height - .5) * np.pi
>
> def coord2img(x, y):
>     # Try to stay inside ellipse to minimize background effects
>     xround = yround = np.floor
>     if x<0: xround = lambda z: -np.floor(-z)
>     if y<0: yround = lambda z: -np.floor(-z)
>     return (xround((x / (2 * np.sqrt(8)) + .5) * width),
>             yround((y / (2 * np.sqrt(2)) + .5) * height))
>
> for a in range(width):
>     sys.stdout.write("\r%i%%" % ( 100 * a / width ))
>     sys.stdout.flush()
>     for b in range(height):
>         x, y = img2coord(a, b)
>         x, y = mollweide(x, y)
>         x, y = coord2img(x, y)
>         result[b, a] = arr[y-1, x-1]
>
> misc.imsave(outfile, result)
> sys.stdout.write("\nDone\n")

Post a reply to this message