Блог

Мои проблемы состоят в следующем: мне даны две пары углов (в сферических координатах), которые состоят из двух частей — азимута и угла широты. Если мы бесконечно удлиним оба угла (тем самым увеличивая их соответствующие радиусы), чтобы создать длинную линию, указывающую в направлении, заданном парой углов, то моя цель — определить

1. если они пересекаются или находятся очень близко друг к другу и
2. где именно они пересекаются.

В настоящее время я попробовал несколько методов:

1. Наиболее очевидным является итеративное сравнение каждого радиуса до тех пор, пока между ними не будет либо совпадения, либо достаточно малого расстояния. (Когда я говорю сравнить каждый радиус, я имею в виду преобразование каждой сферической координаты в декартову, а затем нахождение евклидова расстояния между ними). Однако эта среда выполнения равна $ O (n ^ {2}) $, что чрезвычайно медленно, если я пытаюсь масштабировать эту программу

2. Второй наиболее очевидный метод — использовать пакет оптимизации для определения этого расстояния. К сожалению, я не могу использовать пакет оптимизации итеративно, и после одного экземпляра алгоритм оптимизации повторяет один и тот же ответ, что бесполезно.

3. Наименее очевидным методом является прямое вычисление (с использованием исчисления) точных радиусов от углов. Хотя это быстрый метод, он не очень точен.

Примечание: хотя может показаться простым, что пересечение всегда находится в нулевом начале координат (0,0,0), это не ВСЕГДА так. Некоторые точки никогда не пересекаются.

Код для метода (1)

 def match1(azimuth_recon_1,colatitude_recon_1,azimuth_recon_2, colatitude_recon_2,centroid_1,centroid_2 ):


   # Constants: tolerance factor and extremely large distance
    tol = 3e-2
    prevDist = 99999999

    # Initialize a list of radii to loop through
    # Checking iteravely for a solution
    for r1 in list(np.arange(0,5,tol)):
        for r2 in list(np.arange(0,5,tol)):

            # Get the estimates
            estimate_1 = np.array(spher2cart(r1,azimuth_recon_1,colatitude_recon_1))   np.array(centroid_1)
            estimate_2 = np.array(spher2cart(r2,azimuth_recon_2,colatitude_recon_2))  np.array(centroid_2)

            # Calculate the euclidean distance between them
            dist = np.array(np.sqrt(np.einsum('i...,i...', (estimate_1 - estimate_2), (estimate_1 - estimate_2)))[:,np.newaxis])

        # Compare the distance to this tolerance
            if dist < tol: 
                if dist == 0:
                    return estimate_1, [], True
                else:
                    return estimate_1, estimate_2, False

            ## If the distance is too big break out of the loop
            if dist > prevDist:
                prevDist = 9999999
                break
            prevDist = dist

    return [], [], False

  

Код для метода (3)

 def match2(azimuth_recon_1,colatitude_recon_1,azimuth_recon_2, colatitude_recon_2,centriod_1,centroid_2):

   # Set a Tolerance factor
    tol = 3e-2

    def calculate_radius_2(azimuth_1,colatitude_1,azimuth_2,colatitude_2):

        """Return radius 2 using both pairs of angles (azimuth and colatitude). Equation is provided in the document"""

        return 1/((1-(math.sin(azimuth_1)*math.sin(azimuth_2)*math.cos(colatitude_1-colatitude_2))
         math.cos(azimuth_1)*math.cos(azimuth_2))**2)

    def calculate_radius_1(radius_2,azimuth_1,colatitude_1,azimuth_2,colatitude_2):

        """Returns radius 1 using both pairs of angles (azimuth and colatitude) and radius 2. 
    Equation provided in document"""

        return (radius_2)*((math.sin(azimuth_1)*math.sin(azimuth_2)*math.cos(colatitude_1-colatitude_2))
         math.cos(azimuth_1)*math.cos(azimuth_2))

    # Compute radius 2
    radius_2 = calculate_radius_2(azimuth_recon_1,colatitude_recon_1,azimuth_recon_2,colatitude_recon_2)

    #Compute radius 1
    radius_1 = calculate_radius_1(radius_2,azimuth_recon_1,colatitude_recon_1,azimuth_recon_2,colatitude_recon_2)

    # Get the estimates
    estimate_1 = np.array(spher2cart(radius_1,azimuth_recon_1,colatitude_recon_1))  np.array(centroid_1)
    estimate_2 = np.array(spher2cart(radius_2,azimuth_recon_2,colatitude_recon_2))  np.array(centroid_2) 

    # Calculate the euclidean distance between them
    dist = np.array(np.sqrt(np.einsum('i...,i...', (estimate_1 - estimate_2), (estimate_1 - estimate_2)))[:,np.newaxis])

    # Compare the distance to this tolerance
    if dist < tol: 
        if dist == 0:
            return estimate_1, [], True
        else:
            return estimate_1, estimate_2, False
    else:
        return [], [], False
  

Мой вопрос двоякий:

1. Есть ли более быстрый и точный способ найти радиусы для обеих точек?

2. Если да, то как мне это сделать?

РЕДАКТИРОВАТЬ: я думаю о том, чтобы просто создать два numpy-массива из двух радиусов, а затем сравнить их с помощью numpy boolean logic. Тем не менее, я бы все равно сравнивал их итеративно. Есть ли более быстрый способ выполнить это сравнение?

Используйте kd-tree для таких ситуаций. Он легко найдет минимальное расстояние:

 def match(azimuth_recon_1,colatitude_recon_1,azimuth_recon_2, colatitude_recon_2,centriod_1,centroid_2):

    cartesian_1 = np.array([np.cos(azimuth_recon_1)*np.sin(colatitude_recon_1),np.sin(azimuth_recon_1)*np.sin(colatitude_recon_1),np.cos(colatitude_recon_1)]) #[np.newaxis,:]
    cartesian_2 = np.array([np.cos(azimuth_recon_2)*np.sin(colatitude_recon_2),np.sin(azimuth_recon_2)*np.sin(colatitude_recon_2),np.cos(colatitude_recon_2)]) #[np.newaxis,:]

    # Re-center them via adding the centroid
    estimate_1 = r1*cartesian_1.T   np.array(centroid_1)[np.newaxis,:] 
    estimate_2 = r2*cartesian_2.T   np.array(centroid_2)[np.newaxis,:]

    # Add them to the output list
    n = estimate_1.shape[0]
    outputs_list_1.append(estimate_1)
    outputs_list_2.append(estimate_2)

    # Reshape them so that they are in proper format
    a = np.array(outputs_list_1).reshape(len(two_pair_mic_list)*n,3)
    b = np.array(outputs_list_2).reshape(len(two_pair_mic_list)*n,3)

    # Get the difference
    c = a - b

    # Put into a KDtree
    tree = spatial.KDTree(c)

    # Find the indices where the radius (distance between the points) is 3e-3 or less 
    indices = tree.query_ball_tree(3e-3)
  

Это выведет список индексов, расстояние между которыми составляет 3e-3 или меньше. Теперь все, что вам нужно будет сделать, это использовать список индексов со списком оценок, чтобы найти точные точки. И вот оно, это сэкономит вам много времени и пространства!