Registro de conjuntos de puntos

Dados dos conjuntos de puntos, el registro rígido produce una transformación rígida que mapea un conjunto de puntos al otro. Una transformación rígida se define como una transformación que no cambia la distancia entre dos puntos cualesquiera. Típicamente, dicha transformación consiste en traslación y rotación . ^{[ 12 ]} En raras ocasiones, el conjunto de puntos también puede reflejarse. En robótica y visión por computadora, el registro rígido tiene la mayor cantidad de aplicaciones.

Dados dos conjuntos de puntos, el registro no rígido produce una transformación no rígida que mapea un conjunto de puntos al otro. Las transformaciones no rígidas incluyen transformaciones afines como el escalado y el mapeo de cizallamiento . Sin embargo, en el contexto del registro de conjuntos de puntos, el registro no rígido generalmente implica una transformación no lineal. Si se conocen los modos propios de variación del conjunto de puntos, la transformación no lineal puede parametrizarse mediante los valores propios. ^{[ 13 ]} Una transformación no lineal también puede parametrizarse como una spline de placa delgada . ^{[ 14 ] [ 13 ]}

Algunos métodos de registro de conjuntos de puntos emplean algoritmos que resuelven el problema más general de correspondencia de grafos . ^{[ 11 ]} Sin embargo, la complejidad computacional de estos métodos suele ser elevada y se limitan a registros rígidos. En este artículo, solo consideraremos algoritmos para el registro rígido, donde se supone que la transformación incluye rotaciones y traslaciones 3D (posiblemente también un escalado uniforme).

La PCL (Point Cloud Library) es un marco de código abierto para el procesamiento de nubes de puntos n-dimensionales y geometría 3D . Incluye varios algoritmos de registro de puntos. ^{[ 15 ]}

En el caso más simple, se puede suponer que todas las correspondencias son correctas, lo que significa que los puntos${\displaystyle m_{i},s_{i}\in \mathbb {R} ^{3}}$se generan de la siguiente manera:

dónde${\displaystyle l>0}$es un factor de escala uniforme (en muchos casos${\displaystyle l=1}$se supone),${\displaystyle R\in {\text{SO}}(3)}$es una matriz de rotación 3D adecuada (${\displaystyle {\text{SO}}(d)}$es el grupo ortogonal especial de grado${\displaystyle d}$),${\displaystyle t\in \mathbb {R} ^{3}}$es un vector de traslación 3D y${\displaystyle \epsilon _{i}\in \mathbb {R} ^{3}}$modela el ruido aditivo desconocido ( por ejemplo, ruido gaussiano ). Específicamente, si el ruido${\displaystyle \epsilon _{i}}$Se supone que sigue una distribución gaussiana isotrópica de media cero con desviación estándar${\displaystyle \sigma _{i}}$, es decir,${\displaystyle \epsilon _{i}\sim {\mathcal {N}}(0,\sigma _{i}^{2}I_{3})}$, entonces se puede demostrar que la siguiente optimización produce la estimación de máxima verosimilitud para la escala, rotación y traslación desconocidas:

Nótese que cuando el factor de escala es 1 y el vector de traslación es cero, la optimización recupera la formulación del problema de Wahba . A pesar de la no convexidad de la optimización ( cb.2 ) debido a la no convexidad del conjunto${\displaystyle {\text{SO}}(3)}$El trabajo fundamental de Berthold KP Horn demostró que ( cb.2 ) admite en realidad una solución de forma cerrada, al desacoplar la estimación de escala, rotación y traslación. ^{[ 16 ]} Resultados similares fueron descubiertos por Arun et al . ^{[ 17 ]} Además, para encontrar una transformación única${\displaystyle (l,R,t)}$, al menos${\displaystyle N=3}$Se requieren puntos no colineales en cada conjunto de puntos.

Más recientemente, Briales y González-Jiménez han desarrollado una relajación semidefinida utilizando la dualidad lagrangiana , para el caso en que el conjunto de modelos${\displaystyle {\mathcal {M}}}$contiene diferentes primitivas 3D como puntos, líneas y planos (que es el caso cuando el modelo${\displaystyle {\mathcal {M}}}$es una malla 3D). ^{[ 18 ]} Curiosamente, la relajación semidefinida es empíricamente ajustada, es decir, se puede extraer una solución globalmente óptima certificable de la solución de la relajación semidefinida.

Se sabe que la formulación de mínimos cuadrados ( cb.2 ) tiene un rendimiento arbitrariamente malo en presencia de valores atípicos . Una correspondencia de valores atípicos es un par de mediciones${\displaystyle s_{i}\leftrightarrow m_{i}}$que se aparta del modelo generativo ( cb.1 ). En este caso, se puede considerar un modelo generativo diferente como sigue: ^{[ 19 ]}

donde si el${\displaystyle i-}$el par${\displaystyle s_{i}\leftrightarrow m_{i}}$es un valor interno, entonces obedece al modelo libre de valores atípicos ( cb.1 ), es decir,${\displaystyle s_{i}}$se obtiene de${\displaystyle m_{i}}$mediante una transformación espacial más algo de ruido pequeño; sin embargo, si la${\displaystyle i-}$el par${\displaystyle s_{i}\leftrightarrow m_{i}}$es un valor atípico, entonces${\displaystyle s_{i}}$puede ser cualquier vector arbitrario${\displaystyle o_{i}}$Dado que no se sabe de antemano qué correspondencias son valores atípicos, el registro robusto bajo el modelo generativo ( cb.3 ) es de suma importancia para la visión por computadora y la robótica implementadas en el mundo real, porque las técnicas actuales de coincidencia de características tienden a producir correspondencias altamente corruptas donde más${\displaystyle 95\%}$de las correspondencias pueden ser valores atípicos. ^{[ 20 ]}

A continuación, describimos varios paradigmas comunes para un registro robusto.

El algoritmo iterativo de punto más cercano (ICP) fue introducido por Besl y McKay. ^{[ 36 ]} El algoritmo realiza un registro rígido de forma iterativa alternando en (i) dada la transformación, encontrando el punto más cercano en${\displaystyle {\mathcal {S}}}$por cada punto en${\displaystyle {\mathcal {M}}}$; y (ii) dadas las correspondencias, encontrar la mejor transformación rígida resolviendo el problema de mínimos cuadrados ( cb.2 ). Como tal, funciona mejor si la pose inicial de${\displaystyle {\mathcal {M}}}$está suficientemente cerca de${\displaystyle {\mathcal {S}}}$En pseudocódigo , el algoritmo básico se implementa de la siguiente manera:

algoritmo ICP( M , S ) θ := θ 0 mientras no esté registrado: X := ∅ para m i ∊ T ( M , θ ): ŝ i := punto más cercano en S a m i X := X + ⟨ m i , ŝ i ⟩ θ := least_squares( X ) devolver θ

Aquí, la función least_squaresrealiza una optimización de mínimos cuadrados para minimizar la distancia en cada uno de los${\displaystyle \langle m_{i},{\hat {s}}_{i}\rangle }$pares, utilizando las soluciones de forma cerrada de Horn ^{[ 16 ]} y Arun. ^{[ 17 ]}

Debido a que la función de costo del registro depende de encontrar el punto más cercano en${\displaystyle {\mathcal {S}}}$a cada punto en${\displaystyle {\mathcal {M}}}$, puede cambiar mientras el algoritmo se ejecuta. Por lo tanto, es difícil demostrar que ICP convergerá exactamente al óptimo local. ^{[ 37 ]} De hecho, empíricamente, ICP y EM-ICP no convergen al mínimo local de la función de costo. ^{[ 37 ]} No obstante, debido a que ICP es intuitivo de entender y sencillo de implementar, sigue siendo el algoritmo de registro de conjuntos de puntos más utilizado. ^{[ 37 ]} Se han propuesto muchas variantes de ICP, que afectan todas las fases del algoritmo, desde la selección y el emparejamiento de puntos hasta la estrategia de minimización. ^{[ 13 ] [ 38 ]} Por ejemplo, el algoritmo de maximización de la esperanza se aplica al algoritmo ICP para formar el método EM-ICP, y el algoritmo de Levenberg-Marquardt se aplica al algoritmo ICP para formar el método LM-ICP . ^{[ 12 ]}

El método de correspondencia de puntos robusto (RPM) fue introducido por Gold et al. ^{[ 39 ]} . Este método realiza el registro utilizando recocido determinista y asignación suave de correspondencias entre conjuntos de puntos. Mientras que en ICP la correspondencia generada por la heurística del vecino más cercano es binaria, RPM utiliza una correspondencia suave donde la correspondencia entre dos puntos cualesquiera puede ser cualquier valor entre 0 y 1, aunque finalmente converge a 0 o 1. Las correspondencias encontradas en RPM son siempre uno a uno, lo cual no siempre ocurre en ICP. ^{[ 14 ]} Sea${\displaystyle m_{i}}$ser el${\displaystyle i}$punto en${\displaystyle {\mathcal {M}}}$y${\displaystyle s_{j}}$ser el${\displaystyle j}$punto en${\displaystyle {\mathcal {S}}}$La matriz de coincidencia${\displaystyle \mathbf {\mu } }$se define de la siguiente manera:

El problema se define entonces como: Dados dos conjuntos de puntos${\displaystyle {\mathcal {M}}}$y${\displaystyle {\mathcal {S}}}$encontrar la transformación afín${\displaystyle T}$y la matriz de coincidencia${\displaystyle \mathbf {\mu } }$que mejor los relaciona. ^{[ 39 ]} Conocer la transformación óptima facilita la determinación de la matriz de correspondencia, y viceversa. Sin embargo, el algoritmo RPM determina ambas simultáneamente. La transformación puede descomponerse en un vector de traslación y una matriz de transformación :

${\displaystyle T(m)=\mathbf {A} m+\mathbf {t} }$

La matriz${\displaystyle \mathbf {A} }$en 2D se compone de cuatro parámetros separados${\displaystyle \lbrace a,\theta ,b,c\rbrace }$, que son la escala, la rotación y los componentes de corte vertical y horizontal respectivamente. La función de coste es entonces:

sujeto a${\textstyle \forall j~\sum _{i=1}^{M}\mu _{ij}\leq 1}$,${\textstyle \forall i~\sum _{j=1}^{N}\mu _{ij}\leq 1}$,${\textstyle \forall ij~\mu _{ij}\in \lbrace 0,1\rbrace }$. El${\displaystyle \alpha }$El término sesga el objetivo hacia una correlación más fuerte al disminuir el costo si la matriz de coincidencia tiene más unos en ella. La función${\displaystyle g(\mathbf {A} )}$Sirve para regularizar la transformación afín penalizando los valores grandes de los componentes de escala y cizallamiento:

${\displaystyle g(\mathbf {A} (a,\theta ,b,c))=\gamma (a^{2}+b^{2}+c^{2})}$

para algún parámetro de regularización${\displaystyle \gamma }$.

El método RPM optimiza la función de costo utilizando el algoritmo Softassign . Aquí se derivará el caso 1D. Dado un conjunto de variables${\displaystyle \lbrace Q_{j}\rbrace }$dónde${\displaystyle Q_{j}\in \mathbb {R} ^{1}}$. Una variable${\displaystyle \mu _{j}}$está asociado con cada uno${\displaystyle Q_{j}}$de tal manera que${\textstyle \sum _{j=1}^{J}\mu _{j}=1}$El objetivo es encontrar${\displaystyle \mathbf {\mu } }$que maximiza${\textstyle \sum _{j=1}^{J}\mu _{j}Q_{j}}$Esto puede formularse como un problema continuo introduciendo un parámetro de control.${\displaystyle \beta >0}$En el método de recocido determinista , el parámetro de control${\displaystyle \beta }$aumenta lentamente a medida que se ejecuta el algoritmo.${\displaystyle \mathbf {\mu } }$ser:

Esto se conoce como la función softmax .${\displaystyle \beta }$aumenta, se aproxima a un valor binario como se desea en la ecuación ( rpm.1 ). El problema ahora puede generalizarse al caso 2D, donde en lugar de maximizar${\textstyle \sum _{j=1}^{J}\mu _{j}Q_{j}}$Se maximiza lo siguiente:

dónde

${\displaystyle Q_{ij}=-(\lVert s_{j}-\mathbf {t} -\mathbf {A} m_{i}\rVert ^{2}-\alpha )=-{\frac {\partial \operatorname {cost} }{\partial \mu _{ij}}}}$

Esto es sencillo, excepto que ahora las restricciones en${\displaystyle \mu }$son restricciones matriciales doblemente estocásticas :${\textstyle \forall j~\sum _{i=1}^{M}\mu _{ij}=1}$y${\textstyle \forall i~\sum _{j=1}^{N}\mu _{ij}=1}$. Por lo tanto, el denominador de la ecuación ( rpm.3 ) no puede expresarse simplemente para el caso 2D. Para satisfacer las restricciones, es posible utilizar un resultado debido a Sinkhorn, ^{[ 39 ]} que establece que una matriz doblemente estocástica se obtiene a partir de cualquier matriz cuadrada con todas las entradas positivas mediante el proceso iterativo de normalizaciones alternas de filas y columnas. Así, el algoritmo se escribe de la siguiente manera: ^{[ 39 ]}

algoritmo RPM2D${\displaystyle ({\mathcal {M}},{\mathcal {S}})}$t := 0 a , θ segundo , c := 0 β := β 0${\displaystyle {\hat {\mu }}_{ij}:=1+\epsilon }$mientras β < β f : mientras μ no haya convergido:  // actualizar los parámetros de correspondencia mediante asignación suave${\displaystyle Q_{ij}:=-{\frac {\partial \operatorname {cost} }{\partial \mu _{ij}}}}$${\displaystyle \mu _{ij}^{0}:=\exp(\beta Q_{ij})}$// aplicar el método de Sinkhorn mientras${\displaystyle {\hat {\mu }}}$no ha convergido: // actualización${\displaystyle {\hat {\mu }}}$normalizando en todas las filas:${\displaystyle {\hat {\mu }}_{ij}^{1}:={\frac {{\hat {\mu }}_{ij}^{0}}{\sum _{i=1}^{M+1}{\hat {\mu }}_{ij}^{0}}}}$// actualizar${\displaystyle {\hat {\mu }}}$normalizando todas las columnas:${\displaystyle {\hat {\mu }}_{ij}^{0}:={\frac {{\hat {\mu }}_{ij}^{1}}{\sum _{j=1}^{N+1}{\hat {\mu }}_{ij}^{1}}}}$// Actualizar los parámetros de pose mediante descenso de coordenadas. Actualizar θ usando una solución analítica. actualizar t usando solución analítica Actualizar a, b, c usando el método de Newton.${\displaystyle \beta :=\beta _{r}\beta }$${\displaystyle \gamma :={\frac {\gamma }{\beta _{r}}}}$devolver a, b, c, θ y t

donde el parámetro de control de recocido determinista${\displaystyle \beta }$inicialmente está configurado para${\displaystyle \beta _{0}}$y aumenta por factor${\displaystyle \beta _{r}}$hasta que alcance el valor máximo${\displaystyle \beta _{f}}$. Las sumas en los pasos de normalización suman a${\displaystyle M+1}$y${\displaystyle N+1}$en lugar de simplemente${\displaystyle M}$y${\displaystyle N}$porque las restricciones en${\displaystyle \mu }$son desigualdades. Como tales,${\displaystyle M+1}$y${\displaystyle N+1}$Los elementos th son variables de holgura .

El algoritmo también puede extenderse a conjuntos de puntos en 3D o dimensiones superiores. Las restricciones en la matriz de correspondencia${\displaystyle \mathbf {\mu } }$son los mismos en el caso 3D que en el caso 2D. Por lo tanto, la estructura del algoritmo permanece inalterada, siendo la principal diferencia la forma en que se resuelven las matrices de rotación y traslación. ^{[ 39 ]}

Coincidencia de puntos robusta de estrías de placa delgada

El algoritmo de coincidencia de puntos robustos de spline de placa delgada (TPS-RPM) de Chui y Rangarajan amplía el método RPM para realizar un registro no rígido al parametrizar la transformación como un spline de placa delgada . ^{[ 14 ]} Sin embargo, debido a que la parametrización del spline de placa delgada solo existe en tres dimensiones, el método no se puede extender a problemas que involucren cuatro o más dimensiones.

El método de correlación de kernel (KC) para el registro de conjuntos de puntos fue introducido por Tsin y Kanade. ^{[ 37 ]} En comparación con ICP, el algoritmo KC es más robusto frente a datos ruidosos. A diferencia de ICP, donde, para cada punto del modelo, solo se considera el punto de escena más cercano, aquí cada punto de escena afecta a cada punto del modelo. ^{[ 37 ]} Por lo tanto, este es un algoritmo de registro de múltiples enlaces . Para alguna función kernel${\displaystyle K}$, la correlación del núcleo${\displaystyle KC}$de dos puntos${\displaystyle x_{i},x_{j}}$se define así: ^{[ 37 ]}

La función kernel${\displaystyle K}$El núcleo elegido para el registro de conjuntos de puntos suele ser simétrico y no negativo, similar a los utilizados en la estimación de densidad de la ventana de Parzen . El núcleo gaussiano se utiliza habitualmente por su simplicidad, aunque se pueden sustituir otros como el núcleo de Epanechnikov y el núcleo tricubo. ^{[ 37 ]} La correlación del núcleo de un conjunto completo de puntos${\displaystyle {\mathcal {\chi }}}$se define como la suma de las correlaciones del núcleo de cada punto del conjunto con todos los demás puntos del conjunto: ^{[ 37 ]}

El logaritmo de KC de un conjunto de puntos es proporcional, con un factor constante, a la entropía de la información . Obsérvese que KC es una medida de la "compacidad" del conjunto de puntos; obviamente, si todos los puntos del conjunto estuvieran en la misma ubicación, KC tendría un valor elevado. La función de coste del algoritmo de registro de conjuntos de puntos para algún parámetro de transformación${\displaystyle \theta }$se define de la siguiente manera:

Algunas manipulaciones algebraicas dan como resultado:

La expresión se simplifica al observar que${\displaystyle KC({\mathcal {S}})}$es independiente de${\displaystyle \theta }$. Además, suponiendo un registro rígido,${\displaystyle KC(T({\mathcal {M}},\theta ))}$es invariante cuando${\displaystyle \theta }$cambia porque la distancia euclidiana entre cada par de puntos permanece igual bajo una transformación rígida . Por lo tanto, la ecuación anterior se puede reescribir como:

Las estimaciones de densidad del núcleo se definen como:

${\displaystyle P_{\mathcal {M}}(x,\theta )={\frac {1}{M}}\sum _{m\in {\mathcal {M}}}K(x,T(m,\theta ))}$

${\displaystyle P_{\mathcal {S}}(x)={\frac {1}{N}}\sum _{s\in {\mathcal {S}}}K(x,s)}$

Se puede demostrar entonces que la función de coste es la correlación de las dos estimaciones de densidad del núcleo:

Una vez establecida la función de coste , el algoritmo simplemente utiliza el descenso de gradiente para encontrar la transformación óptima. Calcular la función de coste desde cero en cada iteración es computacionalmente costoso, por lo que se utiliza una versión discreta de la función de coste de la ecuación ( kc.6 ). Las estimaciones de densidad del núcleo${\displaystyle P_{\mathcal {M}},P_{\mathcal {S}}}$Se puede evaluar en puntos de la cuadrícula y almacenar en una tabla de consulta . A diferencia del ICP y métodos relacionados, no es necesario encontrar el vecino más cercano, lo que permite que el algoritmo KC sea relativamente sencillo en su implementación.

En comparación con ICP y EM-ICP para conjuntos de puntos 2D y 3D ruidosos, el algoritmo KC es menos sensible al ruido y produce un registro correcto con mayor frecuencia. ^{[ 37 ]}

La deriva coherente de puntos (CPD) fue introducida por Myronenko y Song. ^{[ 13 ] [ 41 ]} El algoritmo adopta un enfoque probabilístico para alinear conjuntos de puntos, similar al método GMM KC. A diferencia de los enfoques anteriores para el registro no rígido que asumen un modelo de transformación de spline de placa delgada , CPD es independiente del modelo de transformación utilizado. El conjunto de puntos${\displaystyle {\mathcal {M}}}$representa los centroides del modelo de mezcla gaussiana (GMM). Cuando los dos conjuntos de puntos están alineados de forma óptima, la correspondencia es el máximo de la probabilidad posterior del GMM para un punto de datos dado. Para preservar la estructura topológica de los conjuntos de puntos, se fuerza a los centroides del GMM a moverse de forma coherente como un grupo. Se utiliza el algoritmo de maximización de la esperanza para optimizar la función de coste. ^{[ 13 ]}

Sea que haya M puntos en${\displaystyle {\mathcal {M}}}$y N puntos en${\displaystyle {\mathcal {S}}}$La función de densidad de probabilidad GMM para un punto s es:

donde, en D dimensiones,${\displaystyle p(s|i)}$es la distribución gaussiana centrada en el punto${\displaystyle m_{i}\in {\mathcal {M}}}$.

${\displaystyle p(s|i)={\frac {1}{(2\pi \sigma ^{2})^{D/2}}}\exp {\left(-{\frac {\lVert s-m_{i}\rVert ^{2}}{2\sigma ^{2}}}\right)}}$

Las probabilidades de pertenencia${\displaystyle P(i)={\frac {1}{M}}}$es igual para todos los componentes GMM. El peso de la distribución uniforme se denota como${\displaystyle w\in [0,1]}$El modelo de mezcla es entonces:

Los centroides GMM se reparametrizan mediante un conjunto de parámetros.${\displaystyle \theta }$estimado maximizando la verosimilitud. Esto es equivalente a minimizar la función de log-verosimilitud negativa :

donde se supone que los datos son independientes e idénticamente distribuidos . La probabilidad de correspondencia entre dos puntos${\displaystyle m_{i}}$y${\displaystyle s_{j}}$se define como la probabilidad posterior del centroide GMM dado el punto de datos:

${\displaystyle P(i|s_{j})={\frac {P(i)p(s_{j}|i)}{p(s_{j})}}}$

El algoritmo de maximización de la esperanza (EM) se utiliza para encontrar${\displaystyle \theta }$y${\displaystyle \sigma ^{2}}$El algoritmo EM consta de dos pasos. Primero, en el paso E o paso de estimación , adivina los valores de los parámetros (valores de parámetros "antiguos") y luego utiliza el teorema de Bayes para calcular las distribuciones de probabilidad posteriores.${\displaystyle P^{\text{old}}(i,s_{j})}$de los componentes de la mezcla. En segundo lugar, en el paso M o paso de maximización , los "nuevos" valores de los parámetros se encuentran minimizando la esperanza de la función de log-verosimilitud negativa completa, es decir, la función de costo:

Ignorando constantes independientes de${\displaystyle \theta }$y${\displaystyle \sigma }$La ecuación ( cpd.4 ) se puede expresar de la siguiente manera:

dónde

${\displaystyle N_{\mathbf {P} }=\sum _{j=0}^{N}\sum _{i=0}^{M}P^{\text{old}}(i|s_{j})\leq N}$

con${\displaystyle N=N_{\mathbf {P} }}$solo si${\displaystyle w=0}$Las probabilidades posteriores de los componentes GMM calculadas utilizando valores de parámetros anteriores${\displaystyle P^{\text{old}}}$es:

Minimizar la función de costo en la ecuación ( cpd.5 ) necesariamente disminuye la función de log-verosimilitud negativa E en la ecuación ( cpd.3 ) a menos que ya se encuentre en un mínimo local. ^{[ 13 ]} Por lo tanto, el algoritmo se puede expresar utilizando el siguiente pseudocódigo, donde los conjuntos de puntos${\displaystyle {\mathcal {M}}}$y${\displaystyle {\mathcal {S}}}$están representados como${\displaystyle M\times D}$y${\displaystyle N\times D}$matrices${\displaystyle \mathbf {M} }$y${\displaystyle \mathbf {S} }$respectivamente: ^{[ 13 ]}

algoritmo CPD${\displaystyle ({\mathcal {M}},{\mathcal {S}})}$θ := θ 0 inicializar 0 ≤ w ≤ 1${\displaystyle \sigma ^{2}:={\frac {1}{DNM}}\sum _{j=1}^{N}\sum _{i=1}^{M}\lVert s_{j}-m_{i}\rVert ^{2}}$mientras no esté registrado: // Paso E, calcular P para i ∊ [1, M ] y j ∊ [1, N ]:${\displaystyle p_{ij}:={\frac {\exp \left(-{\frac {1}{2\sigma ^{2}}}\lVert s_{j}-T(m_{i},\theta )\rVert ^{2}\right)}{\sum _{k=1}^{M}\exp \left(-{\frac {1}{2\sigma ^{2}}}\lVert s_{j}-T(m_{k},\theta )\rVert ^{2}\right)+(2\pi \sigma ^{2})^{\frac {D}{2}}{\frac {w}{1-w}}{\frac {M}{N}}}}}$// Paso M, resolver para la transformación óptima { θ , σ 2 } := solve ( S , M , P ) return θ

donde el vector${\displaystyle \mathbf {1} }$es un vector columna de unos. La solvefunción difiere según el tipo de registro realizado. Por ejemplo, en el registro rígido, la salida es una escala a , una matriz de rotación.${\displaystyle \mathbf {R} }$y un vector de traslación${\displaystyle \mathbf {t} }$. El parámetro${\displaystyle \theta }$se puede escribir como una tupla de estos:

${\displaystyle \theta =\lbrace a,\mathbf {R} ,\mathbf {t} \rbrace }$

que se inicializa a uno, la matriz identidad y un vector columna de ceros:

${\displaystyle \theta _{0}=\lbrace 1,\mathbf {I} ,\mathbf {0} \rbrace }$

El conjunto de puntos alineados es:

${\displaystyle T(\mathbf {M} )=a\mathbf {M} \mathbf {R} ^{T}+\mathbf {1} \mathbf {t} ^{T}}$

La solve_rigidfunción para el registro rígido se puede escribir entonces de la siguiente manera, con la derivación del álgebra explicada en el artículo de Myronenko de 2010. ^{[ 13 ]}

solve_rigid ( S , M , P ) N P := 1 T P1${\displaystyle \mu _{s}:={\frac {1}{N_{\mathbf {P} }}}\mathbf {S} ^{T}\mathbf {P} ^{T}\mathbf {1} }$${\displaystyle \mu _{m}:={\frac {1}{N_{\mathbf {P} }}}\mathbf {M} ^{T}\mathbf {P} \mathbf {1} }$${\displaystyle {\hat {\mathbf {S} }}:=\mathbf {S} -\mathbf {1} \mu _{s}^{T}}$${\displaystyle {\hat {\mathbf {M} }}:=\mathbf {M} -\mathbf {1} \mu _{m}^{T}}$${\displaystyle \mathbf {A} :={\hat {\mathbf {S} ^{T}}}\mathbf {P} ^{T}{\hat {\mathbf {M} }}}$U , V := svd ( A ) // la descomposición en valores singulares de A = UΣV T C := diag(1, …, 1, det( UV T )) // diag( ξ ) es la matriz diagonal formada a partir del vector ξ R := UCV T${\displaystyle a:={\frac {\operatorname {tr} (\mathbf {A} ^{T}\mathbf {R} )}{\operatorname {tr} (\mathbf {{\hat {\mathbf {M} }}^{T}\operatorname {diag} (\mathbf {P} \mathbf {1} ){\hat {\mathbf {M} }}} )}}}$// tr es la traza de una matriz t := μ s − a R μ m${\displaystyle \sigma ^{2}:={\frac {1}{N_{\mathbf {P} }D}}(\operatorname {tr} (\mathbf {{\hat {\mathbf {S} }}^{T}\operatorname {diag} (\mathbf {P} ^{T}\mathbf {1} ){\hat {\mathbf {S} }}} )-a\operatorname {tr} (\mathbf {A} ^{T}\mathbf {R} ))}$devolver { a , R , t }, σ 2

Para el registro afín, donde el objetivo es encontrar una transformación afín en lugar de una rígida, la salida es una matriz de transformación afín.${\displaystyle \mathbf {B} }$y una traducción${\displaystyle \mathbf {t} }$de tal manera que el conjunto de puntos alineados sea:

${\displaystyle T(\mathbf {M} )=\mathbf {M} \mathbf {B} ^{T}+\mathbf {1} \mathbf {t} ^{T}}$

La solve_affinefunción para el registro rígido se puede escribir entonces de la siguiente manera, con la derivación del álgebra explicada en el artículo de Myronenko de 2010. ^{[ 13 ]}

resolver_afín ( S , M , P ) N P := 1 T P1${\displaystyle \mu _{s}:={\frac {1}{N_{\mathbf {P} }}}\mathbf {S} ^{T}\mathbf {P} ^{T}\mathbf {1} }$${\displaystyle \mu _{m}:={\frac {1}{N_{\mathbf {P} }}}\mathbf {M} ^{T}\mathbf {P} \mathbf {1} }$${\displaystyle {\hat {\mathbf {S} }}:=\mathbf {S} -\mathbf {1} \mu _{s}^{T}}$${\displaystyle {\hat {\mathbf {M} }}:=\mathbf {M} -\mathbf {1} \mu _{m}^{T}}$${\displaystyle \mathbf {B} :=({\hat {\mathbf {S} }}^{T}\mathbf {P} ^{T}{\hat {\mathbf {M} }})({\hat {\mathbf {M} }}^{T}\operatorname {diag} (\mathbf {P} \mathbf {1} ){\hat {\mathbf {M} }})^{-1}}$t := μ s − B μ m${\displaystyle \sigma ^{2}:={\frac {1}{N_{\mathbf {P} }D}}(\operatorname {tr} ({\hat {\mathbf {S} }}^{T}\operatorname {diag} (\mathbf {P} ^{T}\mathbf {1} ){\hat {\mathbf {S} }})-\operatorname {tr} ({\hat {\mathbf {S} }}^{T}\mathbf {P} ^{T}{\hat {\mathbf {M} }}\mathbf {B} ^{T}))}$devolver { B , t }, σ 2

También es posible utilizar CPD con registro no rígido mediante una parametrización derivada mediante cálculo de variaciones . ^{[ 13 ]}

Las sumas de distribuciones gaussianas se pueden calcular en tiempo lineal utilizando la transformada rápida de Gauss (FGT). ^{[ 13 ]} En consecuencia, la complejidad temporal de CPD es${\displaystyle O(M+N)}$, que es asintóticamente mucho más rápido que${\displaystyle O(MN)}$métodos. ^{[ 13 ]}

Este algoritmo fue introducido en 2013 por H. Assalih para facilitar el registro de imágenes de sonar. ^{[ 47 ]} Este tipo de imágenes tienden a tener altos niveles de ruido, por lo que se espera que haya muchos valores atípicos en los conjuntos de puntos a coincidir. SCS ofrece una alta robustez frente a valores atípicos y puede superar el rendimiento de ICP y CPD en presencia de valores atípicos. SCS no utiliza optimización iterativa en espacios de alta dimensión y no es ni probabilístico ni espectral. SCS puede coincidir con transformaciones rígidas y no rígidas, y funciona mejor cuando la transformación objetivo está entre tres y seis grados de libertad .