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