In [7]:
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
import cv2
from time import time
import skimage
import skimage.transform
from skimage import measure
from skimage import morphology
from skimage.filters import threshold_otsu
IMAGE_FILE = '/home/jason/Downloads/train/trichodesmium_bowtie/104373.jpg'
# IMAGE_FILE = '/home/jason/Downloads/train/tunicate_salp_chains/111699.jpg'
# IMAGE_FILE = '/home/jason/Downloads/train/protist_star/101029.jpg'
# IMAGE_FILE = '/home/jason/Downloads/train/artifacts_edge/106668.jpg'
im = 255 - cv2.imread(IMAGE_FILE, cv2.IMREAD_GRAYSCALE)
OUT_SHAPE = (64, 64)
imshow = lambda im: plt.imshow(im, cmap='gray', interpolation='none')
In [8]:
def get_largest_region(im, show_plots=False):
# find the largest nonzero region
imthr = np.array(im)
thresh_fn = lambda x: x > np.mean(x)
# thresh_fn = lambda x: x > threshold_otsu(x)*.9
imthr = np.where(thresh_fn(im),1.,0.)
imdilated = morphology.dilation(imthr, np.ones((4,4)))
labels = measure.label(imdilated)
labels = imthr * labels
labels = labels.astype(int)
if show_plots:
f = plt.figure(figsize=(12,3))
sub1 = plt.subplot(1,4,1)
imshow(im)
sub2 = plt.subplot(1,4,2)
imshow(imthr)
sub3 = plt.subplot(1, 4, 3)
imshow(imdilated)
sub4 = plt.subplot(1, 4, 4)
plt.imshow(labels)
sub1.set_title("Original Image")
sub2.set_title("Thresholded Image")
sub3.set_title("Dilated Image")
sub4.set_title("Labeled Image")
regions = measure.regionprops(labels)
regionmaxprop = None
for regionprop in regions:
# check to see if the region is at least 50% nonzero
if sum(imthr[labels == regionprop.label])*1.0/regionprop.area < 0.50:
continue
if regionmaxprop is None:
regionmaxprop = regionprop
if regionmaxprop.filled_area < regionprop.filled_area:
regionmaxprop = regionprop
return regionmaxprop, labels
In [9]:
tic = time()
regionmax, labels = get_largest_region(im)
bw = np.where(labels == regionmax.label,True,False)
print time() - tic
plt.figure()
imshow(bw)
Out[9]:
In [10]:
OUT_SHAPE = (64, 64)
center_in = np.array(regionmax.bbox).reshape((2, 2)).mean(axis=0)[::-1]
center_out = np.array(OUT_SHAPE) / 2. - 0.5
tform_center = skimage.transform.SimilarityTransform(translation=-center_out)
tform_uncenter = skimage.transform.SimilarityTransform(translation=center_in)
def build_tform(tform_center, tform_uncenter,
zoom=1.0, rotation=0, shear=0, translation=(0, 0)):
tform_augment = skimage.transform.AffineTransform(scale=(1/zoom, 1/zoom), rotation=np.deg2rad(rotation), shear=np.deg2rad(shear), translation=translation)
tform = tform_center + tform_augment + tform_uncenter # shift to center, augment, shift back (for the rotation/shearing)
return tform
def fast_warp(img, tf, output_shape=OUT_SHAPE, mode='constant'):
"""
This wrapper function is about five times faster than skimage.transform.warp, for our use case.
"""
# m = tf._matrix
m = tf.params
img_wf = skimage.transform._warps_cy._warp_fast(img, m, output_shape=output_shape, mode=mode)
return img_wf
In [11]:
theta = regionmax.orientation
print 'Rotate:', np.rad2deg(theta)
max_side = np.diff(np.array(regionmax.bbox).reshape((2, 2)), axis=0).max()
max_side
alpha = OUT_SHAPE[0] / float(max_side)
print 'Zoom:', alpha
tform = build_tform(tform_center, tform_uncenter,
zoom=alpha, rotation=np.rad2deg(theta))
im_w = fast_warp(img=im, tf=tform)
imshow(im_w)
Out[11]:
In [15]:
def get_features(im_file, out_shape=OUT_SHAPE, verbose=False, show_plots=False):
im = 255 - cv2.imread(im_file, cv2.IMREAD_GRAYSCALE) # subject is white, bckg is black
regionmax, labels = get_largest_region(im, show_plots=show_plots)
h, w = sorted(im.shape)
f_size = os.stat(im_file).st_size
if regionmax:
ext = regionmax.extent
hu = regionmax.moments_hu
sol = regionmax.solidity
bw = np.where(labels == regionmax.label,True,False)
if show_plots:
plt.figure()
imshow(bw)
center_in = np.array(regionmax.bbox).reshape((2, 2)).mean(axis=0)[::-1]
center_out = np.array(out_shape) / 2. - 0.5
tf_cent = skimage.transform.SimilarityTransform(translation=-center_out)
tf_uncent = skimage.transform.SimilarityTransform(translation=center_in)
theta = regionmax.orientation
max_side = np.diff(np.array(regionmax.bbox).reshape((2, 2)), axis=0).max()
alpha = out_shape[0] / float(max_side)
if verbose:
print 'Rotate:', np.rad2deg(theta)
print 'Zoom:', alpha
tform = build_tform(tform_center=tf_cent, tform_uncenter=tf_uncent,
zoom=alpha, rotation=np.rad2deg(theta))
im_w = fast_warp(img=im, tf=tform)
if show_plots:
plt.figure()
imshow(im_w)
else:
# Could not find a valid region
# Just scale the image to the appropriate size
# and return empty features
im_w = skimage.transform.resize(im, out_shape, mode='constant')
ext, hu, sol = 0, np.zeros(7), 0
return (im_w, (f_size, h, w, ext, hu, sol))
In [17]:
im_norm, feats = get_features(IMAGE_FILE, out_shape=OUT_SHAPE, verbose=False, show_plots=False)
imshow(im_norm)
print feats
In [193]:
imshow(q)
Out[193]:
In [96]:
import plyvel
import caffe_pb2
import os
from le import le
import taxonomy as tax
from sklearn.preprocessing import StandardScaler
from time import time
import glob
import os
TRAIN_PATH = '/home/jason/Downloads/train/'
DB_PATH = '/home/jason/documents/sandbox/orientation_normed_lvl/'
OUT_SHAPE = (64, 64)
In [59]:
datum = caffe_pb2.Datum()
db = plyvel.DB(DB_PATH, create_if_missing=True)
im_files = glob.glob(os.path.join(TRAIN_PATH, '*', '*.jpg'))
invalids = []
feats_all = []
tic = time()
for ii, im_file in enumerate(im_files):
if ii % 5000 == 0:
print '%d images processed in %d s' % (ii, time() - tic)
y_str = os.path.basename(os.path.dirname(im_file))
y = le.transform(y_str)
im_w, (f_size, h, w, ext, hu, sol) = get_features(im_file, out_shape=OUT_SHAPE)
if not sol:
invalids.append((y_str, im_file))
feats_all.append(np.r_[f_size, h, w, ext, hu, sol]) # Used to make std scaler
# Fill in Datum
datum.channels = 1
datum.height = OUT_SHAPE[0]
datum.width = OUT_SHAPE[1]
datum.data = im_w.astype('uint8').tobytes()
datum.label = y
#float_data = 6;
#encoded = 7 [default = false];
datum.label0 = tax.encode_parent(y_str, 0)
datum.label1 = tax.encode_parent(y_str, 1)
datum.label2 = tax.encode_parent(y_str, 2)
datum.label3 = tax.encode_parent(y_str, 3)
datum.label4 = tax.encode_parent(y_str, 4)
datum.label5 = tax.encode_parent(y_str, 5)
datum.label6 = tax.encode_parent(y_str, 6)
# Normalized Features (Normalize on second pass)
datum.orig_space = f_size
datum.orig_height = h
datum.orig_width = w
datum.extent = ext
datum.hu1 = hu[0]
datum.hu2 = hu[1]
datum.hu3 = hu[2]
datum.hu4 = hu[3]
datum.hu5 = hu[4]
datum.hu6 = hu[5]
datum.hu7 = hu[6]
datum.solidity = sol
k = os.path.basename(im_file)
v = datum.SerializeToString()
db.put(k, v)
db.close()
In [77]:
feats_arr = np.array(feats_all)
scaler = StandardScaler().fit(feats_arr)
In [82]:
# Scaling features
# yeah this could be less ugly, but whatever
db = plyvel.DB(DB_PATH)
wb = db.write_batch()
for key, value in db:
datum.ParseFromString(value)
extra_feats = np.array([
datum.orig_space,
datum.orig_height,
datum.orig_width,
datum.extent,
datum.hu1,
datum.hu2,
datum.hu3,
datum.hu4,
datum.hu5,
datum.hu6,
datum.hu7,
datum.solidity,
])
scaled_feats = scaler.transform(extra_feats)
datum.orig_space = scaled_feats[0]
datum.orig_height = scaled_feats[1]
datum.orig_width = scaled_feats[2]
datum.extent = scaled_feats[3]
datum.hu1 = scaled_feats[4]
datum.hu2 = scaled_feats[5]
datum.hu3 = scaled_feats[6]
datum.hu4 = scaled_feats[7]
datum.hu5 = scaled_feats[8]
datum.hu6 = scaled_feats[9]
datum.hu7 = scaled_feats[10]
datum.solidity = scaled_feats[11]
v = datum.SerializeToString()
wb.put(key, v)
wb.write()
db.close()
In [111]:
ims = []
db = plyvel.DB(DB_PATH)
for key, value in db:
datum.ParseFromString(value)
im = np.array(Image.frombytes('L', OUT_SHAPE, datum.data))
ims.append(im)
In [101]:
from PIL import Image
datum.ParseFromString(v)
im = np.array(Image.frombytes('L', OUT_SHAPE, datum.data))
In [112]:
ims_arr = np.array(ims)
In [115]:
m = ims_arr.mean(axis=0)
In [117]:
m.mean()
Out[117]:
In [48]:
import preproc as pp
reload(pp)
from time import time
import cv2
import cPickle
import shelve
from contextlib import closing
import glob
import os
TRAIN_PATH = '/media/raid_arr/data/ndsb/train/'
train_files = glob.glob(os.path.join(TRAIN_PATH, '*', '*.jpg'))
IM_FILE = '/media/raid_arr/data/ndsb/train/acantharia_protist/100224.jpg'
OUT_SHAPE = (64, 64)
In [50]:
feats_d = {}
tic = time()
for ii, im_file in enumerate(train_files):
toc = time() - tic
if ii % 5000 == 0:
print ii, toc
im_name = os.path.basename(im_file)
im = 255 - cv2.imread(im_file, cv2.IMREAD_GRAYSCALE)
rp, l = pp.get_largest_region(im)
feats_d[im_name] = pp.get_rp_features(rp)
In [53]:
cPickle.dump(feats_d, open('/media/raid_arr/data/ndsb/region_feats.p', 'wb'), protocol=2)
In [54]:
cPickle.load(open('/media/raid_arr/data/ndsb/region_feats.q', 'rb'), protocol=2)
In [95]:
import preproc as pp
reload(pp)
from time import time
import glob
import os
import cv2
import matplotlib.pyplot as plt
%matplotlib inline
IMAGE_FILE = '/home/jason/Downloads/train/trichodesmium_bowtie/104373.jpg'
# IMAGE_FILE = '/home/jason/Downloads/train/tunicate_salp_chains/111699.jpg'
# IMAGE_FILE = '/home/jason/Downloads/train/protist_star/101029.jpg'
# IMAGE_FILE = '/home/jason/Downloads/train/artifacts_edge/106668.jpg'
im = 255 - cv2.imread(IMAGE_FILE, cv2.IMREAD_GRAYSCALE)
OUT_SHAPE = (64, 64)
imshow = lambda im: plt.imshow(im, cmap='gray', interpolation='none')
In [162]:
im_n, feat = pp.get_features(IMAGE_FILE, OUT_SHAPE, verbose=True,
zoom_musig=(1, 0.02),
rotation_musig=(0, 10),
shear_musig=(0, 10),
translation_musig=(0, 2))
imshow(im_n)
Out[162]:
In [ ]: