As hydrodynamic models of planetary nebulae advance from 1D to 2D calculations, it becomes desirable to make the same step for reconstruction techniques, which aim at deriving from the 2D intensity distributions of PNe images their 3D structure. A basic step is the determination of a symmetry axis and its orientation in space which can be described by two angles, one measured in the plane of the sky, and one measured with respect to the tangential plane. While the first one can be determined from the images, e.g. by applying the criterion of maximum normalized correlation between the object halves, which we found to yield the best results, the second angle is treated as a free parameter (see below). The steps of the iterative reconstruction algorithm are the following (with image(x, y) = input image):

1. 1. Fill a 3D kartesian density grid(xyz) (random or continuous), observing the constraint ∫ grid(xyz)dz = image(x, y), where z is along the line of sight.

2. 2. Transform to cylindrical symmetry system (r, z', ϕ) and read out density (r, z') averaging over ϕ.

3. 3. Return density(r, z') into grid(xyz) and normalize the ∫ grid(xyz)dz to match the input image.

4. 4. Return to (2.) as long as two subsequent density(r, z') estimations differ by more than a specified limit.