#Load data
import numpy as n
import logging
logging.basicConfig(level=logging.DEBUG,format='%(levelname)s %(message)s')
import variables,file_io
from xfab import tools,detector
from ImageD11 import blobcorrector
from pylab import *


A_id = variables.refarray().A_id


class uvmap:
    def __init__(self,param=None,grain=None,reconst=None):
        self.param = param
        self.grain = grain
        self.K = -2*n.pi/self.param['wavelength']
        self.S = n.array([[1, 0, 0],[0, 1, 0],[0, 0, 1]])
        self.reconst = reconst   # reconst is a dictinary of grain numbers with a list 
                                 # reflections for which to make uv_maps
        self.R = tools.detect_tilt(self.param['tilt_x'],
                                   self.param['tilt_y'],
                                   self.param['tilt_z'])
        self.setno = 0  # so far always 0
        if self.param['spatial'] != None:
            self.spatial = blobcorrector.correctorclass(self.param['spatial'])
        else:
            self.spatial = None

        #print self.grain[0].refs
        self.horspacing = 1
        self.vertspacing = 1

    def makemap(self):
        print 'Calculating and saving uvmaps. This might take a while, please be patient...'

        z = n.array([0,0,1],dtype=n.float)
        cnt=0

        
        #Loop over all grains
        for grainno in self.reconst.keys():
            
            # Open .hst file of the grain 
            hst = file_io.readhst(self.param['direc'],self.param['stem'])
            hst._open(grainno = grainno, setno = self.setno)
            hst.get_file_pos()

            # Open .uvm file for outputting of uv-maps
            filename = '%s/%s_gr%0.4d_set%0.4d.uvm' \
                %(self.param['direc'],self.param['stem'],grainno,self.setno)
            uvmfile = open(filename,'w')

            # write types to file
            self.param['u_range']= 2*self.param['u_max']+1
            self.param['v_range']= 2*self.param['v_max']+1
            param_out = ['u_range','v_range']
            write_param = '# '
            for p in param_out:
                write_param = write_param + '%s ' %(p)
            uvmfile.write(write_param+'\n')
            write_param = ''
            for p in param_out:
                write_param = write_param + '%i ' %(self.param[p])
            uvmfile.write(write_param+'\n')


            gr_pos = n.array(self.param['pos_grains'][grainno])

            for reflno in self.reconst[grainno]:
                cnt = cnt+1
                print "################################################"
                print "################# REFL %4i    #################" %reflno
                print "################################################"
                print self.grain[grainno].refs[0,1]
                omega_center = self.grain[grainno].refs[reflno,A_id['omega']]*n.pi/180.0
#                self.grain[grainno].refs[reflno,A_id['ref_id']]
                #[refl_id hkl dety detz detyc detzc twotheta eta omega Gw AOI framectr border] = loadref(output_dir, output_fileprefix, grainno, reflno, ceil(reflno/max_refl_in_set));

                # REMEMBER THAT frame_end start can be negative if part of box is outside omega_range
                if self.grain[grainno].refs[reflno,A_id['frame_start']] < 0 or \
                        self.grain[grainno].refs[reflno,A_id['frame_end']] < 0:
                    continue
                                                                                                           
                                               
                frames = n.arange(self.grain[grainno].refs[reflno,A_id['frame_start']],
                                  self.grain[grainno].refs[reflno,A_id['frame_end']]+1)
                #print 'frames ', frames
                omega_limits = n.array([self.param['frame_omega'][frames[0]],
                                self.param['frame_omega'][frames[-1]]+self.param['omega_step']])
                #print 'omega_limits', omega_limits
                omega_limits = omega_limits*n.pi/180.

                hst.read(grainno = grainno, setno = self.setno,file_pos=hst.filepos[reflno])
                #self.hst.read(self.grainno,self.setno,self.hst.filepos[reflno])
                map = hst.shoebox.astype(n.int16)-self.param['darkoffset']
                #map = hst.shoebox.astype(n.uint16)
                
                # Determine eta range over central pixel 
                dety = self.grain[grainno].refs[reflno,A_id['dety']]
                detz = self.grain[grainno].refs[reflno,A_id['detz']]
                hkl = [self.grain[grainno].refs[reflno,A_id['h']],
                       self.grain[grainno].refs[reflno,A_id['k']],
                       self.grain[grainno].refs[reflno,A_id['l']]]

                # Determine Area-of-Interst 
                sbox_yh = (self.param['sbox_y']-1)/2
                sbox_zh = (self.param['sbox_z']-1)/2

                aoi_y = n.int(round(self.grain[grainno].refs[reflno,A_id['detyd']]))+\
                    n.array([-sbox_yh,sbox_yh+1])
                aoi_z = n.int(round(self.grain[grainno].refs[reflno,A_id['detzd']]))+\
                    n.array([-sbox_zh,sbox_zh+1])




                corner_vector = n.zeros((4,2))
                corner_vector[0,:] = [(n.ceil(dety-self.param['dety_center']))*self.param['y_size'], 
                                      (n.ceil(detz-self.param['detz_center']))*self.param['z_size']]
                corner_vector[1,:] = [(n.floor(dety-self.param['dety_center']))*self.param['y_size'],
                                      (n.floor(detz-self.param['detz_center']))*self.param['z_size']]
                corner_vector[2,:] = [(n.ceil(dety-self.param['dety_center']))*self.param['y_size'],
                                      (n.floor(detz-self.param['detz_center']))*self.param['z_size']]
                corner_vector[3,:] = [(n.floor(dety-self.param['dety_center']))*self.param['y_size'],
                                      (n.ceil(detz-self.param['detz_center']))*self.param['z_size']]
                #print corner_vector
                eta_pix = n.zeros((4,4))
                for i in range(4):
                    for j in range(i+1,4):
                        eta_pix[i,j] = n.arccos(n.dot(corner_vector[i],corner_vector[j])/\
                                               (n.linalg.norm(corner_vector[i])*n.linalg.norm(corner_vector[j])))
                #print eta_pix
                eta_pix= eta_pix.max()    #eta_pix is in radians
                #print eta_pix
                eta_step = eta_pix/2.
                #print eta_step

                origo = n.array([-aoi_y[0], -aoi_z[0]]) + 0.5 # adding 0.5 in order to
                                                              # make the corner of the AOI (0,0)
                u_max = self.param['u_max']
                v_max = self.param['v_max']
                u_scale = self.param['u_scale']
                v_scale = self.param['v_scale']
                uvm = n.zeros((u_max*2+1,v_max*2+1))
                etamap = uvm.copy()
                ommap = uvm.copy()
#                 sbmap = uvm.copy()

                B = tools.FormB(self.param['unit_cell'])
                U = self.grain[grainno].U
                Gc0 = n.dot(B,hkl)
                nGc0 = n.linalg.norm(Gc0)
                y0 = Gc0/nGc0
                u = n.cross(y0, z)
                u = u/n.linalg.norm(u)
                v = n.cross(u, y0)
                u_count = -1
                for iu in range(-u_max,u_max+1):
                     u_point = u_scale*iu
                     u_count = u_count +1
                     v_count = -1
                     for iv in range(-v_max,v_max+1):
                         v_point = v_scale*iv
                         v_count = v_count +1
        
                         dGc = nGc0*(u_point*u + v_point*v )
                         Gc = Gc0 + dGc
                         Gw = n.dot(U,Gc)
                         Glen = n.linalg.norm(Gw)
                         tth = 2*n.arcsin(Glen/(2*n.abs(self.K)))
                         (Omega,eta) = tools.find_omega_wedge(Gw,tth,self.param['wedge'])
                         if len(Omega) > 0:
                             solution = n.argmin(n.abs(n.array(Omega)-omega_center))
                             eta = eta[solution]
                             omega = Omega[solution]
                             Om = tools.OMEGA(omega)
                             Gt = n.dot(Om,Gw)
                             #print eta*180.0/n.pi
                         else:
                             print 'No solution found for Omega'
                             continue
                         etamap[u_count,v_count] = eta*180./n.pi
                         ommap[u_count,v_count] = omega*180./n.pi
                         
                         #NEED TO FIND OUT WHAT THIS IS?!
#                         FrameInShoeBox = ceil((my_omega-wmin)/wstep)-frames(1)+1
                         #print "WHY the HELL AM I in HERE", omega,Omega*180./n.pi,(n.array(Omega)-omega_center*n.pi/180.0),(omega*n.pi/180.0)
                         
                         FrameInShoeBox = n.floor((omega-omega_limits[0])/(n.pi/180.0*self.param['omega_step']))
 
                         

                         # If uv point refers to region outside shoebox, set integral to 0 and exit
                         if FrameInShoeBox < 0 or FrameInShoeBox > self.param['sbox_omega']-1:

                             bin_int = 0
                             uvm[u_count,v_count] = bin_int
                             continue
                         #print 'DO I EVER COME HERE!!!!!!!!!!!!!!!!!!!!!',FrameInShoeBox

                         [tx,ty]= n.dot(Om[:2,:2],gr_pos[:2])
                         tz = gr_pos[2]

                         tth_low  = tth - self.param['tth_range']*n.pi/180.
                         tth_high = tth + self.param['tth_range']*n.pi/180.   
                         eta_low  = eta - eta_pix
                         eta_high = eta + eta_pix

                         r1 = n.zeros((2,2))
                         r2 = n.zeros((2,2))
                         r1[0] =  detector.det_coor2(tth_low, eta_low,
                                                     self.param['distance'],
                                                     self.param['y_size'], self.param['z_size'],
                                                     self.param['dety_center'], self.param['detz_center'],
                                                     self.R,
                                                     tx, ty, tz,)
                         r1[1] =  detector.det_coor2(tth_high, eta_low,
                                                     self.param['distance'],
                                                     self.param['y_size'], self.param['z_size'],
                                                     self.param['dety_center'], self.param['detz_center'],
                                                     self.R, 
                                                     tx, ty, tz,)
                         r2[0] =  detector.det_coor2(tth_low, eta_high,
                                                     self.param['distance'],
                                                     self.param['y_size'],
                                                     self.param['z_size'],
                                                     self.param['dety_center'],
                                                     self.param['detz_center'],
                                                     self.R, 
                                                     tx, ty, tz,)
                         r2[1] =  detector.det_coor2(tth_high, eta_high,
                                                     self.param['distance'],
                                                     self.param['y_size'],
                                                     self.param['z_size'],
                                                     self.param['dety_center'],
                                                     self.param['detz_center'],
                                                     self.R, 
                                                     tx, ty, tz,)
                         #print r1
                         #print r2

                         # Correct dety, detz according to spatial file
                         if self.spatial != None:
                             r1[0] = detector.distort(r1[0],
                                                      self.param['o11'],
                                                      self.param['o12'],
                                                      self.param['o21'],
                                                      self.param['o22'],
                                                      self.param['dety_size'],
                                                      self.param['detz_size'],
                                                      self.spatial)
                             r1[1] = detector.distort(r1[1],
                                                      self.param['o11'],
                                                      self.param['o12'],
                                                      self.param['o21'],
                                                      self.param['o22'],
                                                      self.param['dety_size'],
                                                      self.param['detz_size'],
                                                      self.spatial)
                             r2[0] = detector.distort(r2[0],
                                                      self.param['o11'],
                                                      self.param['o12'],
                                                      self.param['o21'],
                                                      self.param['o22'],
                                                      self.param['dety_size'],
                                                      self.param['detz_size'],
                                                      self.spatial)
                             r1[1] = detector.distort(r1[1],
                                                      self.param['o11'],
                                                      self.param['o12'],
                                                      self.param['o21'],
                                                      self.param['o22'],
                                                      self.param['dety_size'],
                                                      self.param['detz_size'],
                                                      self.spatial)

                         #print 'u_point,v_point',u_point,v_point
                         #print 'dGc', dGc
                         #print 'Gc', Gc
                         #print r1,'\n',r2

                         r1[0] = r1[0] + origo
                         r1[1] = r1[1] + origo
                         r2[0] = r2[0] + origo
                         r2[1] = r2[1] + origo
                         #print r1
                         #print r2

                         sizeoflayer = [self.param['sbox_y'],self.param['sbox_z']]
                         propmatrix = radial_int(r1, r2, sizeoflayer, self.horspacing, self.vertspacing)
                         #print propmatrix 
                         #print map[FrameInShoeBox]
                         finalmatrix = propmatrix*map[FrameInShoeBox]
                         bin_int = finalmatrix.sum()
                         uvm[u_count,v_count] = bin_int
                
                print 'Done.'
                uvm = uvm/uvm.sum()
                self.write_uvm(uvmfile,uvm,grainno,reflno)
            uvmfile.close()
        
        return (etamap,ommap,uvm)
                            
#                 finalmatrix = propmatrix.*double(map(:,:,FrameInShoeBox))
#                 bin_int = sum(sum(finalmatrix))
                        
#                         #end
#                         sbmap(u_count,v_count) = FrameInShoeBox
#                     #end
#                 #end
   # uvm = uvm/sum(sum(uvm)) #Normalisation 


    def write_uvm(self, file = None, peak = None, grainno = None, reflno = None):
    #write shoebox as one string of number separated by a space
    #no of elements is sbox_omega*sbox_y*sbox_z
    #fastest zdet_relative, 
    #medium  ydet_relative,
    #slow    omega layer
    #To read shoebox 
    #1d_shoebox =n.fromstring(file.readline,sep=' ',dtype=dtype)
    #shoebox = 1d_shoebox.reshape(sbox_omega,sbox_y,sbox_z)
        file.write('REFL_ID = %i\n' %self.grain[grainno].refs[reflno,A_id['ref_id']])
        file.write('SPOT_ID = %i\n' %self.grain[grainno].refs[reflno,A_id['spot_id']])
        peak.tofile(file,sep=' ')
        file.write(' ')
        file.write('\n')

def radial_int(r1, r2, sizeoflayer, horspacing, vertspacing):
    #print '>>>>>>>>>>>>>>>>>>>> IN RADIAL_INT ROUTINE >>>>>>>>>>>>>>>>>>'

    reflecta = 1
    reflectb = 1

    #print 'r1 - orig', r1
    #print 'r2 - orig', r2
    # TRY TO REORDER POINTs
    #print 'ORDER POINTS'
    points = n.vstack((r1,r2))
    order_points = points[points[:,0].argsort()]
    r1 = order_points[:2]
    r2 = order_points[2:]



    #print 'sort r1', r1
    #print 'sort r2',r2

    # We always want translations in the positive direction horizontally, take
    # this into account when creating the lines... I.e. if the translation on
    # the horizontal axis isn't in the positive direction then rearrange the
    # points such that it will be.

    if r1[1,0]-r1[0,0] < 0 or (r1[1,0]-r1[0,0] == 0 and r1[1,1]-r1[0,1] < 0):
        r1 = n.flipud(r1)
        #print 'RADIAL_INT: FLIPS r1 UD'
    if r2[1,0]-r2[0,0] < 0  or (r2[1,0]-r2[0,0] == 0 and r2[1,1]-r2[0,1] < 0):
        #print 'RADIAL_INT: FLIPS r2 UD'
        r2 = n.flipud(r2)


    if r1[1,1]-r1[0,1] < 0 or (r1[1,0]-r1[0,0] == 0 and r1[1,1]-r1[0,1] < 0):
        r1 = n.flipud(r1)
        #print 'RADIAL_INT: FLIPS r1 UD'
    if r2[1,1]-r2[0,1] < 0  or (r2[1,0]-r2[0,0] == 0 and r2[1,1]-r2[0,1] < 0):
        #print 'RADIAL_INT: FLIPS r2 UD'
        r2 = n.flipud(r2)

    #Just to make sure we have the original values for lines 3 and 4...

    p1 = r1.copy()
    p2 = r2.copy()

    #print 'r1',r1
    #print 'r2',r2

    # Make primary line equations into the form y = mz + c.
    # If gradient would be infinity (at eta = 0, 180 or 360), set it to a very large
    # number (shouldn't happen when reflections are filtered close to the
    # spindle axis!!)

    diff1 = r1[1,0]-r1[0,0]
    diff2 = r2[1,0]-r2[0,0]

    if diff1==0:
        m1 = n.inf
        c1 = r1[1,0]
    else:
        m1 = (r1[1,1]-r1[0,1])/diff1
        c1 = r1[1,1]-m1*r1[1,0]
    if diff2 == 0:
        m2 = n.inf
        c2 =  r2[1,0]
    else:
        m2 = (r2[1,1]-r2[0,1])/diff2
        c2 = r2[1,1]-m2*r2[1,0]

    #print 'r1',r1
    #print 'r2',r2

    #print 'm1, c1',m1,c1
    #print 'm2, c2',m2,c2
    
    x = n.mean(r2[:,0])
    #print 'x,',x
    yline1 = m1*x+c1
    yline2 = m2*x+c2
    #print 'yline1', yline1
    #print 'yline2', yline2
    # Make sure primary line 1 is under line 2
    #if r1[0,1] > r2[0,1] and n.sign(m1) == n.sign(m2)  or \
    if  yline1 > yline2 and n.sign(m1) == n.sign(m2)  or \
            r1[0,0] > r2[0,0] and m1 == n.inf and m2 == n.inf  or \
            n.sign(m1) == -1 and n.sign(m2) == 1 and r1[1,0] < r2[0,0] or \
            n.sign(m1) ==  1 and n.sign(m2) ==-1 and r1[0,0] > r2[1,0]:
        #print '>>>>>>>>>> interchanging lines', r1[0,1], r2[0,1] 

        tmp = m1
        m1 = m2
        m2 = tmp
        tmp = c1
        c1 = c2
        c2 = tmp
        tmp = r1.copy()
        r1 = r2.copy()
        r2 = tmp.copy()


    # Initialize propmatrices
    propmatrix1 = n.zeros(sizeoflayer)
    propmatrix2 = n.zeros(sizeoflayer)

    # Find rectangle of interest - will set everything outside this rectangle to
    # zero in the propmatrix because those areas fall outside the area of
    # integration
    # This rectangle will be used for lines 3 and 4 as well

    #print 'r1',r1
    #print 'r2',r2
    
    rectangle_lowpoint = n.array([n.floor(n.min([n.min(r1[:,0]), n.min(r2[:,0])])),
                                  n.floor(n.min([n.min(r1[:,1]), n.min(r2[:,1])]))])
    rectangle_highpoint = n.array([n.ceil(n.max([n.max(r1[:,0]), n.max(r2[:,0])])),
                                   n.ceil(n.max([n.max(r1[:,1]), n.max(r2[:,1])]))])

    
    #print 'rectangle_lowpoint ',    rectangle_lowpoint
    #print 'rectangle_higpoint ',    rectangle_highpoint


    # If parts of the rectangle lie outside the shoebox, then rectangle needs to
    # be resized to fit the shoebox.

    rectangle_lowpoint = n.clip(rectangle_lowpoint,0,n.inf)
    rectangle_highpoint = n.clip(rectangle_highpoint,0,n.inf)

    if rectangle_lowpoint[0] >= sizeoflayer[0]:
        rectangle_lowpoint[0] = sizeoflayer[0]
    if rectangle_lowpoint[1] >= sizeoflayer[1]:
        rectangle_lowpoint[1] = sizeoflayer[1]
    if rectangle_highpoint[0] >= sizeoflayer[0]:
        rectangle_highpoint[0] = sizeoflayer[0]
    if rectangle_highpoint[1] >= sizeoflayer[1]:
        rectangle_highpoint[1] = sizeoflayer[1]

    #print 'rectangle_lowpoint ',    rectangle_lowpoint
    #print 'rectangle_higpoint ',    rectangle_highpoint

    # If line is vertical rotate AOI for integration and rotate back afterwards
    if m1 == n.inf:
        m1 = 0.0
        propmatrix1 = integrate(c1, m1,
                                [rectangle_lowpoint[1],rectangle_lowpoint[0]],
                                [rectangle_highpoint[1],rectangle_highpoint[0]],
                                n.transpose(propmatrix1),
                                vertspacing,
                                horspacing)
        propmatrix1 = n.transpose(propmatrix1)
    else:        
        propmatrix1 = integrate(c1, m1,
                                rectangle_lowpoint, 
                                rectangle_highpoint, 
                                propmatrix1, 
                                horspacing, 
                                vertspacing)

    # If line is vertical rotate AOI for integration and rotate back afterwards
    if m2 == n.inf:
        m2 = 0.0
        propmatrix2 = integrate(c2, m2, 
                                [rectangle_lowpoint[1],rectangle_lowpoint[0]],
                                [rectangle_highpoint[1],rectangle_highpoint[0]],
                                n.transpose(propmatrix2), 
                                vertspacing,
                                horspacing)
        propmatrix2 = n.transpose(propmatrix2)
    else:        
        propmatrix2 = integrate(c2, m2, 
                                rectangle_lowpoint, 
                                rectangle_highpoint, 
                                propmatrix2, 
                                horspacing, 
                                vertspacing)

    # Specific case: at around 90 deg, the lines may go in opposite directions,
    # leading to wrong integrals.
    # Check if one of the lines goes in the opposite direction as the other

    # HOS NOTE - I do not think this part is necessary in the Python implemention 

    if not ( n.sign(m1) == n.sign(m2) ):
        #print "DDDDDDDDDDDDDDDDDD IN DDDDDDDDDDDDDDDDDDDDDD"
        propmatrix2[rectangle_lowpoint[0]:rectangle_highpoint[0],
                    rectangle_lowpoint[1]:rectangle_highpoint[1]] = \
            n.ones(propmatrix2[rectangle_lowpoint[0]:rectangle_highpoint[0],
                             rectangle_lowpoint[1]:rectangle_highpoint[1]].shape) - \
                 propmatrix2[rectangle_lowpoint[0]:rectangle_highpoint[0],
                             rectangle_lowpoint[1]:rectangle_highpoint[1]]

    #if n.sign(m1) != n.sign(m2):   #This corrects for "reverse integral" in case line 1 has 
    #        reflecta = -1  #a negative gradient but line 2 has a positive gradient

    #print 'propmatrix1: \n',propmatrix1[19:23,6:9]


    
    #print 'propmatrix2: \n',propmatrix2[19:23,6:9]
    #print 'reflecta: ',reflecta

    propmatrix_no1 = reflecta*(propmatrix2 - propmatrix1)

#     figure(1)
#     imshow(propmatrix_no1,interpolation='nearest')
#     show()
#     close(1)

#    print 'propmatrix_no1: \n',propmatrix_no1

    r1 = p1.copy()
    r2 = p2.copy()

    # Exchange r21 with r12 (for the definition of the new lines)
    tmp = r2[0].copy()
    r2[0] = r1[1].copy()
    r1[1] = tmp.copy()
    #print 'r1 reset ex',r1
    #print 'r2 reset ex',r2


    # We always want translations in the positive direction horizontally, take
    # this into account when creating the lines... I.e. if the translation on
    # the horizontal axis isn't in the positive direction then rearrange the
    # points such that it will be.
    if r1[1,0]-r1[0,0] < 0:
        #print 'RADIAL_INT: FLIPS r1 UD for 3/4'
        r1 = n.flipud(r1)
    if r2[1,0]-r2[0,0] < 0:
        #print 'RADIAL_INT: FLIPS r2 UD for 3/4'
        r2 = n.flipud(r2)


    # Make line equations into the form y = mz + c.
    # if the gradient would be infinity (that is when the reflection is exactly at
    # eta=90 or eta=270 - it doesn't happen at the vertical spindle axis because
    # reflections are filtered there), then set it to a very large number

    diff1 = r1[1,0]-r1[0,0]
    diff2 = r2[1,0]-r2[0,0]

    if diff1==0:
        diff1 = 10e-30
    if diff2 == 0:
        diff2 = 10e-30


    m3 = (r1[1,1]-r1[0,1])/diff1
    c3 = r1[1,1]-m3*r1[1,0]
    m4 = (r2[1,1]-r2[0,1])/diff2
    c4 = r2[1,1]-m4*r2[1,0]


    # Make sure line 3 is under line 4

    if r1[0,1] > r2[0,1] and n.sign(m3) == n.sign(m4) or \
         n.sign(m3) == -1 and n.sign(m4) ==  1 and r1[1,0] < r2[0,0] or \
         n.sign(m3) ==  1 and n.sign(m4) == -1 and r1[0,0] > r2[1,0]:
        #print 'permute line 3 and 4'
        tmp = m3
        m3 = m4
        m4 = tmp
        tmp = c3
        c3 = c4
        c4 = tmp


    # Initialize propmatrices
    propmatrix3 = n.zeros(sizeoflayer)
    propmatrix4 = n.zeros(sizeoflayer)

    propmatrix3 = integrate(c3, m3, rectangle_lowpoint, rectangle_highpoint, propmatrix3, horspacing, vertspacing)
    propmatrix4 = integrate(c4, m4, rectangle_lowpoint, rectangle_highpoint, propmatrix4, horspacing, vertspacing)


    # Specific case: at around 90 deg, the lines may go in opposite directions,
    # leading to wrong integrals.
    # Check if one of the lines goes in the opposite direction as the other
#     if not (n.sign(m3) == n.sign(m4)):
#         propmatrix4[rectangle_lowpoint[0]:rectangle_highpoint[0],rectangle_lowpoint[1]:rectangle_highpoint[1]] = \
#             n.ones(propmatrix4[rectangle_lowpoint[0]:rectangle_highpoint[0],
#                                rectangle_lowpoint[1]:rectangle_highpoint[1]].shape) -\
#                    propmatrix4[rectangle_lowpoint[0]:rectangle_highpoint[0],
#                                rectangle_lowpoint[1]:rectangle_highpoint[1]]
#         if n.sign(m3)<0:  #This corrects for "reverse integral" in case line 1 has a
#             reflectb=-1   # negative gradient but line 2 has a positive gradient


#    print 'propmatrix3: \n',propmatrix3
#    print 'propmatrix4: \n',propmatrix4

    propmatrix_no2 = reflectb*(propmatrix4 - propmatrix3)
#     figure(1)
#     imshow(propmatrix4,interpolation='nearest')
#     show()
#     close(1)

#    print 'propmatrix_no2:\n',propmatrix_no2

    # Finally find the propmatrix
    propmatrix = propmatrix_no1*propmatrix_no2
#    if propmatrix.min() < 0.0:
#        raise ValueError, 'The lines calculated from the the points r1, r2 cross'
    return propmatrix

    # end of RADIAL_INT






def integrate(c, m, rectangle_lowpoint, rectangle_highpoint, propmatrix_n, horspacing, vertspacing):
 #Find integral under line
 #Loop over columns in propmatrix_n starting from rectangle_lowpoint
 pixel_hit1 = n.zeros((2))
 pixel_hit2 = n.zeros((2))
 #DEBUG print 'UVM-INT rect lowpoint',rectangle_lowpoint
 #DEBUG print 'UVM-INT rect highpoint',rectangle_highpoint
 #DEBUG print 'UVM-INT columno range',range(int(rectangle_lowpoint[0]),int(rectangle_highpoint[0]))

 for columno in range(int(rectangle_lowpoint[0]),int(rectangle_highpoint[0])):
     
    #Find maximum value of line in that column : left or right side of column
    
    if m == n.inf:
        maxvalue = c#rectangle_highpoint[1]
    else:
        maxvalue = n.max([c+m*(columno+1), c+m*(columno)])
    #print 'maxvalue', maxvalue

    #If maxvalue is smaller than the 2nd coordinate of rectangle_lowpoint
    #then set entire column to zero. Else find minvalue and integrate
    #proportionally 

    if maxvalue < rectangle_lowpoint[1]:
        #propmatrix_n(:,columno) = zeros(size(map(:,1,floor(my_omega))));
        # ^ not necessary since already zero by definition ^
        #print 'DROPS OUT OF integrate'
        pass
    else: 
        #minimum value of the line in that column: left or right side
        if m == n.inf:
            minvalue =  c#rectangle_lowpoint[1]
        else:            
            minvalue = n.min([c+m*(columno), c+m*(columno+1)])
        #print 'minvalue', minvalue

        #Set all the pixels below floor(minvalue) to 1
        if rectangle_lowpoint[1] < n.floor(minvalue):
            propmatrix_n[columno,
                         rectangle_lowpoint[1]:n.min([n.floor(minvalue),
                                                      rectangle_highpoint[1]])]\
                         = n.ones(propmatrix_n[columno,
                                               rectangle_lowpoint[1]:n.min([n.floor(minvalue),
                                                                            rectangle_highpoint[1]])].shape)
        
        #Set the remaining pixels (between floor(minvalue) and
        #ceil(maxvalue) proportionately... 
            
        #Area calculation method: calculate where the line intersects
        # the pixel in pixel coordinates, then find path length, then
        # set up triangle and rectangle areas. 
            
        # If sign(m)>0 then gradient is positive, if sign(m)<0 then
        # gradient is negative. 
            
        #OBS subtracted 1 from rectangle_lowpoint[1] to make rowno not extend outside matrix
        # Probably should do this in making rectangle_lowpoint and highpoint

        #print rectangle_lowpoint[1], n.floor(minvalue), rectangle_highpoint[1],n.ceil(maxvalue)

        for rowno in range(n.max([int(rectangle_lowpoint[1]),  int(n.floor(minvalue))]), 
                           n.min([int(rectangle_highpoint[1]), int(n.ceil(maxvalue))])):
                #DEBUG print "WE ARE IN ROWNO", rowno
                #Pixel intersections. 
                #By design, pixel_hit1 has the lower horizontal axis value
                pixel_hit1[0] = n.max([columno,
                                       n.min([(rowno-c)/m, (rowno+1-c)/m])])
                pixel_hit1[1] = m*pixel_hit1[0]+c
                
                pixel_hit2[0] = n.min([columno+1, 
                                     n.max([(rowno-c)/m, (rowno+1-c)/m])])
                pixel_hit2[1] = m*pixel_hit2[0]+c
                
                #Convert to pixel coordinates, ie normalize
                col_row = n.array([columno, rowno])
                hor_vert_spacing = n.array([horspacing, vertspacing])
                pixel_hit1 = (pixel_hit1-col_row)/hor_vert_spacing
                pixel_hit2 = (pixel_hit2-col_row)/hor_vert_spacing
                
                #Triangle area
                triangle_area = 0.5*n.abs((pixel_hit1[0]-pixel_hit2[0])*(pixel_hit1[1]-pixel_hit2[1]))
                
                #Rectangle areas
                
                rec_area1 = pixel_hit1[0]*pixel_hit1[1]   #area to the bottom left side of hit1 
                rec_area2 = (1-pixel_hit2[0])*pixel_hit2[1]   #area to the bottom right side of hit2
                rec_area3 = (pixel_hit2[0]-pixel_hit1[0])*n.min([pixel_hit1[1], pixel_hit2[1]])             
                proportion = triangle_area + rec_area1 + rec_area2 + rec_area3
               
                propmatrix_n[columno,rowno] = proportion         
                #print 'propmatrix_n',propmatrix_n

 return propmatrix_n
