Changeset f549e37 in sasmodels for sasmodels/data.py
 Timestamp:
 Mar 4, 2018 8:36:10 PM (5 years ago)
 Branches:
 master, core_shell_microgels, magnetic_model, ticket1257vesicleproduct, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
 Children:
 d18d6dd
 Parents:
 fc3ae1b
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

sasmodels/data.py
r2d81cfe rf549e37 706 706 else: 707 707 vmin_scaled, vmax_scaled = vmin, vmax 708 plt.imshow(plottable.reshape(len(data.x_bins), len(data.y_bins)), 708 nx, ny = len(data.x_bins), len(data.y_bins) 709 x_bins, y_bins, image = _build_matrix(data, plottable) 710 plt.imshow(image, 709 711 interpolation='nearest', aspect=1, origin='lower', 710 712 extent=[xmin, xmax, ymin, ymax], … … 714 716 return vmin, vmax 715 717 718 719 # === The following is modified from sas.sasgui.plottools.PlotPanel 720 def _build_matrix(self, plottable): 721 """ 722 Build a matrix for 2d plot from a vector 723 Returns a matrix (image) with ~ square binning 724 Requirement: need 1d array formats of 725 self.data, self.qx_data, and self.qy_data 726 where each one corresponds to z, x, or y axis values 727 728 """ 729 # No qx or qy given in a vector format 730 if self.qx_data is None or self.qy_data is None \ 731 or self.qx_data.ndim != 1 or self.qy_data.ndim != 1: 732 return self.x_bins, self.y_bins, plottable 733 734 # maximum # of loops to fillup_pixels 735 # otherwise, loop could never stop depending on data 736 max_loop = 1 737 # get the x and y_bin arrays. 738 x_bins, y_bins = _get_bins(self) 739 # set zero to None 740 741 #Note: Can not use scipy.interpolate.Rbf: 742 # 'cause too many data points (>10000)<=JHC. 743 # 1d array to use for weighting the data point averaging 744 #when they fall into a same bin. 745 weights_data = np.ones([self.data.size]) 746 # get histogram of ones w/len(data); this will provide 747 #the weights of data on each bins 748 weights, xedges, yedges = np.histogram2d(x=self.qy_data, 749 y=self.qx_data, 750 bins=[y_bins, x_bins], 751 weights=weights_data) 752 # get histogram of data, all points into a bin in a way of summing 753 image, xedges, yedges = np.histogram2d(x=self.qy_data, 754 y=self.qx_data, 755 bins=[y_bins, x_bins], 756 weights=plottable) 757 # Now, normalize the image by weights only for weights>1: 758 # If weight == 1, there is only one data point in the bin so 759 # that no normalization is required. 760 image[weights > 1] = image[weights > 1] / weights[weights > 1] 761 # Set image bins w/o a data point (weight==0) as None (was set to zero 762 # by histogram2d.) 763 image[weights == 0] = None 764 765 # Fill empty bins with 8 nearest neighbors only when at least 766 #one None point exists 767 loop = 0 768 769 # do while loop until all vacant bins are filled up up 770 #to loop = max_loop 771 while (weights == 0).any(): 772 if loop >= max_loop: # this protects neverending loop 773 break 774 image = _fillup_pixels(self, image=image, weights=weights) 775 loop += 1 776 777 return x_bins, y_bins, image 778 779 def _get_bins(self): 780 """ 781 get bins 782 set x_bins and y_bins into self, 1d arrays of the index with 783 ~ square binning 784 Requirement: need 1d array formats of 785 self.qx_data, and self.qy_data 786 where each one corresponds to x, or y axis values 787 """ 788 # find max and min values of qx and qy 789 xmax = self.qx_data.max() 790 xmin = self.qx_data.min() 791 ymax = self.qy_data.max() 792 ymin = self.qy_data.min() 793 794 # calculate the range of qx and qy: this way, it is a little 795 # more independent 796 x_size = xmax  xmin 797 y_size = ymax  ymin 798 799 # estimate the # of pixels on each axes 800 npix_y = int(np.floor(np.sqrt(len(self.qy_data)))) 801 npix_x = int(np.floor(len(self.qy_data) / npix_y)) 802 803 # bin size: x & ydirections 804 xstep = x_size / (npix_x  1) 805 ystep = y_size / (npix_y  1) 806 807 # max and min taking account of the bin sizes 808 xmax = xmax + xstep / 2.0 809 xmin = xmin  xstep / 2.0 810 ymax = ymax + ystep / 2.0 811 ymin = ymin  ystep / 2.0 812 813 # store x and y bin centers in q space 814 x_bins = np.linspace(xmin, xmax, npix_x) 815 y_bins = np.linspace(ymin, ymax, npix_y) 816 817 return x_bins, y_bins 818 819 def _fillup_pixels(self, image=None, weights=None): 820 """ 821 Fill z values of the empty cells of 2d image matrix 822 with the average over upto next nearest neighbor points 823 824 :param image: (2d matrix with some zi = None) 825 826 :return: image (2d array ) 827 828 :TODO: Find better way to do forloop below 829 830 """ 831 # No image matrix given 832 if image is None or np.ndim(image) != 2 \ 833 or np.isfinite(image).all() \ 834 or weights is None: 835 return image 836 # Get bin size in y and x directions 837 len_y = len(image) 838 len_x = len(image[1]) 839 temp_image = np.zeros([len_y, len_x]) 840 weit = np.zeros([len_y, len_x]) 841 # do forloop for all pixels 842 for n_y in range(len(image)): 843 for n_x in range(len(image[1])): 844 # find only null pixels 845 if weights[n_y][n_x] > 0 or np.isfinite(image[n_y][n_x]): 846 continue 847 else: 848 # find 4 nearest neighbors 849 # check where or not it is at the corner 850 if n_y != 0 and np.isfinite(image[n_y  1][n_x]): 851 temp_image[n_y][n_x] += image[n_y  1][n_x] 852 weit[n_y][n_x] += 1 853 if n_x != 0 and np.isfinite(image[n_y][n_x  1]): 854 temp_image[n_y][n_x] += image[n_y][n_x  1] 855 weit[n_y][n_x] += 1 856 if n_y != len_y  1 and np.isfinite(image[n_y + 1][n_x]): 857 temp_image[n_y][n_x] += image[n_y + 1][n_x] 858 weit[n_y][n_x] += 1 859 if n_x != len_x  1 and np.isfinite(image[n_y][n_x + 1]): 860 temp_image[n_y][n_x] += image[n_y][n_x + 1] 861 weit[n_y][n_x] += 1 862 # go 4 next nearest neighbors when no nonzero 863 # neighbor exists 864 if n_y != 0 and n_x != 0 and \ 865 np.isfinite(image[n_y  1][n_x  1]): 866 temp_image[n_y][n_x] += image[n_y  1][n_x  1] 867 weit[n_y][n_x] += 1 868 if n_y != len_y  1 and n_x != 0 and \ 869 np.isfinite(image[n_y + 1][n_x  1]): 870 temp_image[n_y][n_x] += image[n_y + 1][n_x  1] 871 weit[n_y][n_x] += 1 872 if n_y != len_y and n_x != len_x  1 and \ 873 np.isfinite(image[n_y  1][n_x + 1]): 874 temp_image[n_y][n_x] += image[n_y  1][n_x + 1] 875 weit[n_y][n_x] += 1 876 if n_y != len_y  1 and n_x != len_x  1 and \ 877 np.isfinite(image[n_y + 1][n_x + 1]): 878 temp_image[n_y][n_x] += image[n_y + 1][n_x + 1] 879 weit[n_y][n_x] += 1 880 881 # get it normalized 882 ind = (weit > 0) 883 image[ind] = temp_image[ind] / weit[ind] 884 885 return image 886 887 716 888 def demo(): 717 889 # type: () > None
Note: See TracChangeset
for help on using the changeset viewer.