New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

possible isofill issue #519

Closed
doutriaux1 opened this Issue Jul 25, 2014 · 19 comments

Comments

Projects
None yet
3 participants
@doutriaux1
Member

doutriaux1 commented Jul 25, 2014

Celine gave me some data that plot fine in boxfill but seem to be way off in isofill

@doutriaux1 doutriaux1 added this to the 2.0.0 milestone Jul 25, 2014

@doutriaux1 doutriaux1 added the VCS label Jul 25, 2014

@aashish24

This comment has been minimized.

Contributor

aashish24 commented Jul 30, 2014

@doutriaux1 Can you post the data somewhere?

@doutriaux1

This comment has been minimized.

Member

doutriaux1 commented Jul 30, 2014

will do.

@doutriaux1

This comment has been minimized.

Member

doutriaux1 commented Jul 30, 2014

@aashish24 I think this is a VTK/Kitware one, I'm doing something wrong but I don't know what.

branch: vcs_celine_isof_bug
test: Test/test_celine_iso_bug.py
On the attached picture, notice how the data around florida is "White" the actual values near there are 1.5 or so, so it should be red-ish
I left the comments on in this branch that shows the lookuptable colors are set properly as well as the levels. Not sure why it is white. "sFilter" issue? Bandedcontour issue?
bad

@doutriaux1

This comment has been minimized.

Member

doutriaux1 commented Jul 30, 2014

Also on this picture you can see why jerry thinks it looks "horrible" everythig is indeed "blurry" (look at the red colors). bug #512

@aashish24

This comment has been minimized.

Contributor

aashish24 commented Jul 30, 2014

@doutriaux1 are you converting the data to the image somehow? The polydata cannot be blurry. That would be really odd. I guess we have to look at your code together?

@doutriaux1

This comment has been minimized.

Member

doutriaux1 commented Jul 30, 2014

using a BUNCH of filter
sFilter = vtk.vtkDataSetSurfaceFilter()
sFilter.SetInputData(ug) # ug here is a structured grid
mapper = vtk.vtkPolyDataMapper()
cot = vtk.vtkBandedPolyDataContourFilter()
cot.ClippingOn()
cot.SetInputData(sFilter.GetOutput())
cot.SetNumberOfContours(len(l))
for j,v in enumerate(l):
cot.SetValue(j,l[j])
cot.Update()
lut = vtk.vtkLookupTable()
print "NCOLORS:",len(COLS[i])
lut.SetNumberOfTableValues(len(COLS[i]))
for j,color in enumerate(COLS[i]):
r,g,b = cmap.index[color]
print l[j],vcs.colors.rgb2str(r_2.55,g_2.55,b*2.55),l[j+1]
lut.SetTableValue(j,r/100.,g/100.,b/100.)
mapper.SetLookupTable(lut)
mapper.SetScalarRange(lmn,lmx)

@doutriaux1

This comment has been minimized.

Member

doutriaux1 commented Jul 30, 2014

will post here a pure vtk code that does what I do. it wil be easier to debug

@aashish24

This comment has been minimized.

Contributor

aashish24 commented Jul 31, 2014

@doutriaux1 thanks.

@doutriaux1

This comment has been minimized.

Member

doutriaux1 commented Jul 31, 2014

here you go:

from vtk.util import numpy_support as VN
f=cdms2.open("celine.nc")
data=f("data")
ug = vtk.vtkStructuredGrid()
lat = data.getLatitude()
lon = data.getLongitude()
lat = lat[:,numpy.newaxis]_numpy.ones(lon.shape)[numpy.newaxis,:]
lon = lon[numpy.newaxis,:]_numpy.ones(lat.shape)[:,numpy.newaxis]
ug.SetDimensions(lat.shape[1],lat.shape[0],1)
lon = numpy.ma.ravel(lon)
lat = numpy.ma.ravel(lat)
sh = list(lat.shape)
sh.append(1)
lon = numpy.ma.reshape(lon,sh)
lat = numpy.ma.reshape(lat,sh)
z = numpy.zeros(lon.shape)
m3 = numpy.concatenate((lon,lat),axis=1)
m3 = numpy.concatenate((m3,z),axis=1)
xm=lon.min()
xM=lon.max()
ym=lat.min()
yM=lat.max()

# First create the points/vertices (in vcs terms)

deep = True
pts = vtk.vtkPoints()

## Convert nupmy array to vtk ones

ppV = VN.numpy_to_vtk(m3,deep=deep)
pts.SetData(ppV)
ug.SetPoints(pts)
data = VN.numpy_to_vtk(data.filled(0.).flat,deep=True)
ug.GetPointData().SetScalars(data)
sFilter = vtk.vtkDataSetSurfaceFilter()
sFilter.SetInputData(ug)
sFilter.Update()
l = [-4.0584776401519775, -2, -1.5, -1.0, -0.4, -0.1, 0, 0.1, 0.4, 1.0, 1.5, 2, 7.875967502593994]
cols = [[0, 57, 100],[0, 100, 33],[27, 100, 0],[88, 100, 0], [100, 78, 0], [100, 100, 100], [100, 100, 100], [100, 32, 0], [100, 17, 0], [75, 0, 0], [45, 0, 0], [49, 0, 24],]
mapper = vtk.vtkPolyDataMapper()
cot = vtk.vtkBandedPolyDataContourFilter()
cot.ClippingOn()
cot.SetInputData(sFilter.GetOutput())
cot.SetNumberOfContours(len(l))
for j,v in enumerate(l):
  cot.SetValue(j,l[j])
cot.Update()
mapper.SetInputConnection(cot.GetOutputPort())
lut = vtk.vtkLookupTable()
lut.SetNumberOfTableValues(len(cols))
for j,color in enumerate(cols):
  r,g,b = cols[j]
  lut.SetTableValue(j,r/100.,g/100.,b/100.)
mapper.SetLookupTable(lut)
mapper.SetScalarRange(-4.05847764015,7.87596750259)
act = vtk.vtkActor()
act.SetMapper(mapper)
ren=vtk.vtkRenderer()
ren.AddActor(act)
renWin=vtk.vtkRenderWindow()
renWin.SetSize(800,600)
renWin.AddRenderer(ren)
renWin.Render()
raw_input("ok press enter now")

@aashish24

This comment has been minimized.

Contributor

aashish24 commented Jul 31, 2014

Ok. @jbeezley and I are looking to the code now. We will report back once we know what's going on.

@doutriaux1

This comment has been minimized.

Member

doutriaux1 commented Jul 31, 2014

Ok thank you, thank you . thank you! I suspect (but no idea) a celltodata or datatocell call needs (not) to be made somewhere. BTW Plotting in boxfill the data comes out correctly

@jbeezley

This comment has been minimized.

Contributor

jbeezley commented Jul 31, 2014

I'm far from a VTK expert, but I made a few changes that seem to help.

import vtk,cdms2,numpy

from vtk.util import numpy_support as VN

f=cdms2.open("celine.nc")
data=f("data")
ug = vtk.vtkStructuredGrid()
lat = data.getLatitude()
lon = data.getLongitude()
lat = lat[:,numpy.newaxis]*numpy.ones(lon.shape)[numpy.newaxis,:]
lon = lon[numpy.newaxis,:]*numpy.ones(lat.shape)[:,numpy.newaxis]
ug.SetDimensions(lat.shape[1],lat.shape[0],1)
lon = numpy.ma.ravel(lon)
lat = numpy.ma.ravel(lat)
sh = list(lat.shape)
sh.append(1)
lon = numpy.ma.reshape(lon,sh)
lat = numpy.ma.reshape(lat,sh)
z = numpy.zeros(lon.shape)
m3 = numpy.concatenate((lon,lat),axis=1)
m3 = numpy.concatenate((m3,z),axis=1)
xm=lon.min()
xM=lon.max()
ym=lat.min()
yM=lat.max()
# First create the points/vertices (in vcs terms)
deep = True
pts = vtk.vtkPoints()
## Convert nupmy array to vtk ones
ppV = VN.numpy_to_vtk(m3,deep=deep)
pts.SetData(ppV)
ug.SetPoints(pts)
data = VN.numpy_to_vtk(data.filled(0.).flat,deep=True)
ug.GetPointData().SetScalars(data)

sFilter = vtk.vtkDataSetSurfaceFilter()
sFilter.SetInputData(ug)
sFilter.Update()
l = [-4.0584776401519775, -2, -1.5, -1.0, -0.4, -0.1, 0, 0.1, 0.4, 1.0, 1.5, 2, 7.875967502593994]
cols = [[0, 57, 100],[0, 100, 33],[27, 100, 0],[88, 100, 0], [100, 78, 0], [100, 100, 100], [100, 100, 100], [100, 32, 0], [100, 17, 0], [75, 0, 0], [45, 0, 0], [49, 0, 24],]

mapper = vtk.vtkPolyDataMapper()
cot = vtk.vtkBandedPolyDataContourFilter()
cot.ClippingOn()
cot.SetInputData(sFilter.GetOutput())
cot.SetNumberOfContours(len(l))
for j,v in enumerate(l):
  cot.SetValue(j,l[j])
cot.SetScalarModeToIndex()
cot.Update()

mapper.SetInputConnection(cot.GetOutputPort())
lut = vtk.vtkLookupTable()
lut.SetNumberOfTableValues(len(cols))
for j,color in enumerate(cols):
  r,g,b = cols[j]
  lut.SetTableValue(j,r/100.,g/100.,b/100.)
mapper.SetLookupTable(lut)
mapper.SetScalarRange(0, len(l) - 1)
mapper.SetScalarModeToUseCellData()

act = vtk.vtkActor()
act.SetMapper(mapper)

ren=vtk.vtkRenderer()
ren.AddActor(act)
renWin=vtk.vtkRenderWindow()
renWin.SetSize(800,600)
renWin.AddRenderer(ren)
renWin.Render()

iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)
ren.SetBackground(1, 1, 1)
iren = vtk.vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)
iren.Start()
raw_input("ok press enter now")

The contours need to be adjusted, but here is the result.
image_new

@doutriaux1

This comment has been minimized.

Member

doutriaux1 commented Jul 31, 2014

yes! I confirmed it works. AS I suspected it had something to do with celldata vs pointdata.

@aashish24 since we set the data on point, why do we need to use:
mapper.SetScalarModeToUseCellData()

@doutriaux1

This comment has been minimized.

Member

doutriaux1 commented Jul 31, 2014

and that also seems to fix the bluriness that jerry was complaining about.
@aashish24 @jbeezley so should i use:
ug.GetPointData().SetScalars(data)
or would it be better to use:
ug.GetCellData().SetScalars(data)

@doutriaux1

This comment has been minimized.

Member

doutriaux1 commented Jul 31, 2014

oohh. i my "masking code" which I don't need here i see that I do:
p2c = vtk.vtkPointDataToCellData()
p2c.SetInputData(grid2)
geoFilter.SetInputConnection(p2c.GetOutputPort())

@aashish24

This comment has been minimized.

Contributor

aashish24 commented Jul 31, 2014

@doutriaux1 I think you want to use cell data since you want them filled. If you use point data then it will interpolate colors based on the scalar which I don't think you want. Also, you may also have to duplicate points if you use the cell data (may be). In VCS old code, did you have data on the cell or on the grid point? @jbeezley: cool!!

@aashish24

This comment has been minimized.

Contributor

aashish24 commented Aug 6, 2014

@doutriaux1 can we close this one now?

@doutriaux1

This comment has been minimized.

Member

doutriaux1 commented Aug 8, 2014

when it's in master yes branch is vcs_celine_bug_isof

@aashish24

This comment has been minimized.

Contributor

aashish24 commented Aug 8, 2014

In b040125 now.

@aashish24 aashish24 closed this Aug 8, 2014

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment