Changeset f549e37 in sasmodels for sasmodels/data.py


Ignore:
Timestamp:
Mar 4, 2018 8:36:10 PM (7 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
d18d6dd
Parents:
fc3ae1b
Message:

allow plotting of non-gridded data in sasmodels

File:
1 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/data.py

    r2d81cfe rf549e37  
    706706    else: 
    707707        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, 
    709711               interpolation='nearest', aspect=1, origin='lower', 
    710712               extent=[xmin, xmax, ymin, ymax], 
     
    714716    return vmin, vmax 
    715717 
     718 
     719# === The following is modified from sas.sasgui.plottools.PlotPanel 
     720def _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 never-ending loop 
     773            break 
     774        image = _fillup_pixels(self, image=image, weights=weights) 
     775        loop += 1 
     776 
     777    return x_bins, y_bins, image 
     778 
     779def _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- & y-directions 
     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 
     819def _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 up-to 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 for-loop 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 for-loop 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 non-zero 
     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 
    716888def demo(): 
    717889    # type: () -> None 
Note: See TracChangeset for help on using the changeset viewer.