三维的热力图模型 导出3D热力图(

# -*- coding: utf-8 -*- """ Created on Thu Oct 2 01:49:21 2017 @author: Pasztor Balazs """ import numpy as np import pyopencl as cl import math import time import visvis as vv from operator import itemgetter from itertools import cycle from vispy import app as vispyapp from vispy import scene from vispy.color import get_colormaps, BaseColormap from vispy.visuals.transforms import STTransform class mytimer: def __init__(self): self.last_time = time.perf_counter() #self.last_time = time.clock() self.times = [['process', 'timed', 'started', 'ended']] def measure(self, str): cur_time = time.perf_counter() #cur_time = time.clock() elapsed_time = cur_time - self.last_time self.times += [[str, elapsed_time, self.last_time, cur_time]] self.last_time = cur_time def print_times(self): print('{0[0]:40}: {0[1]:7.7}, {0[2]:7.7}, {0[3]:7.7}'.format(self.times[0])) for t in self.times[1:]: print('{0:40}: {1:7.4f}, {2:7.4f}, {3:7.4f}'.format(t[0], t[1], t[2], t[3])) ctx = cl.create_some_context() queue = cl.CommandQueue(ctx) mf = cl.mem_flags def convert_np_data(rawdata): datadim = len(rawdata) datalength = len(rawdata[0]) data = np.zeros((datadim, datalength)) for d in range(datadim): data[d] = rawdata[d] return data def bincount_CPU(data, bincounts = [10, 10]): data = convert_np_data(data) datadim = len(data) datalength = len(data[0]) #data = np.zeros((datadim, datalength)) #for d in range(datadim): # data[d] = rawdata[d] #Min/Max of data per dimensions dataMins = np.zeros(datadim) dataMaxs = np.zeros(datadim) for d in range(datadim): dataMins[d] = np.amin(data[d]) dataMaxs[d] = np.amax(data[d]) # BINNING start bindim = len(bincounts) if(datadim < bindim): print('The specified bins have more dimensions than the data.', ' Extra dimenions will be ignored during binning') bincounts = bincounts[:datadim] bindim = datadim #for d in range (bindim): # #if less bincounts are specified than the dimension of the data # if len(binCounts) < (d + 1): # binCounts.append(1) #calculating the widths of the bins binwidths = np.zeros(datadim) for d in range (bindim): binwidths[d] = (dataMaxs[d]-dataMins[d])/(bincounts[d]-2) #creating the bin-summarize array, with increased size for guards, NaNs, missing bins = np.zeros(tuple(v+0 for v in bincounts)) for e in range(datalength): currbinnum = list() for d in range(bindim): currbinnum += [int((data[d, e]-dataMins[d])/binwidths[d]) + 1 ] # SUMMERIZE -- currently count bins[tuple(currbinnum)] += 1 return bins def show_visvis(bins_h, axislabels=None): app = vv.use() # Load volume vol = bins_h vv.figure(1); vv.clf() RS = ['ray'] tt = [] a = vv.subplot(1,1,1) t = vv.volshow(vol) t.interpolate = False #vv.title('Volume renderer') t.colormap = vv.CM_HOT t.renderStyle = 'ray' tt.append(t) if axislabels is not None: a.axis.xLabel = axislabels[0].replace('_', ' ') a.axis.yLabel = axislabels[1].replace('_', ' ') a.axis.zLabel = axislabels[2].replace('_', ' ') # Create colormap editor in first axes cme = vv.ColormapEditor(vv.gcf(), *tt[:]) # Run app app.Create() app.Run() # create colormaps that work well for translucent and additive volume rendering class TransFire(BaseColormap): glsl_map = """ vec4 translucent_fire(float t) { return vec4(pow(t, 0.5), t, t*t, max(0, t*1.05 - 0.05)); } """ class TransGrays(BaseColormap): glsl_map = """ vec4 translucent_grays(float t) { return vec4(t, t, t, t*0.05); } """ class TransRare(BaseColormap): glsl_map = """ vec4 translucent_grays(float t) { if (t >= 0.0001){ float c = 0.3 + 0.7 * t; return vec4(c, c, c, 0.1); } return vec4(0.0, 0.0, 0.0, 0.0); } """ class TransFireRare(BaseColormap): glsl_map = """ vec4 translucent_fire(float t) { if (t >= 0.0001){ return vec4(pow(t, 0.5), t, t*t, 0.1); } return vec4(0.0, 0.0, 0.0, 0.0); } """ class TransFireRareOld(BaseColormap): glsl_map = """ vec4 translucent_fire(float t) { if (t >= 0.0001){ return vec4(0.3 + 0.7* pow(t, 0.5), 0.3 + 0.7* t, 0.3 + 0.7* t*t, 0.1); } return vec4(0.0, 0.0, 0.0, 0.0); } """ # Setup colormap iterators opaque_cmaps = cycle(get_colormaps()) translucent_cmaps = cycle([TransFire(), TransGrays(), TransRare(), TransFireRare()]) opaque_cmap = next(opaque_cmaps) translucent_cmap = next(translucent_cmaps) def show_vispy(bins_h): vol = bins_h # Prepare canvas canvas = scene.SceneCanvas(keys='interactive', app='pyqt5', size=(800, 600), show=True) canvas.measure_fps() # Set up a viewbox to display the image with interactive pan/zoom view = canvas.central_widget.add_view() # Set whether we are emulating a 3D texture emulate_texture = False #OpenGL ES miatt kell a true csak # Create the volume visuals, only one is visible volume1 = scene.visuals.Volume(vol, parent=view.scene, threshold=0.225, emulate_texture=emulate_texture) # Create three cameras (Fly, Turntable and Arcball) fov = 60. cam1 = scene.cameras.FlyCamera(parent=view.scene, fov=fov, name='Fly') cam2 = scene.cameras.TurntableCamera(parent=view.scene, fov=fov, name='Turntable') cam3 = scene.cameras.ArcballCamera(parent=view.scene, fov=fov, name='Arcball') view.camera = cam1 # Select turntable at first # Create an XYZAxis visual axis = scene.visuals.XYZAxis(parent=view) s = STTransform(translate=(250, 250), scale=(550, 550, 550, 1)) affine = s.as_matrix() axis.transform = affine # Implement axis connection with cam2 @canvas.events.mouse_move.connect def on_mouse_move(event): if event.button == 1 and event.is_dragging: axis.transform.reset() axis.transform.rotate(cam2.roll, (0, 0, 1)) axis.transform.rotate(cam2.elevation, (1, 0, 0)) axis.transform.rotate(cam2.azimuth, (0, 1, 0)) axis.transform.scale((50, 50, 0.001)) axis.transform.translate((50., 50.)) axis.update() # Implement key presses @canvas.events.key_press.connect def on_key_press(event): global opaque_cmap, translucent_cmap if event.text == '1': cam_toggle = {cam1: cam2, cam2: cam3, cam3: cam1} view.camera = cam_toggle.get(view.camera, cam2) print(view.camera.name + ' camera') if view.camera is cam2: axis.visible = True else: axis.visible = False elif event.text == '2': print('2-pressed') methods = ['mip', 'translucent', 'iso', 'additive'] method = methods[(methods.index(volume1.method) + 1) % 4] print("Volume render method: %s" % method) cmap = opaque_cmap if method in ['mip', 'iso'] else translucent_cmap volume1.method = method volume1.cmap = cmap elif event.text == '3': if volume1.method in ['mip', 'iso']: cmap = opaque_cmap = next(opaque_cmaps) else: cmap = translucent_cmap = next(translucent_cmaps) volume1.cmap = cmap elif event.text == '0': cam1.set_range() cam3.set_range() elif event.text != '' and event.text in '-+': s = -0.025 if event.text == '+' else 0.025 volume1.threshold += s th = volume1.threshold print("Isosurface thresh