started to add curl calculation to postprocessingMath.f90

restructured colormap calculation. Algorithm is working now, but input/output still has to be done
This commit is contained in:
Martin Diehl 2011-11-23 20:08:48 +00:00
parent 81e9607316
commit 267b8ac30c
8 changed files with 335 additions and 275 deletions

View File

@ -0,0 +1,36 @@
#!/usr/bin/python
# -*- coding: iso-8859-1 -*-
def write_gsmh(RGB_vector,name):
colormap = open(str(name) + '.map',"w")
colormap.write('View.ColorTable = {\n')
for i in range(len(RGB_vector)-1):
colormap.write('{'+str((RGB_vector[0][i])*255.0)+','+str((RGB_vector[0][i])*255.0)+','+str((RGB_vector[0][i])*255.0)+'},\n')
colormap.write('{'+str((RGB_vector[0][-1])*255.0)+','+str((RGB_vector[0][-1])*255.0)+','+str((RGB_vector[0][-1])*255.0)+'}}')
file.close(colormap)
def write_paraview(RGB_vector,name):
colormap = open(str(name) + '.xml',"w")
colormap.write('<ColorMap name = "'+ str(name)+ '" space = "RGB">\n')
for i in range(len(RGB_vector)):
colormap.write('<Point x="'+str(i)+'" o="1" r="'+str(RGB_vector[0][i])+'" g="'+str(RGB_vector[1][i])+'" b="'+str(RGB_vector[2][i])+'"/>\n')
colormap.write('</ColorMap>')
file.close(colormap)
def write_paraview2(RGB_vector,name):
colormap = open(str(name) + '.xml',"w")
colormap.write('<ColorMap name = "'+ str(name)+ '" space = "RGB">\n')
for i in range(len(RGB_vector)/3):
colormap.write('<Point x="'+str(i)+'" o="1" r="'+str(RGB_vector[i*3])+'" g="'+str(RGB_vector[i*3+1])+'" b="'+str(RGB_vector[i*3+2])+'"/>\n')
colormap.write('</ColorMap>')
file.close(colormap)
def write_raw(RGB_vector,name):
colormap = open(str(name) + '.colormap',"w")
colormap.write('ColorMap name = ' + str(name)+'\n')
for i in range(len(RGB_vector)):
colormap.write(str(RGB_vector[0][i])+'\t'+str(RGB_vector[1][i])+'\t'+str(RGB_vector[2][i])+'\n')
file.close(colormap)
def read_raw(filename):
print 'void'

View File

@ -0,0 +1,173 @@
#!/usr/bin/python
# -*- coding: iso-8859-1 -*-
import math
# from http://code.activestate.com/recipes/121574-matrix-vector-multiplication/
def matmult(m, v):
nrows = len(m)
w = [None] * nrows
for row in range(nrows):
w[row] = reduce(lambda x,y: x+y, map(lambda x,y: x*y, m[row], v))
return w
# convert H(ue) S(aturation) L(uminance) to R(ot) G(elb) B(lau)
# with S,L,R,G,B running from 0 to 1, H running from 0 to 360
# from http://en.wikipedia.org/wiki/HSL_and_HSV
def HSL2RGB(HSL):
RGB = [0.0,0.0,0.0]
H_strich = HSL[0]/60.0
c = (1.0- abs(2.0 * HSL[2] - 1.0))*HSL[1]
x = c*(1.0- abs(H_strich%2-1.0))
m = HSL[2] -.5*c
if (0.0 <= H_strich)and(H_strich<1.0):
RGB[0] = c + m
RGB[1] = x + m
RGB[2] = 0.0 + m
elif (1.0 <= H_strich)and(H_strich<2.0):
RGB[0] = x + m
RGB[1] = c + m
RGB[2] = 0.0 + m
elif (2.0 <= H_strich)and(H_strich<3.0):
RGB[0] = 0.0 + m
RGB[1] = c + m
RGB[2] = x + m
elif (3.0 <= H_strich)and(H_strich<4.0):
RGB[0] = 0.0 + m
RGB[1] = x + m
RGB[2] = c + m
elif (4.0 <= H_strich)and(H_strich<5.0):
RGB[0] = x + m
RGB[1] = 0.0 + m
RGB[2] = c + m
elif (5.0 <= H_strich)and(H_strich<=6.0):
RGB[0] = c + m
RGB[1] = 0.0 + m
RGB[2] = x + m
for i in range(3):
RGB[i] = min(RGB[i],1.0)
RGB[i] = max(RGB[i],0.0)
return RGB
# convert R(ot) G(elb) B(lau) to H(ue) S(aturation) L(uminance)
# with S,L,R,G,B running from 0 to 1, H running from 0 to 360
# from http://130.113.54.154/~monger/hsl-rgb.html
def RGB2HSL(RGB):
HSL = [0.0,0.0,0.0]
maxcolor = max(RGB)
mincolor = min(RGB)
HSL[2] = (maxcolor + mincolor)/2.0
if(mincolor == maxcolor):
HSL[0] = 0.0
HSL[1] = 0.0
else:
if (HSL[2]<0.5):
HSL[1] = (maxcolor - mincolor)/(maxcolor + mincolor)
else:
HSL[1] = (maxcolor - mincolor)/(2.0 -maxcolor -mincolor)
if (maxcolor == RGB[0]):
HSL[0] = 0.0 + (RGB[1] - RGB[2])/(maxcolor - mincolor)
elif (maxcolor == RGB[1]):
HSL[0] = 2.0 + (RGB[2] - RGB[0])/(maxcolor - mincolor)
elif (maxcolor == RGB[2]):
HSL[0] = 4.0 + (RGB[0] - RGB[1])/(maxcolor - mincolor)
HSL[0] = HSL[0]*60.0
if (HSL[0] < 0.0):
HSL[0] = HSL[0] + 360.0
for i in range(2):
HSL[i+1] = min(HSL[i+1],1.0)
HSL[i+1] = max(HSL[i+1],0.0)
return HSL
# convert R(ot) G(elb) B(lau) to CIE XYZ
# with all values in the range of 0 to 1
# from http://www.cs.rit.edu/~ncs/color/t_convert.html
def RGB2XYZ(RGB):
XYZ = [0.0,0.0,0.0]
RGB_lin = [0.0,0.0,0.0]
for i in range(3):
if (RGB[i] > 0.04045):
RGB_lin[i] = ((RGB[i]+0.0555)/1.0555)**2.4
else:
RGB_lin[i] = RGB[i]/12.92
convert =[[0.412453,0.357580,0.180423],[0.212671,0.715160,0.072169],[0.019334,0.119193,0.950227]]
XYZ = matmult(convert,RGB_lin)
for i in range(3):
XYZ[i] = min(XYZ[i],1.0)
XYZ[i] = max(XYZ[i],0.0)
return XYZ
# convert CIE XYZ R(ot) G(elb) B(lau)
# with all values in the range of 0 to 1
# from http://www.cs.rit.edu/~ncs/color/t_convert.html
def XYZ2RGB(XYZ):
RGB_lin = [0.0,0.0,0.0]
RGB = [0.0,0.0,0.0]
convert =[[3.240479,-1.537150,-0.498535],[-0.969256,1.875992,0.041556],[0.055648,-0.204043,1.057311]]
RGB_lin = matmult(convert,XYZ)
for i in range(3):
if (RGB_lin[i] > 0.0031308):
RGB[i] = ((RGB_lin[i])**(1.0/2.4))*1.0555-0.0555
else:
RGB[i] = RGB_lin[i]*12.92
for i in range(3):
RGB[i] = min(RGB[i],1.0)
RGB[i] = max(RGB[i],0.0)
return RGB
# convert CIE Lab to CIE XYZ
# with XYZ in the range of 0 to 1
# from http://en.wikipedia.org/wiki/Lab_color_space, http://www.cs.rit.edu/~ncs/color/t_convert.html
def CIELab2XYZ(Lab,white):
XYZ = [0.0,0.0,0.0]
temp = [0.0,0.0,0.0]
temp[0] = 1.0/116.0 *(Lab[0] + 16.0) + 1/500.0 * Lab[1]
temp[1] = 1.0/116.0 *(Lab[0] + 16.0)
temp[2] = 1.0/116.0 *(Lab[0] + 16.0) - 1/200.0 * Lab[2]
for i in range(3):
if (temp[i] > 6.0/29.0):
temp[i] = temp[i]**(3.0)
else:
temp[i] =3 * (6.0/29.0)**2 * (temp[i]- 4.0/29.0)
XYZ[i] = white[i] * temp[i]
for i in range(3):
XYZ[i] = min(XYZ[i],1.0)
XYZ[i] = max(XYZ[i],0.0)
return XYZ
# convert CIE XYZ to CIE Lab
# with XYZ in the range of 0 to 1
# from http://en.wikipedia.org/wiki/Lab_color_space, http://www.cs.rit.edu/~ncs/color/t_convert.html
def XYZ2CIELab(XYZ,white):
temp = [0.0,0.0,0.0]
Lab = [0.0,0.0,0.0]
for i in range(3):
temp[i] = XYZ[i]/white[i]
if (temp[i] > (6.0/29.0)**3.0):
temp[i] = temp[i]**(1.0/3.0)
else:
temp[i] = 1.0/3.0 * (29.0/6.0)**2.0 * temp[i] + 4.0/29.0
Lab[0] = 116.0 * temp[1] - 16.0
Lab[1] = 500.0 *(temp[0] - temp[1])
Lab[2] = 200.0 *(temp[1] - temp[2])
return Lab
# convert Cie Lab to msh colorspace
# from http://www.cs.unm.edu/~kmorel/documents/ColorMaps/DivergingColorMapWorkshop.xls
def CIELab2Msh(Lab):
Msh = [0.0,0.0,0.0]
Msh[0] = math.sqrt(Lab[0]**2.0 + Lab[1]**2.0 + Lab[2]**2.0)
if (Msh[0] != 0.0):
Msh[1] = math.acos(Lab[0]/Msh[0])
if (Lab[1] != 0.0):
Msh[2] = math.atan2(Lab[2],Lab[1])
return Msh
# convert msh colorspace to Cie Lab
# from http://www.cs.unm.edu/~kmorel/documents/ColorMaps/DivergingColorMapWorkshop.xls
def Msh2CIELab(Msh):
Lab = [0.0,0.0,0.0]
Lab[0] = Msh[0] * math.cos(Msh[1])
Lab[1] = Msh[0] * math.sin(Msh[1]) * math.cos(Msh[2])
Lab[2] = Msh[0] * math.sin(Msh[1]) * math.sin(Msh[2])
return Lab

View File

@ -0,0 +1,66 @@
#!/usr/bin/python
# -*- coding: iso-8859-1 -*-
import math, convert_colormodels
def rad_dif(Msh1,Msh2,white):
HSL1 = convert_colormodels.RGB2HSL(convert_colormodels.XYZ2RGB(convert_colormodels.CIELab2XYZ(convert_colormodels.Msh2CIELab(Msh1),white)))
HSL2 = convert_colormodels.RGB2HSL(convert_colormodels.XYZ2RGB(convert_colormodels.CIELab2XYZ(convert_colormodels.Msh2CIELab(Msh2),white)))
return abs(HSL1[0]*math.pi/180.0-HSL2[0]*math.pi/180.0)
def adjust_hue(Msh_sat,M_unsat):
if Msh_sat[0] >= M_unsat:
return Msh_sat[2]
else:
hSpin = Msh_sat[1]*math.sqrt((M_unsat)**2.0-(Msh_sat[0])**2)/(Msh_sat[0]*math.sin(Msh_sat[1]))
if Msh_sat[2] > - math.pi/3.0:
return Msh_sat[2] + hSpin
else:
return Msh_sat[2] - hSpin
def interpolate_color(RGB1,RGB2,white,interp):
Msh1 = convert_colormodels.CIELab2Msh(convert_colormodels.XYZ2CIELab(convert_colormodels.RGB2XYZ(RGB1),white))
Msh2 = convert_colormodels.CIELab2Msh(convert_colormodels.XYZ2CIELab(convert_colormodels.RGB2XYZ(RGB2),white))
Msh_mid = [0.0,0.0,0.0]
if ((Msh1[1] > 0.05 and Msh2[1] > 0.05) and rad_dif(Msh1,Msh2,white) > math.pi/3.0):
Msh_mid[0] = max(Msh1[0],Msh2[0],88.0)
if interp < 0.5:
Msh2[0] = Msh_mid[0]
Msh2[1] = 0.0
Msh2[2] = 0.0
interp = 2.0*interp
else:
Msh1[0] = Msh_mid[0]
Msh1[1] = 0.0
Msh1[2] = 0.0
interp = 2.0*interp - 1.0
if (Msh1[1] < 0.05) and (Msh2[1] > 0.05):
Msh1[2] = adjust_hue(Msh2,Msh1[0])
elif (Msh2[1] < 0.05) and (Msh1[1] > 0.05):
Msh2[2] = adjust_hue(Msh1,Msh2[0])
for i in range(3):
Msh_mid[i] = (1.0-interp)*Msh1[i] + interp* Msh2[i]
return convert_colormodels.XYZ2RGB(convert_colormodels.CIELab2XYZ(convert_colormodels.Msh2CIELab(Msh_mid),white))
test1 = [0.231372549,0.298039216,0.752941176]
test2 = [0.705882353,0.015686275,0.149019608]
#test1 = [0.0,1.0,24.0/255.0]
#test1 = [0.5,0.5,0.5]
#test2 = [0.0,151.0/255.0,21.0/255.0]
x=0.950456
y=1.0
z=1.088754
white = [x, y, z]
test = test1
iteration = 33
delta = 1.0/float(iteration-1)
f = -delta
for i in range(iteration):
f = f + delta
test = test + interpolate_color(test1,test2,white,f)
test = test + test2
for i in range(iteration):
print i
print test[3*i+3]*255.0
print test[3*i+1+3]*255.0
print test[3*i+2+3]*255.0, '\n'

View File

@ -1,65 +0,0 @@
#!/usr/bin/python
# -*- coding: iso-8859-1 -*-
# This script is used to generate colormaps for gmsh (http://geuz.org/gmsh/)
# The script writes 360 files. Each file contains one colormap.
# More information on the used colors space can be found at http://en.wikipedia.org/wiki/HSL_and_HSV
# written by M. Diehl, m.diehl@mpie.de
import math
print '******************************************************************************'
print ' Write colormaps for gmsh'
print ''
print 'Suitable for datasets that have only positive/negative values'
print 'The colors are described using the HSL model.'
print 'Each of the 360 generated colormaps uses one value of (H)ue.'
print 'The colormaps runs at constant H from given (L)ightness and (S)aturation'
print 'to given L and S'
print 'L is distribute linearly, S changes as the square root of a linear list.'
print 'Suitable values: L_start = S_start = 1 , L_end = S_end = 0.'
print '******************************************************************************'
print ''
startL = float(raw_input('Please enter start value for (L)ightness: '))
endL = float(raw_input('Please enter end value for L: '))
startS = float(raw_input('Please enter start value for (S)aturation: '))
endS = float(raw_input('Please enter end value for S: '))
steps = int(raw_input('Please enter steps/resolution: '))
for h in range(0,360):
colormap = open('colormap_' + str(h).zfill(3) + '.map',"w")
colormap.write('View.ColorTable = {\n')
for i in range(0,steps):
h_strich = h/60.0
if(h_strich>6.0):
h_strich = h_strich-6.0
c = (1- abs(2*(startL + i*(endL-startL)/steps)-1))*math.sqrt(startS + i*(endS-startS)/steps)
x = c*(1- abs(h_strich%2-1))
m = (startL + i*(endL-startL)/steps) -.5*c
if (0.0 <= h_strich)and(h_strich<1.0):
colormap.write('{'+str((c+m)*255.0)+','+str((x+m)*255.0)+','+str((0.0+m)*255.0)+'},\n')
elif (1.0 <= h_strich)and(h_strich<2.0):
colormap.write('{'+str((x+m)*255.0)+','+str((c+m)*255.0)+','+str((0.0+m)*255.0)+'},\n')
elif (2.0 <= h_strich)and(h_strich<3.0):
colormap.write('{'+str((0.0+m)*255.0)+','+str((c+m)*255.0)+','+str((x+m)*255.0)+'},\n')
elif (3.0 <= h_strich)and(h_strich<4.0):
colormap.write('{'+str((0.0+m)*255.0)+','+str((x+m)*255.0)+','+str((c+m)*255.0)+'},\n')
elif (4.0 <= h_strich)and(h_strich<5.0):
colormap.write('{'+str((x+m)*255.0)+','+str((0.0+m)*255.0)+','+str((c+m)*255.0)+'},\n')
elif (5.0 <= h_strich)and(h_strich<=6.0):
colormap.write('{'+str((c+m)*255.0)+','+str((0.0+m)*255.0)+','+str((x+m)*255.0)+'},\n')
c = (1- abs(2*(startL)-1))*(startS)
x = c*(1- abs(h_strich%2-1))
m = (startL) -.5*c
if (0.0 <= h_strich)and(h_strich<1.0):
colormap.write('{'+str((c+m)*255.0)+','+str((x+m)*255.0)+','+str((0.0+m)*255.0)+'}};')
elif (1.0 <= h_strich)and(h_strich<2.0):
colormap.write('{'+str((x+m)*255.0)+','+str((c+m)*255.0)+','+str((0.0+m)*255.0)+'}};')
elif (2.0 <= h_strich)and(h_strich<3.0):
colormap.write('{'+str((0.0+m)*255.0)+','+str((c+m)*255.0)+','+str((x+m)*255.0)+'}};')
elif (3.0 <= h_strich)and(h_strich<4.0):
colormap.write('{'+str((0.0+m)*255.0)+','+str((x+m)*255.0)+','+str((c+m)*255.0)+'}};')
elif (4.0 <= h_strich)and(h_strich<5.0):
colormap.write('{'+str((x+m)*255.0)+','+str((0.0+m)*255.0)+','+str((c+m)*255.0)+'}};')
elif (5.0 <= h_strich)and(h_strich<=6.0):
colormap.write('{'+str((c+m)*255.0)+','+str((0.0+m)*255.0)+','+str((x+m)*255.0)+'}};')

View File

@ -1,87 +0,0 @@
#!/usr/bin/python
# -*- coding: iso-8859-1 -*-
# This script is used to generate colormaps for gmsh (http://geuz.org/gmsh/)
# The script writes 360 files. Each file contains one colormap.
# More information on the used colors space can be found at http://en.wikipedia.org/wiki/HSL_and_HSV
# written by M. Diehl, m.diehl@mpie.de
import math
print '******************************************************************************'
print ' Write colormaps for gmsh'
print ''
print 'Suitable for datasets running from negative to positive values.'
print 'The colors are described using the HSL model.'
print 'Each of the 360 generated colormaps uses two values of (H)ue.'
print 'The colormaps runs at constant H from given (L)ightness and (S)aturation'
print 'to given L_min and S_min and goes to H+180° (with given L and S)'
print 'Suitable values: L = L_min =.5, S = 1, and S_min=0,'
print '******************************************************************************'
print ''
startL = float(raw_input('Please enter start value for (L)ightness: '))
endL = float(raw_input('Please enter minimum value for L: '))
startS = float(raw_input('Please enter start value for (S)aturation: '))
endS = float(raw_input('Please enter minimum value for S: '))
steps = int(raw_input('Please enter steps/resolution: '))
steps = steps/2
for h in range(0,360):
colormap = open('colormap_' + str(h).zfill(3) + '.map',"w")
colormap.write('View.ColorTable = {\n')
i=0
h_strich = h/60.0
if(h_strich>6.0):
h_strich = h_strich-6.0
for j in range(0,steps+1):
# let L run linearly from 0 to 1, let S be the square root of a linear list from 0 to 1
c = (1- abs(2*(startL - j*(startL-endL)/steps)-1))*(startS - j*(startS-endS)/steps)
x = c*(1- abs(h_strich%2-1))
m = (startL - j*(startL-endL)/steps) -.5*c
if (0.0 <= h_strich)and(h_strich<1.0):
colormap.write('{'+str((c+m)*255.0)+','+str((x+m)*255.0)+','+str((0.0+m)*255.0)+'},\n')
elif (1.0 <= h_strich)and(h_strich<2.0):
colormap.write('{'+str((x+m)*255.0)+','+str((c+m)*255.0)+','+str((0.0+m)*255.0)+'},\n')
elif (2.0 <= h_strich)and(h_strich<3.0):
colormap.write('{'+str((0.0+m)*255.0)+','+str((c+m)*255.0)+','+str((x+m)*255.0)+'},\n')
elif (3.0 <= h_strich)and(h_strich<4.0):
colormap.write('{'+str((0.0+m)*255.0)+','+str((x+m)*255.0)+','+str((c+m)*255.0)+'},\n')
elif (4.0 <= h_strich)and(h_strich<5.0):
colormap.write('{'+str((x+m)*255.0)+','+str((0.0+m)*255.0)+','+str((c+m)*255.0)+'},\n')
elif (5.0 <= h_strich)and(h_strich<=6.0):
colormap.write('{'+str((c+m)*255.0)+','+str((0.0+m)*255.0)+','+str((x+m)*255.0)+'},\n')
i = i+1
h_strich = (h+180.0)/60.0
if(h_strich>6.0):
h_strich = h_strich-6.0
for j in range(1,steps):
c = (1- abs(2*(endL+j*(startL-endL)/steps)-1))*(endS+j*(startS-endS)/steps)
x = c*(1- abs(h_strich%2-1))
m = (endL+j*(startL-endL)/steps) -.5*c
if (0.0 <= h_strich)and(h_strich<1.0):
colormap.write('{'+str((c+m)*255.0)+','+str((x+m)*255.0)+','+str((0.0+m)*255.0)+'},\n')
elif (1.0 <= h_strich)and(h_strich<2.0):
colormap.write('{'+str((x+m)*255.0)+','+str((c+m)*255.0)+','+str((0.0+m)*255.0)+'},\n')
elif (2.0 <= h_strich)and(h_strich<3.0):
colormap.write('{'+str((0.0+m)*255.0)+','+str((c+m)*255.0)+','+str((x+m)*255.0)+'},\n')
elif (3.0 <= h_strich)and(h_strich<4.0):
colormap.write('{'+str((0.0+m)*255.0)+','+str((x+m)*255.0)+','+str((c+m)*255.0)+'},\n')
elif (4.0 <= h_strich)and(h_strich<5.0):
colormap.write('{'+str((x+m)*255.0)+','+str((0.0+m)*255.0)+','+str((c+m)*255.0)+'},\n')
elif (5.0 <= h_strich)and(h_strich<=6.0):
colormap.write('{'+str((c+m)*255.0)+','+str((0.0+m)*255.0)+','+str((x+m)*255.0)+'},\n')
i=i+1
c = (1- abs(2*(startL)-1))*(startS)
x = c*(1- abs(h_strich%2-1))
m = (startL) -.5*c
if (0.0 <= h_strich)and(h_strich<1.0):
colormap.write('{'+str((c+m)*255.0)+','+str((x+m)*255.0)+','+str((0.0+m)*255.0)+'}};')
elif (1.0 <= h_strich)and(h_strich<2.0):
colormap.write('{'+str((x+m)*255.0)+','+str((c+m)*255.0)+','+str((0.0+m)*255.0)+'}};')
elif (2.0 <= h_strich)and(h_strich<3.0):
colormap.write('{'+str((0.0+m)*255.0)+','+str((c+m)*255.0)+','+str((x+m)*255.0)+'}};')
elif (3.0 <= h_strich)and(h_strich<4.0):
colormap.write('{'+str((0.0+m)*255.0)+','+str((x+m)*255.0)+','+str((c+m)*255.0)+'}};')
elif (4.0 <= h_strich)and(h_strich<5.0):
colormap.write('{'+str((x+m)*255.0)+','+str((0.0+m)*255.0)+','+str((c+m)*255.0)+'}};')
elif (5.0 <= h_strich)and(h_strich<=6.0):
colormap.write('{'+str((c+m)*255.0)+','+str((0.0+m)*255.0)+','+str((x+m)*255.0)+'}};')

View File

@ -1,50 +0,0 @@
#!/usr/bin/python
# -*- coding: iso-8859-1 -*-
# This script is used to generate colormaps for paraview (www.paraview.org)
# The script writes 360 files. Each file contains one colormap.
# More information on the used colors space can be found at http://en.wikipedia.org/wiki/HSL_and_HSV
# written by M. Diehl, m.diehl@mpie.de
import math
print '******************************************************************************'
print ' Write colormaps for paraview'
print ''
print 'Suitable for datasets that have only positive/negative values'
print 'The colors are described using the HSL model.'
print 'Each of the 360 generated colormaps uses one value of (H)ue.'
print 'The colormaps runs at constant H from given (L)ightness and (S)aturation'
print 'to given L and S'
print 'L is distribute linearly, S changes as the square root of a linear list.'
print 'Suitable values: L_start = S_start = 1 , L_end = S_end = 0.'
print '******************************************************************************'
print ''
startL = float(raw_input('Please enter start value for (L)ightness: '))
endL = float(raw_input('Please enter end value for L: '))
startS = float(raw_input('Please enter start value for (S)aturation: '))
endS = float(raw_input('Please enter end value for S: '))
steps = int(raw_input('Please enter steps/resolution: '))
for h in range(0,360):
colormap = open('colormap_' + str(h) + '.xml',"w")
colormap.write('<ColorMap name = "' + str(h) + '" space = "RGB">\n')
for i in range(0,steps+1):
h_strich = h/60.0
if(h_strich>6.0):
h_strich = h_strich-6.0
c = (1- abs(2*(startL + i*(endL-startL)/steps)-1))*math.sqrt(startS + i*(endS-startS)/steps)
x = c*(1- abs(h_strich%2-1))
m = (startL + i*(endL-startL)/steps) -.5*c
if (0.0 <= h_strich)and(h_strich<1.0):
colormap.write('<Point x="'+str(i)+'" o="1" r="'+str(c+m)+'" g="'+str(x+m)+'" b="'+str(0.0+m)+'"/>\n')
elif (1.0 <= h_strich)and(h_strich<2.0):
colormap.write('<Point x="'+str(i)+'" o="1" r="'+str(x+m)+'" g="'+str(c+m)+'" b="'+str(0.0+m)+'"/>\n')
elif (2.0 <= h_strich)and(h_strich<3.0):
colormap.write('<Point x="'+str(i)+'" o="1" r="'+str(0.0+m)+'" g="'+str(c+m)+'" b="'+str(x+m)+'"/>\n')
elif (3.0 <= h_strich)and(h_strich<4.0):
colormap.write('<Point x="'+str(i)+'" o="1" r="'+str(0.0+m)+'" g="'+str(x+m)+'" b="'+str(c+m)+'"/>\n')
elif (4.0 <= h_strich)and(h_strich<5.0):
colormap.write('<Point x="'+str(i)+'" o="1" r="'+str(x+m)+'" g="'+str(0.0+m)+'" b="'+str(c+m)+'"/>\n')
elif (5.0 <= h_strich)and(h_strich<=6.0):
colormap.write('<Point x="'+str(i)+'" o="1" r="'+str(c+m)+'" g="'+str(0.0+m)+'" b="'+str(x+m)+'"/>\n')
colormap.write('</ColorMap>')

View File

@ -1,73 +0,0 @@
#!/usr/bin/python
# -*- coding: iso-8859-1 -*-
# This script is used to generate colormaps for paraview (www.paraview.org)
# The script writes 360 files. Each file contains one colormap.
# More information on the used colors space can be found at http://en.wikipedia.org/wiki/HSL_and_HSV
# written by M. Diehl, m.diehl@mpie.de
import math
print '******************************************************************************'
print ' Write colormaps for paraview'
print ''
print 'Suitable for datasets running from negative to positive values.'
print 'The colors are described using the HSL model.'
print 'Each of the 360 generated colormaps uses two values of (H)ue.'
print 'The colormaps runs at constant H from given (L)ightness and (S)aturation'
print 'to given L_min and S_min and goes to H+180° (with given L and S)'
print 'Suitable values: L = L_min =.5, S = 1, and S_min=0,'
print '******************************************************************************'
print ''
startL = float(raw_input('Please enter start value for (L)ightness: '))
endL = float(raw_input('Please enter minimum value for L: '))
startS = float(raw_input('Please enter start value for (S)aturation: '))
endS = float(raw_input('Please enter minimum value for S: '))
steps = int(raw_input('Please enter steps/resolution: '))
steps = steps/2
for h in range(0,360):
colormap = open('colormap_' + str(h).zfill(3) + '.xml',"w")
colormap.write('<ColorMap name = "Complementary Color ' + str(h).zfill(3) + '" space = "RGB">\n')
i=0
h_strich = h/60.0
if(h_strich>6.0):
h_strich = h_strich-6.0
for j in range(0,steps+1):
# let L run linearly from 0 to 1, let S be the square root of a linear list from 0 to 1
c = (1- abs(2*(startL - j*(startL-endL)/steps)-1))*(startS - j*(startS-endS)/steps)
x = c*(1- abs(h_strich%2-1))
m = (startL - j*(startL-endL)/steps) -.5*c
if (0.0 <= h_strich)and(h_strich<1.0):
colormap.write('<Point x="'+str(i)+'" o="1" r="'+str(c+m)+'" g="'+str(x+m)+'" b="'+str(0.0+m)+'"/>\n')
elif (1.0 <= h_strich)and(h_strich<2.0):
colormap.write('<Point x="'+str(i)+'" o="1" r="'+str(x+m)+'" g="'+str(c+m)+'" b="'+str(0.0+m)+'"/>\n')
elif (2.0 <= h_strich)and(h_strich<3.0):
colormap.write('<Point x="'+str(i)+'" o="1" r="'+str(0.0+m)+'" g="'+str(c+m)+'" b="'+str(x+m)+'"/>\n')
elif (3.0 <= h_strich)and(h_strich<4.0):
colormap.write('<Point x="'+str(i)+'" o="1" r="'+str(0.0+m)+'" g="'+str(x+m)+'" b="'+str(c+m)+'"/>\n')
elif (4.0 <= h_strich)and(h_strich<5.0):
colormap.write('<Point x="'+str(i)+'" o="1" r="'+str(x+m)+'" g="'+str(0.0+m)+'" b="'+str(c+m)+'"/>\n')
elif (5.0 <= h_strich)and(h_strich<=6.0):
colormap.write('<Point x="'+str(i)+'" o="1" r="'+str(c+m)+'" g="'+str(0.0+m)+'" b="'+str(x+m)+'"/>\n')
i = i+1
h_strich = (h+180.0)/60.0
if(h_strich>6.0):
h_strich = h_strich-6.0
for j in range(1,steps+1):
c = (1- abs(2*(endL+j*(startL-endL)/steps)-1))*(endS+j*(startS-endS)/steps)
x = c*(1- abs(h_strich%2-1))
m = (endL+j*(startL-endL)/steps) -.5*c
if (0.0 <= h_strich)and(h_strich<1.0):
colormap.write('<Point x="'+str(i)+'" o="1" r="'+str(c+m)+'" g="'+str(x+m)+'" b="'+str(0.0+m)+'"/>\n')
elif (1.0 <= h_strich)and(h_strich<2.0):
colormap.write('<Point x="'+str(i)+'" o="1" r="'+str(x+m)+'" g="'+str(c+m)+'" b="'+str(0.0+m)+'"/>\n')
elif (2.0 <= h_strich)and(h_strich<3.0):
colormap.write('<Point x="'+str(i)+'" o="1" r="'+str(0.0+m)+'" g="'+str(c+m)+'" b="'+str(x+m)+'"/>\n')
elif (3.0 <= h_strich)and(h_strich<4.0):
colormap.write('<Point x="'+str(i)+'" o="1" r="'+str(0.0+m)+'" g="'+str(x+m)+'" b="'+str(c+m)+'"/>\n')
elif (4.0 <= h_strich)and(h_strich<5.0):
colormap.write('<Point x="'+str(i)+'" o="1" r="'+str(x+m)+'" g="'+str(0.0+m)+'" b="'+str(c+m)+'"/>\n')
elif (5.0 <= h_strich)and(h_strich<=6.0):
colormap.write('<Point x="'+str(i)+'" o="1" r="'+str(c+m)+'" g="'+str(0.0+m)+'" b="'+str(x+m)+'"/>\n')
i=i+1
colormap.write('</ColorMap>')

View File

@ -869,6 +869,66 @@ subroutine calculate_mises(res_x,res_y,res_z,tensor,vm)
enddo; enddo; enddo
end subroutine calculate_mises
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine curl_fft(res_x,res_y,res_z,vec_tens,geomdim,field,divergence_field)
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! calculates divergence field using integration in Fourier space
!use vec_tens to decide if tensor (3) or vector (1)
implicit none
integer res_x, res_y, res_z, vec_tens
real*8 geomdim(3)
real*8 field(res_x,res_y,res_z,vec_tens,3)
real*8 curl_fft(res_x,res_y,res_z,vec_tens,3)
complex*8 field_fft(res_x/2_pInt+1,res_y,res_z,vec_tens,3)
real*8 xi(res_x,res_y,res_z,3)
complex*16 img
integer i, j, k
real*8, parameter :: pi = 3.14159265358979323846264338327950288419716939937510
integer*8 :: plan_fft(2)
img = cmplx(0.0,1.0)
call dfftw_plan_many_dft_r2c(plan_fft(1),3,(/res_x,res_y,res_z/),vec_tens*3,&
curl_fft,(/res_x,res_y,res_z/),1,res_x*res_y*res_z,&
field_fft,(/res_x/2+1,res_y,res_z/),1,(res_x/2+1)*res_y*res_z,32) ! 32 =FFTW_PATIENT
call dfftw_plan_many_dft_c2r(plan_fft(2),3,(/res_x,res_y,res_z/),vec_tens*3,&
field_fft,(/res_x/2+1,res_y,res_z/),1,(res_x/2+1)*res_y*res_z,&
curl_fft,(/res_x,res_y,res_z/),1,res_x*res_y*res_z,32) ! 32 = FFTW_PATIENT
! field_copy is destroyed during plan creation
curl_fft = field
call dfftw_execute_dft_r2c(plan_fft(1), field_copy, field_fft)
xi = 0.0
do k = 0, res_z-1
do j = 0, res_y-1
do i = 0, res_x/2
xi(i+1,j+1,k+1,:) = (/real(i),real(j),real(k)/)/geomdim
if(k==res_z/2) xi(i+1,j+1,k+1,3)= 0.0 ! set highest frequencies to zero
if(j==res_y/2) xi(i+1,j+1,k+1,2)= 0.0
if(i==res_x/2) xi(i+1,j+1,k+1,1)= 0.0
enddo; enddo; enddo
do k = 1, res_z
do j = 1, res_y
do i = 1, res_x/2+1
divergence_field_fft(i,j,k,1) = sum(field_fft(i,j,k,1,:)*xi(i,j,k,:))
if(vec_tens == 3) then
divergence_field_fft(i,j,k,2) = sum(field_fft(i,j,k,2,:)*xi(i,j,k,:))
divergence_field_fft(i,j,k,3) = sum(field_fft(i,j,k,3,:)*xi(i,j,k,:))
endif
enddo; enddo; enddo
divergence_field_fft = divergence_field_fft*img*2.0*pi
call dfftw_execute_dft_c2r(plan_fft(2), divergence_field_fft, divergence_field)
end subroutine curl_fft
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
subroutine divergence_fft(res_x,res_y,res_z,vec_tens,geomdim,field,divergence_field)
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++