disk_vmag(mesh) yields the visual magnification based on the projection of disks on the cortical
      surface into the visual field.

    All options accepted by mag_data() are accepted by disk_vmag().
12455
def disk_vmag(hemi, retinotopy='any', to=None, **kw):
    '''
    
    '''
    mdat = mag_data(hemi, retinotopy=retinotopy, **kw)
    if pimms.is_vector(mdat): return tuple([face_vmag(m, to=to) for m in mdat])
    elif pimms.is_vector(mdat.keys(), 'int'):
        return pimms.lazy_map({k: curry(lambda k: face_vmag(mdat[k], to=to), k)
                               for k in six.iterkeys(mdat)})
    #TODO: implement the disk_vmag calculation using mdat
    # convert to the appropriate type according to the to param
    raise NotImplementedError()
294565
def backgroundMMD(self, catalog, method='cloud-in-cells', weights=None):
        

        # Select objects in annulus
        cut_annulus = self.roi.inAnnulus(catalog.lon,catalog.lat)
        mag_1 = catalog.mag_1[cut_annulus]
        mag_2 = catalog.mag_2[cut_annulus]

        # Units are (deg^2)
        solid_angle = ugali.utils.binning.take2D(self.solid_angle_mmd, mag_2, mag_1,
                                                 self.roi.bins_mag, self.roi.bins_mag)

        # Weight each object before binning
        # Divide by solid angle and bin size in magnitudes to get number density 
        # [objs / deg^2 / mag^2]
        if weights is None:
            number_density = (solid_angle*self.roi.delta_mag**2)**(-1)
        else:
            number_density = weights*(solid_angle*self.roi.delta_mag**2)**(-1)

        method = str(method).lower()
        if method == 'cloud-in-cells':
            # Apply cloud-in-cells algorithm
            mmd_background = ugali.utils.binning.cloudInCells(mag_2,mag_1,
                                                              [self.roi.bins_mag,self.roi.bins_mag],
                                                              weights=number_density)[0]
        elif method == 'bootstrap':
            # Not implemented
            raise ValueError("Bootstrap method not implemented")
            mag_1 + (mag_1_err * np.random.normal(0, 1., len(mag_1)))
            mag_2 + (mag_2_err * np.random.normal(0, 1., len(mag_2)))

        elif method == 'histogram':
            # Apply raw histogram
            mmd_background = np.histogram2d(mag_1,mag_2,bins=[self.roi.bins_mag,self.roi.bins_mag],
                                            weights=number_density)[0]

        elif method == 'kde':
            # Gridded kernel density estimator
            logger.warning("### KDE not implemented properly")
            mmd_background = ugali.utils.binning.kernelDensity(mag_2,mag_1,
                                                               [self.roi.bins_mag,self.roi.bins_mag],
                                                               weights=number_density)[0]
        elif method == 'uniform':
            logger.warning("### WARNING: Uniform CMD")
            hist = np.histogram2d(mag_1,mag_2,bins=[self.roi.bins_mag,self.roi.bins_mag], 
                                  weights=number_density)[0]
            mmd_background = np.mean(hist)*np.ones(hist.shape)
            observable = (self.solid_angle_mmd > self.minimum_solid_angle)
            mmd_background *= observable
            return mmd_background
        else:
            raise ValueError("Unrecognized method: %s"%method)
        ## Account for the objects that spill out of the observable space
        ## But what about the objects that spill out to red colors??
        #for index_color in range(0, len(self.roi.centers_color)):
        #    for index_mag in range(0, len(self.roi.centers_mag)):
        #        if self.solid_angle_cmd[index_mag][index_color] < self.minimum_solid_angle:
        #            cmd_background[index_mag - 1][index_color] += cmd_background[index_mag][index_color]
        #            cmd_background[index_mag][index_color] = 0.
        #            break

        mmd_area = self.solid_angle_mmd*self.roi.delta_mag**2 # [deg^2 * mag^2]

        # ADW: This accounts for leakage to faint magnitudes
        # But what about the objects that spill out to red colors??
        # Maximum obsevable magnitude index for each color (uses the fact that
        # np.argmin returns first minimum (zero) instance found.
        # NOTE: More complicated maps may have holes causing problems

        observable = (self.solid_angle_mmd > self.minimum_solid_angle)
        index_mag_1 = observable.argmin(axis=0) - 1
        index_mag_2 = np.arange(len(self.roi.centers_mag))
        # Add the cumulative leakage back into the last bin of the CMD
        leakage = (mmd_background * ~observable).sum(axis=0)
        ### mmd_background[[index_mag_1,index_mag_2]] += leakage
        # Zero out all non-observable bins
        ### mmd_background *= observable

        # Avoid dividing by zero by setting empty bins to the value of the 
        # minimum filled bin of the CMD. This choice is arbitrary and 
        # could be replaced by a static minimum, some fraction of the 
        # CMD maximum, some median clipped minimum, etc. However, should 
        # be robust against outliers with very small values.
        min_mmd_background = max(mmd_background[mmd_background > 0.].min(),
                                 1e-4*mmd_background.max())
        mmd_background[observable] = mmd_background[observable].clip(min_mmd_background)

        ### # ADW: This is a fudge factor introduced to renormalize the CMD
        ### # to the number of input stars in the annulus. While leakage
        ### # will still smooth the distribution, it shouldn't result in 
        ### fudge_factor = len(mag) / float((cmd_background*cmd_area).sum())
        ### cmd_background *= fudge_factor

        return mmd_background
756334
def _createMagConversionDict():
    
    magnitude_conversion_filepath = resource_stream(__name__, 'data/magnitude_conversion.dat')
    raw_table = np.loadtxt(magnitude_conversion_filepath, '|S5')

    magDict = {}
    for row in raw_table:
        if sys.hexversion >= 0x03000000:
            starClass = row[1].decode("utf-8")  # otherwise we get byte ints or b' caused by 2to3
            tableData = [x.decode("utf-8") for x in row[3:]]
        else:
            starClass = row[1]
            tableData = row[3:]
        magDict[starClass] = tableData

    return magDict
967365
function computeMirrorDimensions({source, ...args}) {
  return withPromise((resolve) => {
    const sourceRect = source.getBoundingClientRect();
    resolve({source, sourceRect, ...args});
  });
}
1066259
function model(M, prop) {
  const bmiLens = L(
    prop,
    // if we don't have <prop> property yet, then use these default values
    L.defaults({weight: 80, height: 180}),
    // add "read-only" bmi property to our BMI model that is derived from weight and height
    L.augment({bmi: ({weight: w, height: h}) => Math.round(w / (h * h * 0.0001))})
  )
  return M.lens(bmiLens)
}
1181809
function () {

    /**
     * @type {number}
     */
    this.motionBlurFactor = 0.99;
    /**
     * Vector field lookup image
     * @type {clay.Texture2D}
     */
    this.vectorFieldTexture = new Texture2D({
        type: Texture.FLOAT,
        // minFilter: Texture.NEAREST,
        // magFilter: Texture.NEAREST,
        flipY: false
    });

    /**
     * Particle life range
     * @type {Array.<number>}
     */
    this.particleLife = [5, 20];

    this._particleType = 'point';

    /**
     * @type {number}
     */
    this._particleSize = 1;

    /**
     * @type {Array.<number>}
     */
    this.particleColor = [1, 1, 1, 1];

    /**
     * @type {number}
     */
    this.particleSpeedScaling = 1.0;

    /**
     * @type {clay.Texture2D}
     */
    this._thisFrameTexture = null;

    this._particlePass = null;
    this._spawnTexture = null;
    this._particleTexture0 = null;
    this._particleTexture1 = null;

    this._particlePointsMesh = null;

    this._surfaceFrameBuffer = null;

    this._elapsedTime = 0.0;

    this._scene = null;
    this._camera = null;

    this._lastFrameTexture = null;

    // this._temporalSS = new TemporalSS(50);

    // this._antialising = false;

    this._supersampling = 1;

    this._downsampleTextures = [];

    this._width = 512;
    this._height = 512;

    this.init();
}
835649
def _create_affine(x_axis, y_axis, z_axis, image_pos, voxel_sizes):
    

    # Create affine matrix (http://nipy.sourceforge.net/nibabel/dicom/dicom_orientation.html#dicom-slice-affine)

    affine = numpy.array(
        [[x_axis[0] * voxel_sizes[0], y_axis[0] * voxel_sizes[1], z_axis[0] * voxel_sizes[2], image_pos[0]],
         [x_axis[1] * voxel_sizes[0], y_axis[1] * voxel_sizes[1], z_axis[1] * voxel_sizes[2], image_pos[1]],
         [x_axis[2] * voxel_sizes[0], y_axis[2] * voxel_sizes[1], z_axis[2] * voxel_sizes[2], image_pos[2]],
         [0, 0, 0, 1]])
    return affine
1063385
function createVoxelMesh(gl, voxels, voxelSideTextureIDs, voxelSideTextureSizes, position, pad, that) {
  //Create mesh
  var vert_data = createAOMesh(voxels, voxelSideTextureIDs, voxelSideTextureSizes)
  var vertexArrayObjects = {}

  if (vert_data === null) {
    // no vertices allocated
  } else {
    //Upload triangle mesh to WebGL
    var vert_buf = createBuffer(gl, vert_data)
    var triangleVAO = createVAO(gl, [
      { "buffer": vert_buf,
        "type": gl.UNSIGNED_BYTE,
        "size": 4,
        "offset": 0,
        "stride": 8,
        "normalized": false
      },
      { "buffer": vert_buf,
        "type": gl.UNSIGNED_BYTE,
        "size": 4,
        "offset": 4,
        "stride": 8,
        "normalized": false
      }
    ])
    triangleVAO.length = Math.floor(vert_data.length/8)

    vertexArrayObjects.surface = triangleVAO
  }

  // move the chunk into place
  var modelMatrix = mat4.create()
  var w = voxels.shape[2] - pad  // =[1]=[0]=game.chunkSize
  var translateVector = [
    position[0] * w,
    position[1] * w,
    position[2] * w]

  mat4.translate(modelMatrix, modelMatrix, translateVector)

  //Bundle result and return
  var result = {
    vertexArrayObjects: vertexArrayObjects, // other plugins can add their own VAOs
    center: [w>>1, w>>1, w>>1],
    radius: w,
    modelMatrix: modelMatrix
  }

  if (that) that.emit('meshed', result, gl, vert_data, voxels)

  return result
}
799169
def getVoxelsScalar(self, vmin=None, vmax=None):
                
        nx, ny, nz = self.image.GetDimensions()
        
        voxdata = np.zeros([nx, ny, nz])
        lsx = range(nx)
        lsy = range(ny)
        lsz = range(nz)
        renorm = vmin is not None and vmax is not None
        for i in lsx:
            for j in lsy:
                for k in lsz:
                    s = self.image.GetScalarComponentAsFloat(i, j, k, 0)
                    if renorm:
                        s = (s - vmin) / (vmax - vmin)
                    voxdata[i, j, k] = s
        return voxdata
671552
def showSegmentation(
        segmentation,
        voxelsize_mm=np.ones([3, 1]),
        degrad=4,
        label=1,
        smoothing=True
        ):
    
    labels = []

    segmentation = segmentation[::degrad, ::degrad, ::degrad]

    # import pdb; pdb.set_trace()
    mesh_data = seg2fem.gen_mesh_from_voxels_mc(segmentation, voxelsize_mm*degrad)
    if smoothing:
        mesh_data.coors = seg2fem.smooth_mesh(mesh_data)
    else:
        mesh_data = seg2fem.gen_mesh_from_voxels_mc(segmentation, voxelsize_mm * 1.0e-2)
        # mesh_data.coors +=
    vtk_file = "mesh_geom.vtk"
    mesh_data.write(vtk_file)
    QApplication(sys.argv)
    view = viewer.QVTKViewer(vtk_file)
    view.exec_()

    return labels
769664
def recompute_mag(mpc_in, skip_centroids=False):
    
    # TODO this really shouldn't need to build a 'reading' to get the cutout...

    from ossos.downloads.cutouts import downloader
    dlm = downloader.ImageCutoutDownloader()

    mpc_obs = deepcopy(mpc_in)
    assert isinstance(mpc_obs, mp_ephem.ephem.Observation)
    assert isinstance(mpc_obs.comment, mp_ephem.ephem.OSSOSComment)

    if mpc_obs.null_observation:
        return mpc_obs

    parts = re.search('(?P<expnum>\d{7})(?P<type>\S)(?P<ccd>\d\d)', mpc_obs.comment.frame)
    if parts is None:
        return mpc_obs
    expnum = parts.group('expnum')
    ccd = parts.group('ccd')
    file_type = parts.group('type')

    observation = astrom.Observation(expnum, file_type, ccd)
    assert isinstance(observation, astrom.Observation)

    ast_header = _connection_error_wrapper(storage._get_sghead, int(expnum))[int(ccd)+1]

    filter_value = None
    for keyword in ['FILTER', 'FILT1 NAME']:
       filter_value = ast_header.get(keyword, None)
       if filter_value is not None:
         if filter_value.startswith('gri'):
             filter_value = 'w'
         else:
             filter_value = filter_value[0]
         break
    # The ZP for the current astrometric lines is the pipeline one.  The new ZP is in the astheader file.
    new_zp = ast_header.get('PHOTZP')

    # The .zeropoint.used value is likely the one used for the original photometry.
    old_zp = _connection_error_wrapper(storage.get_zeropoint, int(expnum), int(ccd))

    reading = astrom.SourceReading(float(mpc_obs.comment.x), float(mpc_obs.comment.y), float(mpc_obs.comment.x),
                                   float(mpc_obs.comment.y), mpc_obs.coordinate.ra.degree,
                                   mpc_obs.coordinate.dec.degree, float(mpc_obs.comment.x), float(mpc_obs.comment.y),
                                   observation, ssos=True, from_input_file=True, null_observation=False,
                                   discovery=mpc_obs.discovery)

    cutout = _connection_error_wrapper(dlm.download_cutout, reading, needs_apcor=True)
    cutout._zmag = new_zp

    if math.fabs(cutout.zmag - old_zp) > 0.3:
        logging.warn("Large change in zeropoint detected: {}  -> {}".format(old_zp, cutout.zmag))

    try:
        PHOT = cutout.get_observed_magnitude(centroid=not skip_centroids and mpc_obs.note1 != "H")
        x = PHOT['XCENTER']
        y = PHOT['YCENTER']
        mag = PHOT['MAG']
        merr = PHOT['MERR']
        cutout.update_pixel_location((x, y), hdu_index=cutout.extno)
        x, y = cutout.observed_source_point
    except Exception as ex:
        logging.error("ERROR: {}".format(str(ex)))
        return mpc_obs

    try:
        if mpc_obs.comment.mag_uncertainty is not None and mpc_obs.comment.mag is not None and math.fabs(mpc_obs.comment.mag - mag) > 3.5 * mpc_obs.comment.mag_uncertainty:
           logging.warn("recomputed magnitude shift large: {} --> {}".format(mpc_obs.mag, mag[0]))
        if math.sqrt((x.value - mpc_obs.comment.x) ** 2 + (y.value - mpc_obs.comment.y) ** 2) > 1.9:
            logging.warn("Centroid shifted ({},{}) -> ({},{})".format(mpc_obs.comment.x,
                                                                      mpc_obs.comment.y,
                                                                      x.value,
                                                                      y.value))
    except Exception as ex:
        logging.error(str(ex))

    # Don't use the new X/Y for Hand measured entries.  (although call to get_observed_magnitude should have changed)
    if str(mpc_obs.note1) != "H" and not skip_centroids:
        mpc_obs.comment.x = x.value
        mpc_obs.comment.y = y.value

    try:
        mag = float(mag)
    except:
        return mpc_obs

    if math.isnan(mag):
        return mpc_obs

    if mag > 10:
        mpc_obs._band = filter_value
        mpc_obs.comment.mag = mag
        mpc_obs.comment.mag_uncertainty = merr

    # Update the mpc record magnitude if previous value existed here.
    if (mpc_obs.mag is not None or (mpc_obs.mag is None and mpc_in.comment.photometry_note[0] == "Z")) and mag > 10:
        mpc_obs.mag = mag

    return mpc_obs
769838
def _get_supercell(obj, rCut=5.0):
    rCutHard = rCut + 5  # Giving extra space for hard cutOff
    

    cell_vectors = obj.get_cell()
    a1, a2, a3 = cell_vectors[0], cell_vectors[1], cell_vectors[2]

    # vectors perpendicular to two cell vectors
    b1 = np.cross(a2, a3, axis=0)
    b2 = np.cross(a3, a1, axis=0)
    b3 = np.cross(a1, a2, axis=0)
    # projections onto perpendicular vectors
    p1 = np.dot(a1, b1) / np.dot(b1, b1) * b1
    p2 = np.dot(a2, b2) / np.dot(b2, b2) * b2
    p3 = np.dot(a3, b3) / np.dot(b3, b3) * b3
    xyz_arr = np.linalg.norm(np.array([p1, p2, p3]), axis=1)
    cell_images = np.ceil(rCutHard/xyz_arr)
    nx = int(cell_images[0])
    ny = int(cell_images[1])
    nz = int(cell_images[2])
    suce = obj * (1+2*nx, 1+2*ny, 1+2*nz)
    shift = obj.get_cell()

    shifted_suce = suce.copy()
    shifted_suce.translate(-shift[0]*nx - shift[1]*ny - shift[2]*nz)

    return shifted_suce
445743
private static void convolve(double[] xreal, double[] ximag, double[] yreal, double[] yimag, double[] outreal, double[] outimag) {
//        if (xreal.length != ximag.length || xreal.length != yreal.length || yreal.length != yimag.length || xreal.length != outreal.length || outreal.length != outimag.length)
//            throw new IllegalArgumentException("Mismatched lengths");

        int n = xreal.length;

        FFT(xreal, ximag);
        FFT(yreal, yimag);
        for (int i = 0; i < n; i++) {
            double temp = xreal[i] * yreal[i] - ximag[i] * yimag[i];
            ximag[i] = ximag[i] * yreal[i] + xreal[i] * yimag[i];
            xreal[i] = temp;
        }
        inverseTransform(xreal, ximag);

        // Scaling (because this FFT implementation omits it)
        for (int i = 0; i < n; i++) {
            outreal[i] = xreal[i] / n;
            outimag[i] = ximag[i] / n;
        }
    }
618504
def _extract_gaussian_gradient_magnitude(image, mask = slice(None), sigma = 1, voxelspacing = None):
    
    # set voxel spacing
    if voxelspacing is None:
        voxelspacing = [1.] * image.ndim
        
    # determine gaussian kernel size in voxel units
    sigma = _create_structure_array(sigma, voxelspacing)
        
    return _extract_intensities(scipy_gaussian_gradient_magnitude(image, sigma), mask)
1024153
def magnification(self, x, y, kwargs, diff=diff):
        
        f_xx, f_xy, f_yx, f_yy = self.hessian(x, y, kwargs, diff=diff)
        det_A = (1 - f_xx) * (1 - f_yy) - f_xy*f_yx
        return 1/det_A
1219553
public static void convolveHV(Kernel kernel, int[] inPixels, int[] outPixels, int width, int height, boolean alpha, int edgeAction) {
		int index = 0;
		float[] matrix = kernel.getKernelData( null );
		int rows = kernel.getHeight();
		int cols = kernel.getWidth();
		int rows2 = rows/2;
		int cols2 = cols/2;

		for (int y = 0; y < height; y++) {
			for (int x = 0; x < width; x++) {
				float r = 0, g = 0, b = 0, a = 0;

				for (int row = -rows2; row <= rows2; row++) {
					int iy = y+row;
					int ioffset;
					if (0 <= iy && iy < height)
						ioffset = iy*width;
					else if ( edgeAction == CLAMP_EDGES )
						ioffset = y*width;
					else if ( edgeAction == WRAP_EDGES )
						ioffset = ((iy+height) % height) * width;
					else
						continue;
					int moffset = cols*(row+rows2)+cols2;
					for (int col = -cols2; col <= cols2; col++) {
						float f = matrix[moffset+col];

						if (f != 0) {
							int ix = x+col;
							if (!(0 <= ix && ix < width)) {
								if ( edgeAction == CLAMP_EDGES )
									ix = x;
								else if ( edgeAction == WRAP_EDGES )
									ix = (x+width) % width;
								else
									continue;
							}
							int rgb = inPixels[ioffset+ix];
							a += f * ((rgb >> 24) & 0xff);
							r += f * ((rgb >> 16) & 0xff);
							g += f * ((rgb >> 8) & 0xff);
							b += f * (rgb & 0xff);
						}
					}
				}
				int ia = alpha ? PixelUtils.clamp((int)(a+0.5)) : 0xff;
				int ir = PixelUtils.clamp((int)(r+0.5));
				int ig = PixelUtils.clamp((int)(g+0.5));
				int ib = PixelUtils.clamp((int)(b+0.5));
				outPixels[index++] = (ia << 24) | (ir << 16) | (ig << 8) | ib;
			}
		}
	}
705441
def bInit( self, vrModel, vrDiffuseTexture ):
        ""
        # create and bind a VAO to hold state for this model
        self.m_glVertArray = glGenVertexArrays(1)
        glBindVertexArray( self.m_glVertArray )
        # Populate a vertex buffer
        self.m_glVertBuffer = glGenBuffers(1)
        glBindBuffer( GL_ARRAY_BUFFER, self.m_glVertBuffer )
        glBufferData( GL_ARRAY_BUFFER, sizeof( openvr.RenderModel_Vertex_t ) * vrModel.unVertexCount, vrModel.rVertexData, GL_STATIC_DRAW )
        # Identify the components in the vertex buffer
        glEnableVertexAttribArray( 0 )
        glVertexAttribPointer( 0, 3, GL_FLOAT, False, sizeof( openvr.RenderModel_Vertex_t ), openvr.RenderModel_Vertex_t.vPosition.offset )
        glEnableVertexAttribArray( 1 )
        glVertexAttribPointer( 1, 3, GL_FLOAT, False, sizeof( openvr.RenderModel_Vertex_t ), openvr.RenderModel_Vertex_t.vNormal.offset )
        glEnableVertexAttribArray( 2 )
        glVertexAttribPointer( 2, 2, GL_FLOAT, False, sizeof( openvr.RenderModel_Vertex_t ), openvr.RenderModel_Vertex_t.rfTextureCoord.offset )
        # Create and populate the index buffer
        self.m_glIndexBuffer = glGenBuffers(1)
        glBindBuffer( GL_ELEMENT_ARRAY_BUFFER, self.m_glIndexBuffer )
        glBufferData( GL_ELEMENT_ARRAY_BUFFER, sizeof( ctypes.c_uint16 ) * vrModel.unTriangleCount * 3, vrModel.rIndexData, GL_STATIC_DRAW )
        glBindVertexArray( 0 )
        # create and populate the texture
        self.m_glTexture = glGenTextures(1)
        glBindTexture( GL_TEXTURE_2D, self.m_glTexture )
        glTexImage2D( GL_TEXTURE_2D, 0, GL_RGBA, vrDiffuseTexture.unWidth, vrDiffuseTexture.unHeight,
            0, GL_RGBA, GL_UNSIGNED_BYTE, vrDiffuseTexture.rubTextureMapData )
        # If this renders black ask McJohn what's wrong.
        glGenerateMipmap(GL_TEXTURE_2D)
        glTexParameteri( GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE )
        glTexParameteri( GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE )
        glTexParameteri( GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR )
        glTexParameteri( GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR_MIPMAP_LINEAR )
        fLargest = glGetFloatv( GL_MAX_TEXTURE_MAX_ANISOTROPY_EXT )
        glTexParameterf( GL_TEXTURE_2D, GL_TEXTURE_MAX_ANISOTROPY_EXT, fLargest )
        glBindTexture( GL_TEXTURE_2D, 0 )
        self.m_unVertexCount = vrModel.unTriangleCount * 3
        return True