def load_spc(fname):
    """Load data from Becker&Hickl SPC files.

    Returns:
        3 numpy arrays: timestamps, detector, nanotime
    """
    spc_dtype = np.dtype([("field0", "<u2"), ("b", "<u1"), ("c", "<u1"), ("a", "<u2")])
    data = np.fromfile(fname, dtype=spc_dtype)

    nanotime = 4095 - np.bitwise_and(data["field0"], 0x0FFF)
    detector = data["c"]

    # Build the macrotime (timestamps) using in-place operation for efficiency
    timestamps = data["b"].astype("int64")
    np.left_shift(timestamps, 16, out=timestamps)
    timestamps += data["a"]

    # extract the 13-th bit from data['field0']
    overflow = np.bitwise_and(np.right_shift(data["field0"], 13), 1)
    overflow = np.cumsum(overflow, dtype="int64")

    # Add the overflow bits
    timestamps += np.left_shift(overflow, 24)

    return timestamps, detector, nanotime

    def _make_DO_array(self, data, channels):
        """ Get the port ordering in the final integer
        Parameters
        ----------
        data: dict
            Mapping from channel names to per-channel data arrays. Each array is a series of
            integer samples. Each sample is an integer representation of the output state of the
            corresponding channel.
        """
        ports = []
        for ch in channels:
            for port_name, line_name in ch.line_pairs:
                if port_name not in ports:
                    ports.append(port_name)

        # Final int array
        out = np.zeros(len(data.values()[0]), dtype=np.uint32)

        for ch in channels:
            arr = data[ch.name]
            for i, (port_name, line_name) in enumerate(ch.line_pairs):
                line_num = int(line_name[4:])
                bits = np.bitwise_and(arr, (1 << i))  # Mask out the user-input bits

                left_shift_amount = line_num - i
                if left_shift_amount > 0:
                    byts = np.left_shift(bits, left_shift_amount)
                elif left_shift_amount < 0:
                    byts = np.right_shift(bits, -left_shift_amount)
                else:
                    byts = bits

                byte_num = ports.index(port_name)
                out += np.left_shift(byts, 8*byte_num)
        return out

    def __invertibleToRGB(self, rgbVarr, varrIdx, colorLutStruct):
        '''
        Decode an RGB image rendered in INVERTIBLE_LUT mode. The image encodes
        float values as colors, so the RGB value is first decoded into its
        represented float value and then the color table is applied.
        '''
        w0 = np.left_shift(
            rgbVarr[varrIdx[0], varrIdx[1], 0].astype(np.uint32), 16)
        w1 = np.left_shift(
            rgbVarr[varrIdx[0], varrIdx[1], 1].astype(np.uint32), 8)
        w2 = rgbVarr[varrIdx[0], varrIdx[1], 2]

        value = np.bitwise_or(w0, w1)
        value = np.bitwise_or(value, w2)
        value = np.subtract(value.astype(np.int32), 1)
        normalized_val = np.divide(value.astype(float), 0xFFFFFE)

        # Map float value to color lut (use a histogram to support non-uniform
        # colormaps (fetch bin indices))
        colorLut = colorLutStruct.lut
        bins = colorLutStruct.adjustedBins

        idx = np.digitize(normalized_val, bins)
        idx = np.subtract(idx, 1)

        valueImage = np.zeros([rgbVarr.shape[0], rgbVarr.shape[1]],
                              dtype=np.uint32)
        valueImage[varrIdx[0], varrIdx[1]] = idx

        return colorLut[valueImage]

def _dec10216(inbuf):
    inbuf = np.fromstring(inbuf, dtype=np.uint8)
    arr10 = inbuf.astype(np.uint16)
    arr16 = np.zeros((len(arr10) / 5 * 4,), dtype=np.uint16)
    arr10_len = (len(arr16) * 5) / 4
    arr10 = arr10[:arr10_len] # adjust size
    """
    /*
     * pack 4 10-bit words in 5 bytes into 4 16-bit words
     * 
     * 0       1       2       3       4       5
     * 01234567890123456789012345678901234567890
     * 0         1         2         3         4
     */      
    ip = &in_buffer[i];
    op = &out_buffer[j];
    op[0] = ip[0]*4 + ip[1]/64;
    op[1] = (ip[1] & 0x3F)*16 + ip[2]/16;
    op[2] = (ip[2] & 0x0F)*64 + ip[3]/4;
    op[3] = (ip[3] & 0x03)*256 +ip[4];
    """
    arr16.flat[::4] = np.left_shift(arr10[::5], 2) + \
        np.right_shift((arr10[1::5]), 6)
    arr16.flat[1::4] = np.left_shift((arr10[1::5] & 63), 4) + \
        np.right_shift((arr10[2::5]), 4)
    arr16.flat[2::4] = np.left_shift(arr10[2::5] & 15, 6) + \
        np.right_shift((arr10[3::5]), 2)
    arr16.flat[3::4] = np.left_shift(arr10[3::5] & 3, 8) + \
        arr10[4::5]    
    return arr16.tostring()

    def rearrange_bits(array):
        # Do bit rearrangement for the 10-bit lytro raw format
        # Normalize output to 1.0 as float64
        t0 = array[0::5]
        t1 = array[1::5]
        t2 = array[2::5]
        t3 = array[3::5]
        lsb = array[4::5]

        t0 = np.left_shift(t0, 2) + np.bitwise_and(lsb, 3)
        t1 = np.left_shift(t1, 2) + np.right_shift(np.bitwise_and(lsb, 12), 2)
        t2 = np.left_shift(t2, 2) + np.right_shift(np.bitwise_and(lsb, 48), 4)
        t3 = np.left_shift(t3, 2) + np.right_shift(np.bitwise_and(lsb, 192), 6)

        image = np.zeros(LYTRO_ILLUM_IMAGE_SIZE, dtype=np.uint16)
        image[:, 0::4] = t0.reshape(
            (LYTRO_ILLUM_IMAGE_SIZE[0], LYTRO_ILLUM_IMAGE_SIZE[1] // 4)
        )
        image[:, 1::4] = t1.reshape(
            (LYTRO_ILLUM_IMAGE_SIZE[0], LYTRO_ILLUM_IMAGE_SIZE[1] // 4)
        )
        image[:, 2::4] = t2.reshape(
            (LYTRO_ILLUM_IMAGE_SIZE[0], LYTRO_ILLUM_IMAGE_SIZE[1] // 4)
        )
        image[:, 3::4] = t3.reshape(
            (LYTRO_ILLUM_IMAGE_SIZE[0], LYTRO_ILLUM_IMAGE_SIZE[1] // 4)
        )

        # Normalize data to 1.0 as 64-bit float.
        # Division is by 1023 as the Lytro Illum saves 10-bit raw data.
        return np.divide(image, 1023.0).astype(np.float64)

  def ingest ( self ):
    """Read the stack and ingest"""
    
    with closing (ocpcaproj.OCPCAProjectsDB()) as projdb:
      proj = projdb.loadToken(self.token)
    
    with closing (ocpcadb.OCPCADB(proj)) as db:

      ch = proj.getChannelObj(self.channel)
      # get the dataset configuration
      [[ximagesz, yimagesz, zimagesz],(starttime,endtime)] = proj.datasetcfg.imageSize(self.resolution)
      [xcubedim, ycubedim, zcubedim] = cubedim = proj.datasetcfg.getCubeDims()[self.resolution]
      [xoffset, yoffset, zoffset] = proj.datasetcfg.getOffset()[self.resolution]
    
      # for all specified resolutions
      for resolution in range(0,1,1):

        # extract parameters for iteration
        numxtiles = ximagesz/self.tilesz[0]
        numytiles = yimagesz/self.tilesz[1]

        # Ingest in database aligned slabs in the z dimension
        for slice_number in range(0, zimagesz, zcubedim):

          slab = np.zeros ( [zcubedim,yimagesz,ximagesz], dtype=np.uint32 )
          # over all tiles in that slice
          for b in range(zcubedim):
            for ytile in range(numytiles):
              for xtile in range(numxtiles):

                # if we are at the end of the space, quit
                if slice_number+b <= zimagesz:
                  try:
                    filename = '{}{}/{}/{}/{}.png'.format(self.tilepath, resolution, slice_number+b+zoffset, ytile+17, xtile+16)
                    print "Opening filename {}".format(filename)
                    # add tile to stack
                    imgdata = np.asarray ( Image.open(filename, 'r').convert('RGBA') )
                    imgdata = np.left_shift(imgdata[:,:,3], 24, dtype=np.uint32) | np.left_shift(imgdata[:,:,2], 16, dtype=np.uint32) | np.left_shift(imgdata[:,:,1], 8, dtype=np.uint32) | np.uint32(imgdata[:,:,0])
                    slab [b,ytile*self.tilesz[1]:(ytile+1)*self.tilesz[1],xtile*self.tilesz[0]:(xtile+1)*self.tilesz[0]] = imgdata
                  except IOError, e:
                    print "Failed to open file {}".format(filename)
                    slab [b,ytile*self.tilesz[1]:(ytile+1)*self.tilesz[1],xtile*self.tilesz[0]:(xtile+1)*self.tilesz[0]] = np.zeros([self.tilesz[1], self.tilesz[0]], dtype=np.uint32)

          for y in range (0, yimagesz+1, ycubedim):
            for x in range (0, ximagesz+1, xcubedim):
              
              # getting the cube id and ingesting the data one cube at a time
              zidx = ocplib.XYZMorton ([x/xcubedim, y/ycubedim, (slice_number)/zcubedim])
              cube = Cube.getCube(cubedim, ch.getChannelType(), ch.getDataType())
              cube.zeros()

              xmin, ymin = x, y
              xmax = min (ximagesz, x+xcubedim)
              ymax = min (yimagesz, y+ycubedim)
              zmin = 0
              zmax = min(slice_number+zcubedim, zimagesz+1)
              cube.data[0:zmax-zmin,0:ymax-ymin,0:xmax-xmin] = slab[zmin:zmax, ymin:ymax, xmin:xmax]
              
              if cube.isNotZeros():
                db.putCube(ch, zidx, self.resolution, cube, update=True)

def _wrapx(input, output, nx):
    """
    Wrap the X format column Boolean array into an ``UInt8`` array.

    Parameters
    ----------
    input
        input Boolean array of shape (`s`, `nx`)

    output
        output ``Uint8`` array of shape (`s`, `nbytes`)

    nx
        number of bits
    """

    output[...] = 0  # reset the output
    nbytes = ((nx - 1) // 8) + 1
    unused = nbytes * 8 - nx
    for i in range(nbytes):
        _min = i * 8
        _max = min((i + 1) * 8, nx)
        for j in range(_min, _max):
            if j != _min:
                np.left_shift(output[..., i], 1, output[..., i])
            np.add(output[..., i], input[..., j], output[..., i])

    # shift the unused bits
    np.left_shift(output[..., i], unused, output[..., i])

    def _make_DO_array(self, data, channels):
        """ Get the port ordering in the final integer """
        ports = []
        for ch in channels:
            for port_name, line_name in ch.line_pairs:
                if port_name not in ports:
                    ports.append(port_name)

        # Final int array
        out = np.zeros(len(data.values()[0]), dtype=np.uint32)

        for ch in channels:
            arr = data[ch.name]
            for i, (port_name, line_name) in enumerate(ch.line_pairs):
                line_num = int(line_name[4:])
                bits = np.bitwise_and(arr, (1 << i))  # Mask out the user-input bits

                left_shift_amount = line_num - i
                if left_shift_amount > 0:
                    byts = np.left_shift(bits, left_shift_amount)
                elif left_shift_amount < 0:
                    byts = np.right_shift(bits, -left_shift_amount)
                else:
                    byts = bits

                byte_num = ports.index(port_name)
                out += np.left_shift(byts, 8*byte_num)
        return out

def amagacolor(secret, red, green, blue):
    red = np.left_shift(np.right_shift(red,2),2)
    green = np.left_shift(np.right_shift(green,2),2)
    blue = np.left_shift(np.right_shift(blue,2),2)
    
    secretred,secretgreen,secretblue = separaimatge(secret,2,2,2)
    return red+secretred, green+secretgreen, blue+secretblue

def load_spc(fname):
    """Load data from Becker & Hickl SPC files.

    Returns:
        3 numpy arrays (timestamps, detector, nanotime) and a float
        (timestamps_unit).
    """

    f = open(fname, 'rb')
    # We first decode the first 6 bytes which is a header...
    header = np.fromfile(f, dtype='u2', count=3)
    timestamps_unit = header[1] * 0.1e-9
    num_routing_bits = np.bitwise_and(header[0], 0x000F)  # unused

    # ...and then the remaining records containing the photon data
    spc_dtype = np.dtype([('field0', '<u2'), ('b', '<u1'), ('c', '<u1'),
                          ('a', '<u2')])
    data = np.fromfile(f, dtype=spc_dtype)

    nanotime =  4095 - np.bitwise_and(data['field0'], 0x0FFF)
    detector = data['c']

    # Build the macrotime (timestamps) using in-place operation for efficiency
    timestamps = data['b'].astype('int64')
    np.left_shift(timestamps, 16, out=timestamps)
    timestamps += data['a']

    # extract the 13-th bit from data['field0']
    overflow = np.bitwise_and(np.right_shift(data['field0'], 13), 1)
    overflow = np.cumsum(overflow, dtype='int64')

    # Add the overflow bits
    timestamps += np.left_shift(overflow, 24)

    return timestamps, detector, nanotime, timestamps_unit

def _3dby8toRGB ( indata ):
  """Convert a numpy array of 3d, 8-bit data to 32bit RGB"""

  rgbdata = np.zeros ( indata.shape[0:2], dtype=np.uint32)

  rgbdata = np.uint32(0xFF000000) + np.left_shift(np.uint32(indata[:,:,:,0]),16) + np.left_shift(np.uint32(indata[:,:,:,1]),8) + np.uint32(indata[:,:,:,2])

  return rgbdata

 def _applytwo(self, f1, f2, offset, poses, num, rem):
     bts1 = f1 if self._subset is None else f1[self._subset]
     bts2 = f2 if self._subset is None else f2[self._subset]
     mask = 2**(num-rem)-1
     out = np.bitwise_and(bts2, mask)
     np.left_shift(out, rem, out)
     out += np.right_shift(bts1, offset)
     return out

def _read_3(fid):
    """ Read 3 byte integer from file
    """
    data = np.fromfile(fid, dtype=np.uint8, count=3).astype(np.int32)

    out = np.left_shift(data[0], 16) + np.left_shift(data[1], 8) + data[2]

    return out

def partBy2(x):
    x = np.bitwise_and(x, 0x1f00000000ffff)
    x = np.bitwise_and(np.bitwise_or(x, np.left_shift(x, 32)), 0x1f00000000ffff)
    x = np.bitwise_and(np.bitwise_or(x, np.left_shift(x, 16)), 0x1f0000ff0000ff)
    x = np.bitwise_and(np.bitwise_or(x, np.left_shift(x, 8)), 0x100f00f00f00f00f)
    x = np.bitwise_and(np.bitwise_or(x, np.left_shift(x, 4)), 0x10c30c30c30c30c3)
    x = np.bitwise_and(np.bitwise_or(x, np.left_shift(x, 2)), 0x1249249249249249)
    return x

    def test_shift(self):
        from numpy import left_shift, right_shift, dtype

        assert (left_shift([5, 1], [2, 13]) == [20, 2 ** 13]).all()
        assert (right_shift(10, range(5)) == [10, 5, 2, 1, 0]).all()
        bool_ = dtype("bool").type
        assert left_shift(bool(1), 3) == left_shift(1, 3)
        assert right_shift(bool(1), 3) == right_shift(1, 3)

def fliplr(x,length):
	x1 = array(x)
	x[:] = 0
	for i in xrange(length):
		x2 = array(x1)
		x2 = right_shift(x2,i)
		bitwise_and(x2,1,out=x2)
		left_shift(x2,length-1-i,out=x2)
		x += x2

def main():

  parser = argparse.ArgumentParser(description='Ingest the TIFF data')
  parser.add_argument('token', action="store", type=str, help='Token for the project')
  parser.add_argument('channel', action="store", type=str, help='Channel for the project')
  parser.add_argument('path', action="store", type=str, help='Directory with the image files')
  parser.add_argument('resolution', action="store", type=int, help='Resolution of data')

  result = parser.parse_args()
  
  # Load a database
  with closing (ocpcaproj.OCPCAProjectsDB()) as projdb:
    proj = projdb.loadToken(result.token)

  with closing (ocpcadb.OCPCADB(proj)) as db:

    ch = proj.getChannelObj(result.channel)
    # get the dataset configuration
    [[ximagesz, yimagesz, zimagesz],(starttime,endtime)] = proj.datasetcfg.imageSize(result.resolution)
    [xcubedim,ycubedim,zcubedim] = cubedim = proj.datasetcfg.getCubeDims()[result.resolution]
    [xoffset, yoffset, zoffset] = proj.datasetcfg.getOffset()[result.resolution]

    # Get a list of the files in the directories
    for slice_number in range (zoffset, zimagesz+1, zcubedim):
      slab = np.zeros([zcubedim, yimagesz, ximagesz ], dtype=np.uint32)
      for b in range(zcubedim):
        if (slice_number + b <= zimagesz):
          try:
            # reading the raw data
            file_name = "{}BrainBowExM_{:0>4}.tif".format(result.path, slice_number+b)
            print "Open filename {}".format(file_name)
            img = Image.open(file_name, 'r').convert("RGBA")
            imgdata = np.asarray(img)
            slab[b,:,:] = np.left_shift(imgdata[:,:,3], 24, dtype=np.uint32) | np.left_shift(imgdata[:,:,2], 16, dtype=np.uint32) | np.left_shift(imgdata[:,:,1], 8, dtype=np.uint32) | np.uint32(imgdata[:,:,0])
          except IOError, e:
            print e
            imgdata = np.zeros((yimagesz, ximagesz), dtype=np.uint32)
            slab[b,:,:] = imgdata

      for y in range ( 0, yimagesz+1, ycubedim ):
        for x in range ( 0, ximagesz+1, xcubedim ):

          # Getting a Cube id and ingesting the data one cube at a time
          zidx = ocplib.XYZMorton ( [x/xcubedim, y/ycubedim, (slice_number-zoffset)/zcubedim] )
          cube = Cube.getCube(cubedim, ch.getChannelType(), ch.getDataType())
          cube.zeros()

          xmin = x
          ymin = y
          xmax = min ( ximagesz, x+xcubedim )
          ymax = min ( yimagesz, y+ycubedim )
          zmin = 0
          zmax = min(slice_number+zcubedim, zimagesz+1)

          cube.data[0:zmax-zmin,0:ymax-ymin,0:xmax-xmin] = slab[zmin:zmax, ymin:ymax, xmin:xmax]
          db.putCube(ch, zidx, result.resolution, cube, update=True)

 def _pack_bin(self, pixel_list):
     """
     Internal. Encodes [R,G,B] into 16 bit RGB565
     works on a numpy array (H, W, 3) returns flattened bytes string.
     """
     bits16 = np.zeros(pixel_list.shape[:2], dtype=np.uint16)
     bits16 += np.left_shift(np.bitwise_and(pixel_list[:,:,0], 0xF8), 8)
     bits16 += np.left_shift(np.bitwise_and(pixel_list[:,:,1], 0xFC), 3)
     bits16 += np.right_shift(pixel_list[:,:,2], 3)
     return bits16.tostring()

 def func (src, dest):
     mapped = mapper (src)
     dest.fill (0xFF000000)
     effscratch = (mapped[:,:,0] * 0xFF).astype (np.uint32)
     np.left_shift (effscratch, 16, effscratch)
     np.bitwise_or (dest, effscratch, dest)
     effscratch = (mapped[:,:,1] * 0xFF).astype (np.uint32)
     np.left_shift (effscratch, 8, effscratch)
     np.bitwise_or (dest, effscratch, dest)
     effscratch = (mapped[:,:,2] * 0xFF).astype (np.uint32)
     np.bitwise_or (dest, effscratch, dest)

def dec10to16(data):
    arr10 = data.astype(np.uint16).flat
    new_shape = list(data.shape[:-1]) + [(data.shape[-1] * 8) / 10]
    arr16 = np.zeros(new_shape, dtype=np.uint16)
    arr16.flat[::4] = np.left_shift(arr10[::5], 2) + \
        np.right_shift((arr10[1::5]), 6)
    arr16.flat[1::4] = np.left_shift((arr10[1::5] & 63), 4) + \
        np.right_shift((arr10[2::5]), 4)
    arr16.flat[2::4] = np.left_shift(arr10[2::5] & 15, 6) + \
        np.right_shift((arr10[3::5]), 2)
    arr16.flat[3::4] = np.left_shift(arr10[3::5] & 3, 8) + \
        arr10[4::5]
    return arr16