<s> import math <EOL> import numpy as np <EOL> from scipy import optimize <EOL> from . . _shared . utils import skimage_deprecation , warn <EOL> def _check_data_dim ( data , dim ) : <EOL> if data . ndim != <NUM_LIT:2> or data . shape [ <NUM_LIT:1> ] != dim : <EOL> raise ValueError ( '<STR_LIT>' % dim ) <EOL> def _check_data_atleast_2D ( data ) : <EOL> if data . ndim < <NUM_LIT:2> or data . shape [ <NUM_LIT:1> ] < <NUM_LIT:2> : <EOL> raise ValueError ( '<STR_LIT>' ) <EOL> def _norm_along_axis ( x , axis ) : <EOL> """<STR_LIT>""" <EOL> return np . sqrt ( np . einsum ( '<STR_LIT>' , x , x ) ) <EOL> class BaseModel ( object ) : <EOL> def __init__ ( self ) : <EOL> self . params = None <EOL> class LineModel ( BaseModel ) : <EOL> """<STR_LIT>""" <EOL> def __init__ ( self ) : <EOL> self . params = None <EOL> warn ( skimage_deprecation ( '<STR_LIT>' <EOL> '<STR_LIT>' ) ) <EOL> def estimate ( self , data ) : <EOL> """<STR_LIT>""" <EOL> _check_data_dim ( data , dim = <NUM_LIT:2> ) <EOL> X0 = data . mean ( axis = <NUM_LIT:0> ) <EOL> if data . shape [ <NUM_LIT:0> ] == <NUM_LIT:2> : <EOL> theta = np . arctan2 ( data [ <NUM_LIT:1> , <NUM_LIT:1> ] - data [ <NUM_LIT:0> , <NUM_LIT:1> ] , <EOL> data [ <NUM_LIT:1> , <NUM_LIT:0> ] - data [ <NUM_LIT:0> , <NUM_LIT:0> ] ) <EOL> elif data . shape [ <NUM_LIT:0> ] > <NUM_LIT:2> : <EOL> data = data - X0 <EOL> _ , _ , v = np . linalg . svd ( data ) <EOL> theta = np . arctan2 ( v [ <NUM_LIT:0> , <NUM_LIT:1> ] , v [ <NUM_LIT:0> , <NUM_LIT:0> ] ) <EOL> else : <EOL> raise ValueError ( '<STR_LIT>' ) <EOL> theta = ( theta + np . pi / <NUM_LIT:2> ) % np . pi <EOL> dist = X0 [ <NUM_LIT:0> ] * math . cos ( theta ) + X0 [ <NUM_LIT:1> ] * math . sin ( theta ) <EOL> self . params = ( dist , theta ) <EOL> return True <EOL> def residuals ( self , data ) : <EOL> """<STR_LIT>""" <EOL> _check_data_dim ( data , dim = <NUM_LIT:2> ) <EOL> dist , theta = self . params <EOL> x = data [ : , <NUM_LIT:0> ] <EOL> y = data [ : , <NUM_LIT:1> ] <EOL> return dist - ( x * math . cos ( theta ) + y * math . sin ( theta ) ) <EOL> def predict_x ( self , y , params = None ) : <EOL> """<STR_LIT>""" <EOL> if params is None : <EOL> params = self . params <EOL> dist , theta = params <EOL> return ( dist - y * math . sin ( theta ) ) / math . cos ( theta ) <EOL> def predict_y ( self , x , params = None ) : <EOL> """<STR_LIT>""" <EOL> if params is None : <EOL> params = self . params <EOL> dist , theta = params <EOL> return ( dist - x * math . cos ( theta ) ) / math . sin ( theta ) <EOL> class LineModelND ( BaseModel ) : <EOL> """<STR_LIT>""" <EOL> def estimate ( self , data ) : <EOL> """<STR_LIT>""" <EOL> _check_data_atleast_2D ( data ) <EOL> X0 = data . mean ( axis = <NUM_LIT:0> ) <EOL> if data . shape [ <NUM_LIT:0> ] == <NUM_LIT:2> : <EOL> u = data [ <NUM_LIT:1> ] - data [ <NUM_LIT:0> ] <EOL> norm = np . linalg . norm ( u ) <EOL> if norm > <NUM_LIT:0> : <EOL> u /= norm <EOL> elif data . shape [ <NUM_LIT:0> ] > <NUM_LIT:2> : <EOL> data = data - X0 <EOL> _ , _ , u = np . linalg . svd ( data , full_matrices = False ) <EOL> u = u [ <NUM_LIT:0> ] <EOL> else : <EOL> raise ValueError ( '<STR_LIT>' ) <EOL> self . params = ( X0 , u ) <EOL> return True <EOL> def residuals ( self , data ) : <EOL> """<STR_LIT>""" <EOL> X0 , u = self . params <EOL> return _norm_along_axis ( ( data - X0 ) - <EOL> np . dot ( data - X0 , u ) [ ... , np . newaxis ] * u , <EOL> axis = <NUM_LIT:1> ) <EOL> def predict ( self , x , axis = <NUM_LIT:0> , params = None ) : <EOL> """<STR_LIT>""" <EOL> if params is None : <EOL> params = self . params <EOL> X0 , u = params <EOL> if u [ axis ] == <NUM_LIT:0> : <EOL> raise ValueError ( '<STR_LIT>' % axis ) <EOL> l = ( x - X0 [ axis ] ) / u [ axis ] <EOL> return X0 + l [ ... , np . newaxis ] * u <EOL> def predict_x ( self , y , params = None ) : <EOL> """<STR_LIT>""" <EOL> return self . predict ( y , axis = <NUM_LIT:1> , params = params ) [ : , <NUM_LIT:0> ] <EOL> def predict_y ( self , x , params = None ) : <EOL> """<STR_LIT>""" <EOL> return self . predict ( x , axis = <NUM_LIT:0> , params = params ) [ : , <NUM_LIT:1> ] <EOL> class CircleModel ( BaseModel ) : <EOL> """<STR_LIT>""" <EOL> def estimate ( self , data ) : <EOL> """<STR_LIT>""" <EOL> _check_data_dim ( data , dim = <NUM_LIT:2> ) <EOL> x = data [ : , <NUM_LIT:0> ] <EOL> y = data [ : , <NUM_LIT:1> ] <EOL> A = np . zeros ( ( <NUM_LIT:3> , data . shape [ <NUM_LIT:0> ] ) , dtype = np . double ) <EOL> A [ <NUM_LIT:2> , : ] = - <NUM_LIT:1> <EOL> def dist ( xc , yc ) : <EOL> return np . sqrt ( ( x - xc ) ** <NUM_LIT:2> + ( y - yc ) ** <NUM_LIT:2> ) <EOL> def fun ( params ) : <EOL> xc , yc , r = params <EOL> return dist ( xc , yc ) - r <EOL> def Dfun ( params ) : <EOL> xc , yc , r = params <EOL> d = dist ( xc , yc ) <EOL> A [ <NUM_LIT:0> , : ] = - ( x - xc ) / d <EOL> A [ <NUM_LIT:1> , : ] = - ( y - yc ) / d <EOL> return A <EOL> xc0 = x . mean ( ) <EOL> yc0 = y . mean ( ) <EOL> r0 = dist ( xc0 , yc0 ) . mean ( ) <EOL> params0 = ( xc0 , yc0 , r0 ) <EOL> params , _ = optimize . leastsq ( fun , params0 , Dfun = Dfun , col_deriv = True ) <EOL> self . params = params <EOL> return True <EOL> def residuals ( self , data ) : <EOL> """<STR_LIT>""" <EOL> _check_data_dim ( data , dim = <NUM_LIT:2> ) <EOL> xc , yc , r = self . params <EOL> x = data [ : , <NUM_LIT:0> ] <EOL> y = data [ : , <NUM_LIT:1> ] <EOL> return r - np . sqrt ( ( x - xc ) ** <NUM_LIT:2> + ( y - yc ) ** <NUM_LIT:2> ) <EOL> def predict_xy ( self , t , params = None ) : <EOL> """<STR_LIT>""" <EOL> if params is None : <EOL> params = self . params <EOL> xc , yc , r = params <EOL> x = xc + r * np . cos ( t ) <EOL> y = yc + r * np . sin ( t ) <EOL> return np . concatenate ( ( x [ ... , None ] , y [ ... , None ] ) , axis = t . ndim ) <EOL> class EllipseModel ( BaseModel ) : <EOL> """<STR_LIT>""" <EOL> def estimate ( self , data ) : <EOL> """<STR_LIT>""" <EOL> _check_data_dim ( data , dim = <NUM_LIT:2> ) <EOL> x = data [ : , <NUM_LIT:0> ] <EOL> y = data [ : , <NUM_LIT:1> ] <EOL> N = data . shape [ <NUM_LIT:0> ] <EOL> A = np . zeros ( ( N + <NUM_LIT:5> , <NUM_LIT:2> * N ) , dtype = np . double ) <EOL> A [ <NUM_LIT:0> , : N ] = - <NUM_LIT:1> <EOL> A [ <NUM_LIT:1> , N : ] = - <NUM_LIT:1> <EOL> diag_idxs = np . diag_indices ( N ) <EOL> def fun ( params ) : <EOL> xyt = self . predict_xy ( params [ <NUM_LIT:5> : ] , params [ : <NUM_LIT:5> ] ) <EOL> fx = x - xyt [ : , <NUM_LIT:0> ] <EOL> fy = y - xyt [ : , <NUM_LIT:1> ] <EOL> return np . append ( fx , fy ) <EOL> def Dfun ( params ) : <EOL> xc , yc , a , b , theta = params [ : <NUM_LIT:5> ] <EOL> t = params [ <NUM_LIT:5> : ] <EOL> ct = np . cos ( t ) <EOL> st = np . sin ( t ) <EOL> ctheta = math . cos ( theta ) <EOL> stheta = math . sin ( theta ) <EOL> A [ <NUM_LIT:2> , : N ] = - ctheta * ct <EOL> A [ <NUM_LIT:3> , : N ] = stheta * st <EOL> A [ <NUM_LIT:4> , : N ] = a * stheta * ct + b * ctheta * st <EOL> A [ <NUM_LIT:5> : , : N ] [ diag_idxs ] = a * ctheta * st + b * stheta * ct <EOL> A [ <NUM_LIT:2> , N : ] = - stheta * ct <EOL> A [ <NUM_LIT:3> , N : ] = - ctheta * st <EOL> A [ <NUM_LIT:4> , N : ] = - a * ctheta * ct + b * stheta * st <EOL> A [ <NUM_LIT:5> : , N : ] [ diag_idxs ] = a * stheta * st - b * ctheta * ct <EOL> return A <EOL> params0 = np . empty ( ( N + <NUM_LIT:5> , ) , dtype = np . double ) <EOL> xc0 = x . mean ( ) <EOL> yc0 = y . mean ( ) <EOL> r0 = np . sqrt ( ( x - xc0 ) ** <NUM_LIT:2> + ( y - yc0 ) ** <NUM_LIT:2> ) . mean ( ) <EOL> params0 [ : <NUM_LIT:5> ] = ( xc0 , yc0 , r0 , <NUM_LIT:0> , <NUM_LIT:0> ) <EOL> params0 [ <NUM_LIT:5> : ] = np . arctan2 ( y - yc0 , x - xc0 ) <EOL> params , _ = optimize . leastsq ( fun , params0 , Dfun = Dfun , col_deriv = True ) <EOL> self . params = params [ : <NUM_LIT:5> ] <EOL> return True <EOL> def residuals ( self , data ) : <EOL> """<STR_LIT>""" <EOL> _check_data_dim ( data , dim = <NUM_LIT:2> ) <EOL> xc , yc , a , b , theta = self . params <EOL> ctheta = math . cos ( theta ) <EOL> stheta = math . sin ( theta ) <EOL> x = data [ : , <NUM_LIT:0> ] <EOL> y = data [ : , <NUM_LIT:1> ] <EOL> N = data . shape [ <NUM_LIT:0> ] <EOL> def fun ( t , xi , yi ) : <EOL> ct = math . cos ( t ) <EOL> st = math . sin ( t ) <EOL> xt = xc + a * ctheta * ct - b * stheta * st <EOL> yt = yc + a * stheta * ct + b * ctheta * st <EOL> return ( xi - xt ) ** <NUM_LIT:2> + ( yi - yt ) ** <NUM_LIT:2> <EOL> residuals = np . empty ( ( N , ) , dtype = np . double ) <EOL> t0 = np . arctan2 ( y - yc , x - xc ) - theta <EOL> for i in range ( N ) : <EOL> xi = x [ i ] <EOL> yi = y [ i ] <EOL> t , _ = optimize . leastsq ( fun , t0 [ i ] , args = ( xi , yi ) ) <EOL> residuals [ i ] = np . sqrt ( fun ( t , xi , yi ) ) <EOL> return residuals <EOL> def predict_xy ( self , t , params = None ) : <EOL> """<STR_LIT>""" <EOL> if params is None : <EOL> params = self . params <EOL> xc , yc , a , b , theta = params <EOL> ct = np . cos ( t ) <EOL> st = np . sin ( t ) <EOL> ctheta = math . cos ( theta ) <EOL> stheta = math . sin ( theta ) <EOL> x = xc + a * ctheta * ct - b * stheta * st <EOL> y = yc + a * stheta * ct + b * ctheta * st <EOL> return np . concatenate ( ( x [ ... , None ] , y [ ... , None ] ) , axis = t . ndim ) <EOL> def _dynamic_max_trials ( n_inliers , n_samples , min_samples , probability ) : <EOL> """<STR_LIT>""" <EOL> if n_inliers == <NUM_LIT:0> : <EOL> return np . inf <EOL> nom = <NUM_LIT:1> - probability <EOL> if nom == <NUM_LIT:0> : <EOL> return np . inf <EOL> inlier_ratio = n_inliers / float ( n_samples ) <EOL> denom = <NUM_LIT:1> - inlier_ratio ** min_samples <EOL> if denom == <NUM_LIT:0> : <EOL> return <NUM_LIT:1> <EOL> elif denom == <NUM_LIT:1> : <EOL> return np . inf <EOL> nom = np . log ( nom ) <EOL> denom = np . log ( denom ) <EOL> if denom == <NUM_LIT:0> : <EOL> return <NUM_LIT:0> <EOL> return int ( np . ceil ( nom / denom ) ) <EOL> def ransac ( data , model_class , min_samples , residual_threshold , <EOL> is_data_valid = None , is_model_valid = None , <EOL> max_trials = <NUM_LIT:100> , stop_sample_num = np . inf , stop_residuals_sum = <NUM_LIT:0> , <EOL> stop_probability = <NUM_LIT:1> ) : <EOL> """<STR_LIT>""" <EOL> best_model = None <EOL> best_inlier_num = <NUM_LIT:0> <EOL> best_inlier_residuals_sum = np . inf <EOL> best_inliers = None <EOL> if min_samples < <NUM_LIT:0> : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> if max_trials < <NUM_LIT:0> : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> if stop_probability < <NUM_LIT:0> or stop_probability > <NUM_LIT:1> : <EOL> raise ValueError ( "<STR_LIT>" ) <EOL> if not isinstance ( data , list ) and not isinstance ( data , tuple ) : <EOL> data = [ data ] <EOL> data = list ( data ) <EOL> num_samples = data [ <NUM_LIT:0> ] . shape [ <NUM_LIT:0> ] <EOL> for num_trials in range ( max_trials ) : <EOL> samples = [ ] <EOL> random_idxs = np . random . randint ( <NUM_LIT:0> , num_samples , min_samples ) <EOL> for d in data : <EOL> samples . append ( d [ random_idxs ] ) <EOL> if is_data_valid is not None and not is_data_valid ( * samples ) : <EOL> continue <EOL> sample_model = model_class ( ) <EOL> success = sample_model . estimate ( * samples ) <EOL> if success is not None : <EOL> if not success : <EOL> continue <EOL> if is_model_valid is not None and not is_model_valid ( sample_model , <EOL> * samples ) : <EOL> continue <EOL> sample_model_residuals = np . abs ( sample_model . residuals ( * data ) ) <EOL> sample_model_inliers = sample_model_residuals < residual_threshold <EOL> sample_model_residuals_sum = np . sum ( sample_model_residuals ** <NUM_LIT:2> ) <EOL> sample_inlier_num = np . sum ( sample_model_inliers ) <EOL> if ( <EOL> sample_inlier_num > best_inlier_num <EOL> or ( sample_inlier_num == best_inlier_num <EOL> and sample_model_residuals_sum < best_inlier_residuals_sum ) <EOL> ) : <EOL> best_model = sample_model <EOL> best_inlier_num = sample_inlier_num <EOL> best_inlier_residuals_sum = sample_model_residuals_sum <EOL> best_inliers = sample_model_inliers <EOL> if ( <EOL> best_inlier_num >= stop_sample_num <EOL> or best_inlier_residuals_sum <= stop_residuals_sum <EOL> or num_trials <EOL> >= _dynamic_max_trials ( best_inlier_num , num_samples , <EOL> min_samples , stop_probability ) <EOL> ) : <EOL> break <EOL> if best_inliers is not None : <EOL> for i in range ( len ( data ) ) : <EOL> data [ i ] = data [ i ] [ best_inliers ] <EOL> best_model . estimate ( * data ) <EOL> return best_model , best_inliers </s>
