pro construct_orthogonal_projection, slice_data ; PLOT PROJECTIONS OF QUANTITIES THROUGH ENITRE SIMULATION VOLUME ; Weighting options are none (0), density(1) and density^2(2) so far. @common_blocks.inc ; using grid_info read in data (found in data_dir directory) of all ; necessary grids for projections ; ; slice_ori: is a 3 component vectors with one not equal zero along ; which direction the projection is done. ; ; slice_size has 4 components specifying the left and right edge of ; the projection ; note that slice_ori has to have two elements set to zero since we ; only do orthogonal slices. ; if total(grid_info[0].dim gt 1) ne 3 then begin print, 'orthogonal projections need 3D data' return endif weighted = weight_index slice_data = DBLARR(xy_sl_size,xy_sl_size) if (weighted gt 0) then weight_data = DBLARR(xy_sl_size,xy_sl_size) ; adjust slice_sice ever so slightly to make sure the left edge ; are an integer number of finest resolution pixels slice_size = DOUBLE(slice_size) depthp = DOUBLE(slice_value+depth_value/2.) depthp = nearest_valid_depth_value(depthp, grid_info[0].dim[0]) depthm = (depthp-depth_value) if depthp gt 1.D then begin depthm = 1.D - depth_value depthp = 1.D endif if depthm lt 0. then begin depthm = 0.D depthp = depth_value endif dxdy = DOUBLE(slice_size[2:3]-slice_size[0:1])/DOUBLE([xy_sl_size,xy_sl_size]) slice_size[0] = dxdy[0]*LONG64(slice_size[0]/dxdy[0]) slice_size[1] = dxdy[1]*LONG64(slice_size[1]/dxdy[1]) slice_size[2] = slice_size[0]+dxdy[0]*(xy_sl_size) slice_size[3] = slice_size[1]+dxdy[1]*(xy_sl_size) update_slice_text, slice_text, center, slice_size,image_stack(0).title dxts = DOUBLE(xy_sl_size)/(slice_size[2]-slice_size[0]) dyts = DOUBLE(xy_sl_size)/(slice_size[3]-slice_size[1]) num_dim = 3 ; rough check of input data slice_data[*,*] = 0.D axes = (slice_ori gt 1.e-30) const_sub = where(axes eq max(axes)) const_sub = const_sub[0] other_subs = where(axes eq 0) dom = DomainRightEdge-DomainLeftEdge depthp *= dom[const_sub] depthm *= dom[const_sub] IF (Fix(TOTAL(axes)) gt 1) then BEGIN print, 'construct_orthogonal_projection_data: slice_ori:', $ slice_ori,' not allowed !' STOP ENDIF IF (N_ELEMENTS(grid_info) lt 1) THEN BEGIN print, 'construct_orthogonal_projection: no grid information in grid_info:',$ grid_info STOP ENDIF dim_grid_info = size(grid_info) num_of_grids = dim_grid_info[1] ii_l = 0 ; min left edge and scaling min_left = DOUBLE([slice_size(0),slice_size(1)]) max_right = DOUBLE([slice_size(2),slice_size(3)]) scale_up = min(1.D/(max_right-min_left)) ; use only grids that will show up on the image: ; determine grids that are at least partially in the projection grid_flag = grid_in_projection(grid_info, DOUBLE(slice_ori), DOUBLE(slice_size), [depthm,depthp]) index_range = grid_info.End_index[0]-grid_info.Start_index[0]+1 delta_distance = (DOUBLE(grid_info.Right_edge[0] - grid_info.Left_edge[0])/DOUBLE(index_range)) del_pix = DOUBLE(xy_sl_size)*delta_distance/(max_right(0)-min_left(0)) u_g = grid_info(where((del_pix GE DOUBLE(fastmode)) AND (grid_flag gt 0))) if N_ELEMENTS(where(grid_info[0].dim[*] gt 0)) eq 2 then begin print, 'use slice for 2D data! ' u_g = grid_info endif u_g = u_g[SHIFT((REVERSE(SORT(u_g.level))),20)] dim_grid_info = size(u_g) num_of_grids = dim_grid_info[1] index_range = u_g.End_index-u_g.Start_index+1 delta_distance = (u_g.Right_edge - u_g.Left_edge)/DOUBLE(index_range) n_stuetz = 0L this_level=0 le_o = u_g(0:num_of_grids-1).left_edge(other_subs(0)) le_t = u_g(0:num_of_grids-1).left_edge(other_subs(1)) ri_o = u_g(0:num_of_grids-1).right_edge(other_subs(0)) ri_t = u_g(0:num_of_grids-1).right_edge(other_subs(1)) slice_ori_l = dblarr(3) slice_ori_l(const_sub) = slice_ori(const_sub) slice_ori_l(other_subs) = min_left slice_ori_r = slice_ori_l slice_ori_r(const_sub) = slice_ori(const_sub) slice_ori_r(other_subs) = max_right ii_l0 = (FLOOR((depthm - u_g.Left_edge[const_sub]) $ /delta_distance[const_sub,*]) ) > 0 < (index_range-1) ii_l1 = (FLOOR((slice_ori_l[other_subs[0]] - u_g.Left_edge[other_subs[0]]) $ /delta_distance[other_subs[0],*]) ) > 0 < (index_range-1) ii_l2 = (FLOOR((slice_ori_l[other_subs[1]] - u_g.Left_edge[other_subs[1]]) $ /delta_distance[other_subs[1],*]) ) > 0 < (index_range-1) ii_r0 = (((FLOOR((depthp - u_g.Left_edge[const_sub])/ $ delta_distance[const_sub,*]) ) < (index_range[const_sub,*]-1))) > 0 ii_r1 = ((FLOOR((slice_ori_r[other_subs[0]] - u_g.Left_edge(other_subs[0]))/ $ delta_distance[other_subs[0],*]) ) < (index_range[other_subs[0],*]-1)) > 0 ii_r2 = ((FLOOR((slice_ori_r[other_subs[1]] - u_g.Left_edge(other_subs[1]))/ $ delta_distance[other_subs[1],*]) ) < (index_range[other_subs[1],*]-1)) > 0 start = intarr(num_of_grids,3) & stride = intarr(num_of_grids,3) & righti = intarr(num_of_grids,3) start[*,const_sub] = ii_l0 start[*,other_subs[0]] = ii_l1 start[*,other_subs[1]] = ii_l2 righti[*,const_sub] = ii_r0 righti[*,other_subs[0]] = ii_r1 righti[*,other_subs[1]] = ii_r2 ; read all selected grid data fieldnames = list_str[var_index] if ((weighted gt 0) ) then begin weightfield= list_str[weight_index-1] fieldnames = [fieldnames, weightfield] weightindex = (where(fieldnames eq weightfield))[0] endif vin = 0 npoints = (TOTAL(LONG(righti[*,other_subs[0]] - start[*, other_subs[0]]+1)*(righti[*,other_subs[1]] - start[*, other_subs[1]]+1))) ; print, 'at most will have ', npoints, ' in projection' harr = DBLARR(npoints) x_stuetz = harr y_stuetz = harr z_stue = harr dx_stuetz = harr if (weighted ge 1) then w_stuetz = harr Num_pix_x = xy_sl_size Num_pix_y = xy_sl_size xi = other_subs[0] & yi = other_subs[1] & zi = const_sub if verbose then print, 'construcing image from ', num_of_grids,' grids in '+fieldnames[vin]+' projection' zero = 0. ; loop over all relevant grids ngrids = 401; ( read how many grids in one go?) FOR kk=0L, num_of_grids-1,ngrids DO BEGIN inn = min([kk+ngrids-1, num_of_grids-1]) a = get_data(u_g[kk:inn], fieldnames, const_sub, $ LEFT_INDEX=start[kk:inn,*], RIGHT_INDEX=righti[kk:inn,*]) zero_solution_under_subgrid, u_g[kk:inn], a, ZERO=zero, all_grid_info = u_g print, 'read up to:', inn for i=kk,inn do begin ll = i-kk cg = u_g[i] i_l = [ii_l1[i], ii_l2[i]] i_r = [ii_r1[i], ii_r2[i]] ; project if (weighted eq 0) then begin data2d = (*a[ll].data[vin]) if (zi+1 le N_elements(size(*a[ll].data[vin], /dimension))) then data2d = TOTAL((*a[ll].data[vin]), zi+1, /DOUBLE) endif ELSE BEGIN data2d = TOTAL(fn(a[ll], weightfield)*fn(a[ll],fieldnames[0]), zi+1, /DOUBLE) weight2d = TOTAL(fn(a[ll], weightfield), zi+1, /DOUBLE) ENDELSE ; fast way to construct x and y points dix = i_r[0]-i_l[0]+1 diy = i_r[1]-i_l[1]+1 locs = LINDGEN(dix*diy) ind = ARRAY_INDICES([dix, diy], locs, /DIMENSIONS) ind[0,*] += i_l[0] ind[1,*] += i_l[1] end_i = LONG(dix)*diy - 1 x_stuetz = cg.left_edge[xi] + DOUBLE(ind[0,*])*delta_distance[xi,i] y_stuetz = cg.left_edge[yi] + DOUBLE(ind[1,*])*delta_distance[yi,i] dx_stuetz[0:end_i] = delta_distance[xi,i] z_stue[0:end_i] = DOUBLE(data2d[0:dix-1,0:diy-1])*delta_distance[zi,i] if (weighted ge 1) then w_stuetz[0:end_i] = DOUBLE(weight2d[0:dix-1,0:diy-1])*delta_distance[zi,i] xmi = LONG((x_stuetz[0:end_i] -DOUBLE(slice_size(0)))*dxts) > 0 < (xy_sl_size-1) xma = LONG((x_stuetz[0:end_i]+dx_stuetz[0:end_i]-DOUBLE(slice_size(0)))*dxts) -1 < (xy_sl_size-1) > 0 ymi = LONG((y_stuetz[0:end_i] -DOUBLE(slice_size(1)))*dyts) > 0 < (xy_sl_size-1) yma = LONG((y_stuetz[0:end_i]+dx_stuetz[0:end_i]-DOUBLE(slice_size(1)))*dyts) -1 < (xy_sl_size-1) > 0 ; not elegant or smart but works ... ; just doing a nearest neighbor projection by adding very zones ; contribution to image for k = 0L, end_i do begin xma[k] = xma[k] > xmi[k] yma[k] = yma[k] > ymi[k] slice_data[xmi[k]:xma[k],ymi[k]:yma[k]] += z_stue[k] endfor IF (weighted GE 1) THEN for k = 0L, end_i do weight_data[xmi[k]:xma[k],ymi[k]:yma[k]] += w_stuetz[k] ; now get a list of the central points of the cells that contribute. ; From this grid add only the ones taht are not covered by finer grids ; or have been already covered by grids on the same level painted previously. childs = where( u_g.level ge cg.level and $ ; all overlapping (in projections) grids u_g.num ne cg.num and $ ; on this and higher levels u_g.left_edge[xi] lt cg.right_edge[xi] and $ u_g.left_edge[yi] lt cg.right_edge[yi] and $ u_g.right_edge[xi] gt cg.left_edge[xi] and $ u_g.right_edge[yi] gt cg.left_edge[yi], count) if count gt 0 then begin ; do not consider grids on the same level that we will draw later iind = where(u_g[childs].level eq cg.level and childs gt i, count1) if count1 gt 0 then begin childs[iind] = -1 iind = where(childs ne -1, count1) if iind[0] ne -1 then childs = childs[iind] count = count1 endif endif xc = x_stuetz[0:end_i] + 0.5D * delta_distance[xi,i] yc = y_stuetz[0:end_i] + 0.5D * delta_distance[yi,i] dc = dx_stuetz[0:end_i] undefine, mind for j=0L, count-1 do begin ind = where(xc ge u_g[childs[j]].left_edge[xi] and $ xc le u_g[childs[j]].right_edge[xi] and $ yc ge u_g[childs[j]].left_edge[yi] and $ yc le u_g[childs[j]].right_edge[yi], count2) if count2 gt 0 then mind = (N_elements(mind) gt 0) ? [mind[*],ind[*]] : ind[*] ; all indeces of points covered by finer grid ; print, j, count2, ' are duplicate' endfor if N_elements(mind) gt 0 then begin ; print, N_elements(mind) HH = xc[0:end_i] HH[mind] = -1 ind = where(HH[0:end_i] ne -1, count1) ; print, 'keep:', count1 if count1 gt 0 then begin xc = xc[ind] ; keep only all the other indeces yc = yc[ind] ; keep only all the other indeces dc = dc[ind] ; keep only all the other indeces xs = (N_elements(xs) le 0) ? xc : [xs[*],xc[*]] ys = (N_elements(ys) le 0) ? yc : [ys[*],yc[*]] ds = (N_elements(ds) le 0) ? dc : [ds[*],dc[*]] endif endif else begin xs = (N_elements(xs) le 0) ? xc : [xs[*],xc[*]] ys = (N_elements(ys) le 0) ? yc : [ys[*],yc[*]] ds = (N_elements(ds) le 0) ? dc : [ds[*],dc[*]] endelse ; if (i+1 mod 5) eq 0 then $ if check_cancel() then goto, enough ENDFOR ; over grids ... ENDFOR enough: IF (weighted GE 1) THEN slice_data = slice_data/weight_data n_stuetz = N_elements(xs) x_stuetz = xs y_stuetz = ys dx_stuetz = ds ; now get points from orthogonal slice ; which pixel xmi = LONG((x_stuetz - DOUBLE(slice_size[0]))*dxts) > 0 < (xy_sl_size-1) ymi = LONG((y_stuetz - DOUBLE(slice_size[1]))*dyts) > 0 < (xy_sl_size-1) z_stuetz = x_stuetz ; for i=0L,N_elements(x_stuetz)-1 do z_stuetz[i] = slice_data[xmi[i], ymi[i]] z_stuetz = slice_data[xmi, ymi] slice_data = render_image() ; should I log it after the projection? Yes unless IF shouldbe_logged(list_str[var_index]) THEN BEGIN slice_data = ALOG10(abs(slice_data)) z_stuetz = ALOG10(abs(z_stuetz)) ; print, 'construct_orthogonal_projection: Took log of projected value' ENDIF if verbose then print, 'used ', n_stuetz, '=', SQRT(DOUBLE(n_stuetz)),'^2 points.' if verbose then print, 'projection done' RETURN END ;.compile construct_orthogonal_projection