Hi(gh)!

By request in the process of debugging, here the current version of the 
script file mentioned in the subject:

// begin of code

#version 3.6;

#include "functions.inc"

#declare rds=false;

global_settings
{
   assumed_gamma 2.2
   #if (rds = true)
     #declare amb=0;
     radiosity { } // most basic radiosity, to be detailed in the future!
   #else
     #declare amb=0.05;
   #end
}



#declare path="Rohdaten/"; // path for ASCII file containg elevation data
#declare path_textures="Texturen/";


#declare xdim=1000; // full resolution: 3601;
#declare ydim=1000; // full resolution: 3601;

#declare rd=6378140-422;
#declare rng=9270;
#declare flattening=0;

#declare latstart=41;
#declare latend=42;
#declare longstart=44;
#declare longend=45;
#declare latsize=latend-latstart;
#declare longsize=longend-longstart;
#declare month=7;

#declare res=1000; // full resolution: 3601;

#if (latstart-1 < 0)
   #declare lat="s";
#else
   #declare lat="n";
#end

#if (longstart < 0)
   #declare long="w";
#else
   #declare long="e";
#end

#declare tilename=concat(lat, str(latstart, -2, 0), long, str(longstart, 
-3, 0));

// tiles are named by the coordinates of their southwestern corner!

#declare altdataname = concat(path, tilename, ".txt");

#declare NumVertices = xdim*ydim;
#declare NumFaces = (xdim-1)*(ydim-1)*2;
#declare NumNormals = NumVertices;
#declare NumTextures = NumVertices;
#declare V_vec_Arr = array[NumVertices];
#declare Face_Arr = array[NumFaces];
#declare N_vec_Arr = array[NumNormals];
#declare T_Arr = array[NumTextures];

#declare a = 0;

#declare P_Heightfield=
pigment
{
   image_map
   {
     "Heightfields/n41e044_1000x1000.png"
     interpolate 2
   }
}

#declare P_Texture=
pigment
{
   image_map
   {
     "Texturen/n41e044_2004-07.png"
     interpolate 2
   }
}


#declare a=0;
#declare c=0;
#while (a < ydim)
   #declare b=0;
   #while (b < xdim)
     #declare hcolor = eval_pigment(P_Heightfield, <(0.5+b)*(1/xdim), 
(0.5+a)*(1/ydim), 1>);
     #declare hval = hcolor.red*rng+hcolor.green*(rng/255);
     #declare V_vec_Arr[c] = 
(rd+hval)*<sin(radians(-longstart-(b/xdim)*longsize))*(1-flattening)*cos(radians(latstart+(a/ydim)*latsize)),
		(1-flattening)*sin(radians(latstart+(a/ydim)*latsize)), 
cos(radians(-longstart-(b/xdim)*longsize))*(1-flattening)*cos(radians(latstart+(a/ydim)*latsize))>;
     #declare c=c+1;		
     #declare b=b+1;
   #end
   #declare a=a+1;
#end

#declare a = 0;
#declare c = 0;
#while (a < ydim-1)
   #declare b = 0;
   #while (b < xdim-1)
     #declare Face_Arr[c] = <a*res+b, a*res+b+1, (a+1)*res+b>;
     #declare Face_Arr[c+1] = <a*res+b+1, (a+1)*res+b+1, (a+1)*res+b>;
     #declare c = c+2;
     #declare b = b+1;
   #end
   #declare a = a+1;
#end

#declare a=0; // vcross(C-B,A-B)
#declare c=0;
#while (a < ydim)
   #declare b=0;
   #while (b < xdim)
     #if (a=0 & b=0) // lower left corner: 1 adjacent face
       #declare nno = vcross(V_vec_Arr[xdim]-V_vec_Arr[1], 
V_vec_Arr[0]-V_vec_Arr[1]);
     #else
       #if (a=0 & b>0 & b<xdim-1) // row 0 except lower left and lower 
right corner: 3 adjacent faces per vertex
	#declare nno = vcross(V_vec_Arr[(a+1)*xdim+(b-1)]-V_vec_Arr[b], 
V_vec_Arr[b-1]-V_vec_Arr[b]);
	#declare nno = nno + 
vcross(V_vec_Arr[(a+1)*ydim+(b-1)]-V_vec_Arr[(a+1)*xdim+b], 
V_vec_Arr[b]-V_vec_Arr[(a+1)*xdim+b]);
	#declare nno = nno + 
vcross(V_vec_Arr[(a+1)*ydim+b]-V_vec_Arr[(a+1)*xdim+b+1], 
V_vec_Arr[b]-V_vec_Arr[(a+1)*xdim+b+1]);
	#declare nno = nno/3;
       #else
	#if (a=0 & b=xdim-1) // lower right corner: 2 adjacent faces
	  #declare nno = vcross(V_vec_Arr[xdim*2-2]-V_vec_Arr[xdim-1], 
V_vec_Arr[xdim-2]-V_vec_Arr[xdim-1]);
	  #declare nno = nno + vcross(V_vec_Arr[xdim*2-2]-V_vec_Arr[xdim*2-1], 
V_vec_Arr[xdim-1]-V_vec_Arr[xdim*2-1]);
	  #declare nno = nno/2;
	#else
	  #if (a>0 & a<ydim-1 & b=0) // leftmost column except lower left and 
upper left corner: 3 adjacent faces per vertex
	    #declare nno = vcross(V_vec_Arr[a*ydim]-V_vec_Arr[(a-1)*ydim+1], 
V_vec_Arr[(a-1)*ydim]-V_vec_Arr[(a-1)*ydim+1]);
	    #declare nno = nno + 
vcross(V_vec_Arr[a*ydim+1]-V_vec_Arr[(a-1)*ydim+1], 
V_vec_Arr[a*ydim]-V_vec_Arr[(a-1)*ydim+1]);
	    #declare nno = nno + 
vcross(V_vec_Arr[(a+1)*ydim]-V_vec_Arr[a*ydim+1], 
V_vec_Arr[a*ydim]-V_vec_Arr[a*ydim+1]);
	    #declare nno = nno/3;
	  #else
	    #if (a>0 & a<ydim-1 & b>0 & b<xdim-1) // all rows and columns not 
touching any edge: 6 adjacent faces per vertex
	      #declare nno = 
vcross(V_vec_Arr[(a-1)*xdim+b]-V_vec_Arr[a*xdim+b-1], 
V_vec_Arr[a*xdim+b]-V_vec_Arr[a*xdim+b-1]);
	      #declare nno = nno + 
vcross(V_vec_Arr[(a-1)*xdim+b+1]-V_vec_Arr[(a-1)*xdim+b], 
V_vec_Arr[a*xdim+b]-V_vec_Arr[(a-1)*xdim+b]);
	      #declare nno = nno + 
vcross(V_vec_Arr[a*xdim+b+1]-V_vec_Arr[(a-1)*xdim+b+1], 
V_vec_Arr[a*xdim+b]-V_vec_Arr[(a-1)*xdim+b+1]);
	      #declare nno = nno + 
vcross(V_vec_Arr[(a+1)*xdim+b]-V_vec_Arr[a*xdim+b+1], 
V_vec_Arr[a*xdim+b]-V_vec_Arr[a*xdim+b+1]);
	      #declare nno = nno + 
vcross(V_vec_Arr[(a+1)*xdim+b-1]-V_vec_Arr[(a+1)*xdim+b], 
V_vec_Arr[a*xdim+b]-V_vec_Arr[(a+1)*xdim+b]);
	      #declare nno = nno + 
vcross(V_vec_Arr[a*xdim+b-1]-V_vec_Arr[(a+1)*xdim+b-1], 
V_vec_Arr[a*ydim+b]-V_vec_Arr[(a+1)*xdim+b-1]);
	      #declare nno = nno/6;
	    #else
	      #if (a>0 & a<ydim-1 & b=ydim-1) // rightmost column except lower 
right and upper right corner: 3 adjacent faces per vertex
		#declare nno = 
vcross(V_vec_Arr[(a-1)*xdim+(xdim-1)]-V_vec_Arr[a*xdim+(xdim-2)], 
V_vec_Arr[a*xdim+(xdim-1)]-V_vec_Arr[a*xdim+(xdim-2)]);
		#declare nno = nno + 
vcross(V_vec_Arr[a*xdim+(xdim-2)]-V_vec_Arr[(a+1)*xdim+(xdim-2)], 
V_vec_Arr[a*xdim+(xdim-1)]-V_vec_Arr[(a+1)*xdim+(xdim-2)]);
		#declare nno = nno + 
vcross(V_vec_Arr[(a+1)*xdim+(xdim-2)]-V_vec_Arr[(a+1)*xdim+(xdim-1)], 
V_vec_Arr[a*xdim+(xdim-1)]-V_vec_Arr[(a+1)*xdim+(xdim-1)]);
		#declare nno = nno/3;
	      #else
		#if (a=ydim-1 & b=0) // upper left corner: 2 adjacent faces
		  #declare nno = vcross(V_vec_Arr[(a-1)*xdim+1]-V_vec_Arr[(a-1)*xdim], 
V_vec_Arr[a*xdim]-V_vec_Arr[(a-1)*xdim]);
		  #declare nno = nno + 
vcross(V_vec_Arr[a*xdim+1]-V_vec_Arr[(a-1)*xdim+1], 
V_vec_Arr[a*xdim]-V_vec_Arr[(a-1)*xdim+1]);
		  #declare nno = nno/2;
		#else
		  #if (a=ydim-1 & b>0 & b<xdim-1) // row 3600 except upper left and 
upper right corner: 3 adjacent faces per vertex
		    #declare nno = 
vcross(V_vec_Arr[(a-1)*xdim+b]-V_vec_Arr[a*xdim+b-1], 
V_vec_Arr[a*xdim+b]-V_vec_Arr[a*xdim+b-1]);
		    #declare nno = nno + 
vcross(V_vec_Arr[(a-1)*xdim+b+1]-V_vec_Arr[(a-1)*xdim+b], 
V_vec_Arr[a*xdim+b]-V_vec_Arr[(a-1)*xdim+b]);
		    #declare nno = nno + 
vcross(V_vec_Arr[a*xdim+b+1]-V_vec_Arr[(a-1)*xdim+b+1], 
V_vec_Arr[a*xdim+b]-V_vec_Arr[(a-1)*ydim+b+1]);
		    #declare nno = nno/3;
		  #else
		    #if (a=ydim-1 & b=xdim-1) // upper right corner, 1 adjacent face
		      #declare nno = 
vcross(V_vec_Arr[xdim*ydim-xdim]-V_vec_Arr[xdim*ydim-2], 
V_vec_Arr[xdim*ydim-1]-V_vec_Arr[xdim*ydim-2]);
		    #end
		  #end
		#end
	      #end
	    #end
	  #end
	#end
       #end
     #end
     #warning str(c, 1, 0)
     #declare N_vec_Arr[c] = vnormalize(nno);
     #declare b=b+1;
     #declare c=c+1;
   #end
   #declare a=a+1;
#end



#declare 
lNormVect=vnormalize(<sin(radians(-longstart-0.5))*cos(radians(latstart+0.5)), 
sin(radians(latstart+0.5)),
cos(radians(-longstart-0.5))*cos(radians(latstart+0.5))>);
#declare ApproxLook = lNormVect * rd;

#fopen ES concat(tilename, ".inc") write
#write (ES, concat("#declare ApproxLook = <", vstr(3, 
ApproxLook,",",3,4), ">;\n"))
#write (ES, concat("#declare F_Earthslice = finish {ambient ",str(amb, 
1, 2)," diffuse 1 brilliance 0.25 }\n"))
#write (ES, concat("#declare Earth_Slice_",tilename,"=\nmesh2\n{\n 
vertex_vectors\n     {\n",str(xdim*ydim, 1, 0),"\n"))

#declare sf=1;

#declare i=0;
#while (i<NumVertices)
   #write (ES, concat("    <",vstr(3, V_vec_Arr[i]-ApproxLook, 
",",3,4),">"))
   #if (i=NumVertices-1)
     #write (ES, "\n")
   #else
     #write (ES, ",\n")
   #end
   #declare i=i+1;
#end
#write (ES,"    }\n")
#write (ES,"	normal_vectors\n")
#write (ES,"	{\n")
#write (ES, concat("    ",str(xdim*ydim, 1, 0),"\n"))
#declare i=0;
#while (i<NumNormals)
   #write (ES, concat("    <",vstr(3, N_vec_Arr[i], ",",3,4),">"))
   #if (i=NumNormals-1)
     #write (ES, "\n")
   #else
     #write (ES, ",\n")
   #end
   #declare i=i+1;
#end

#write (ES,"    }\n")
#if (!(tilename="n80e014" | tilename="n80e016" | tilename="n80e017"))
   #declare i=0;
   #write (ES, "  texture_list")
   #write (ES, "  {\n")
   #write (ES, concat("    ",str(NumTextures, 1, 0),"\n"))
   #declare a=0;
   // #declare c=0;
   #while (a<ydim)
     #declare b=0;
     #while (b<xdim)
       // #declare C_Texture = eval_pigment(P_Texture, 
<(0.5+b)*(1/(xdim-1)), (0.5+a)*(1/(ydim-1)), 0>);
       #declare C_Texture = color rgb <0.5, 0.5, 0.5>;
       #warning concat("Iteration ", str(a, 1, 0), "-", str(b, 1, 0))
       #write (ES, "    texture\n")
       #write (ES, "    {\n")
       #write (ES, "      pigment\n")
       #write (ES, "      {\n")
       #write (ES, concat("	   color rgb <",vstr(3, C_Texture, ",", 1, 
7),">\n"))
       #write (ES, "      }\n")
       #write (ES  "	 finish{ F_Earthslice }\n")
       #write (ES, "    }\n"
       // #declare c=c+1;
       #declare b=b+1;
     #end
     #declare a=a+1;
   #end
   #write (ES, "  }\n")
#end


#write (ES, concat("  face_indices\n{\n    ",str(NumFaces, 7, 0),"\n"))
#declare i=0;
#while (i<NumFaces)
   #write (ES, concat("    <",vstr(3, Face_Arr[i], ",",0, 
0),">,",str(div(i,2), 1, 0)))
   #if (i=NumFaces-1)
     #write (ES, "\n")
   #else
     #write (ES, ",\n")
   #end
   #declare i=i+1;
#end
#write (ES, "  }\n")




/* #while (i < NumTextures)
   #write (ES, concat("texture { T_Arr[",str(i, 1, 0)))
   #if (i=NumTextures-1)
     #write (ES, "\n")
   #else
     #write (ES, ",\n")
   #end
   #declare i=i+1;
#end */

#if (tilename="n80e014" | tilename="n80e016" | tilename="n80e017")
   #write (ES, "  texture\n")
   #write (ES, "  {\n")
   #write (ES, "    pigment { color rgb <0.66, 0.65, 0.65> }\n")
   #write (ES, "    finish { ambient 0.1 diffuse 1 brilliance 0.4 }\n")
   #write (ES, "  }\n")
#end
#write (ES, "  double_illuminate\n")
#write (ES, "  translate ApproxLook\n")
#write (ES, "}\n")
#fclose ES

// end of code

See you in Khyberspace!

Yadgar

Post a reply to this message

> By request in the process of debugging, here the current version of the
> script file mentioned in the subject:

A few points of observation:

1. _attach_ the scene as a .pov file to avoid the split-line / wrapping issue

2. You load the same image twice
     #declare P_Heightfield=
     #declare P_Texture=
why not just reuse the same image_map pigment pattern?

3. THEN, you scan the entire image TWICE, which has to take forever
#declare hcolor = eval_pigment(P_Heightfield, <(0.5+b)*(1/xdim),
(0.5+a)*(1/ydim), 1>);
// #declare C_Texture = eval_pigment(P_Texture, <(0.5+b)*(1/(xdim-1)),
(0.5+a)*(1/(ydim-1)), 0>);
and why is one calculation using xdim & ydim, and the other using xdim-1 &
ydim-1 ?

4. it's presently unclear to me how your calculations are working, but it does
seem to me that you are working on the images in the default unit square pigment
pattern.   You should get rid of

#declare xdim=1000; // full resolution: 3601;
#declare ydim=1000; // full resolution: 3601;

and use something like
#declare ImageMap = pigment {image_map {"yourimage.jpg" once} };
#declare Resolution = max_extent (ImageMap);
#declare Resolution = Resolution + <0, 0, 1>;

in order to automatically give you the correct number of pixels in the image.
then if you

scale Resolution

You will have an image_map pigment pattern where you can use integers in your
loops for pixel positions.


5.  I think you're breaking things up into several sequential loops where there
are a lot of redundant things happening, and you would do far better to perform
all of your operations in a single loop.
You store a lot of values in large arrays, but I'm not sure that if it were all
in one loop that you would need to store all of those values.  And of course
that would free up a LOT of memory.

Post a reply to this message

=?UTF-8?Q?J=c3=b6rg_=22Yadgar=22_Bleimann?= <yaz### [at] gmxde> wrote:
> Hi(gh)!
>
> By request in the process of debugging, here the current version of the
> script file mentioned in the subject:
>
> // begin of code
> ...
> // end of code

Yadgar,

Your array V_vec_Arr is a 1 dimensional array.
But it seems like you are using it as if it was a 2 dimensional array.

If so, then you could have defined and used it like this:

#declare V_vec_Arr = array[ydim][xdim];
V_vec_Arr[a][b]

- but instead you seem to define and use it like this:

#declare V_vec_Arr = array[ydim*xdim];
V_vec_Arr[a*xdim + b]

This makes the indexing more awkward and error prone than it needs to be.

E.g.:

You have several lines that look somewhat like this:

#declare nno =
    nno +
    vcross(
        V_vec_Arr[(a + 1)*ydim + (b - 1)] - V_vec_Arr[(a + 1)*xdim + b],
        V_vec_Arr[                b     ] - V_vec_Arr[(a + 1)*xdim + b]
    )
;

- where the "*ydim" part seems wrong to me.

And then you have a line like this:

#if (a > 0 & a < ydim - 1 & b = ydim - 1) // rightmost column except lower

- where I suspect that "b = ydim - 1" should be "b = xdim - 1".

If these are bugs, they will probably not bite you until xdim becomes different
from ydim.

So I suggest that you do your array indexing in a different way.

--
Tor Olav
http://subcube.com
https://github.com/t-o-k

Post a reply to this message

"Tor Olav Kristensen" <tor### [at] TOBEREMOVEDgmailcom> wrote:

> Yadgar,
[... lots of excellent points ...]
> So I suggest that you do your array indexing in a different way.


TOK has spoken.   I can tell you from experience that listening and heeding his
advice is highly recommended.

Post a reply to this message

> Yadgar,
> 
> Your array V_vec_Arr is a 1 dimensional array.
> But it seems like you are using it as if it was a 2 dimensional array.
> 
> If so, then you could have defined and used it like this:
> 
> #declare V_vec_Arr = array[ydim][xdim];
> V_vec_Arr[a][b]
> 
> - but instead you seem to define and use it like this:
> 
> #declare V_vec_Arr = array[ydim*xdim];
> V_vec_Arr[a*xdim + b]
> 
> This makes the indexing more awkward and error prone than it needs to be.
> 
> E.g.:
> 
> You have several lines that look somewhat like this:
> 
> #declare nno =
>      nno +
>      vcross(
>          V_vec_Arr[(a + 1)*ydim + (b - 1)] - V_vec_Arr[(a + 1)*xdim + b],
>          V_vec_Arr[                b     ] - V_vec_Arr[(a + 1)*xdim + b]
>      )
> ;
> 
> - where the "*ydim" part seems wrong to me.
> 
> And then you have a line like this:
> 
> #if (a > 0 & a < ydim - 1 & b = ydim - 1) // rightmost column except lower
> 
> - where I suspect that "b = ydim - 1" should be "b = xdim - 1".
> 
> If these are bugs, they will probably not bite you until xdim becomes different
> from ydim.
> 
> So I suggest that you do your array indexing in a different way.
> 
> --
> Tor Olav
> http://subcube.com
> https://github.com/t-o-k
> 
> 
Arghhh! That DO hurt !
That kind of array indexing gymnastic is only acceptable if you are 
using a language that don't support multidimensional arrays.


Alain

Post a reply to this message