[Из песочницы] Ray Casting Visual Search (RCVS). Простой и быстрый алгоритм поиска схожих по геометрии 3D моделей

u28ifwuqbqnlhhiwwldggmtayp4.png

Для меня эти две модели очень похожи, однако у них нет очевидных характеристик, по которым можно было бы измерить их сходство. У этих моделей разное количество вершин, рёбер и полигонов, они разного размера, к тому же по-разному повёрнуты в пространстве, и у обеих одинаковые трансформации (Положение = [0,0,0], Вращение в радианах = [0,0,0], Масштаб = [1,1,1]). Как определить их подобие?
Если бы эти модели были одинаково ориентированы в пространстве, их можно было бы свести к общему размеру Bounding Box, переместить начало локальных координат в его центр, привести к воксельному виду и сравнить списки вокселей.

hxhkm7cupsn8fts2jgjmpmqmg3o.png
Обезьянки и их воксели.

Если бы эти модели были произвольно ориентированы вдоль осей координат, их можно было бы свести к общему размеру Bounding Box, переместить начало локальных координат в его центр, порендерить их с шести сторон и сравнивать наборы изображений между собой.

1ot1g00gyan9sroahuclx20bbo0.png
Рендеры кривизны поверхности.

Для того, чтобы добиться инвариантности поворота, я решил обратиться к сферической системе координат.

$r = \sqrt{x^2+y^2+z^2}$

Сначала я считал расстояния до вокселей от начала координат, собирал расстояния в упорядоченные списки и сравнивал эти списки. Это дало неплохие результаты, но требовало много времени для приведения к воксельному виду достаточного разрешения и сравнения больших списков.

Краткое описание RCVS


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

2einitowujo6sygimeu5qcsh6ss.gifp3nngs1qtrc92rch2wbfylqgkeu.gif
80 и 1280 лучей соответственно.

Результаты сравнения моделей обезьянки Сюзанны в разных вариантах.

bd2x0mkajcvvj3qlzjeexr1xu6e.png

Suzanne_lo
                          Suzanne_lo | 1             
                 Suzanne_lo_rot_rand | 0.9878376875939173
                   Suzanne_lo_rot_90 | 0.9800766750109137
                 Suzanne_hi_rot_rand | 0.9761078629585936
                          Suzanne_hi | 0.9730722945028376
                   Suzanne_hi_rot_90 | 0.9680325422930756
                               Skull | 0.7697333228143514
             Ceramic_Teapot_rot_rand | 0.6634783235048608
                      Ceramic_Teapot | 0.6496529954470333
                         Grenade_MK2 | 0.5275616759076022


Подробное описание RCVS


Иллюстрации шагов


Перед сравнением 3D моделей (далее объектов), каждая проходит этап подготовки (эту часть я реализовал на Python в Blender):

  1. Начало локальных координат компонентов объекта устанавливается в центр массы объёма.
  2. Объект вписывается в сферу единичного радиуса с центром в начале локальных координат компонентов объекта.
  3. Применяются трансформации размера.
  4. Вокруг объекта строится икосаэдральный сферический многогранник (далее икосфера) единичного радиуса с нормалями, направленными внутрь. Количество полигонов икосферы, на которых проводились тесты: 80, 320, 1280, 5120, 22480. Оптимальным количеством полигонов икосферы по времени и точности в тестах было определено 1280.
  5. Из каждого полигона икосферы по нормали отправляется луч и сохраняется его путь до столкновения с поверхностью объекта, нормаль которой направлена в сторону луча.
  6. Нормали объекта переворачиваются в противоположном направлении.
  7. Из каждого полигона икосферы ещё раз отправляются лучи и сохраняется их путь до встречи с поверхностью объекта, нормали которой направлены внутрь объекта. (Это фиксирует внутреннюю геометрию и/или заполняет слепые зоны предыдущих лучей)
  8. Пути лучей, отправленных с противоположных полигонов компонуются в списки, где первые два пути из пункта 5 и вторые два из пункта 7. Каждый список сортируется по наименьшему значению пути из пункта 5 так, что пути из пунктов 5 и 7 попарно соответствуют друг другу.
  9. Список списков всех путей сортируется по первому значению. Этот список имеет 4 столбца и количество строк, равное половине количества полигонов икосферы.
  10. Список путей лучей записывается в файл.


Список путей лучей выглядит примерно так:

0.00271218840585662  0.2112752957162676 0.00271218840585662  0.2112752957162676
0.004740002405006037 0.210980921765861  0.004740002405006037 0.2109809188911709
0.00660892842734402  0.2066613370130234 0.00660892842734402  0.2066613370130234
…
0.2692299271137439   0.2725085782639611 0.2692299271137439   0.2725085782639611
0.2705489991382905   0.2707842753523402 0.270548999138295    0.2707843056356914
0.2709498404192674   0.271036677664217  0.2709498404192674   0.271036642944409


Когда все объекты получили готовые списки путей лучей их можно сравнивать (эту часть я написал на Rust):

  1. Задаётся пороговое значение допустимой ошибки. В тестах наилучшим было определено значение 0.01, ошибка в пределах 1%.
  2. Создаётся счётчик соответствия.
  3. Для каждой строки в первом списке ищется строка из второго, каждое значение которой отличается от соответствующего значения из первого списка в пределах допустимой ошибки. Если такая строка находится, она не участвует в дальнейшем поиске. Вычисляется попарная абсолютная разница значений из двух списков и вычитается из единицы, среднее арифметическое этих значений прибавляется к счётчику.
  4. Когда пройдены все строки первого списка значение счётчика делится на количество строк. Полученное значение является результатом сравнения.
  5. Если объекты идентичны по форме (независимо от их размера и поворота) результат близок к единице.


Бананы, Витрувианский человек и ограничения RCVS


Алгоритм отлично работает для поиска близких по геометрии объектов, например высоко- и низко-полигональных вариантов одной и той же модели. И иногда неплохо справляется с похожими по смыслу объектами, например с чайниками или гитарами разных форм. Однако, поскольку количество моделей, которым будет соответствовать один и тот же список путей лучей примерно можно оценить в $\frac{n!}{\left(\frac{n}{2}\right)!}$, есть существенная вероятность, что визуально непохожие по геометрии объекты будут иметь счёт совпадения в диапазоне (0.5 — 0.9). В ходе тестов так вышло с бананами, которые на 0.7 совпадают с самолётами и людьми в разных позах, которые в свою очередь на «внушительные» 0.85 совпадают с тиграми и собаками.

Для более точного поиска стоит учитывать размер и/или поворот объектов.

Исходники.


Код на Python для создания словаря противоположных полигонов икосферы в Blender
import bpy
import os

name = bpy.context.active_object.data.name

polys = bpy.data.meshes[name].polygons
print("-----")
length = len(polys)

x = 0
p = 0
dict = []
for i in range(0, length):
    print(i, "/", length)
    for l in range(i, length):
        norm1 = polys[i].normal
        norm2 = polys[l].normal
        if norm1 == -norm2:
            dict.append(str(i)+", "+str(l)) 
            p += 1  


print("Polys", p)
basedir = os.path.dirname(bpy.data.filepath) + "/ico_dicts/"  
if not os.path.exists(basedir):
    os.makedirs(basedir)
file = open(basedir + "/" + str(length) + ".csv","w")
for poly in dict:      
    file.write(poly + "\n")
file.close()
print("----ok----")


Код на Python для подготовки объектов и создания списка путей лучей в Blender.
import bpy
import os
from math import sqrt
import csv


#Icosphere properties
#icosubd:
#5 -> 5120 X 2 rays
#4 -> 1280 X 2 rays
#3 ->  320 X 2 rays
#2 ->   80 X 2 rays
icosubd = 4
size = 1280 #name of directory
radius = 1




def icoRC(icosubd, rad, raydist, normals, target):
    icorad = rad
    bpy.ops.mesh.primitive_ico_sphere_add(radius=icorad, subdivisions=icosubd, enter_editmode=True, location=(0, 0, 0))
    bpy.ops.mesh.normals_make_consistent(inside=normals)
    bpy.ops.object.editmode_toggle()
    

    ico = bpy.context.active_object
    mesh = bpy.data.meshes[ico.data.name]
    rays = []
    size = len(mesh.polygons)
    for poly in mesh.polygons:
        normal = poly.normal
        start = poly.center
        
        ray = target.ray_cast(start, normal, distance=raydist)
      
        if ray[0] == True:
            length = sqrt((start[0]-ray[1][0]) ** 2 + (start[1]-ray[1][1]) ** 2 + (start[2]-ray[1][2]) ** 2)
            rays.append(length / (radius *2))
        else:
            rays.append(-1)
            
        
    result = []
    #Combine opposite rays using CSV dicts
    dictdir = os.path.dirname(bpy.data.filepath) + "/ico_dicts/"
    with open(dictdir+str(size)+'.csv', newline='') as csvfile:
        ico_dict_file = csv.reader(csvfile, delimiter=',')
        for row in ico_dict_file:
            result.append([rays[int(row[0])], rays[int(row[1])]])
            
            
    bpy.ops.object.select_all(action='DESELECT')
    bpy.data.objects[ico.name].select_set(True) 
    bpy.ops.object.delete() 
    
    return result
   
   
    
    
#Batch preparation of objects by scale and origin 
def batch_prepare():   
    for obj in bpy.context.selected_objects:
        bpy.ops.object.origin_set(type='ORIGIN_CENTER_OF_VOLUME', center='MEDIAN')
        max_dist = 0
        bb = bpy.data.objects[obj.name].bound_box
        
        for vert in obj.data.vertices:
            dist = vert.co[0]**2 + vert.co[1]**2 + vert.co[2]**2
            if dist > max_dist:
                max_dist = dist
        
        obj.scale *= (1 / sqrt(max_dist))



def write_rays(dir):  
    for obj in bpy.context.selected_objects:
        bpy.ops.object.transform_apply(location=False, rotation=False, scale=True)
        
        #Cast rays on object with normals turned inside   
        bpy.context.view_layer.objects.active  = obj
        bpy.ops.object.editmode_toggle()
        bpy.ops.mesh.normals_make_consistent(inside=True)
        bpy.ops.object.editmode_toggle()
        result1 = icoRC(icosubd, radius, radius*2, True, obj)
        
        #Cast rays on object with normals turned outside 
        bpy.context.view_layer.objects.active  = obj
        bpy.ops.object.editmode_toggle()
        bpy.ops.mesh.normals_make_consistent(inside=False)
        bpy.ops.object.editmode_toggle()
        result2 = icoRC(icosubd, radius, radius*2, True, obj)
        
        final = []
        #Sort sub-lists and append them to main list
        for i in range(0, len(result1)):
            if result2[i][0] < result2[i][1]:
                final.append([result2[i][0], result2[i][1], result1[i][0], result1[i][1]])
            else:
                final.append([result2[i][1], result2[i][0], result1[i][1], result1[i][0]])
        
        
        
     
        basedir = os.path.dirname(bpy.data.filepath) + "/rays/" + str(size) + "/" 
        if not os.path.exists(basedir):
            os.makedirs(basedir)
        file = open(basedir + "/" + obj.name,"w")
        
        final.sort(key=lambda x: x[0]) #Sort list by first values in sub-lists
     
        #Writing to file   
        for ray in final:
            ray_str = ""
            for dist in ray:   
                ray_str += str(dist) + " "
            file.write(ray_str + "\n")

        file.close()
        
    
    
batch_prepare()
write_rays(size)


Код на Rust для сравнения списков.
use rayon;
use std::collections::HashSet;
use std::fs;
use std::io::{BufRead, BufReader, Error, Read};

//Creates Vector from file
fn read(io: R) -> Result>, Error> {
    let br = BufReader::new(io);
    br.lines()
        .map(|line| {
            line.and_then(|v| {
                Ok(v.split_whitespace()
                    .map(|x| match x.trim().parse::() {
                        Ok(value) => value,
                        Err(_) => -1.0,
                    })
                    .collect::>())
            })
        })
        .collect()
}

//Compares two objects
fn compare(one: &Vec>, two: &Vec>, error: f64) -> f64 {
    let length = one.len();
    let mut count: f64 = 0.0;
    let mut seen: HashSet = HashSet::new();

    for orig in one {
        let size = orig.len();

        for i in 0..length {
            if seen.contains(&i) {
                continue;
            } else {
                let targ = &two[i];

                if (0..size).map(|x| (orig[x] - targ[x]).abs()).sum::() > error * size as f64 {
                    continue;
                }

                let err_sum = (0..size)
                    .map(|x| 1.0 - (orig[x] - targ[x]).abs())
                    .sum::();

                count += err_sum / size as f64;
                seen.insert(i);
                break;
            }
        }
    }
    count / length as f64
}

//Compares one object with others on different threads
fn one_with_others(obj: &(String, Vec>), dict: &Vec<(String, Vec>)>) {
    println!("\n{:}\n", obj.0);
    rayon::scope(|s| {
        for obj2 in dict {
            let obj = obj.clone();
            s.spawn(move |_| {
                let result = compare(&obj.1, &obj2.1, 0.01);
                println!(
                    "{:>36} |{:<100}| {:<14}",
                    obj2.0,
                    "▮".repeat((result * 100.0) as usize),
                    result,
                )
            });
        }
    });
}

fn main() {
    let ray = 1280; //name of directory
    let obj_cur = "Suzanne_lo"; //name of object

    let mut objects = Vec::<(String, Vec>)>::new();
    let path_rays = format!("./rays/{:}/", ray);
    let paths = fs::read_dir(path_rays).unwrap();

    let obj_path = format!("./rays/{:}/{:}", ray, obj_cur);
    let _obj = (
        obj_cur.to_string(),
        read(fs::File::open(obj_path).unwrap()).unwrap(),
    );

    for path in paths {
        let dir = &path.unwrap().path();
        let name = dir.to_str().unwrap().split("/").last().unwrap();

        objects.push((
            name.to_string(),
            read(fs::File::open(dir).unwrap()).unwrap(),
        ));
    }

    one_with_others(&_obj, &objects);
}


Ссылка на исходники и проект Blender с тестовой сценой.

oh4p7sqgz6j9ti6xi8pjmg73li4.png

© Habrahabr.ru