Bildtransformationen

[Yukterez]

zum Hauptartikel ←

[Yukterez]

Der Hauptartikel ist in der Wolfram Sprache geschrieben; man kommt aber, wenn auch etwas umständlicher, auch mit Python, Matlab oder Maple ans Ziel. Beispiel: Mollweide→Equirectangular in

Python 2.7

1: Speichere Original-Bild im .ppm (Portable Pixmap) Format
2: Run Python 2.7.4 plus Scipy, Numpy und Pil Module
3: Run Script:

Code: Alles auswählen

# Python 2.7 (32Bit) Syntax
# Benötigte Module: Numpy, Scipy, PIL

import numpy
from scipy import misc
import sys

infile = "C:\\Users\\Yukterez\\Desktop\\planckanomalies_mollweide.ppm"
outfile = "C:\\Users\\Yukterez\\Desktop\\planckanomalies_equirectangular.ppm"

arr = misc.imread(infile)
result = numpy.empty(arr.shape)
height, width, depth = arr.shape
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
def theta(phi):
    xn = 0
    old = xn+1.
    while abs(xn-old)>1e-6:
        old = xn
        xn = xn-(xn+numpy.sin(xn)-numpy.pi*numpy.sin(phi))/(1+numpy.cos(xn))
    return xn/2.

def mollweide(x, y):
    th = theta(y)
    return numpy.sqrt(8)*x*numpy.cos(th)/numpy.pi, numpy.sqrt(2)*numpy.sin(th)

def img2coord(a, b):
    return (a*1./width-.5)*2*numpy.pi, (b*1./height-.5)*numpy.pi

def coord2img(x, y):
    xround = yround = numpy.floor
    if x<0: xround = lambda z: -numpy.floor(-z)
    if y<0: yround = lambda z: -numpy.floor(-z)
    return (xround((x/(2*numpy.sqrt(8))+.5)*width),
            yround((y/(2*numpy.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")

Codes in Matlab und Maple folgen demnächst.

Mehrsprachig,

[Yukterez]

Anwendungsbeispiel:

Wir haben eine Karte die wir

1) neu zentrieren
2) auf eine Kugel projizieren

wollen. Als Rohmaterial wählen wir die Planck Polarisationskarte. Nachdem das Bild offensichtlich im Hammer-Aitoff-Format vorliegt benötigen wir dazu folgende Formeln:

Code: Alles auswählen

(* Hammer Aitoff >> Equirectangular *)
HA2Eq[{x_, y_}] := {X[x, y], Y[x, y]};
X[x_, y_] := (2 Sqrt[2] Cos[y] Sin[x/2])/Sqrt[1 + Cos[y] Cos[x/2]];
Y[x_, y_] := (Sqrt[2] Sin[y])/Sqrt[1 + Cos[y] Cos[x/2]];

(* Equirectangular >> Hammer Aitoff *)
Eq2HA[{x_, y_}] := {Lha[x, y], Bha[x, y]};
ζ[x_, y_] := Sqrt[1 - x^2/16 - y^2/4];
Bha[x_, y_] := ArcSin[ζ[x, y] y];
Lha[x_, y_] := 2 ArcTan[ζ[x, y] x/(2 (2 ζ[x, y]^2 - 1))];

(* Kugelkoodinaten >> 3D *)
xyz[{x_, y_}] :=
 {Sin[y] Cos[x],
  Sin[y] Sin[x],
  Cos[y]}

(* Rotation auf der x-Achse *)
Xyz[{x_, y_, z_}, α_] :=
 {x Cos[α] - y Sin[α],
  x Sin[α] + y Cos[α], z}

(* Rotation auf der y-Achse *)
xYz[{x_, y_, z_}, β_] :=
 {x Cos[β] + z Sin[β], y, z Cos[β] - x Sin[β]}

(* Rotation auf der z-Achse *)
xyZ[{x_, y_, z_}, ψ_] :=
 {x, y Cos[ψ] - z Sin[ψ], y Sin[ ψ] + z Cos[ψ]}

(* 3D >> Kugelkoordinaten *)
xy[{x_, y_, z_}] :=
 {ArcTan[x, y], ArcCos[z]}

(* Winkeltransformation *)
rm[pic_, α_, β_, ψ_] :=
  xy[xyZ[xYz[Xyz[xyz[pic], α], β], ψ]];
Eq2Eq[{x_, y_}] := rm[{x, y}, α, β, ψ];

(* Equirectangular >> Kugel *)
Kugel[pic_, {X_, Y_, Z_}] := SphericalPlot3D[
   1, {u, 0, π}, {v, 0, 2 π},
   Mesh -> None,
   TextureCoordinateFunction -> ({#5, 1 - #4} &),
   PerformanceGoal -> "Quality",
   PlotStyle -> Directive[Texture[pic]],
   Lighting -> "Neutral",
   Axes -> False,
   RotationAction -> "Clip",
   SphericalRegion -> True,
   Boxed -> False,
   ViewPoint -> {X, Y, Z},
   ViewAngle -> 1/5,
   ImageSize -> 700,
   PlotPoints -> 240];

Jetzt laden wir das Originalbild:

Code: Alles auswählen

(*Rohmaterial laden*)
HA = Import["C:\\Users\\Yukterez\\Desktop\\cmb.polarisation.jpg"]

und verwandeln es in eine Plattkarte:

Code: Alles auswählen

(* Hammer-Aitoff >> Equirectangular *)
Eq = ImageTransformation[HA, HA2Eq,
  DataRange -> {{-2 Sqrt[2], 2 Sqrt[2]}, {-Sqrt[2], Sqrt[2]}},
  PlotRange -> {{-π, π}, {-π/2, π/2}}]

die wir neu zentrieren:

Code: Alles auswählen

(* Neue Zentrierung *)
α = +π; β = +π/4; ψ = +π/4; (* Perspektive *)
Eq2 = ImageTransformation[Eq, Eq2Eq,
  DataRange -> {{-π, π}, {0, π}},
  PlotRange -> {{-π, π}, {0, π}}]

wieder zurückverwandeln:

Code: Alles auswählen

HA2 = ImageTransformation[Eq2, Eq2HA,
  DataRange -> {{-π, π}, {-π/2, π/2}},
  PlotRange -> {{-2 Sqrt[2], 2 Sqrt[2]}, {-Sqrt[2], Sqrt[2]}}]

und auf eine Kugel projizieren:

Code: Alles auswählen

(* Equirectangular >> Kugel *)
α = π; β = 0; ψ = 0; x = 5;
sphere1 = Kugel[Eq,
  {x Cos[α] Cos[β],
   x Cos[ψ] Sin[α] + x Cos[α] Sin[β] Sin[ψ],
   x Sin[α] Sin[ψ] - x Cos[α] Cos[ψ] Sin[β]}]

Anderer Blickwinkel:

Code: Alles auswählen

(* Equirectangular >> Kugel *)
α = π; β = π/4; ψ = -π/4; x = 5;
sphere1 = Kugel[Eq,
  {x Cos[α] Cos[β],
   x Cos[ψ] Sin[α] + x Cos[α] Sin[β] Sin[ψ],
   x Sin[α] Sin[ψ] - x Cos[α] Cos[ψ] Sin[β]}]

,

→ zurück zum Hauptartikel ←