Algorithm and a Method for Characterizing Fractal Volumes

Abstract

A computer implemented method for obtaining an analytical representation of an internal structure and spatial properties distribution of a selected physical domain includes identifying d-dimensional correspondences of measured spatial properties or field distributions; and applying an inverse algorithm to the d-dimensional spatial properties or field distributions to calculate the Weierstrass-Mandelbrot (W-M) fractal model to thereby determine parameters defining an analytical and continuous Weierstrass-Mandelbrot (W-M) representation.

Description

FIELD OF THE INVENTION

The present invention is directed to a method for obtaining an analytical representation of an internal structure and spatial properties distribution of a selected physical domain, and in particular to applying an inverse algorithm to calculate the W-M fractal model and thereby determine parameters defining an analytical and continuous Weierstrass-Mandelbrot (W-M) representation.

BACKGROUND OF THE INVENTION

The need to adequately represent and model the structure and the behavior of physical or virtual volumes containing physical property distributions or physical scalar fields as they are relevant to various applications has motivated the utilization of appropriate fractal functions that can be used to represent such volumes.

A lower dimension approach that is useful for determining the mechanical, thermal, electrical and fluid properties of the static and the sliding contact between two different conductors of heat and electric current involves the full identification of the parameters of the 2D Weirstrass-Mandelbrot (W-M) function, described in (Michopoulos, J. G., and Iliopoulos, A. P., 2011. “Complete High Dimensional Inverse Characterization of Fractal Surfaces”, 2011 ASME International Design Engineering Technical Conferences & Computers and Information in Engineering Conference, Vol. CD-ROM, ASME) (“Michopoulos et al.”) and in U.S. patent application Ser. No. 13/507,968 filed Aug. 8, 2012 and incorporated herein by reference. The present invention represents a generalization of this methodology for higher dimensional fractal volumes.

BRIEF SUMMARY OF THE INVENTION

According to the invention, a method for obtaining an analytical representation of an internal structure and spatial properties distribution of a selected physical domain includes identifying d-dimensional correspondences of measured spatial properties or field distributions; and applying an inverse algorithm to the d-dimensional spatial properties or field distributions to calculate the W-M fractal model to thereby determine parameters defining an analytical and continuous fractal Weierstrass-Mandelbrot (W-M) representation.

Exemplary physical domains for application of the method of the invention are generally any dense or coarse material where the method is applied to its internal structure as enclosed in its defined volume and surface such as a metal alloy, a carbon epoxy composite, a biological tissue, and any physical or virtual entity with a varying distribution of a scalar field. The physical domain may also be a man-made structure, like a microprocessor, an aircraft, a ship, an automobile, a building or similar technological objects. The physical domain alternatively may be a live organism or part of a live organism, or a digital representation like a 1 dimensional, 2 dimensional a 3 dimensional or a d-dimensional image of a real world or of an artificial object—the image can be of any shape. The physical domain may be a landform, like a mountain, or a lake or a part of a landform, a planet, or a part of a planet, a solar system, or a part of a solar system, a cluster of solar systems, or a part of a cluster of solar systems, a galaxy, or a part of a galaxy, a cluster of galaxies, or a part of a cluster of galaxies.

The invention provides an efficient method for determining an analytical representation of the internal structure and spatial properties distribution of various materials such as composite, granular, and porous materials that have been loaded under various conditions as well as any physical field such as electric and/or magnetic fields, temperature, etc. In particular, an example of these materials are carbon epoxy composites for which is identified the parameters of the d-dimensional W-M fractal function from measured data. The invention provides a volume parameterization as a multivariate analytic generalization of its uni-variate version.

The invention provides exact and full characterization of d-dimensional fractal volumes. Applications that require or are facilitated by having such an analytical representation of internal microstructure or spatial distribution of inhomogeneous material properties are enabled by this approach. The invention can be extended to the use of a homogeneous mesh for performing finite element analysis where the properties are represented by a W-M function, thereby bypassing the need for geometric tailoring of the mesh such as it explicitly represents the volumetric structure and replaces it with a generalized W-M fractal function.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a graph showing the uniform sampling of a sphere (dots) and four typical sampling vectors;

FIG. 2 is a a block diagram showing the forward problem perspective and the inverse problem perspective according to the invention;

FIG. 3 is a pseudo-colored isosurface combined with associated density plots of synthetic and inversely identified volumes of a W-M fractal according to the invention, demonstrating the success of the algorithm in this patent;

FIG. 4 is a graph showing the mean error of inverse identification through singular value decomposition vs parameter γ according to the invention;

FIG. 5 compares noisy data with respective identified fractal volumes according to the invention in order to demonstrate the robustness against noise;

FIG. 6 is a graph showing the result of noise in the error of the fractal inverse identification from synthetic data according to the invention;

FIG. 7 shows 8 of the 243 micro tomography slices according to the invention;

FIG. 8 compares original and inversely identified micro-tomographic volumes according to the invention in order to demonstrate the success of the method on real data; and

FIG. 9 is a graph showing the mean absolute error of the inversion as a function of γ according to the invention.

DETAILED DESCRIPTION OF THE INVENTION

d-Dimensional Modeling

It has been established that a fractal surface represents a very rich analytical representation that can model surface topographies of materials at small scales. This seems to be true especially because of the properties of continuity, non-differentiability and self affinity of specific types of fractals, that are also desired properties of surface topographies. The surface power spectra obeys a power-law relationship over a wide range of frequencies, because the surface topographies resemble a random process. Such a surface can be represented in polar coordinates (r, θ) by a complex function W as

$\begin{array}{cc}W\left(r,\theta \right)={L\left(\frac{G}{L}\right)}^{D-2}\sqrt{\frac{\ln \gamma }{M}}\sum _{m=1}^M\sum _{n=0}^N{\gamma }^{\left(D-3\right)n}\left(1-{}^{\frac{2\pi }{L}{\gamma }^nr\cos \left(\theta -a_m\right)}\right){}^{{\phi }_{\mathrm{mn}}},& \left(1\right)\end{array}$

where L is the linear size of the domain under consideration, G is the fractal roughness, D is the fractal dimension, φ_mn is a table of (usually) random phases (with φ_mnε[0,2π]), γ is a parameter that controls the density of the frequencies, M is the number of superimposed ridges and N is a number that is chosen based on the desired highest sampling frequency. The following equation can be used to identify an appropriate value for N,

N=n _max =int[log(L/Lc)/log γ],  (2)

with L_c being a cut-off wavelength, typically defined either by the highest sampling frequency, or by a physical barrier.

It should be noted that the parameter γ must be greater than 1 and it usually takes values in the vicinity of 1.5 because of surface flatness and frequency distribution density considerations. The latter though has been recently debated and only the requirement γ>1 was considered as a valid assumption.

Although the above mentioned representation is very useful in various physical contexts, its generalization to higher dimensions can provide a convenient representation of fields in those dimensions. Eq. (1) can be generalized to d-dimensions by the formula

$\begin{array}{cc}W\left(r\right)==\sqrt{\frac{\gamma }{M}}\sum _{m=1}^M\sum _{n=-\infty }^{\infty }A_m\left(1-{}^{k_0{\gamma }^nn_m·r}\right){{}^{{\phi }_{\mathrm{mn}}}\left(k_0{\gamma }^n\right)}^{D-\left(d+1\right)},& \left(3\right)\end{array}$

where n_m are unit vectors spaced uniformly over a unit d-dimensional hypersphere. k₀ can be substituted by 2 π/L, parameter A can be substituted by 2π(2 π/G)^2-D and n can be bound between 0 and N. Then (3) can be written as:

$\begin{array}{cc}W\left(r\right)={L\left(\frac{G}{L}\right)}^{D-2}\sqrt{\frac{\ln \gamma }{M}}\sum _{m=1}^M\sum _{n=0}^N{\gamma }^{n\left(D-\left(d+1\right)\right)}\left(1-{}^{\frac{2\pi }{L}{\gamma }^nn_m·r}\right){}^{{\phi }_{\mathrm{mn}}}.& \left(4\right)\end{array}$

The construction of the unit vectors n_m can be achieved with the help of the hyperspherical coordinates

$\begin{array}{cc}\begin{array}{c}{x}_1=r\cos {\theta }_1\\ x_2=r\sin {\theta }_1\cos {\theta }_2\\ x_3=r\sin {\theta }_1\sin {\theta }_2\cos {\theta }_3\\ ⋮\\ x_d=r\sin {\theta }_1\sin {\theta }_2\dots \sin {\theta }_{d-2}\sin {\theta }_{d-1}\end{array}\},& \left(5\right)\end{array}$

where θ_i are the angular coordinates. For unit vectors, one can use r=1. It should be noted that in order to avoid traversing the hypersphere twice we can choose to traverse θ₁ . . . θ_d-2 over the range [0,π]. A uniform sampling of the hypersphere can be achieved by setting,

$\begin{array}{cc}\begin{array}{c}{\theta }_1={\theta }_2=\dots ={\theta }_{d-2}=\frac{\left(m-1\right)\pi }{M-1}\\ {\theta }_{d-1}=2p\frac{\left(m-1\right)\pi }{M-1}\end{array}\},m=1\dots M,& \left(6\right)\end{array}$

where pεN−{0} defines the density of the winding around the surface of the hypersphere. For p=1 the hypershpere will be sampled with only a single traversal, while for higher numbers of p, denser windings can be achieved. Using (5) and (6) the unit vector sampling can be written as:

$\begin{array}{cc}\begin{array}{cc}{n}_m=& {\left\{x_1^m,x_2^m,\dots ,x_d^m\right\}}^T=\\ =& {\left\{\begin{array}{c}\cos \left[\frac{\left(m-1\right)\pi }{M-1}\right],\\ ,\sin \left[\frac{\left(m-1\right)\pi }{M-1}\right]\cos \left[\frac{\left(m-1\right)\pi }{M-1}\right],\\ ,{\sin }^2\left[\frac{\left(m-1\right)\pi }{M-1}\right]\cos \left[\frac{\left(m-1\right)\pi }{M-1}\right],\\ ⋮\\ ,{\sin }^{d-2}\left[\frac{\left(m-1\right)\pi }{M-1}\right]\sin \left[2p\frac{\left(m-1\right)\pi }{M-1}\right]\end{array}\right\}}^T\end{array}\}.& \left(7\right)\end{array}$

An example sampling of a 3-sphere is shown in FIG. 1 for p=8 and M=64.

The Inverse Problem

The systemic view of both the forward and the inverse problems are shown in FIG. 2. In the case of the inverse problem, it is shown that its solution produces the parameters of the fractal function provided that field distribution of the volume z(r) has been experimentally determined and therefore is known.

More specifically, given an array of measurements z_ijk^e, i=1 . . . I, j=1 . . . J, k=1 . . . K over a region of size L^d one may then identify the parameters γ and φ_mn of a d-dimensional volume z(r) that best fits those measurements. In previous studies the identification of those parameters was done only for the 2 Dimensional surface. In order to achieve the inverse identification of the volume parameters the approach is similar to that described in Michopoulos et al.

The inner sum in Eq. 4 can be written as,

$\begin{array}{cc}\sum _{n=0}^N{\gamma }^{n\left(D-\left(d+1\right)\right)}\left(1-{}^{\frac{2\pi }{L}{\gamma }^nn_m·r}\right){}^{{\phi }_{\mathrm{mn}}}={c_m\left(r\right)}^T{\phi }_m,& \left(8\right)\end{array}$

where c_m(r) and φ_m are vectors in P^(N+1) given by,

$\begin{array}{cc}{c}_m\left(r\right)={\left\{c_{m0},c_{m1},\dots ,c_{\mathrm{mN}}\right\}}^T={\left\{\begin{array}{c}{\gamma }^{\left(D-\left(d+1\right)\right)0}\left(1-{}^{\frac{2\pi }{L}{\gamma }^0n_m·r}\right),\\ {\gamma }^{\left(D-\left(d+1\right)\right)1}\left(1-{}^{\frac{2\pi }{L}{\gamma }^1n_m·r}\right),...\\ {\gamma }^{\left(D-\left(d+1\right)\right)N}\left(1-{}^{\frac{2\pi }{L}{\gamma }^Nn_m·r}\right)\end{array}\right\}}^T,& \left(9\right)\\ \mathrm{and},{\phi }_m={\left\{{}^{{\mathrm{\phi }}_{m0}},{}^{{\mathrm{\phi }}_{m2}},\dots ,{}^{{\mathrm{\phi }}_{mN}}\right\}}^T.& \left(10\right)\end{array}$

By setting,

$\begin{array}{cc}{\hspace{0.17em}}^{*}Q={L\left(\frac{G}{L}\right)}^{D-2}\sqrt{\frac{\ln \gamma }{M}},& \left(11\right)\end{array}$

Eq. 4 can be written as:

$\begin{array}{cc}W\left(r\right)={\hspace{0.17em}}^{*}Q\sum _{m=1}^Mc_m^T{\phi }_m.& \left(12\right)\end{array}$

By coalescing the vectors c_m and φ_m into larger vectors in P^M(N+1) such as,

c(r)={c₁^T,c₂^T, . . . ,c_M^T}^T,  (13)

and

φ={φ₁^T,φ₂^T, . . . ,φ_M^T}^T,  (14)

Eq. 12 can be written as:

W(r)=*Qc(r)^Tφ.  (15)

For the needs of the inverse problem it is assumed that a number of measurements at points r_t^e exist for a volume represented as z^e(r_t^e), t=1 . . . T, T≧M(N+1). We seek to identify a volume that is described by Eq. 15 and approximates the experimental points z^e(r_t^e).

To solve this problem we first form the following linear system,

$\begin{array}{cc}\left[\begin{array}{c}W\left(r_1\right)\\ W\left(r_2\right)\\ \begin{array}{c}⋮\\ W\left(r_T\right)\end{array}\end{array}\right]=\left[\begin{array}{c}{z}^{}\left(r_1\right)\\ z^{}\left(r_2\right)\\ \begin{array}{c}⋮\\ z^{}\left(r_T\right)\end{array}\end{array}\right],\mathrm{or},& \left(16\right)\\ {\hspace{0.17em}}^{*}Q\left[\begin{array}{c}{c\left(r_1\right)}^T\\ {c\left(r_2\right)}^T\\ \begin{array}{c}\dots \\ {c\left(r_T\right)}^T\end{array}\end{array}\right]\left[\begin{array}{c}{\phi }_1\\ {\phi }_2\\ \dots \\ {\phi }_M\end{array}\right]=\left[\begin{array}{c}{z}^{}\left(r_1\right)\\ z^{}\left(r_2\right)\\ \begin{array}{c}⋮\\ z^{}\left(r_T\right)\end{array}\end{array}\right].& \left(17\right)\end{array}$

By expanding the vectors in Eq. 16, the system can be written as,

$\begin{array}{cc}\mathrm{Cp}=z,\mathrm{with},C={\hspace{0.17em}}^{*}Q\left[\begin{array}{c}{c}_{11}\left(r_1\right),c_{12}\left(r_1\right),\dots ,c_{21}\left(r_1\right),\dots ,c_{\mathrm{MN}}\left(r_1\right)\\ c_{11}\left(r_2\right),c_{12}\left(r_2\right),\dots ,c_{21}\left(r_2\right),\dots ,c_{\mathrm{MN}}\left(r_2\right)\\ ⋮\\ c_{11}\left(r_T\right),c_{12}\left(r_T\right),\dots ,c_{21}\left(r_T\right),\dots ,c_{\mathrm{MN}}\left(r_T\right)\end{array}\right],,p={\left[{\phi }_{11},{\phi }_{12},\dots ,{\phi }_{1N},{\phi }_{21},\dots ,{\phi }_{\mathrm{MN}}\right]}^T\mathrm{and}z=\left[\begin{array}{c}{z}^{}\left(r_1\right)\\ z^{}\left(r_2\right)\\ \begin{array}{c}⋮\\ z^{}\left(r_T\right)\end{array}\end{array}\right].& \left(18\right)\end{array}$

The system of Eq. 17 is an overdetermined system of M(N+1) equations. Since the right hand side vector z contains experimental measurements, it also contains noise; the system cannot, in general, have an exact solution. Nevertheless, one can seek a p, such as Πp-zΠ is minimized, where Π_pΓ. is the vector norm. Such a p is known as the least squares solution to the over-determined system. It should be noted that the left hand side expression of Eq. 16 yields results in the complex domain, but as long as a minimal solution is achieved for real numbers on the right hand side, the imaginary parts will be close to 0. A solution can be given by the following equation,

p=Vy,  (19)

where V is calculated by the Singular Value Decomposition (SVD) of C as,

UDV ^T =C,  (20)

where y is a vector defined as y_t=b_t′/d_t, b=b_t is a vector given by,

b=U ^T p),  (21)

and d_i is the i th entry of the diagonal of D.

The solution of the inverse problem as described by the overdetermined system of Eq. 16 gives the phases φ_mn given known values of the other parameters. In a general volume the only other parameter that is unknown is γ. As was pointed out in Michopoulos et al parameters G and D don't need to be considered as unknowns to be determined in this optimization. This is because for any combination of the phases it is always possible to find new values for φ_mn, that result in generating the same surface.

Numerical Experiments

In order to assess the quality of the high dimensional fractal characterization results of the numerical examples that follow, we define the following error function,

$\begin{array}{cc}\mathrm{error}=\frac{\sum _{i=1}^P\frac{z_i^{}-z_i^d}{{\max \left(z_i\right)}^d-\min \left(z_i^d\right)}}{P},& \left(22\right)\end{array}$

where z_i^e is the value of the scalar field of the experimentally measured (or reference) points, P is the number of those points. In the following examples P is set as P=I×J×K=40×40×40=64000 points. z_i^d is the value of the scalar field of the inversely identified fractal and is equal to the real part of the truncated W-M function of Eq. (4) (z_i^d=Re{W(r)}).

TABLE 1
PARAMETERS OF THE REFERENCE MODEL.
Parameter Value
L         1
G         1 × 106
γ         1.5
D         2.3
M         10
N         16
d         3

To study the feasibility of the invention, a few numerical experiments were designed in the 3D space. The first experiment was based on synthetic data and is aimed at inversely identifying only the phases of a volume constructed by the fractal itself. The original volume is shown in FIG. 3 and was constructed using Eq. (4) with random phases and the parameters as shown in Table 1. The inversely identified volume using the phases resulting from Eq. (19) is presented in FIG. 3. It should be noted that this difference is very small compared to the magnitude of the volume and any discrepancies should be considered as the numerical error of the SVD algorithm.

The second synthetic experiment involved the identification of both the phases and the γ parameter. An exhaustive search approach was adopted in this case, as the sensitivity of the SVD inversion relative to the value of γ is also of interest. For a range of the possible values for parameter γ the inversion of the phases was executed and the value of the error function (Eq. 22) was calculated. The error for various values of γ is presented on FIG. 4. The smallest value for the error is at γ=1.5, which is the one used originally for the generation of the volume (Table 1).

In order to investigate the behavior of the invention under noisy data, we also applied the inversion in the original data of FIG. 3, but for various levels of injected noise. The noise levels varied relative to the range of the scalar field values from 0% to 300%. The results of the inverse identification for noise levels 20%, 40% and 300% are presented in FIG. 5 and show impressively robust response. Even for noise levels of 300% (FIG. 5), where the volume characteristics are visually non-existent, the fractal inversion methodology resembles (FIG. 5) the actual field very well. In FIG. 6 the mean absolute error of the identification for various levels of noise is presented, where it can be observed that as the noise levels increase, the mean absolute error increases linearly.

Although the previous analysis demonstrates the consistency and performance of the invention, it is much more useful when applied to actual data. For this reason, a numerical test is performed based on micro-tomography scans of a failed AS4/3501-6 composite specimen used for 3 degrees of freedom in-plane mechanical loading. AS4/3501-6 is a carbon/epoxy composite used in general purpose structural applications and the associated prepreg is manufactured by Hexcel Corporation. The data set used consisted of 243 tomographic slices of the specimen stored as digital images of resolution 1256×588 pixels each.

Because the data set is very large, the volume characterization was done on a subset of the volumetric data shown in FIG. 7. The subset size was 16×16×16 voxels and was taken from the region depicted by the red square in FIG. 8. The fractal parameters were chosen to be: γ=1.5,D=2.3 μL=1,G=1*10⁶ μM=80 and N=50. The inversion of the volume gives very satisfactory results as shown in FIG. 8. For this test the mean absolute error as calculated by Eq. (22) was 4.74%.

To investigate the effect the choice of γ has on the fractal inversion, a parametric study was performed for various values of γ ranging from 1.02 to 2. For each of those values the fractal phases where identified and the values of the mean absolute error from Eq. (22) were calculated. The results of this study are plotted in FIG. 9. Although the original choice of γ in the previous example resulted in a very satisfactory inverse identification of the failed composite specimen volume, it can be seen from FIG. 9, that even better results can be achieved for a value of γ equal to 1.24. In that particular case the mean absolute error of the identified volume relative to the original was about 0.84% and the resulting volume is depicted in FIG. 7. This result underlines that for optimum characterization purposes on has to consider that the set of parameters to be identified should include both the phases and γ.

It should be noted that the present invention can be accomplished by executing one or more sequences of one or more computer-readable instructions read into a memory of one or more computers from volatile or non-volatile computer-readable media capable of storing and/or transferring computer programs or computer-readable instructions for execution by one or more computers. Volatile computer readable media that can be used can include a compact disk, hard disk, floppy disk, tape, magneto-optical disk, PROM (EPROM, EEPROM, flash EPROM), DRAM, SRAM, SDRAM, or any other magnetic medium; punch card, paper tape, or any other physical medium. Non-volatile media can include a memory such as a dynamic memory in a computer. In addition, computer readable media that can be used to store and/or transmit instructions for carrying out methods described herein can include non-physical media such as an electromagnetic carrier wave, acoustic wave, or light wave such as those generated during radio wave and infrared data communications.

While specific embodiments of the present invention have been shown and described, it should be understood that other modifications, substitutions and alternatives are apparent to one of ordinary skill in the art. Such modifications, substitutions and alternatives can be made without departing from the spirit and scope of the invention, which should be determined from the appended claims.

Claims

1. A computer implemented method for obtaining an analytical representation of an internal structure and spatial properties distribution of a selected physical domain, comprising: identifying d-dimensional correspondences of measured spatial properties or field distributions; and applying an inverse algorithm to the d-dimensional spatial properties or field distributions to calculate the Weierstrass-Mandelbrot (W-M) fractal model to thereby determine parameters defining an analytical and continuous Weierstrass-Mandelbrot (W-M) representation.

2. The method of claim 1, wherein the physical domain is selected from a metal alloy, a carbon epoxy composite, or a biological tissue.

3. The method of claim 1, wherein the physical domain is a man-made structure selected from a microprocessor, an aircraft, a ship, an automobile, a darn, a road layer, a bridge, or a building.

4. The method of claim 1, wherein the physical domain is a live organism or part of a live organism or a in vitro version of them (e.g wood).

5. The method of claim 1, wherein the physical domain is a digital representation of a real or artificial object.

6. The method of claim 5, wherein the digital representation is a spatial distribution of a scalar field within a volume such as magnetic resonance slice images, x-ray tomography data, ultrasonic data trough a volume, or any simulated entity distributed in a volume.

7. The method of claim 1, wherein the physical domain is a landform selected from a mountain, a lake, part of a landform, a planet, part of a planet, a solar system, part of a solar system, a cluster of solar systems, part of a cluster of solar systems, a galaxy, part of a galaxy, a cluster of galaxies, or part of a cluster of galaxies.

8. The method of claim 1, wherein the d-dimensional modeling is obtained by applying a formula

9. The method of claim 8, wherein based on an array of measurements over a selected region of the domain, parameters of the d-dimensional W-M fractal model are identified that best fit the array of measurements by decomposing the problem so as to apply singular value decomposition for determining unknown phases of the W-M fractal model.

10. A computer software product comprising a physical computer-readable medium including stored instructions that, when executed by a computer, cause the computer to: identify d-dimensional correspondences of measured spatial properties or field distributions; andapply an inverse algorithm to the d-dimensional spatial properties or field distributions to calculate the Weierstrass-Mandelbrot (W-M) fractal model to thereby determine parameters defining an analytical and continuous Weierstrass-Mandelbrot (W-M) representation.

11. The computer software product of claim 10, wherein the physical domain is selected from a metal alloy, a carbon epoxy composite, or a biological tissue.

12. The computer software product of claim 10, wherein the physical domain is a man-made structure selected from a microprocessor, an aircraft, a ship, an automobile, a dam, a road layer, a bridge, or a building.

13. The computer software product of claim 10, wherein the physical domain is a live organism or part of a live organism or a in vitro version of them (e.g wood).

14. The computer software product of claim 10, wherein the physical domain is a digital representation of a real or artificial object.

15. The computer software product of claim 14, wherein the digital representation is a spatial distribution of a scalar field within a volume such as magnetic resonance slice images, x-ray tomography data, ultrasonic data trough a volume, or any simulated entity distributed in a volume.

16. The computer software product of claim 10, wherein the physical domain is a landform selected from a mountain, a lake, part of a landform, a planet, part of a planet, a solar system, part of a solar system, a cluster of solar systems, part of a cluster of solar systems, a galaxy, part of a galaxy, a cluster of galaxies, or part of a cluster of galaxies.

17. The computer software product of claim 10, wherein the d-dimensional modeling is obtained by applying a formula

18. The computer software product of claim 17, wherein based on an array of measurements over a selected region of the domain, parameters of the d-dimensional W-M fractal model are identified that best fit the array of measurements by decomposing the problem so as to apply singular value decomposition for determining unknown phases of the W-M fractal model.