diff --git a/es-spatial-plugin/Dockerfile b/es-spatial-plugin/Dockerfile new file mode 100644 index 0000000000..bec0b0cac3 --- /dev/null +++ b/es-spatial-plugin/Dockerfile @@ -0,0 +1,21 @@ +ARG BASE_IMAGE_REPO= + +FROM ${BASE_IMAGE_REPO}docker.elastic.co/elasticsearch/elasticsearch:8.18.7 + +USER root +RUN apt update && apt upgrade -y && apt clean + +COPY target/cmr-es-spatial-plugin-0.1.0-SNAPSHOT.zip /tmp +RUN mkdir -p /var/lib/elasticsearch/tmp &&\ + chown elasticsearch:elasticsearch /var/lib/elasticsearch/tmp &&\ + chown elasticsearch:elasticsearch /usr/share/elasticsearch &&\ + bin/elasticsearch-plugin install discovery-ec2 --batch &&\ + bin/elasticsearch-plugin install file:///tmp/cmr-es-spatial-plugin-0.1.0-SNAPSHOT.zip --batch &&\ + chown -R elasticsearch:elasticsearch /usr/share/elasticsearch/plugins &&\ + echo "root soft memlock unlimited" >> /etc/security/limits.conf &&\ + echo "root hard memlock unlimited" >> /etc/security/limits.conf + +WORKDIR / + +USER 1000 +ENV HOME=/usr/share/elasticsearch diff --git a/es-spatial-plugin/project.clj b/es-spatial-plugin/project.clj index 5e928a7753..837c84237a 100644 --- a/es-spatial-plugin/project.clj +++ b/es-spatial-plugin/project.clj @@ -8,6 +8,9 @@ last (clojure.string/replace "\"" "")))) +(def plugin-jar-name + (str "target/cmr-es-spatial-plugin-" version ".jar")) + (def uberjar-name (str "target/cmr-es-spatial-plugin-" version "-standalone.jar")) @@ -23,30 +26,30 @@ (def es-deps-target-path "es-deps") -(def elastic-version "7.17.25") +(defn get-list-of-dep-jars [] + (let [all-jars (into [] (map #(.getName %) (.listFiles (clojure.java.io/file "target/lib")))) + ;; Minimal set: only what spatial-lib Java code actually needs + allowed-prefixes ["clojure" ; Runtime (spatial-lib has compiled Clojure) + "cmr-spatial-lib" ; Main library + "jafama" ; Math library used by spatial calculations + "primitive-math" ; Math optimizations + "vectorz"]] ; Vector math library + (map #(str "target/lib/" %) (filter (fn [jar-name] (some (fn [prefix] (str/starts-with? jar-name prefix)) allowed-prefixes)) all-jars)))) (defproject nasa-cmr/cmr-es-spatial-plugin "0.1.0-SNAPSHOT" :description "A Elastic Search plugin that enables spatial search entirely within elastic." :url "https://github.com/nasa/Common-Metadata-Repository/tree/master/es-spatial-plugin" + :java-source-paths ["src/java"] + :javac-options ["-target" "11" "-source" "11"] :jvm-opts ^:replace ["-server" "-Dclojure.compiler.direct-linking=true"] :plugins [[lein-shell "0.5.0"]] :profiles {:security {:plugins [[com.livingsocial/lein-dependency-check "1.4.1"]] :dependency-check {:output-format [:all] :suppression-file "resources/security/suppression.xml"}} - :provided {:dependencies [[nasa-cmr/cmr-common-lib "0.1.1-SNAPSHOT" - :exclusions [[com.fasterxml.jackson.core/jackson-core] - [com.fasterxml.jackson.dataformat/jackson-dataformat-cbor] - [com.fasterxml.jackson.dataformat/jackson-dataformat-smile] - [com.fasterxml.jackson.dataformat/jackson-dataformat-yaml]]] - [nasa-cmr/cmr-spatial-lib "0.1.0-SNAPSHOT" - :exclusions [[com.fasterxml.jackson.core/jackson-core] - [com.fasterxml.jackson.dataformat/jackson-dataformat-cbor] - [com.fasterxml.jackson.dataformat/jackson-dataformat-smile] - [com.fasterxml.jackson.dataformat/jackson-dataformat-yaml]]] - [org.elasticsearch/elasticsearch ~elastic-version] - [org.clojure/tools.reader "1.3.2"] - [org.yaml/snakeyaml "1.31"]]} + :provided {:dependencies [[nasa-cmr/cmr-common-lib "0.1.1-SNAPSHOT"] + [nasa-cmr/cmr-spatial-lib "0.1.0-SNAPSHOT"] + [org.elasticsearch/elasticsearch "8.18.7"]]} :es-deps {:dependencies [[nasa-cmr/cmr-spatial-lib "0.1.0-SNAPSHOT" ;; These exclusions will be provided by elasticsearch. :exclusions [[com.dadrox/quiet-slf4j] @@ -54,45 +57,34 @@ [com.fasterxml.jackson.dataformat/jackson-dataformat-cbor] [com.fasterxml.jackson.dataformat/jackson-dataformat-smile] [com.fasterxml.jackson.dataformat/jackson-dataformat-yaml] + [commons-io] [commons-codec] [commons-logging] [joda-time] [org.ow2.asm/asm] [org.ow2.asm/asm-all] - ;; Both lz4 libraries are linked together. yawk - ;; is supposed to be a drop in replacement and - ;; uses the same package name as the original. [net.jpountz.lz4/lz4] - [at.yawk.lz4/lz4-java] [org.locationtech.jts/jts-core] [org.locationtech.jts.JTSVersion] [org.slf4j/slf4j-api]]] - [org.clojure/tools.reader "1.3.2"] + [org.clojure/tools.reader "1.5.0"] [org.clojure/clojure "1.11.2"]] :target-path ~es-deps-target-path :uberjar-name ~es-deps-uberjar-name :jar-name ~es-deps-jar-name :aot []} - :es-plugin {:aot [cmr.elasticsearch.plugins.spatial.script.core - cmr.elasticsearch.plugins.spatial.factory.lfactory - cmr.elasticsearch.plugins.spatial.factory.core - cmr.elasticsearch.plugins.spatial.engine.core - cmr.elasticsearch.plugins.spatial.plugin]} + :jar-deps {:plugins [[org.clojars.jj/copy-deps "1.0.1"]] + :dependencies [[nasa-cmr/cmr-spatial-lib "0.1.0-SNAPSHOT"] + [nasa-cmr/cmr-common-lib "0.1.1-SNAPSHOT"]] + :aot []} :dev {:dependencies [[criterium "0.4.4"] - [cheshire "5.12.0"] - [org.clojure/tools.reader "1.3.2"] + [cheshire "5.13.0"] [nasa-cmr/cmr-common-lib "0.1.1-SNAPSHOT"] [nasa-cmr/cmr-spatial-lib "0.1.0-SNAPSHOT"] - [org.elasticsearch/elasticsearch ~elastic-version] + [org.elasticsearch/elasticsearch "8.18.7"] [org.clojars.gjahad/debug-repl "0.3.3"] - [org.clojure/tools.nrepl "0.2.13"] - [org.clojure/tools.namespace "0.2.11"] - [org.yaml/snakeyaml "1.31"]] - :aot [cmr.elasticsearch.plugins.spatial.script.core - cmr.elasticsearch.plugins.spatial.factory.lfactory - cmr.elasticsearch.plugins.spatial.factory.core - cmr.elasticsearch.plugins.spatial.engine.core - cmr.elasticsearch.plugins.spatial.plugin] + [nrepl/nrepl "1.3.0"] + [org.clojure/tools.namespace "1.2.0"]] :global-vars {*warn-on-reflection* false *assert* false}} :static {} @@ -108,28 +100,27 @@ :kaocha {:dependencies [[lambdaisland/kaocha "1.0.732"] [lambdaisland/kaocha-cloverage "1.0.75"] [lambdaisland/kaocha-junit-xml "0.0.76"]]}} - :aliases {"install-es-deps" ["do" - "with-profile" "es-deps,provided" "clean," - "with-profile" "es-deps,provided" "uberjar," - ;; target-path is being ignored for uberjar. move uberjar to es-deps-target-path. - ["shell" "echo" "inst-es-deps"] - "shell" "mv" ~(str "target/" es-deps-uberjar-name) ~es-deps-target-path] - "install-es-plugin" ["do" - ["shell" "echo" "inst-es-plugin"] - "with-profile" "es-plugin,provided" "clean," - "with-profile" "es-plugin,provided" "uberjar,"] - "package-es-plugin" ["do" + :aliases {"install-es-plugin" ["do" + ["shell" "echo" "Building ES spatial plugin JAR"] + "with-profile" "provided" "clean," + "with-profile" "provided" "jar,"] + "gather-dependencies" ["do" + ["shell" "echo" "Collecting dependent JARs"] + "with-profile" "jar-deps" "copy-deps,"] + "prepare-es-plugin" ["do" "install-es-plugin" - ["shell" "echo" "pack-es-deps"] - "shell" - "zip" - "-j" - ~plugin-zip-name - ~uberjar-name - "resources/plugin/plugin-descriptor.properties"] + "gather-dependencies"] + "package-es-plugin" ~(vec (concat ["do" + ["shell" "echo" "Packaging ES plugin into zip file"] + "shell" + "zip" + "-j" + plugin-zip-name + plugin-jar-name + "resources/plugin/plugin-descriptor.properties"] + (get-list-of-dep-jars))) "build-all" ["do" ["shell" "echo" "build-all"] - "install-es-deps," "install-es-plugin,"] ;; Kaocha test aliases diff --git a/es-spatial-plugin/resources/plugin/plugin-descriptor.properties b/es-spatial-plugin/resources/plugin/plugin-descriptor.properties index 2c9fd2baa0..51d57d8206 100644 --- a/es-spatial-plugin/resources/plugin/plugin-descriptor.properties +++ b/es-spatial-plugin/resources/plugin/plugin-descriptor.properties @@ -13,7 +13,7 @@ classname=cmr.elasticsearch.plugins.SpatialSearchPlugin # use the system property java.specification.version # version string must be a sequence of nonnegative decimal integers # separated by "."'s and may have leading zeros -java.version=1.8 +java.version=11 # # 'elasticsearch.version': version of elasticsearch compiled against elasticsearch.version=7.17.14 diff --git a/es-spatial-plugin/resources/plugin/plugin-security.policy b/es-spatial-plugin/resources/plugin/plugin-security.policy deleted file mode 100644 index 8a04dfe814..0000000000 --- a/es-spatial-plugin/resources/plugin/plugin-security.policy +++ /dev/null @@ -1,6 +0,0 @@ -grant { - permission java.util.PropertyPermission "*", "read,write"; - permission java.lang.RuntimePermission "getClassLoader"; - permission java.lang.RuntimePermission "accessClassInPackage.sun.misc"; - permission java.io.FilePermission ".lein-env", "read"; - }; diff --git a/es-spatial-plugin/src/cmr/elasticsearch/plugins/spatial/engine/core.clj b/es-spatial-plugin/src/cmr/elasticsearch/plugins/spatial/engine/core.clj deleted file mode 100644 index dd7e004621..0000000000 --- a/es-spatial-plugin/src/cmr/elasticsearch/plugins/spatial/engine/core.clj +++ /dev/null @@ -1,41 +0,0 @@ -(ns cmr.elasticsearch.plugins.spatial.engine.core - (:import - (cmr.elasticsearch.plugins SpatialScriptFactory) - (org.elasticsearch.script FilterScript ScriptContext ScriptEngine) - (java.util Map)) - (:gen-class - :name cmr.elasticsearch.plugins.SpatialScriptEngine - :implements [org.elasticsearch.script.ScriptEngine])) - -(import 'cmr.elasticsearch.plugins.SpatialScriptEngine) - -(defn -getType - "Get script lang." - [^SpatialScriptEngine _this] - "cmr_spatial") - -(defn -compile - "Compile script." - [^SpatialScriptEngine this - ^String _script-name - ^String script-source - ^ScriptContext context - ^Map _params] - (cond - (not (.equals context FilterScript/CONTEXT)) - (throw (new IllegalArgumentException - (format "%s scripts cannot be used for context [%s]" - (.getType this) (.name context)))) - - (not (.equals "spatial" script-source)) - (throw (new IllegalArgumentException - (format "Unknown script name %s" script-source))) - - :else - (-> context .factoryClazz (.cast (new SpatialScriptFactory))))) - -#_{:clj-kondo/ignore [:redundant-do]} -(defn -close - "Do nothing." - [_this] - (do)) diff --git a/es-spatial-plugin/src/cmr/elasticsearch/plugins/spatial/factory/core.clj b/es-spatial-plugin/src/cmr/elasticsearch/plugins/spatial/factory/core.clj deleted file mode 100644 index 61e4ba6521..0000000000 --- a/es-spatial-plugin/src/cmr/elasticsearch/plugins/spatial/factory/core.clj +++ /dev/null @@ -1,22 +0,0 @@ -(ns cmr.elasticsearch.plugins.spatial.factory.core - (:import - (cmr.elasticsearch.plugins SpatialScriptLeafFactory) - (java.util Map) - (org.elasticsearch.search.lookup SearchLookup)) - (:gen-class - :name cmr.elasticsearch.plugins.SpatialScriptFactory - :implements [org.elasticsearch.script.FilterScript$Factory - org.elasticsearch.script.FilterScript$LeafFactory] - :state data)) - -(import 'cmr.elasticsearch.plugins.SpatialScriptFactory) - -(defn -isResultDeterministic - "Implies the results are cacheable if true. - See [[https://www.elastic.co/guide/en/elasticsearch/reference/current/modules-scripting-engine.html]]" - [^SpatialScriptFactory _this] - false) - -(defn -newFactory - [^SpatialScriptFactory _this ^Map params ^SearchLookup lookup] - (new SpatialScriptLeafFactory params lookup)) diff --git a/es-spatial-plugin/src/cmr/elasticsearch/plugins/spatial/factory/lfactory.clj b/es-spatial-plugin/src/cmr/elasticsearch/plugins/spatial/factory/lfactory.clj deleted file mode 100644 index 46541f4e80..0000000000 --- a/es-spatial-plugin/src/cmr/elasticsearch/plugins/spatial/factory/lfactory.clj +++ /dev/null @@ -1,90 +0,0 @@ -(ns cmr.elasticsearch.plugins.spatial.factory.lfactory - (:import - (cmr.elasticsearch.plugins SpatialScript) - (java.util Map) - (org.apache.logging.log4j LogManager) - (org.elasticsearch.script DocReader) - (org.elasticsearch.common.xcontent.support XContentMapValues) - (org.elasticsearch.search.lookup SearchLookup)) - (:require - [clojure.string :as string] - [cmr.spatial.relations :as relations] - [cmr.spatial.serialize :as srl]) - (:gen-class - :name cmr.elasticsearch.plugins.SpatialScriptLeafFactory - :implements [org.elasticsearch.script.FilterScript$LeafFactory] - :constructors {[java.util.Map org.elasticsearch.search.lookup.SearchLookup] []} - :init init - :state data)) - -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -;; Begin factory helper functions ;; -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; - -(def parameters - "The parameters to the Spatial script" - [:ords :ords-info]) - -(defn- extract-params - "Extracts the parameters from the params map given in the script." - [script-params] - (when script-params - (into {} (for [param parameters] - [param (XContentMapValues/nodeStringValue - (get script-params (name param)) nil)])))) - -(defn- assert-required-parameters - "Asserts that all the parameters are supplied or it throws an exception." - [params] - (when-not (every? params parameters) - (throw (IllegalArgumentException. - (str "Missing one or more of required parameters: " - (clojure.string/join parameters ", ")))))) - -(defn- convert-params - "Convert the comma separated string into a vector of integers." - [[k v]] - [k (map #(Integer. ^String %) (string/split v #","))]) - -(defn- params->spatial-shape - [params] - (let [{:keys [ords-info ords]} (->> params - (map convert-params) - (into {}))] - (first (srl/ords-info->shapes ords-info ords)))) - -(defn get-intersects-fn [script-params] - (let [params (extract-params script-params) - shape (params->spatial-shape params)] - (assert-required-parameters params) - (try - (relations/shape->intersects-fn shape) - (catch Exception e - (.error (LogManager/getLogger "cmr_spatial_lfactory") (format "Unable to create intersects function for shapes [%s]" (pr-str shape)) e) - (throw (ex-info "An exception occurred creating intersects-fn" {:shape shape :params params} e)))))) - -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -;; End factory helper functions ;; -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; - -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -;; Begin leaf factory functions ;; -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; - -(import 'cmr.elasticsearch.plugins.SpatialScriptLeafFactory) - -(defn- -init [^Map params ^SearchLookup lookup] - [[] {:params params - :lookup lookup}]) - -(defn -newInstance [^SpatialScriptLeafFactory this ^DocReader doc-reader] - (let [^Map params (-> this .data :params)] - (SpatialScript. - ^Object (get-intersects-fn params) - ^Map params - ^SearchLookup (-> this .data :lookup) - ^DocReader doc-reader))) - -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -;; End leaf factory functions ;; -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; diff --git a/es-spatial-plugin/src/cmr/elasticsearch/plugins/spatial/plugin.clj b/es-spatial-plugin/src/cmr/elasticsearch/plugins/spatial/plugin.clj deleted file mode 100644 index 579a603511..0000000000 --- a/es-spatial-plugin/src/cmr/elasticsearch/plugins/spatial/plugin.clj +++ /dev/null @@ -1,22 +0,0 @@ -(ns cmr.elasticsearch.plugins.spatial.plugin - (:import - (org.elasticsearch.script ScriptModule) - (org.elasticsearch.plugins ScriptPlugin) - (org.elasticsearch.common.settings Settings) - (org.elasticsearch.plugins Plugin) - (cmr.elasticsearch.plugins SpatialScriptEngine) - (java.util Collection Collections)) - (:gen-class - :name cmr.elasticsearch.plugins.SpatialSearchPlugin - :extends org.elasticsearch.plugins.Plugin - :implements [org.elasticsearch.plugins.ScriptPlugin])) - -(defn -getScriptEngine - "Spatial script engine." - [_this ^Settings _settings ^Collection _contexts] - (new SpatialScriptEngine)) - -(defn -getContexts - "Return script contexts." - [_this] - (Collections/emptyList)) diff --git a/es-spatial-plugin/src/cmr/elasticsearch/plugins/spatial/script/core.clj b/es-spatial-plugin/src/cmr/elasticsearch/plugins/spatial/script/core.clj deleted file mode 100644 index d16d2bb2c5..0000000000 --- a/es-spatial-plugin/src/cmr/elasticsearch/plugins/spatial/script/core.clj +++ /dev/null @@ -1,88 +0,0 @@ -(ns cmr.elasticsearch.plugins.spatial.script.core - (:require - [cmr.common.util :as u] - [cmr.spatial.serialize :as srl]) - (:import - (java.util Map) - (org.apache.logging.log4j LogManager) - (org.elasticsearch.script DocReader) - (org.elasticsearch.search.lookup FieldLookup - LeafDocLookup - LeafStoredFieldsLookup - LeafSearchLookup - SearchLookup)) - (:gen-class - :name cmr.elasticsearch.plugins.SpatialScript - :extends org.elasticsearch.script.FilterScript - :constructors {[java.lang.Object - java.util.Map - org.elasticsearch.search.lookup.SearchLookup - org.elasticsearch.script.DocReader] - [java.util.Map - org.elasticsearch.search.lookup.SearchLookup - org.elasticsearch.script.DocReader]} - :methods [[getFields [] org.elasticsearch.search.lookup.LeafStoredFieldsLookup]] - :init init - :state data)) - -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -;; Begin script helper functions ;; -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; - -(defn- get-from-fields - [^LeafStoredFieldsLookup lookup key] - (when (and lookup key (.containsKey lookup key)) - (when-let [^FieldLookup field-lookup (.get lookup key)] - (seq (.getValues field-lookup))))) - -(defn doc-intersects? - "Returns true if the doc contains a ring that intersects the ring passed in." - [^LeafStoredFieldsLookup lookup intersects-fn] - - (if-let [ords-info (get-from-fields lookup "ords-info")] - (let [ords (get-from-fields lookup "ords") - shapes (srl/ords-info->shapes ords-info ords)] - (try - ;; Must explicitly return true or false or elastic search will complain - (if (u/any-true? intersects-fn shapes) - true - false) - (catch Throwable t - (.error (LogManager/getLogger "cmr_spatial_script") t) - (throw (ex-info "An exception occurred checking for intersections" {:shapes shapes} t))))) - false)) - -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -;; End script helper functions ;; -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; - -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -;; Begin script functions ;; -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; - -(import 'cmr.elasticsearch.plugins.SpatialScript) - -(defn ^LeafStoredFieldsLookup -getFields - [^SpatialScript this] - (-> this .data :search-lookup .fields)) - -(defn ^Map -getDoc - [^SpatialScript this] - (-> this .data :search-lookup .doc)) - -;; Need to override setDocument for more control over lookup -(defn -setDocument - [^SpatialScript this doc-id] - (-> this .data :search-lookup (.setDocument doc-id))) - -(defn- -init [^Object intersects-fn ^Map params ^SearchLookup lookup ^DocReader doc-reader] - [[params lookup doc-reader] {:intersects-fn intersects-fn - :search-lookup (.getLeafSearchLookup lookup (.getLeafReaderContext doc-reader))}]) - -(defn -execute [^SpatialScript this] - (doc-intersects? (.getFields this) - (-> this .data :intersects-fn))) - -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; -;; End script functions ;; -;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; diff --git a/es-spatial-plugin/src/java/cmr/elasticsearch/plugins/SpatialSearchPlugin.java b/es-spatial-plugin/src/java/cmr/elasticsearch/plugins/SpatialSearchPlugin.java new file mode 100644 index 0000000000..68c144550a --- /dev/null +++ b/es-spatial-plugin/src/java/cmr/elasticsearch/plugins/SpatialSearchPlugin.java @@ -0,0 +1,267 @@ +package cmr.elasticsearch.plugins; + +import org.elasticsearch.plugins.Plugin; +import org.elasticsearch.plugins.ScriptPlugin; +import org.elasticsearch.script.DocReader; +import org.elasticsearch.script.ScriptEngine; +import org.elasticsearch.script.ScriptContext; +import org.elasticsearch.script.FilterScript; +import org.elasticsearch.common.settings.Settings; +import org.elasticsearch.search.lookup.SearchLookup; +import org.elasticsearch.search.lookup.LeafSearchLookup; +import org.elasticsearch.search.lookup.LeafStoredFieldsLookup; +import org.elasticsearch.common.xcontent.support.XContentMapValues; +import org.apache.logging.log4j.LogManager; +import org.apache.logging.log4j.Logger; + +import cmr.spatial.serialize.OrdsInfoShapes; +import cmr.spatial.relations.ShapeIntersections; +import cmr.spatial.relations.ShapePredicate; +import cmr.spatial.shape.SpatialShape; + +import java.io.IOException; +import java.util.Collection; +import java.util.List; +import java.util.Map; +import java.util.Set; +import java.util.stream.Collectors; + +/** + * Pure Java Elasticsearch plugin for CMR spatial searches. + * Uses the Java spatial library to avoid Clojure runtime security issues in ES8+. + */ +public class SpatialSearchPlugin extends Plugin implements ScriptPlugin { + + public SpatialSearchPlugin(Settings settings) {} + + @Override + public ScriptEngine getScriptEngine(Settings settings, Collection> contexts) { + return new SpatialScriptEngine(); + } + + /** + * Script engine implementation for spatial searches. + */ + private static class SpatialScriptEngine implements ScriptEngine { + private static final Logger logger = LogManager.getLogger(SpatialScriptEngine.class); + + @Override + public String getType() { + return "cmr_spatial"; + } + + @Override + public T compile(String scriptName, String source, ScriptContext context, Map params) { + if (!context.instanceClazz.equals(FilterScript.class)) { + throw new IllegalArgumentException( + getType() + " scripts cannot be used for context [" + context.name + "]"); + } + + if (!source.equals("spatial")) { + throw new IllegalArgumentException("Unknown script name [" + source + "]"); + } + + return context.factoryClazz.cast(new SpatialScriptFactory()); + } + + @Override + public Set> getSupportedContexts() { + return Set.of(FilterScript.CONTEXT); + } + } + + /** + * Factory for creating leaf factories with parsed spatial parameters. + */ + private static class SpatialScriptFactory implements FilterScript.Factory { + private static final Logger logger = LogManager.getLogger(SpatialScriptFactory.class); + + @Override + public boolean isResultDeterministic() { + return false; + } + + @Override + public FilterScript.LeafFactory newFactory(Map params, SearchLookup lookup) { + // Extract and validate parameters + String ordsParam = XContentMapValues.nodeStringValue(params.get("ords"), null); + String ordsInfoParam = XContentMapValues.nodeStringValue(params.get("ords-info"), null); + + if (ordsParam == null || ordsInfoParam == null) { + throw new IllegalArgumentException( + "Missing required parameters: ords and ords-info must be provided"); + } + + // Parse comma-separated integers + List ords = parseIntList(ordsParam); + List ordsInfo = parseIntList(ordsInfoParam); + + // Convert to spatial shape and create intersects predicate + try { + List shapes = OrdsInfoShapes.ordsInfoToShapes(ordsInfo, ords); + if (shapes.isEmpty()) { + throw new IllegalArgumentException("No shapes could be parsed from ords-info/ords"); + } + + SpatialShape queryShape = shapes.get(0); + ShapePredicate intersectsFn = ShapeIntersections.createIntersectsFn(queryShape); + + return new SpatialLeafFactory(intersectsFn, lookup); + + } catch (Exception e) { + logger.error("Failed to create intersects function from shape parameters", e); + throw new IllegalArgumentException("Unable to create spatial intersects function", e); + } + } + + private static List parseIntList(String csvString) { + String[] parts = csvString.split(","); + return java.util.Arrays.stream(parts) + .map(String::trim) + .map(Integer::parseInt) + .collect(Collectors.toList()); + } + } + + /** + * Leaf factory that creates FilterScript instances per segment. + */ + private static class SpatialLeafFactory implements FilterScript.LeafFactory { + private final ShapePredicate intersectsFn; + private final SearchLookup lookup; + + public SpatialLeafFactory(ShapePredicate intersectsFn, SearchLookup lookup) { + this.intersectsFn = intersectsFn; + this.lookup = lookup; + } + + @Override + public FilterScript newInstance(DocReader reader) throws IOException { + return new SpatialScript(intersectsFn, lookup, reader); + } + } + + /** + * FilterScript that executes spatial intersection tests on documents. + */ + private static class SpatialScript extends FilterScript { + private static final Logger logger = LogManager.getLogger(SpatialScript.class); + + private final ShapePredicate intersectsFn; + private final SearchLookup searchLookup; + private LeafSearchLookup leafLookup; + + public SpatialScript(ShapePredicate intersectsFn, SearchLookup lookup, DocReader reader) { + super(Map.of(), lookup, reader); + this.intersectsFn = intersectsFn; + this.searchLookup = lookup; + // Try to get leaf lookup - DocReader should be associated with a LeafReaderContext + try { + // In ES 8, DocReader provides the context + Object ctx = reader.getClass().getMethod("getLeafReaderContext").invoke(reader); + if (ctx instanceof org.apache.lucene.index.LeafReaderContext) { + this.leafLookup = lookup.getLeafSearchLookup((org.apache.lucene.index.LeafReaderContext) ctx); + } + } catch (Exception e) { + logger.warn("Could not get leaf lookup from reader", e); + } + } + + @Override + public void setDocument(int docId) { + super.setDocument(docId); + // Update the leaf lookup for this document + if (leafLookup != null) { + leafLookup.setDocument(docId); + } + } + + @Override + public boolean execute() { + try { + if (leafLookup == null) { + logger.warn("Leaf lookup not initialized"); + return false; + } + + // Access fields through the leaf lookup + LeafStoredFieldsLookup fields = leafLookup.fields(); + + if (fields == null) { + logger.debug("Could not access fields lookup"); + return false; + } + + // Extract ords-info and ords from document fields + List ordsInfo = getFieldValues(fields, "ords-info"); + List ords = getFieldValues(fields, "ords"); + + if (ordsInfo == null || ords == null) { + return false; + } + + // Deserialize shapes from document + List docShapes = OrdsInfoShapes.ordsInfoToShapes(ordsInfo, ords); + + // Test if any document shape intersects the query shape + for (SpatialShape docShape : docShapes) { + if (intersectsFn.intersects(docShape)) { + return true; + } + } + + return false; + + } catch (Exception e) { + logger.error("Error executing spatial script", e); + throw new RuntimeException("Spatial intersection test failed", e); + } + } + + private LeafStoredFieldsLookup getFieldsLookup() { + try { + // Access the SearchLookup to get fields + if (searchLookup != null) { + // ES 7/8 API: searchLookup should give us access to leaf lookup + java.lang.reflect.Method getLeafSearchLookup = + searchLookup.getClass().getMethod("getLeafSearchLookup", org.apache.lucene.index.LeafReaderContext.class); + // Get the reader context from our DocReader (which wraps it) + // For now, try to get it from the FilterScript base + return null; // Will fix this properly + } + } catch (Exception e) { + logger.warn("Could not access leaf lookup", e); + } + return null; + } + + @SuppressWarnings("unchecked") + private List getFieldValues(LeafStoredFieldsLookup fields, String key) { + if (fields == null || key == null) { + return null; + } + + try { + if (!fields.containsKey(key)) { + return null; + } + + Object fieldLookup = fields.get(key); + if (fieldLookup == null) { + return null; + } + + // FieldLookup should have a getValues() method + java.lang.reflect.Method getValues = fieldLookup.getClass().getMethod("getValues"); + Object values = getValues.invoke(fieldLookup); + if (values instanceof List) { + return (List) values; + } + } catch (Exception e) { + logger.debug("Could not access field values for key: " + key, e); + } + + return null; + } + } +} diff --git a/spatial-lib/project.clj b/spatial-lib/project.clj index d7998479ec..f520278abc 100644 --- a/spatial-lib/project.clj +++ b/spatial-lib/project.clj @@ -1,6 +1,7 @@ (defproject nasa-cmr/cmr-spatial-lib "0.1.0-SNAPSHOT" :description "A spatial library for the CMR." :url "https://github.com/nasa/Common-Metadata-Repository/tree/master/spatial-lib" + :java-source-paths ["src/java"] :dependencies [[nasa-cmr/cmr-common-lib "0.1.1-SNAPSHOT"] [net.jafama/jafama "2.3.1"] [net.mikera/core.matrix "0.54.0"] diff --git a/spatial-lib/src/cmr/spatial/arc.clj b/spatial-lib/src/cmr/spatial/arc.clj index c1e0cfde15..cf893aff17 100644 --- a/spatial-lib/src/cmr/spatial/arc.clj +++ b/spatial-lib/src/cmr/spatial/arc.clj @@ -467,5 +467,5 @@ (extend-protocol d/DerivedCalculator - cmr.spatial.arc.Arc + cmr.spatial.internal.arc.Arc (calculate-derived ^Arc [^Arc a] a)) diff --git a/spatial-lib/src/cmr/spatial/arc_line_segment_intersections.clj b/spatial-lib/src/cmr/spatial/arc_line_segment_intersections.clj index 55841d515a..e54b66fa47 100644 --- a/spatial-lib/src/cmr/spatial/arc_line_segment_intersections.clj +++ b/spatial-lib/src/cmr/spatial/arc_line_segment_intersections.clj @@ -7,8 +7,7 @@ [cmr.common.util :as u] [cmr.spatial.line-segment :as s] [primitive-math]) - (:import cmr.spatial.arc.Arc - cmr.spatial.line_segment.LineSegment + (:import cmr.spatial.line_segment.LineSegment cmr.spatial.point.Point)) (primitive-math/use-primitive-operators) @@ -69,11 +68,11 @@ (defn line-segment-arc-intersections "Returns a list of the points where the line segment intersects the arc." - [^LineSegment ls ^Arc arc] + [^LineSegment ls arc] (let [ls-mbr (.mbr ls) - arc-mbr1 (.mbr1 arc) - arc-mbr2 (.mbr2 arc) + arc-mbr1 (:mbr1 arc) + arc-mbr2 (:mbr2 arc) mbr1-intersects (m/intersects-br? :geodetic ls-mbr arc-mbr1) mbr2-intersects (and arc-mbr2 (m/intersects-br? :geodetic ls-mbr arc-mbr2)) arc-points (a/arc->points arc) @@ -132,7 +131,7 @@ "Determines if line 1 and 2 intersect. A line can be an arc or a line segment." [line1 line2] - (if (= (type line2) Arc) + (if (= (type line2) cmr.spatial.arc.Arc) (intersections-with-arc line1 line2) (intersections-with-line-segment line1 line2))) @@ -149,7 +148,8 @@ [ls point] (s/point-on-segment? ls point)) - Arc + ;; Clojure Arc record type + cmr.spatial.arc.Arc (intersections-with-arc [arc1 arc2] (a/intersections arc1 arc2)) diff --git a/spatial-lib/src/cmr/spatial/cartesian_ring.clj b/spatial-lib/src/cmr/spatial/cartesian_ring.clj index 9d8bbf9b6a..36423c7e12 100644 --- a/spatial-lib/src/cmr/spatial/cartesian_ring.clj +++ b/spatial-lib/src/cmr/spatial/cartesian_ring.clj @@ -6,7 +6,8 @@ [cmr.spatial.line-segment :as s] [cmr.spatial.derived :as d] [cmr.common.dev.record-pretty-printer :as record-pretty-printer]) - (:import cmr.spatial.line_segment.LineSegment)) + (:import + cmr.spatial.line_segment.LineSegment)) (primitive-math/use-primitive-operators) (def external-point @@ -30,7 +31,10 @@ line-segments ;; the minimum bounding rectangle - mbr]) + mbr + + ;; Cached Java ring for intersection operations + java-ring]) (record-pretty-printer/enable-record-pretty-printing CartesianRing) @@ -49,29 +53,19 @@ lines))) (defn covers-point? - "Determines if a ring covers the given point. The algorithm works by counting the number of times - an arc between the point and a known external point crosses the ring. An even count means the point - is external. An odd count means the point is inside the ring." - [^CartesianRing ring point] - - ;; Only do real intersection if the mbr covers the point. - (when (mbr/cartesian-covers-point? (.mbr ring) point) - (or - ;; The point is actually one of the rings points - (contains? (.point_set ring) point) - ;; otherwise we'll do the real intersection algorithm - (let [;; Create the test segment - crossing-line (s/line-segment point external-point) - intersections (lines-and-line-intersections (.line_segments ring) crossing-line)] - (or (odd? (count intersections)) - ;; if the point itself is one of the intersections then the ring covers it - (intersections point)))))) + "Determines if a ring covers the given point." + [^cmr.spatial.cartesian_ring.CartesianRing ring ^cmr.spatial.point.Point point] + ;; Use cached Java ring if available, otherwise create on-the-fly + (let [java-ring (or (.java_ring ring) + (cmr.spatial.internal.ring.CartesianRing/createRing (vec (.points ring)))) + java-point (cmr.spatial.shape.Point. (.lon point) (.lat point))] + (.coversPoint java-ring java-point))) (defn ring "Creates a new ring with the given points. If the other fields of a ring are needed. The calculate-derived function should be used to populate it." [points] - (->CartesianRing (mapv p/with-cartesian-equality points) nil nil nil)) + (->CartesianRing (mapv p/with-cartesian-equality points) nil nil nil nil)) (defn ring->line-segments "Determines the line-segments from the points in the ring." @@ -111,5 +105,6 @@ (if (.line_segments ring) ring (let [^CartesianRing ring (assoc ring :point-set (set (.points ring))) - ^CartesianRing ring (assoc ring :line-segments (ring->line-segments ring))] - (assoc ring :mbr (ring->mbr ring)))))) + ^CartesianRing ring (assoc ring :line-segments (ring->line-segments ring)) + ^CartesianRing ring (assoc ring :mbr (ring->mbr ring))] + (assoc ring :java-ring (cmr.spatial.internal.ring.CartesianRing/createRing (vec (.points ring)))))))) diff --git a/spatial-lib/src/cmr/spatial/geodetic_ring.clj b/spatial-lib/src/cmr/spatial/geodetic_ring.clj index ea4d931834..421d923fb0 100644 --- a/spatial-lib/src/cmr/spatial/geodetic_ring.clj +++ b/spatial-lib/src/cmr/spatial/geodetic_ring.clj @@ -7,7 +7,7 @@ [cmr.spatial.mbr :as mbr] [cmr.spatial.point :as p] [primitive-math]) - (:import cmr.spatial.arc.Arc)) + (:import)) (primitive-math/use-primitive-operators) @@ -44,7 +44,9 @@ ;; Three points that are not within the ring. These are used to test if a point is inside or ;; outside a ring. We generate multiple external points so that we have a backup if one external ;; point is antipodal to a point we're checking is inside a ring. - external-points]) + external-points + ;; Cached Java ring for intersection operations + java-ring]) (record-pretty-printer/enable-record-pretty-printing GeodeticRing) @@ -81,27 +83,13 @@ " point. Ring: " (pr-str ring) " point: " (pr-str point))))))) (defn covers-point? - "Determines if a ring covers the given point. The algorithm works by counting the number of times - an arc between the point and a known external point crosses the ring. An even count means the point - is external. An odd count means the point is inside the ring." - [ring point] - ;; The pre check is necessary for rings which might contain both north and south poles - {:pre [(> (count (:external-points ring)) 0)]} - - (or (and (:contains-north-pole ring) (p/is-north-pole? point)) - (and (:contains-south-pole ring) (p/is-south-pole? point)) - ;; Only do real intersection if the mbr covers the point. - (when (mbr/geodetic-covers-point? (:mbr ring) point) - (if ((:point-set ring) point) - true ; The point is actually one of the rings points - ;; otherwise we'll do the real intersection algorithm - (let [external-point (choose-external-point ring point) - ;; Create the test arc - crossing-arc (a/arc point external-point) - intersections (arcs-and-arc-intersections (:arcs ring) crossing-arc)] - (or (odd-long? (count intersections)) - ;; if the point itself is one of the intersections then the ring covers it - (intersections point))))))) + "Determines if a ring covers the given point." + [ring ^cmr.spatial.point.Point point] + ;; Use cached Java ring if available, otherwise create on-the-fly + (let [java-ring (or (:java-ring ring) + (cmr.spatial.internal.ring.GeodeticRing/createRing (vec (:points ring)))) + java-point (cmr.spatial.shape.Point. (.lon point) (.lat point))] + (.coversPoint java-ring java-point))) (defn- arcs->course-rotation-direction "Calculates the rotation direction of the arcs of a ring. Will be one of :clockwise, @@ -121,10 +109,10 @@ (let [courses (loop [courses (transient []) arcs arcs] (if (empty? arcs) (persistent! courses) - (let [^Arc a (first arcs)] + (let [a (first arcs)] (recur (-> courses - (conj! (.initial_course a)) - (conj! (.ending_course a))) + (conj! (:initial-course a)) + (conj! (:ending-course a))) (rest arcs))))) ;; Add the first turn angle on again to complete the turn courses (conj courses (first courses))] @@ -134,7 +122,7 @@ "Creates a new ring with the given points. If the other fields of a ring are needed. The calculate-derived function should be used to populate it." [points] - (->GeodeticRing (mapv p/with-geodetic-equality points) nil nil nil nil nil nil nil)) + (->GeodeticRing (mapv p/with-geodetic-equality points) nil nil nil nil nil nil nil nil)) (defn contains-both-poles? "Returns true if a ring contains both the north pole and the south pole" @@ -190,12 +178,12 @@ (or (.mbr ring) (let [arcs (ring->arcs ring) {:keys [contains-north-pole contains-south-pole]} (ring->pole-containment ring) - br (reduce (fn [br ^Arc arc] - (let [mbr1 (.mbr1 arc) + br (reduce (fn [br arc] + (let [mbr1 (:mbr1 arc) br (if br (mbr/union br mbr1) mbr1) - mbr2 (.mbr2 arc)] + mbr2 (:mbr2 arc)] (if mbr2 (mbr/union br mbr2) br))) @@ -235,4 +223,5 @@ (assoc ring :arcs (ring->arcs ring)) (ring->pole-containment ring) (assoc ring :mbr (ring->mbr ring)) - (assoc ring :external-points (ring->external-points ring)))))) + (assoc ring :external-points (ring->external-points ring)) + (assoc ring :java-ring (cmr.spatial.internal.ring.GeodeticRing/createRing (vec (:points ring)))))))) diff --git a/spatial-lib/src/cmr/spatial/line_string.clj b/spatial-lib/src/cmr/spatial/line_string.clj index 995052b760..7b5e02217f 100644 --- a/spatial-lib/src/cmr/spatial/line_string.clj +++ b/spatial-lib/src/cmr/spatial/line_string.clj @@ -18,7 +18,8 @@ (:import (cmr.spatial.arc Arc) (cmr.spatial.line_segment LineSegment) - (cmr.spatial.mbr Mbr))) + (cmr.spatial.mbr Mbr) + (cmr.spatial.geometry LineStringIntersections))) (primitive-math/use-primitive-operators) @@ -151,44 +152,34 @@ (defn covers-point? "Returns true if the line covers the point" - [^LineString line point] - (let [point-set (.point_set line) - segments (.segments line)] - (or (contains? point-set point) - (u/any-true? #(point-on-segment? % point) segments)))) + [^LineString line ^cmr.spatial.point.Point point] + ;; Delegate to Java implementation + (let [java-line (cmr.spatial.shape.LineString. + (name (.coordinate_system line)) + (vec (p/points->ords (.points line)))) + java-point (cmr.spatial.shape.Point. (.lon point) (.lat point))] + (LineStringIntersections/coversPoint java-line java-point))) (defn intersects-br? "Returns true if the line intersects the br" [^LineString line ^Mbr br] - (when (m/intersects-br? (.coordinate_system line) (.mbr line) br) - (if (m/single-point? br) - (covers-point? line (p/point (.west br) (.north br))) - - (let [coord-sys (.coordinate_system line)] - (or - ;; Does the br cover any points of the line? - (u/any-true? #(m/covers-point? coord-sys br %) (.points line)) - ;; Does the line contain any points of the br? - (u/any-true? #(covers-point? line %) (m/corner-points br)) - ;; Do any of the sides intersect? - (let [segments (.segments line) - mbr-segments (s/mbr->line-segments br)] - (loop [segments segments] - (when-let [segment (first segments)] - (let [intersects? (loop [mbr-segments mbr-segments] - (when-let [mbr-segment (first mbr-segments)] - (or (seq (asi/intersections segment mbr-segment)) - (recur (rest mbr-segments)))))] - (or intersects? (recur (rest segments)))))))))))) - + ;; Delegate to Java implementation + (let [java-line (cmr.spatial.shape.LineString. + (name (.coordinate_system line)) + (vec (p/points->ords (.points line)))) + java-mbr (cmr.spatial.shape.Mbr. (.west br) (.north br) (.east br) (.south br))] + (LineStringIntersections/intersectsMbr java-line java-mbr))) (defn intersects-line-string? - "Returns true if the line string instersects the other line string" + "Returns true if the line string intersects the other line string" [line1 line2] - (u/any-true? (fn [[s1 s2]] - (seq (asi/intersections s1 s2))) - (for [segment1 (:segments line1) - segment2 (:segments line2)] - [segment1 segment2]))) + ;; Delegate to Java implementation + (let [java-line1 (cmr.spatial.shape.LineString. + (name (:coordinate-system line1)) + (vec (p/points->ords (:points line1)))) + java-line2 (cmr.spatial.shape.LineString. + (name (:coordinate-system line2)) + (vec (p/points->ords (:points line2))))] + (LineStringIntersections/intersectsLineString java-line1 java-line2))) (extend-protocol v/SpatialValidation cmr.spatial.line_string.LineString diff --git a/spatial-lib/src/cmr/spatial/lr_binary_search.clj b/spatial-lib/src/cmr/spatial/lr_binary_search.clj index e79d803f03..5b78d87a05 100644 --- a/spatial-lib/src/cmr/spatial/lr_binary_search.clj +++ b/spatial-lib/src/cmr/spatial/lr_binary_search.clj @@ -264,7 +264,13 @@ (do (warn "Use mbr from one of the points in the polygon because lr is not found " (pr-str polygon)) - (m/point->mbr (-> polygon :rings first :points first)))))))) + ;; Get first point - handle both polygon and ring + (let [first-point (if (:rings polygon) + ;; It's a polygon - get first point of first ring + (-> polygon :rings first :points first) + ;; It's a ring - get first point directly + (-> polygon :points first))] + (m/point->mbr first-point)))))))) (comment (require '[cmr.spatial.kml :as kml]) diff --git a/spatial-lib/src/cmr/spatial/mbr.clj b/spatial-lib/src/cmr/spatial/mbr.clj index 84a72738c4..82f69b2104 100644 --- a/spatial-lib/src/cmr/spatial/mbr.clj +++ b/spatial-lib/src/cmr/spatial/mbr.clj @@ -10,7 +10,9 @@ [cmr.spatial.validation :as sv] [cmr.spatial.messages :as msg] [cmr.common.dev.record-pretty-printer :as record-pretty-printer]) - (:import cmr.spatial.point.Point)) + (:import + cmr.spatial.point.Point + [cmr.spatial.geometry MbrIntersections])) (primitive-math/use-primitive-operators) @@ -299,40 +301,23 @@ (defn non-crossing-intersects-br? "Specialized version of intersects-br? for two mbrs that don't cross the antimeridian. - Returns true if the mbr intersects the other bounding rectangle." - [coord-sys ^Mbr m1 ^Mbr m2] + Returns true if the mbr intersects the other bounding rectangle." + [coord-sys ^cmr.spatial.mbr.Mbr m1 ^cmr.spatial.mbr.Mbr m2] (pj/assert (not (or (crosses-antimeridian? m1) (crosses-antimeridian? m2)))) - (let [w1 (.west m1) - n1 (.north m1) - e1 (.east m1) - s1 (.south m1) - w2 (.west m2) - n2 (.north m2) - e2 (.east m2) - s2 (.south m2) - m1-touches-north? (double-approx= n1 90.0 0.0000001) - m1-touches-south? (double-approx= s1 -90.0 0.0000001) - m2-touches-north? (double-approx= n2 90.0 0.0000001) - m2-touches-south? (double-approx= s2 -90.0 0.0000001)] - (or (and (range-intersects? w1 e1 w2 e2) - (range-intersects? s1 n1 s2 n2)) - (and (= coord-sys :geodetic) - (or (and m1-touches-north? m2-touches-north?) - (and m1-touches-south? m2-touches-south?)))))) + ;; Delegate to Java implementation + (let [java-mbr1 (cmr.spatial.shape.Mbr. (.west m1) (.north m1) (.east m1) (.south m1)) + java-mbr2 (cmr.spatial.shape.Mbr. (.west m2) (.north m2) (.east m2) (.south m2))] + (MbrIntersections/nonCrossingIntersects (name coord-sys) java-mbr1 java-mbr2))) (defn intersects-br? "Returns true if the mbr intersects the other bounding rectangle" - [coord-sys ^Mbr mbr ^Mbr other-br] - (if (and (not (crosses-antimeridian? mbr)) (not (crosses-antimeridian? other-br))) - ;; optimized case for mbrs that don't cross the antimeridian - (non-crossing-intersects-br? coord-sys mbr other-br) - (let [[m1-east m1-west] (split-across-antimeridian mbr) - [m2-east m2-west] (split-across-antimeridian other-br)] - (or (non-crossing-intersects-br? coord-sys m1-east m2-east) - (and m2-west (non-crossing-intersects-br? coord-sys m1-east m2-west)) - (and m1-west (non-crossing-intersects-br? coord-sys m1-west m2-east)) - (and m1-west m2-west (non-crossing-intersects-br? coord-sys m1-west m2-west)))))) + [coord-sys ^cmr.spatial.mbr.Mbr mbr ^cmr.spatial.mbr.Mbr other-br] + ;; Delegate to Java implementation which handles all edge cases + (let [java-mbr1 (cmr.spatial.shape.Mbr. (.west mbr) (.north mbr) (.east mbr) (.south mbr)) + java-mbr2 (cmr.spatial.shape.Mbr. (.west other-br) (.north other-br) + (.east other-br) (.south other-br))] + (MbrIntersections/mbrsIntersect (name coord-sys) java-mbr1 java-mbr2))) (defn intersections "Returns the intersection of the two minimum bounding rectangles. This could return multiple mbrs diff --git a/spatial-lib/src/cmr/spatial/point.clj b/spatial-lib/src/cmr/spatial/point.clj index 336fc212b6..3d3bf8ef0e 100644 --- a/spatial-lib/src/cmr/spatial/point.clj +++ b/spatial-lib/src/cmr/spatial/point.clj @@ -12,7 +12,9 @@ [pjstadig.assertions :as pj] [cmr.spatial.derived :as d] [cmr.spatial.validation :as v] - [cmr.spatial.messages :as msg])) + [cmr.spatial.messages :as msg]) + (:import + [cmr.spatial.geometry PointIntersections])) (primitive-math/use-primitive-operators) @@ -394,19 +396,11 @@ (approx= ([expected n] (approx= expected n DELTA)) - ([^Point p1 ^Point p2 ^double delta] - (let [lon1 (.lon p1) - lat1 (.lat p1) - lon2 (.lon p2) - lat2 (.lat p2) - np? #(double-approx= ^double % 90.0 delta) - sp? #(double-approx= ^double % -90.0 delta) - am? #(double-approx= ^double (abs ^double %) 180.0 delta)] - (and (double-approx= lat1 lat2 delta) - (or (double-approx= lon1 lon2 delta) - (and (am? lon1) (am? lon2)) - (and (np? lat1) (np? lat2)) - (and (sp? lat1) (sp? lat2)))))))) + ([^cmr.spatial.point.Point p1 ^cmr.spatial.point.Point p2 ^double delta] + ;; Delegate to Java implementation + (let [java-p1 (cmr.spatial.shape.Point. (.lon p1) (.lat p1)) + java-p2 (cmr.spatial.shape.Point. (.lon p2) (.lat p2))] + (PointIntersections/pointsApproxEqual java-p1 java-p2 delta))))) (extend-protocol d/DerivedCalculator diff --git a/spatial-lib/src/cmr/spatial/points_validation_helpers.clj b/spatial-lib/src/cmr/spatial/points_validation_helpers.clj index a35b49020a..b1ead39029 100644 --- a/spatial-lib/src/cmr/spatial/points_validation_helpers.clj +++ b/spatial-lib/src/cmr/spatial/points_validation_helpers.clj @@ -20,7 +20,7 @@ (points->pairs points identity)) ([points xform] (let [modified-points (map xform points)] - (conj (for [idx (range (dec (count points)))] + (conj (for [^long idx (range (dec (count points)))] (vector (nth modified-points idx) (nth modified-points (inc idx)))))))) diff --git a/spatial-lib/src/cmr/spatial/relations.clj b/spatial-lib/src/cmr/spatial/relations.clj index e479853505..04a2753547 100644 --- a/spatial-lib/src/cmr/spatial/relations.clj +++ b/spatial-lib/src/cmr/spatial/relations.clj @@ -17,7 +17,9 @@ cmr.spatial.line_string.LineString cmr.spatial.mbr.Mbr cmr.spatial.point.Point - cmr.spatial.polygon.Polygon)) + cmr.spatial.polygon.Polygon + cmr.spatial.relations.ShapeIntersections + cmr.spatial.relations.ShapePredicate)) (defprotocol SpatialRelations "Defines functions for determining relations between different spatial areas." @@ -318,3 +320,23 @@ ;; Shape is the second argument so that the polymorphic protocol dispatch can be used ;; on the first argument. (f (d/calculate-derived other-shape) shape)))) + +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;; Java Interop API +;; Provides access to Java intersection implementations for use in ES plugin +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; + +(defn java-shape->intersects-fn + "Creates an intersection function from a Java shape (Point, Mbr, etc.). + Returns a Java ShapePredicate that tests if another Java shape intersects. + + This is intended for use in the ES spatial plugin where working with Java + shapes directly avoids the overhead of converting to/from Clojure records. + + Example: + (let [java-point (cmr.spatial.shape.Point. 10.0 20.0) + predicate (java-shape->intersects-fn java-point) + other-point (cmr.spatial.shape.Point. 10.0 20.0)] + (.intersects predicate other-point)) ;; => true" + [java-shape] + (ShapeIntersections/createIntersectsFn java-shape)) diff --git a/spatial-lib/src/cmr/spatial/ring_relations.clj b/spatial-lib/src/cmr/spatial/ring_relations.clj index a8bfb40da1..0db9767aac 100644 --- a/spatial-lib/src/cmr/spatial/ring_relations.clj +++ b/spatial-lib/src/cmr/spatial/ring_relations.clj @@ -14,7 +14,8 @@ cmr.spatial.cartesian_ring.CartesianRing cmr.spatial.geodetic_ring.GeodeticRing cmr.spatial.point.Point - cmr.spatial.mbr.Mbr)) + cmr.spatial.mbr.Mbr + [cmr.spatial.geometry RingIntersections])) (primitive-math/use-primitive-operators) (defn ring @@ -85,23 +86,31 @@ [_] :geodetic)) +(defn- ->java-ring + "Converts a Clojure ring (GeodeticRing or CartesianRing) to its Java equivalent. + Returns Java rings unchanged." + [ring] + (cond + ;; Already Java rings - return unchanged + (instance? cmr.spatial.internal.ring.GeodeticRing ring) ring + (instance? cmr.spatial.internal.ring.CartesianRing ring) ring + ;; Clojure rings - convert to Java + (instance? GeodeticRing ring) (cmr.spatial.internal.ring.GeodeticRing/createRing (:points ring)) + (instance? CartesianRing ring) (cmr.spatial.internal.ring.CartesianRing/createRing (:points ring)) + :else (throw (IllegalArgumentException. (str "Unsupported ring type: " (class ring)))))) + +(defn- ->java-mbr + "Converts a Clojure MBR to a Java MBR. Returns Java MBRs unchanged." + [mbr] + (if (instance? cmr.spatial.shape.Mbr mbr) + mbr + (cmr.spatial.shape.Mbr. (:west mbr) (:north mbr) (:east mbr) (:south mbr)))) + (defn intersects-ring? "Returns true if the rings intersect each other." [r1 r2] - (or - ;; Do any of the line-segments intersect? - ;; Performance enhancement: this should use the multiple arc intersection algorithm to avoid O(N^2) intersections - (some (fn [[line1 line2]] - (seq (asi/intersections line1 line2))) - (for [ls1 (segments r1) - ls2 (segments r2)] - [ls1 ls2])) - - ;; Is ring 2 inside ring 1? Only one point check is required - (covers-point? r1 (first (:points r2))) - - ;; Is ring 1 inside ring 2? Only one point check is required - (covers-point? r2 (first (:points r1))))) + ;; Delegate to Java implementation + (RingIntersections/intersectsRing (->java-ring r1) (->java-ring r2))) (defn covers-ring? "Returns true if the ring covers the other ring." @@ -162,21 +171,8 @@ (defn intersects-br? "Returns true if the ring intersects the br" [ring ^Mbr br] - (when (m/intersects-br? (coordinate-system ring) (:mbr ring) br) - (if (m/single-point? br) - (covers-point? ring (p/point (.west br) (.north br))) - - (or - ;; Does the br cover any points of the ring? - (if (= :geodetic (coordinate-system ring)) - (some #(m/geodetic-covers-point? br %) (:points ring)) - (some #(m/cartesian-covers-point? br %) (:points ring))) - - ;; Does the ring completely contain the br? We only need to check one point of the br - (covers-point? ring (p/point (.west br) (.north br))) - - ;; Do any of the sides intersect? - (lines-intersects-br-sides? (segments ring) br))))) + ;; Delegate to Java implementation + (RingIntersections/intersectsBr (->java-ring ring) (->java-mbr br))) (defn intersects-line-string? "Returns true if the ring intersects the line" diff --git a/spatial-lib/src/cmr/spatial/ring_validations.clj b/spatial-lib/src/cmr/spatial/ring_validations.clj index bcff1d8a20..e71d89a440 100644 --- a/spatial-lib/src/cmr/spatial/ring_validations.clj +++ b/spatial-lib/src/cmr/spatial/ring_validations.clj @@ -8,7 +8,7 @@ [cmr.spatial.validation :as v] [primitive-math]) #_{:clj-kondo/ignore [:unused-import]} - (:import cmr.spatial.arc.Arc)) + (:import cmr.spatial.internal.arc.Arc)) (primitive-math/use-primitive-operators) (defn- ring-closed-validation diff --git a/spatial-lib/src/cmr/spatial/serialize.clj b/spatial-lib/src/cmr/spatial/serialize.clj index 92f2f6fd4e..a6b7c167cd 100644 --- a/spatial-lib/src/cmr/spatial/serialize.clj +++ b/spatial-lib/src/cmr/spatial/serialize.clj @@ -20,7 +20,10 @@ [cmr.spatial.mbr :as m] [cmr.spatial.point :as p] [cmr.spatial.polygon :as poly] - [cmr.spatial.ring-relations :as rr])) + [cmr.spatial.ring-relations :as rr]) + (:import + [cmr.spatial.serialize OrdsInfoShapes] + [cmr.spatial.shape SpatialShape Point Mbr Polygon Ring LineString])) ;; Some thoughts about how to store the elasticsearch data in a way that preserves space and accuracy. (comment @@ -215,29 +218,39 @@ ;; Create a combined sequence of all the shape ordinates :ords (u/mapcatv :ords infos)})) +(defn- java-shape->clojure-shape + "Converts a Java shape object to its Clojure record equivalent." + [java-shape] + (cond + (instance? Point java-shape) + (p/point (.getLon ^Point java-shape) (.getLat ^Point java-shape)) + + (instance? Mbr java-shape) + (m/mbr (.getWest ^Mbr java-shape) (.getNorth ^Mbr java-shape) + (.getEast ^Mbr java-shape) (.getSouth ^Mbr java-shape)) + + (instance? Polygon java-shape) + (let [^Polygon poly-obj java-shape + rings (.getRings poly-obj) + coord-sys (keyword (.getCoordinateSystem poly-obj)) + clj-rings (mapv (fn [^Ring ring] + (rr/ords->ring (keyword (.getCoordinateSystem ring)) (.getOrdinates ring))) + rings)] + (poly/polygon coord-sys clj-rings)) + + (instance? LineString java-shape) + (let [^LineString ls-obj java-shape] + (l/ords->line-string (keyword (.getCoordinateSystem ls-obj)) (.getOrdinates ls-obj))) + + :else + (throw (Exception. (str "Unknown Java shape type: " (class java-shape)))))) + (defn ords-info->shapes - "Converts the ords-info data and ordinates into a sequence of shapes" + "Converts the ords-info data and ordinates into a sequence of shapes. + Uses the Java implementation for deserialization, then converts back to Clojure records." [ords-info ords] - (loop [ords-info-pairs (partition 2 ords-info) - ords (vec ords) - shapes []] - (if (empty? ords-info-pairs) - shapes - (let [[int-type size] (first ords-info-pairs) - type (integer->shape-type int-type) - _ (when-not type - (errors/internal-error! - (format "Could not get a shape type from integer [%s]. Ords info: %s" - int-type (pr-str ords-info)))) - shape-ords (subvec ords 0 size) - shape (stored-ords->shape type shape-ords) - ;; If shape is a hole add it to the last shape - shapes (if (or (= type :geodetic-hole) (= type :cartesian-hole)) - (update-in shapes [(dec (count shapes)) :rings] conj shape) - (conj shapes shape))] - (recur (rest ords-info-pairs) - (subvec ords size) - shapes))))) + (let [java-shapes (OrdsInfoShapes/ordsInfoToShapes ords-info ords)] + (mapv java-shape->clojure-shape java-shapes))) ;; This comment left in to show how to test the performance of a spatial search intersection between ;; a search area and one possible area. diff --git a/spatial-lib/src/java/cmr/spatial/geometry/CircleIntersections.java b/spatial-lib/src/java/cmr/spatial/geometry/CircleIntersections.java new file mode 100644 index 0000000000..aed7ce324b --- /dev/null +++ b/spatial-lib/src/java/cmr/spatial/geometry/CircleIntersections.java @@ -0,0 +1,222 @@ +package cmr.spatial.geometry; + +import cmr.spatial.shape.Circle; +import cmr.spatial.shape.LineString; +import cmr.spatial.shape.Mbr; +import cmr.spatial.shape.Point; +import cmr.spatial.shape.Polygon; +import cmr.spatial.shape.Ring; +import cmr.spatial.internal.ring.GeodeticRing; +import cmr.spatial.math.MathUtils; + +import java.util.ArrayList; +import java.util.Collections; +import java.util.List; + +/** + * Implements Circle intersection testing. + * Translates covers-point? and circle->polygon from circle.clj. + * + * A circle is defined by a center point and radius (in meters). + * For complex intersections, circles are converted to polygon approximations. + */ +public class CircleIntersections { + + /** + * Minimum radius in meters. + */ + public static final double MIN_RADIUS = 10.0; + + /** + * Maximum radius in meters. + */ + public static final double MAX_RADIUS = 6000000.0; + + /** + * Radius of the earth in meters for polygon approximation. + * Matches EARTH_RADIUS_APPROX in circle.clj. + */ + private static final double EARTH_RADIUS_APPROX = 6378137.0; + + /** + * Default number of points for polygon approximation. + * More points = more accurate but slower. + */ + private static final int DEFAULT_NUM_POINTS = 64; + + /** + * Validates circle parameters. + * + * @param circle The circle to validate + * @param numPoints Number of points (if applicable) + * @throws IllegalArgumentException if parameters are invalid + */ + private static void validateCircleParams(Circle circle, int numPoints) { + double radius = circle.getRadius(); + + if (radius < MIN_RADIUS || radius > MAX_RADIUS) { + throw new IllegalArgumentException( + String.format("Circle radius must be between %f and %f meters, got: %f", + MIN_RADIUS, MAX_RADIUS, radius)); + } + + if (numPoints < 3) { + throw new IllegalArgumentException( + String.format("numPoints must be >= 3 for a valid polygon, got: %d", numPoints)); + } + } + + /** + * Tests if a circle covers a point. + * A circle covers a point if the distance from the center to the point + * is less than or equal to the radius. + * + * @param circle The circle to test + * @param point The point to check + * @return true if the circle covers the point + */ + public static boolean coversPoint(Circle circle, Point point) { + validateCircleParams(circle, 3); // numPoints not used here, but validate radius + double distance = MathUtils.distance(circle.getCenter(), point); + return distance <= circle.getRadius(); + } + + /** + * Converts a circle to a polygon approximation with the default number of points. + * + * @param circle The circle to convert + * @return A geodetic polygon approximating the circle + */ + public static Polygon circleToPolygon(Circle circle) { + return circleToPolygon(circle, DEFAULT_NUM_POINTS); + } + + /** + * Converts a circle to a polygon approximation. + * + * Algorithm from circle.clj:circle->polygon, which is based on: + * https://github.com/gabzim/circle-to-polygon + * + * Creates a polygon by calculating n points around the circle's perimeter + * using spherical geometry. Points are calculated by moving along great circles + * at the specified radius from the center. + * + * @param circle The circle to convert + * @param numPoints Number of points in the polygon (more = more accurate) + * @return A geodetic polygon approximating the circle + */ + public static Polygon circleToPolygon(Circle circle, int numPoints) { + validateCircleParams(circle, numPoints); + + Point center = circle.getCenter(); + double radius = circle.getRadius(); + + double lon1Rad = Math.toRadians(center.getLon()); + double lat1Rad = Math.toRadians(center.getLat()); + double rRad = radius / EARTH_RADIUS_APPROX; + + List points = new ArrayList<>(); + + for (int i = 0; i <= numPoints; i++) { + double theta = -2.0 * Math.PI * (i / (double) numPoints); + + // Calculate point at this angle around the circle + double pLatRad = Math.asin( + Math.sin(lat1Rad) * Math.cos(rRad) + + Math.cos(lat1Rad) * Math.sin(rRad) * Math.cos(theta) + ); + + double pLonDelta = Math.atan2( + Math.sin(theta) * Math.sin(rRad) * Math.cos(lat1Rad), + Math.cos(rRad) - Math.sin(lat1Rad) * Math.sin(pLatRad) + ); + + double pLon = sanitizeLongitude(Math.toDegrees(lon1Rad + pLonDelta)); + double pLat = Math.toDegrees(pLatRad); + + points.add(new Point(pLon, pLat)); + } + + // Replace last point with first to ensure closed ring + // (handles precision issues) + points.set(points.size() - 1, points.get(0)); + + // Convert points to ordinates + List ordinates = new ArrayList<>(); + for (Point p : points) { + ordinates.add(p.getLon()); + ordinates.add(p.getLat()); + } + + // Create a Ring and wrap in Polygon + Ring ring = new Ring("geodetic", ordinates); + return new Polygon("geodetic", Collections.singletonList(ring)); + } + + /** + * Tests if a circle intersects a polygon. + * Converts the circle to a polygon approximation and delegates to polygon intersection. + * + * @param circle The circle to test + * @param polygon The polygon to check + * @return true if the circle intersects the polygon + */ + public static boolean intersectsPolygon(Circle circle, Polygon polygon) { + Polygon circlePoly = circleToPolygon(circle); + return PolygonIntersections.intersectsPolygon(circlePoly, polygon); + } + + /** + * Tests if a circle intersects a ring. + * Converts the circle to a polygon approximation and delegates to polygon-ring intersection. + * + * @param circle The circle to test + * @param ring The ring to check + * @return true if the circle intersects the ring + */ + public static boolean intersectsRing(Circle circle, Ring ring) { + Polygon circlePoly = circleToPolygon(circle); + return PolygonIntersections.intersectsRing(circlePoly, ring); + } + + /** + * Tests if a circle intersects a bounding rectangle. + * Converts the circle to a polygon approximation and delegates to polygon-BR intersection. + * + * @param circle The circle to test + * @param mbr The bounding rectangle to check + * @return true if the circle intersects the bounding rectangle + */ + public static boolean intersectsBr(Circle circle, Mbr mbr) { + Polygon circlePoly = circleToPolygon(circle); + return PolygonIntersections.intersectsBr(circlePoly, mbr); + } + + /** + * Tests if a circle intersects a line string. + * Converts the circle to a polygon approximation and delegates to polygon-linestring intersection. + * + * @param circle The circle to test + * @param lineString The line string to check + * @return true if the circle intersects the line string + */ + public static boolean intersectsLineString(Circle circle, LineString lineString) { + Polygon circlePoly = circleToPolygon(circle); + return PolygonIntersections.intersectsLineString(circlePoly, lineString); + } + + /** + * Sanitizes longitude to be within -180 to 180 range. + * + * @param lon Longitude in degrees + * @return Sanitized longitude + */ + private static double sanitizeLongitude(double lon) { + if (lon > 180.0) { + return lon - 360.0; + } else if (lon < -180.0) { + return lon + 360.0; + } + return lon; + } +} diff --git a/spatial-lib/src/java/cmr/spatial/geometry/LineStringIntersections.java b/spatial-lib/src/java/cmr/spatial/geometry/LineStringIntersections.java new file mode 100644 index 0000000000..56dc313aef --- /dev/null +++ b/spatial-lib/src/java/cmr/spatial/geometry/LineStringIntersections.java @@ -0,0 +1,348 @@ +package cmr.spatial.geometry; + +import cmr.spatial.internal.arc.Arc; +import cmr.spatial.internal.arc.ArcLineSegmentIntersections; +import cmr.spatial.shape.LineString; +import cmr.spatial.shape.Mbr; +import cmr.spatial.shape.Point; +import cmr.spatial.internal.segment.LineSegment; +import cmr.spatial.math.MathUtils; + +import java.util.ArrayList; +import java.util.HashSet; +import java.util.List; +import java.util.Set; + +/** + * Implements LineString intersection testing. + * Translates covers-point?, intersects-br?, and intersects-line-string? from Clojure. + */ +public class LineStringIntersections { + + /** + * Tests if a line string covers a point. + * A line string covers a point if the point is one of the line's points + * or if the point lies on one of the line's segments. + * + * @param lineString The line string to test + * @param point The point to check + * @return true if the line string covers the point + */ + public static boolean coversPoint(LineString lineString, Point point) { + String coordSystem = lineString.getCoordinateSystem(); + List ordinates = lineString.getOrdinates(); + + // Convert ordinates to points + List points = ordinatesToPoints(ordinates); + + // Check if point is in the point set + if (points.contains(point)) { + return true; + } + + // Check if point lies on any segment + if ("geodetic".equals(coordSystem)) { + List arcs = pointsToArcs(points); + for (Arc arc : arcs) { + if (arc.pointOnArc(point)) { + return true; + } + } + } else { // cartesian + List segments = pointsToLineSegments(points); + for (LineSegment segment : segments) { + if (segment.pointOnSegment(point)) { + return true; + } + } + } + + return false; + } + + /** + * Tests if a line string intersects a minimal bounding rectangle (MBR). + */ + public static boolean intersectsMbr(LineString lineString, Mbr mbr) { + String coordSystem = lineString.getCoordinateSystem(); + List ordinates = lineString.getOrdinates(); + List points = ordinatesToPoints(ordinates); + + // Calculate line MBR + Mbr lineMbr = calculateLineMbr(coordSystem, points); + + // Quick reject if MBRs don't intersect + if (!MbrIntersections.mbrsIntersect(coordSystem, lineMbr, mbr)) { + return false; + } + + // Check if MBR is single point + if (isSinglePoint(mbr)) { + Point testPoint = new Point(mbr.getWest(), mbr.getNorth()); + return coversPoint(lineString, testPoint); + } + + // Check if MBR covers any line points + for (Point p : points) { + if (mbrCoversPoint(coordSystem, mbr, p)) { + return true; + } + } + + // Check if line covers any MBR corner points + List corners = mbrCornerPoints(mbr); + for (Point corner : corners) { + if (coversPoint(lineString, corner)) { + return true; + } + } + + // Check if any segments intersect MBR edges + List mbrSegments = LineSegment.mbrToLineSegments(mbr); + + if ("geodetic".equals(coordSystem)) { + List arcs = pointsToArcs(points); + for (Arc arc : arcs) { + for (LineSegment mbrSeg : mbrSegments) { + if (ArcLineSegmentIntersections.intersects(arc, mbrSeg)) { + return true; + } + } + } + } else { // cartesian + List segments = pointsToLineSegments(points); + for (LineSegment seg : segments) { + for (LineSegment mbrSeg : mbrSegments) { + if (ArcLineSegmentIntersections.intersects(seg, mbrSeg)) { + return true; + } + } + } + } + + return false; + } + + /** + * Tests if two line strings intersect. + */ + public static boolean intersectsLineString(LineString lineString1, LineString lineString2) { + List points1 = ordinatesToPoints(lineString1.getOrdinates()); + List points2 = ordinatesToPoints(lineString2.getOrdinates()); + + String coordSys1 = lineString1.getCoordinateSystem(); + String coordSys2 = lineString2.getCoordinateSystem(); + + // Get segments for both lines + Object[] segments1 = getSegments(coordSys1, points1); + Object[] segments2 = getSegments(coordSys2, points2); + + // Check all segment pairs + for (Object seg1 : segments1) { + for (Object seg2 : segments2) { + if (ArcLineSegmentIntersections.intersects(seg1, seg2)) { + return true; + } + } + } + + return false; + } + + // Helper methods + + private static List ordinatesToPoints(List ordinates) { + List points = new ArrayList<>(); + for (int i = 0; i < ordinates.size(); i += 2) { + points.add(new Point(ordinates.get(i), ordinates.get(i + 1))); + } + return points; + } + + private static List pointsToArcs(List points) { + List arcs = new ArrayList<>(); + for (int i = 0; i < points.size() - 1; i++) { + try { + arcs.add(Arc.createArc(points.get(i), points.get(i + 1))); + } catch (IllegalArgumentException e) { + // Skip degenerate or antipodal point pairs + } + } + return arcs; + } + + private static List pointsToLineSegments(List points) { + return LineSegment.pointsToLineSegments(points); + } + + private static Object[] getSegments(String coordSystem, List points) { + if ("geodetic".equals(coordSystem)) { + return pointsToArcs(points).toArray(); + } else { + return pointsToLineSegments(points).toArray(); + } + } + + private static Mbr calculateLineMbr(String coordSystem, List points) { + if (points.isEmpty()) { + return new Mbr(0, 0, 0, 0); + } + + if ("geodetic".equals(coordSystem)) { + // For geodetic coordinates, compute per-segment arc MBRs and union them + Mbr result = null; + + for (int i = 0; i < points.size() - 1; i++) { + Point p1 = points.get(i); + Point p2 = points.get(i + 1); + + try { + // Create arc and get its MBRs (may be 1 or 2 depending on geometry) + Arc arc = Arc.createArc(p1, p2); + Mbr[] arcMbrs = arc.getMbrs(); + + for (Mbr arcMbr : arcMbrs) { + if (result == null) { + result = arcMbr; + } else { + result = unionMbrs(result, arcMbr); + } + } + } catch (IllegalArgumentException e) { + // If arc creation fails (duplicate or antipodal points), + // fall back to treating as a single point + if (result == null) { + result = new Mbr(p1.getLon(), p1.getLat(), p1.getLon(), p1.getLat()); + } + } + } + + return result != null ? result : new Mbr(0, 0, 0, 0); + } else { + // For cartesian coordinates, use simple min/max + double west = points.get(0).getLon(); + double east = west; + double north = points.get(0).getLat(); + double south = north; + + for (Point p : points) { + double lon = p.getLon(); + double lat = p.getLat(); + if (lon < west) west = lon; + if (lon > east) east = lon; + if (lat > north) north = lat; + if (lat < south) south = lat; + } + + return new Mbr(west, north, east, south); + } + } + + /** + * Unions two MBRs into a single MBR that contains both. + */ + private static Mbr unionMbrs(Mbr mbr1, Mbr mbr2) { + double w1 = mbr1.getWest(); + double e1 = mbr1.getEast(); + double w2 = mbr2.getWest(); + double e2 = mbr2.getEast(); + boolean crosses1 = crossesAntimeridian(mbr1); + boolean crosses2 = crossesAntimeridian(mbr2); + double west; + double east; + + if (crosses1 && crosses2) { + west = Math.min(w1, w2); + east = Math.max(e1, e2); + if (west <= east) { + west = -180.0; + east = 180.0; + } + } else if (crosses1 || crosses2) { + Mbr crossing = crosses1 ? mbr1 : mbr2; + Mbr other = crosses1 ? mbr2 : mbr1; + w1 = crossing.getWest(); + e1 = crossing.getEast(); + w2 = other.getWest(); + e2 = other.getEast(); + + double westDist = w1 - w2; + double eastDist = e2 - e1; + + if (westDist <= 0.0 || eastDist <= 0.0) { + west = w1; + east = e1; + } else if (eastDist < westDist) { + west = w1; + east = e2; + } else { + west = w2; + east = e1; + } + + if (west <= east) { + west = -180.0; + east = 180.0; + } + } else { + if (w1 > w2) { + double tmp = w1; + w1 = w2; + w2 = tmp; + tmp = e1; + e1 = e2; + e2 = tmp; + } + + west = Math.min(w1, w2); + east = Math.max(e1, e2); + + double dist = east - west; + double altWest = w2; + double altEast = e1; + double altDist = (180.0 - altWest) + (altEast - -180.0); + + if (altDist < dist) { + west = altWest; + east = altEast; + } + } + + double north = Math.max(mbr1.getNorth(), mbr2.getNorth()); + double south = Math.min(mbr1.getSouth(), mbr2.getSouth()); + return new Mbr(west, north, east, south); + } + + private static boolean crossesAntimeridian(Mbr mbr) { + return mbr.getWest() > mbr.getEast(); + } + + private static boolean isSinglePoint(Mbr mbr) { + return MathUtils.doubleApprox(mbr.getWest(), mbr.getEast(), MathUtils.DELTA) + && MathUtils.doubleApprox(mbr.getNorth(), mbr.getSouth(), MathUtils.DELTA); + } + + private static boolean mbrCoversPoint(String coordSystem, Mbr mbr, Point point) { + double lon = point.getLon(); + double lat = point.getLat(); + + boolean latInRange = lat >= mbr.getSouth() && lat <= mbr.getNorth(); + if (!latInRange) return false; + + if ("geodetic".equals(coordSystem) && mbr.getWest() > mbr.getEast()) { + // Crosses antimeridian + return lon >= mbr.getWest() || lon <= mbr.getEast(); + } else { + return lon >= mbr.getWest() && lon <= mbr.getEast(); + } + } + + private static List mbrCornerPoints(Mbr mbr) { + List corners = new ArrayList<>(); + corners.add(new Point(mbr.getWest(), mbr.getNorth())); // NW + corners.add(new Point(mbr.getEast(), mbr.getNorth())); // NE + corners.add(new Point(mbr.getEast(), mbr.getSouth())); // SE + corners.add(new Point(mbr.getWest(), mbr.getSouth())); // SW + return corners; + } +} diff --git a/spatial-lib/src/java/cmr/spatial/geometry/MbrIntersections.java b/spatial-lib/src/java/cmr/spatial/geometry/MbrIntersections.java new file mode 100644 index 0000000000..adb7fe9164 --- /dev/null +++ b/spatial-lib/src/java/cmr/spatial/geometry/MbrIntersections.java @@ -0,0 +1,137 @@ +package cmr.spatial.geometry; + +import cmr.spatial.shape.Mbr; +import cmr.spatial.math.MathUtils; +import java.util.ArrayList; +import java.util.List; + +/** + * Implements MBR-to-MBR intersection testing. + * Translates the intersects-br? protocol method from Clojure. + */ +public class MbrIntersections { + + /** + * Tests if two MBRs intersect using geodetic (spherical) coordinate system. + * Handles both simple cases (neither crosses antimeridian) and complex cases + * (one or both cross the antimeridian). + * + * @param mbr1 First minimum bounding rectangle + * @param mbr2 Second minimum bounding rectangle + * @return true if the MBRs intersect + */ + public static boolean mbrsIntersect(Mbr mbr1, Mbr mbr2) { + return mbrsIntersect("geodetic", mbr1, mbr2); + } + + /** + * Tests if two MBRs intersect using the specified coordinate system. + * + * @param coordSys Coordinate system: "geodetic" or "cartesian" + * @param mbr1 First minimum bounding rectangle + * @param mbr2 Second minimum bounding rectangle + * @return true if the MBRs intersect + */ + public static boolean mbrsIntersect(String coordSys, Mbr mbr1, Mbr mbr2) { + // Optimized path: if neither crosses antimeridian, use simpler logic + if (!crossesAntimeridian(mbr1) && !crossesAntimeridian(mbr2)) { + return nonCrossingIntersects(coordSys, mbr1, mbr2); + } + + // Complex path: handle antimeridian crossing + List mbr1_parts = splitAcrossAntimeridian(mbr1); + List mbr2_parts = splitAcrossAntimeridian(mbr2); + + // Check all 4 combinations of split parts + for (Mbr m1 : mbr1_parts) { + for (Mbr m2 : mbr2_parts) { + if (nonCrossingIntersects(coordSys, m1, m2)) { + return true; + } + } + } + return false; + } + + /** + * Tests if two non-antimeridian-crossing MBRs intersect. + * For two normal (non-wrapping) rectangles, intersection requires both: + * 1. Longitude ranges to overlap + * 2. Latitude ranges to overlap + * Special case: In geodetic coordinates, rectangles touching at poles count as intersecting. + * + * @param coordSys Coordinate system + * @param mbr1 First MBR (must not cross antimeridian) + * @param mbr2 Second MBR (must not cross antimeridian) + * @return true if the MBRs intersect + */ + public static boolean nonCrossingIntersects(String coordSys, Mbr mbr1, Mbr mbr2) { + double w1 = mbr1.getWest(); + double n1 = mbr1.getNorth(); + double e1 = mbr1.getEast(); + double s1 = mbr1.getSouth(); + + double w2 = mbr2.getWest(); + double n2 = mbr2.getNorth(); + double e2 = mbr2.getEast(); + double s2 = mbr2.getSouth(); + + // Check if poles are touched + boolean m1_touches_north = MathUtils.doubleApprox(n1, 90.0, 0.0000001); + boolean m1_touches_south = MathUtils.doubleApprox(s1, -90.0, 0.0000001); + boolean m2_touches_north = MathUtils.doubleApprox(n2, 90.0, 0.0000001); + boolean m2_touches_south = MathUtils.doubleApprox(s2, -90.0, 0.0000001); + + // Main intersection check: both longitude and latitude ranges must intersect + boolean ranges_intersect = MathUtils.rangeIntersects(w1, e1, w2, e2) + && MathUtils.rangeIntersects(s1, n1, s2, n2); + + if (ranges_intersect) { + return true; + } + + // Geodetic special case: touching at poles counts as intersection + if ("geodetic".equals(coordSys)) { + if ((m1_touches_north && m2_touches_north) + || (m1_touches_south && m2_touches_south)) { + return true; + } + } + + return false; + } + + /** + * Returns true if an MBR crosses the antimeridian (±180° boundary). + * An MBR crosses if its western bound is greater than its eastern bound. + * + * @param mbr The minimum bounding rectangle + * @return true if the MBR crosses the antimeridian + */ + public static boolean crossesAntimeridian(Mbr mbr) { + return mbr.getWest() > mbr.getEast(); + } + + /** + * Splits an MBR across the antimeridian if it crosses. + * If the MBR crosses the antimeridian, returns a list of 2 MBRs (east and west portions). + * If the MBR doesn't cross, returns a list containing just the original MBR. + * An MBR that crosses is split at ±180° into two non-crossing parts. + * + * @param mbr The minimum bounding rectangle + * @return List of 1 or 2 MBRs + */ + public static List splitAcrossAntimeridian(Mbr mbr) { + List result = new ArrayList<>(); + + if (crossesAntimeridian(mbr)) { + // Split into east and west portions + result.add(new Mbr(mbr.getWest(), mbr.getNorth(), 180.0, mbr.getSouth())); + result.add(new Mbr(-180.0, mbr.getNorth(), mbr.getEast(), mbr.getSouth())); + } else { + result.add(mbr); + } + + return result; + } +} diff --git a/spatial-lib/src/java/cmr/spatial/geometry/PointIntersections.java b/spatial-lib/src/java/cmr/spatial/geometry/PointIntersections.java new file mode 100644 index 0000000000..4f0015a8ab --- /dev/null +++ b/spatial-lib/src/java/cmr/spatial/geometry/PointIntersections.java @@ -0,0 +1,92 @@ +package cmr.spatial.geometry; + +import cmr.spatial.shape.Point; +import cmr.spatial.math.MathUtils; + +/** + * Implements Point-to-Point intersection testing. + * Translates the covers-point? protocol method for Point from Clojure. + */ +public class PointIntersections { + + /** + * Tests if two points are equal (intersect). + * Points are considered equal if they are approximately equal within DELTA tolerance. + * + * @param p1 First point + * @param p2 Second point + * @return true if points are approximately equal + */ + public static boolean pointsIntersect(Point p1, Point p2) { + return pointsApproxEqual(p1, p2, MathUtils.DELTA); + } + + /** + * Tests if two points are approximately equal within the given delta. + * Handles special cases for geodetic (spherical) equality: + * - Latitudes must match within tolerance + * - Longitudes can match directly, or both be on antimeridian, or both at poles + * - All longitudes are equivalent at either pole (±90°) + * - -180 and 180 are equivalent longitudes + * + * @param p1 First point + * @param p2 Second point + * @param delta Tolerance for approximate equality + * @return true if points are approximately equal within tolerance + */ + public static boolean pointsApproxEqual(Point p1, Point p2, double delta) { + double lat1 = p1.getLat(); + double lat2 = p2.getLat(); + double lon1 = p1.getLon(); + double lon2 = p2.getLon(); + + // Latitudes must match + if (!MathUtils.doubleApprox(lat1, lat2, delta)) { + return false; + } + + // Check longitude equivalence + if (MathUtils.doubleApprox(lon1, lon2, delta)) { + return true; + } + + // Check if both on antimeridian (-180 or 180) + if (isOnAntimeridian(lon1) && isOnAntimeridian(lon2)) { + return true; + } + + // Check if both at north pole + if (isNorthPole(lat1, delta) && isNorthPole(lat2, delta)) { + return true; + } + + // Check if both at south pole + if (isSouthPole(lat1, delta) && isSouthPole(lat2, delta)) { + return true; + } + + return false; + } + + /** + * Returns true if latitude is approximately the north pole (90.0). + */ + public static boolean isNorthPole(double lat, double delta) { + return MathUtils.doubleApprox(lat, 90.0, delta); + } + + /** + * Returns true if latitude is approximately the south pole (-90.0). + */ + public static boolean isSouthPole(double lat, double delta) { + return MathUtils.doubleApprox(lat, -90.0, delta); + } + + /** + * Returns true if longitude is on the antimeridian (-180 or 180). + * Uses tolerance-based comparison to match Point.onAntimeridian() logic. + */ + private static boolean isOnAntimeridian(double lon) { + return Math.abs(Math.abs(lon) - 180.0) < MathUtils.DELTA; + } +} diff --git a/spatial-lib/src/java/cmr/spatial/geometry/PointMbrIntersections.java b/spatial-lib/src/java/cmr/spatial/geometry/PointMbrIntersections.java new file mode 100644 index 0000000000..1655817a9b --- /dev/null +++ b/spatial-lib/src/java/cmr/spatial/geometry/PointMbrIntersections.java @@ -0,0 +1,119 @@ +package cmr.spatial.geometry; + +import cmr.spatial.shape.Point; +import cmr.spatial.shape.Mbr; +import cmr.spatial.math.MathUtils; + +/** + * Implements Point-to-Mbr intersection testing. + * Translates the covers-point? and intersects-br? protocol methods from Clojure. + */ +public class PointMbrIntersections { + + /** + * Tests if a point intersects with an MBR. + * A point intersects an MBR if the MBR covers the point. + * + * @param point The point to test + * @param mbr The minimum bounding rectangle + * @return true if the MBR covers the point + */ + public static boolean pointIntersectsMbr(Point point, Mbr mbr) { + return geodetic_covers_point(mbr, point); + } + + /** + * Tests if an MBR covers a point using geodetic (spherical) coordinate system. + * Handles special cases for poles and normal latitude/longitude coverage. + * + * @param mbr The minimum bounding rectangle + * @param point The point to test + * @return true if the MBR covers the point + */ + public static boolean geodetic_covers_point(Mbr mbr, Point point) { + return geodetic_covers_point(mbr, point, MathUtils.COVERS_TOLERANCE); + } + + /** + * Tests if an MBR covers a point using geodetic coordinate system with custom delta. + * + * @param mbr The minimum bounding rectangle + * @param point The point to test + * @param delta Tolerance for coverage checks + * @return true if the MBR covers the point + */ + public static boolean geodetic_covers_point(Mbr mbr, Point point, double delta) { + double lat = point.getLat(); + double lon = point.getLon(); + + // Special case: north pole + if (PointIntersections.isNorthPole(lat, delta)) { + return covers_lat(mbr, 90.0, delta); + } + + // Special case: south pole + if (PointIntersections.isSouthPole(lat, delta)) { + return covers_lat(mbr, -90.0, delta); + } + + // Normal case: check latitude and longitude coverage + return covers_lat(mbr, lat, delta) + && geodetic_lon_range_covers_lon(mbr.getWest(), mbr.getEast(), lon, delta); + } + + /** + * Returns true if the MBR covers the given latitude. + * A latitude is covered if it falls within the MBR's north and south bounds. + * + * @param mbr The minimum bounding rectangle + * @param lat The latitude to check + * @param delta Tolerance + * @return true if the MBR covers the latitude + */ + public static boolean covers_lat(Mbr mbr, double lat, double delta) { + double north = mbr.getNorth() + delta; + double south = mbr.getSouth() - delta; + return lat >= south && lat <= north; + } + + /** + * Returns true if a longitude is within the geodetic range. + * Handles antimeridian crossing: if the range crosses ±180°, it wraps around. + * For a normal range (west < east), the longitude must be between them. + * For a crossing range (west > east), the longitude must be either >= west or <= east. + * + * @param west Western bound of the longitude range + * @param east Eastern bound of the longitude range + * @param lon The longitude to check + * @param delta Tolerance + * @return true if longitude is in range + */ + public static boolean geodetic_lon_range_covers_lon(double west, double east, + double lon, double delta) { + double west_adjusted = west - delta; + double east_adjusted = east + delta; + boolean crosses_antimeridian = west_adjusted > east_adjusted; + + if (crosses_antimeridian) { + // Range wraps around antimeridian + return lon >= west_adjusted || lon <= east_adjusted; + } else { + // Normal range + return lon >= west_adjusted && lon <= east_adjusted; + } + } + + /** + * Helper method: returns true if latitude is north pole. + */ + private static boolean isNorthPole(double lat, double delta) { + return PointIntersections.isNorthPole(lat, delta); + } + + /** + * Helper method: returns true if latitude is south pole. + */ + private static boolean isSouthPole(double lat, double delta) { + return PointIntersections.isSouthPole(lat, delta); + } +} diff --git a/spatial-lib/src/java/cmr/spatial/geometry/PolygonIntersections.java b/spatial-lib/src/java/cmr/spatial/geometry/PolygonIntersections.java new file mode 100644 index 0000000000..c93e89c461 --- /dev/null +++ b/spatial-lib/src/java/cmr/spatial/geometry/PolygonIntersections.java @@ -0,0 +1,243 @@ +package cmr.spatial.geometry; + +import cmr.spatial.shape.LineString; +import cmr.spatial.shape.Mbr; +import cmr.spatial.shape.Point; +import cmr.spatial.shape.Polygon; +import cmr.spatial.shape.Ring; +import cmr.spatial.math.MathUtils; +import java.util.List; + +/** + * Implements Polygon intersection testing. + * Translates covers-point?, covers-br?, intersects-br?, intersects-ring?, intersects-polygon?, + * and intersects-line-string? from polygon.clj. + * + * A polygon consists of an outer boundary ring and optional inner rings (holes). + * Intersection and coverage tests delegate to the ring operations. + */ +public class PolygonIntersections { + + /** + * Tests if a polygon covers a point. + * A polygon covers a point if: + * - The outer boundary ring covers the point, AND + * - NONE of the holes cover the point + * + * @param polygon The polygon to test + * @param point The point to check + * @return true if the polygon covers the point + */ + public static boolean coversPoint(Polygon polygon, Point point) { + Ring boundary = polygon.getBoundary(); + if (boundary == null) { + return false; + } + + // Convert Ring data holder to Java ring implementation + Object javaRing = RingIntersections.createJavaRing(boundary); + + // Outer ring must cover the point + if (!RingIntersections.coversPoint(javaRing, point)) { + return false; + } + + // None of the holes should cover the point + for (Ring hole : polygon.getHoles()) { + Object javaHole = RingIntersections.createJavaRing(hole); + if (RingIntersections.coversPoint(javaHole, point)) { + return false; + } + } + + return true; + } + + /** + * Tests if a polygon covers a minimal bounding rectangle (MBR). + * A polygon covers an MBR if: + * - The outer boundary ring covers the MBR, AND + * - NONE of the holes intersect the MBR + * + * @param polygon The polygon to test + * @param mbr The bounding rectangle + * @return true if the polygon covers the MBR + */ + public static boolean coversBr(Polygon polygon, Mbr mbr) { + Ring boundary = polygon.getBoundary(); + if (boundary == null) { + return false; + } + + // Convert Ring data holder to Java ring implementation + Object javaRing = RingIntersections.createJavaRing(boundary); + + // Outer ring must cover the MBR + if (!RingIntersections.coversBr(javaRing, mbr)) { + return false; + } + + // None of the holes should intersect the MBR + for (Ring hole : polygon.getHoles()) { + Object javaHole = RingIntersections.createJavaRing(hole); + if (RingIntersections.intersectsBr(javaHole, mbr)) { + return false; + } + } + + return true; + } + + /** + * Tests if a polygon intersects a minimal bounding rectangle (MBR). + * A polygon intersects an MBR if: + * - The outer boundary ring intersects the MBR, AND + * - NONE of the holes cover the MBR + * + * @param polygon The polygon to test + * @param mbr The bounding rectangle + * @return true if the polygon intersects the MBR + */ + public static boolean intersectsBr(Polygon polygon, Mbr mbr) { + Ring boundary = polygon.getBoundary(); + if (boundary == null) { + return false; + } + + // Convert Ring data holder to Java ring implementation + Object javaRing = RingIntersections.createJavaRing(boundary); + + // Outer ring must intersect the MBR + if (!RingIntersections.intersectsBr(javaRing, mbr)) { + return false; + } + + // None of the holes should cover the MBR + for (Ring hole : polygon.getHoles()) { + Object javaHole = RingIntersections.createJavaRing(hole); + if (RingIntersections.coversBr(javaHole, mbr)) { + return false; + } + } + + return true; + } + + /** + * Tests if a polygon intersects a ring. + * A polygon intersects a ring if: + * - The outer boundary ring intersects the ring, AND + * - NONE of the holes cover the ring + * + * @param polygon The polygon to test + * @param ring The ring to check + * @return true if the polygon intersects the ring + */ + public static boolean intersectsRing(Polygon polygon, Ring ring) { + Ring boundary = polygon.getBoundary(); + if (boundary == null) { + return false; + } + + // Convert Ring data holders to Java ring implementations + Object javaBoundary = RingIntersections.createJavaRing(boundary); + Object javaRing = RingIntersections.createJavaRing(ring); + + // Outer ring must intersect the other ring + if (!RingIntersections.intersectsRing(javaBoundary, javaRing)) { + return false; + } + + // None of the holes should cover the ring + for (Ring hole : polygon.getHoles()) { + Object javaHole = RingIntersections.createJavaRing(hole); + if (RingIntersections.coversRing(javaHole, javaRing)) { + return false; + } + } + + return true; + } + + /** + * Tests if a polygon intersects a line string. + * A polygon intersects a line string if: + * - The outer boundary ring intersects the line string, AND + * - The line string is not completely covered by any of the holes + * + * @param polygon The polygon to test + * @param lineString The line string to check + * @return true if the polygon intersects the line string + */ + public static boolean intersectsLineString(Polygon polygon, LineString lineString) { + Ring boundary = polygon.getBoundary(); + if (boundary == null) { + return false; + } + + // Convert Ring data holder to Java ring implementation + Object javaRing = RingIntersections.createJavaRing(boundary); + + // Outer ring must intersect the line string + if (!RingIntersections.intersectsLineString(javaRing, lineString)) { + return false; + } + + // None of the holes should cover the line string + for (Ring hole : polygon.getHoles()) { + Object javaHole = RingIntersections.createJavaRing(hole); + if (RingIntersections.coversLineString(javaHole, lineString)) { + return false; + } + } + + return true; + } + + /** + * Tests if a polygon intersects another polygon. + * A polygon intersects another polygon if: + * - The outer boundary ring of poly1 intersects the outer boundary ring of poly2, AND + * - NONE of the holes of poly1 cover the outer boundary of poly2, AND + * - NONE of the holes of poly2 cover the outer boundary of poly1 + * + * @param polygon1 First polygon + * @param polygon2 Second polygon + * @return true if the polygons intersect + */ + public static boolean intersectsPolygon(Polygon polygon1, Polygon polygon2) { + Ring boundary1 = polygon1.getBoundary(); + Ring boundary2 = polygon2.getBoundary(); + + if (boundary1 == null || boundary2 == null) { + return false; + } + + // Convert Ring data holders to Java ring implementations + Object javaBoundary1 = RingIntersections.createJavaRing(boundary1); + Object javaBoundary2 = RingIntersections.createJavaRing(boundary2); + + // Outer rings must intersect + if (!RingIntersections.intersectsRing(javaBoundary1, javaBoundary2)) { + return false; + } + + // None of the holes of poly1 should cover the boundary of poly2 + for (Ring hole : polygon1.getHoles()) { + Object javaHole = RingIntersections.createJavaRing(hole); + if (RingIntersections.coversRing(javaHole, javaBoundary2)) { + return false; + } + } + + // None of the holes of poly2 should cover the boundary of poly1 + for (Ring hole : polygon2.getHoles()) { + Object javaHole = RingIntersections.createJavaRing(hole); + if (RingIntersections.coversRing(javaHole, javaBoundary1)) { + return false; + } + } + + return true; + } +} diff --git a/spatial-lib/src/java/cmr/spatial/geometry/RingIntersections.java b/spatial-lib/src/java/cmr/spatial/geometry/RingIntersections.java new file mode 100644 index 0000000000..075917e24a --- /dev/null +++ b/spatial-lib/src/java/cmr/spatial/geometry/RingIntersections.java @@ -0,0 +1,805 @@ +package cmr.spatial.geometry; + +import cmr.spatial.internal.arc.Arc; +import cmr.spatial.internal.arc.ArcLineSegmentIntersections; +import cmr.spatial.shape.Mbr; +import cmr.spatial.shape.Point; +import cmr.spatial.internal.segment.LineSegment; +import cmr.spatial.math.MathUtils; +import cmr.spatial.internal.ring.CartesianRing; +import cmr.spatial.internal.ring.GeodeticRing; + +import java.util.ArrayList; +import java.util.HashSet; +import java.util.List; +import java.util.Set; + +/** + * Provides high-level ring intersection APIs that work with both GeodeticRing and CartesianRing. + * + * This class implements the core ring intersection algorithms from ring_relations.clj: + * - Ring-to-ring intersections and covers + * - Ring-to-bounding-rectangle intersections and covers + * - Type-safe dispatch for geodetic vs cartesian coordinate systems + * + * Key algorithms: + * - intersectsRing: O(n²) segment intersection check with point containment fallback + * - coversRing: All points covered + no segment intersections + * - intersectsBr: MBR pre-check + point containment + segment-side intersections + * - coversBr: MBR covers + corner points inside + no improper intersections + * + * Implementation based on cmr.spatial.ring-relations Clojure namespace. + */ +public class RingIntersections { + + /** + * Tolerance for approximate equality in MBR covers calculations. + */ + private static final double COVERS_TOLERANCE = MathUtils.COVERS_TOLERANCE; + + // ==================== Ring-to-Ring Intersections ==================== + + /** + * Tests if two rings intersect each other. + * + * Two rings intersect if: + * 1. Any of their segments intersect (checks all pairs - O(n²)), OR + * 2. Ring2's first point is inside ring1, OR + * 3. Ring1's first point is inside ring2 + * + * Type dispatch handles both GeodeticRing (using Arc segments) and CartesianRing + * (using LineSegment segments). + * + * Algorithm from ring_relations.clj:intersects-ring? + * + * @param ring1 First ring (GeodeticRing or CartesianRing) + * @param ring2 Second ring (GeodeticRing or CartesianRing) + * @return true if the rings intersect + * @throws IllegalArgumentException if rings are not GeodeticRing or CartesianRing + */ + public static boolean intersectsRing(Object ring1, Object ring2) { + Object[] segments1 = getSegments(ring1); + Object[] segments2 = getSegments(ring2); + Point[] points1 = getRingPoints(ring1); + Point[] points2 = getRingPoints(ring2); + + // Do any segments intersect? O(n²) check + for (Object seg1 : segments1) { + for (Object seg2 : segments2) { + if (ArcLineSegmentIntersections.intersects(seg1, seg2)) { + return true; + } + } + } + + // Is ring2 inside ring1? Only need to check first point + if (points2.length > 0 && ringCoversPoint(ring1, points2[0])) { + return true; + } + + // Is ring1 inside ring2? Only need to check first point + if (points1.length > 0 && ringCoversPoint(ring2, points1[0])) { + return true; + } + + return false; + } + + /** + * Tests if ring1 completely covers ring2. + * + * A ring covers another ring if: + * 1. All points of ring2 are inside ring1 (uses covers-point? for each), AND + * 2. None of ring2's segments intersect ring1's segments + * + * Algorithm from ring_relations.clj:covers-ring? + * + * @param ring1 Outer ring (GeodeticRing or CartesianRing) + * @param ring2 Inner ring to test (GeodeticRing or CartesianRing) + * @return true if ring1 completely covers ring2 + * @throws IllegalArgumentException if rings are not GeodeticRing or CartesianRing + */ + public static boolean coversRing(Object ring1, Object ring2) { + Point[] points2 = getRingPoints(ring2); + Object[] segments1 = getSegments(ring1); + Object[] segments2 = getSegments(ring2); + + // All points of ring2 must be inside ring1 + for (Point p : points2) { + if (!ringCoversPoint(ring1, p)) { + return false; + } + } + + // No segments of ring2 should intersect segments of ring1 + for (Object seg2 : segments2) { + for (Object seg1 : segments1) { + if (ArcLineSegmentIntersections.intersects(seg1, seg2)) { + return false; + } + } + } + + return true; + } + + // ==================== Ring-to-Point Intersections ==================== + + /** + * Tests if a ring covers (contains) a point. + * + * @param ring The ring (GeodeticRing or CartesianRing) + * @param point The point to test + * @return true if the ring covers the point + * @throws IllegalArgumentException if ring type is not supported + */ + public static boolean coversPoint(Object ring, Point point) { + return ringCoversPoint(ring, point); + } + + // ==================== Ring-to-LineString Intersections ==================== + + /** + * Tests if a ring intersects a line string. + * A ring intersects a line string if: + * - Any line string point is covered by the ring, OR + * - Any line string segment intersects any ring segment + * + * @param ring The ring (GeodeticRing or CartesianRing) + * @param lineString The line string + * @return true if the ring intersects the line string + */ + public static boolean intersectsLineString(Object ring, cmr.spatial.shape.LineString lineString) { + // Get line string points and segments + List linePoints = ordinatesToPoints(lineString.getOrdinates()); + + // Check if any line points are covered by the ring + for (cmr.spatial.shape.Point point : linePoints) { + if (coversPoint(ring, point)) { + return true; + } + } + + // Get segments from both shapes + Object[] ringSegments = getSegments(ring); + Object[] lineSegments = getLineStringSegments(lineString, linePoints); + + // Check if any line segment intersects any ring segment + for (Object lineSeg : lineSegments) { + for (Object ringSeg : ringSegments) { + if (cmr.spatial.internal.arc.ArcLineSegmentIntersections.intersects(lineSeg, ringSeg)) { + return true; + } + } + } + + return false; + } + + /** + * Tests if a ring covers a line string. + * A ring covers a line string if: + * - All line string points are covered by the ring, AND + * - None of the line string segments intersect the ring segments + * (no crossings - the line must stay entirely inside the ring) + * + * @param ring The ring (GeodeticRing or CartesianRing) + * @param lineString The line string + * @return true if the ring covers the line string + */ + public static boolean coversLineString(Object ring, cmr.spatial.shape.LineString lineString) { + // Get line string points and segments + List linePoints = ordinatesToPoints(lineString.getOrdinates()); + + // All line points must be covered by the ring + for (cmr.spatial.shape.Point point : linePoints) { + if (!coversPoint(ring, point)) { + return false; + } + } + + // Get segments from both shapes + Object[] ringSegments = getSegments(ring); + Object[] lineSegments = getLineStringSegments(lineString, linePoints); + + // No line segment should intersect any ring segment + // (If they intersect, the line crosses the ring boundary) + for (Object lineSeg : lineSegments) { + for (Object ringSeg : ringSegments) { + if (cmr.spatial.internal.arc.ArcLineSegmentIntersections.intersects(lineSeg, ringSeg)) { + return false; + } + } + } + + return true; + } + + // ==================== Ring-to-BR Intersections ==================== + + /** + * Tests if a ring intersects a bounding rectangle. + * + * Algorithm from ring_relations.clj:intersects-br? + * + * Pre-check: Ring's MBR must intersect BR (early exit if not) + * + * Special case: If BR is a single point, delegate to coversPoint + * + * General case - ring intersects BR if any of: + * 1. BR covers any points of the ring (geodetic vs cartesian check) + * 2. Ring contains any corner point of the BR (check just one corner) + * 3. Any ring segment intersects any BR side (optimized loop) + * + * @param ring The ring (GeodeticRing or CartesianRing) + * @param br The bounding rectangle + * @return true if the ring intersects the BR + * @throws IllegalArgumentException if ring is not GeodeticRing or CartesianRing + */ + public static boolean intersectsBr(Object ring, Mbr br) { + String coordSys = getCoordinateSystem(ring); + Mbr ringMbr = getRingMbr(ring); + + // Pre-check: do MBRs intersect? + if (!MbrIntersections.mbrsIntersect(coordSys, ringMbr, br)) { + return false; + } + + // Special case: single-point BR + if (isSinglePoint(br)) { + Point point = new Point(br.getWest(), br.getNorth()); + return ringCoversPoint(ring, point); + } + + Point[] ringPoints = getRingPoints(ring); + boolean isGeodetic = "geodetic".equals(coordSys); + + // Does BR cover any ring points? + for (Point p : ringPoints) { + if (isGeodetic) { + if (geodeticCoversPoint(br, p)) { + return true; + } + } else { + if (cartesianCoversPoint(br, p)) { + return true; + } + } + } + + // Does ring contain any BR corner point? (just check one) + Point brCorner = new Point(br.getWest(), br.getNorth()); + if (ringCoversPoint(ring, brCorner)) { + return true; + } + + // Do any ring segments intersect BR sides? + Object[] segments = getSegments(ring); + return linesIntersectsBrSides(segments, br); + } + + /** + * Tests if a ring completely covers a bounding rectangle. + * + * Algorithm from ring_relations.clj:covers-br? + * + * A ring covers a BR if: + * 1. The ring's MBR covers the BR + * 2. All corner points of the BR are inside the ring + * 3. The ring segments do not intersect the BR except at corner/ring points + * + * @param ring The ring (GeodeticRing or CartesianRing) + * @param br The bounding rectangle + * @return true if the ring completely covers the BR + * @throws IllegalArgumentException if ring is not GeodeticRing or CartesianRing + */ + public static boolean coversBr(Object ring, Mbr br) { + String coordSys = getCoordinateSystem(ring); + Mbr ringMbr = getRingMbr(ring); + + // Ring's MBR must cover the BR + if (!coversMbr(coordSys, ringMbr, br)) { + return false; + } + + // All corner points must be inside the ring + List cornerPoints = getCornerPoints(br); + for (Point corner : cornerPoints) { + if (!ringCoversPoint(ring, corner)) { + return false; + } + } + + // Ring segments should not intersect BR except at acceptable points + Set acceptablePoints = new HashSet<>(); + for (Point p : getRingPoints(ring)) { + acceptablePoints.add(p); + } + for (Point p : cornerPoints) { + acceptablePoints.add(p); + } + + List intersections = brIntersections(ring, br); + + // Either no intersections, or all intersections are acceptable points + if (intersections.isEmpty()) { + return true; + } + + for (Point intersection : intersections) { + if (!containsApproximatePoint(acceptablePoints, intersection)) { + return false; + } + } + + return true; + } + + // ==================== Helper Methods ==================== + + /** + * Returns intersection points where ring segments intersect BR sides. + * + * Algorithm from ring_relations.clj:br-intersections + * + * Special case: If BR is a single point, check if point lies on any ring segment + * General case: Find all intersections between ring segments and BR line segments + * + * @param ring The ring (GeodeticRing or CartesianRing) + * @param br The bounding rectangle + * @return List of intersection points + */ + public static List brIntersections(Object ring, Mbr br) { + List result = new ArrayList<>(); + + String coordSys = getCoordinateSystem(ring); + Mbr ringMbr = getRingMbr(ring); + + // Pre-check: do MBRs intersect? + if (!MbrIntersections.mbrsIntersect(coordSys, ringMbr, br)) { + return result; + } + + // Special case: single-point BR + if (isSinglePoint(br)) { + Point point = new Point(br.getWest(), br.getNorth()); + Object[] segments = getSegments(ring); + + for (Object seg : segments) { + if (intersectsPoint(seg, point)) { + result.add(point); + return result; + } + } + return result; + } + + // General case: find all segment-side intersections + Object[] ringSegments = getSegments(ring); + List brSides = LineSegment.mbrToLineSegments(br); + + for (Object ringSeg : ringSegments) { + for (LineSegment brSide : brSides) { + List segmentIntersections = + ArcLineSegmentIntersections.intersections(ringSeg, brSide); + result.addAll(segmentIntersections); + } + } + + return result; + } + + /** + * Optimized check if any ring segment intersects any BR side. + * + * Algorithm from ring_relations.clj:lines-intersects-br-sides? + * + * Avoids full cartesian product iteration by checking each segment against + * all sides with early exit. BR can have up to 6 sides if it crosses the + * antimeridian. + * + * @param segments Array of ring segments (Arc[] or LineSegment[]) + * @param br The bounding rectangle + * @return true if any segment intersects any BR side + */ + public static boolean linesIntersectsBrSides(Object[] segments, Mbr br) { + List sides = LineSegment.mbrToLineSegments(br); + + // Optimized: check each ring segment against all sides with early exit + for (Object segment : segments) { + for (LineSegment side : sides) { + if (ArcLineSegmentIntersections.intersects(segment, side)) { + return true; + } + } + } + + return false; + } + + // ==================== Type Dispatch Helpers ==================== + + /** + * Gets the segments (arcs or line segments) from a ring. + * + * @param ring GeodeticRing or CartesianRing (already converted from Clojure by caller) + * @return Arc[] for geodetic rings, LineSegment[] for cartesian rings (Java objects) + * @throws IllegalArgumentException if ring type is not supported + */ + private static Object[] getSegments(Object ring) { + if (ring instanceof GeodeticRing) { + return ((GeodeticRing) ring).getArcs(); + } else if (ring instanceof CartesianRing) { + return ((CartesianRing) ring).getLineSegments(); + } else { + throw new IllegalArgumentException( + "Ring must be GeodeticRing or CartesianRing (Clojure rings should be converted by caller), got: " + + ring.getClass().getName()); + } + } + + /** + * Gets the coordinate system of a ring. + * + * @param ring GeodeticRing or CartesianRing + * @return "geodetic" or "cartesian" + * @throws IllegalArgumentException if ring type is not supported + */ + private static String getCoordinateSystem(Object ring) { + if (ring instanceof GeodeticRing) { + return "geodetic"; + } else if (ring instanceof CartesianRing) { + return "cartesian"; + } else { + throw new IllegalArgumentException( + "Ring must be GeodeticRing or CartesianRing (Clojure rings should be converted by caller), got: " + + ring.getClass().getName()); + } + } + + /** + * Gets the points array from a ring. + * + * @param ring GeodeticRing or CartesianRing + * @return Array of points defining the ring + * @throws IllegalArgumentException if ring type is not supported + */ + private static Point[] getRingPoints(Object ring) { + if (ring instanceof GeodeticRing) { + return ((GeodeticRing) ring).getPoints(); + } else if (ring instanceof CartesianRing) { + return ((CartesianRing) ring).getPoints(); + } else { + throw new IllegalArgumentException( + "Ring must be GeodeticRing or CartesianRing (Clojure rings should be converted by caller), got: " + + ring.getClass().getName()); + } + } + + /** + * Gets the MBR from a ring. + * + * @param ring GeodeticRing or CartesianRing + * @return The ring's minimum bounding rectangle + * @throws IllegalArgumentException if ring type is not supported + */ + private static Mbr getRingMbr(Object ring) { + if (ring instanceof GeodeticRing) { + return ((GeodeticRing) ring).getMbr(); + } else if (ring instanceof CartesianRing) { + return ((CartesianRing) ring).getMbr(); + } else { + throw new IllegalArgumentException( + "Ring must be GeodeticRing or CartesianRing (Clojure rings should be converted by caller), got: " + + ring.getClass().getName()); + } + } + + /** + * Tests if a ring covers a point using the appropriate algorithm. + * + * @param ring GeodeticRing or CartesianRing + * @param point The point to test + * @return true if the ring covers the point + * @throws IllegalArgumentException if ring type is not supported + */ + private static boolean ringCoversPoint(Object ring, Point point) { + if (ring instanceof GeodeticRing) { + return ((GeodeticRing) ring).coversPoint(point); + } else if (ring instanceof CartesianRing) { + return ((CartesianRing) ring).coversPoint(point); + } else { + throw new IllegalArgumentException( + "Ring must be GeodeticRing or CartesianRing (Clojure rings should be converted by caller), got: " + + ring.getClass().getName()); + } + } + + // ==================== MBR Utility Methods ==================== + + /** + * Tests if an MBR is a single point. + * + * @param mbr The bounding rectangle + * @return true if west==east and north==south + */ + private static boolean isSinglePoint(Mbr mbr) { + return MathUtils.doubleApprox(mbr.getWest(), mbr.getEast(), MathUtils.DELTA) && + MathUtils.doubleApprox(mbr.getNorth(), mbr.getSouth(), MathUtils.DELTA); + } + + /** + * Returns the corner points of an MBR. + * Order: upper-left, upper-right, lower-right, lower-left + * + * @param br The bounding rectangle + * @return List of 4 corner points + */ + private static List getCornerPoints(Mbr br) { + List corners = new ArrayList<>(); + double west = br.getWest(); + double east = br.getEast(); + double north = br.getNorth(); + double south = br.getSouth(); + + corners.add(new Point(west, north)); // Upper left + corners.add(new Point(east, north)); // Upper right + corners.add(new Point(east, south)); // Lower right + corners.add(new Point(west, south)); // Lower left + + return corners; + } + + /** + * Tests if MBR1 covers MBR2 using the specified coordinate system. + * + * @param coordSys "geodetic" or "cartesian" + * @param mbr1 First MBR (should cover mbr2) + * @param mbr2 Second MBR (should be covered) + * @return true if mbr1 covers mbr2 + */ + private static boolean coversMbr(String coordSys, Mbr mbr1, Mbr mbr2) { + List corners = getCornerPoints(mbr2); + + if ("geodetic".equals(coordSys)) { + for (Point corner : corners) { + if (!geodeticCoversPoint(mbr1, corner)) { + return false; + } + } + } else { + for (Point corner : corners) { + if (!cartesianCoversPoint(mbr1, corner)) { + return false; + } + } + } + + return true; + } + + /** + * Tests if an MBR covers a point using geodetic coordinate system. + * Handles poles and antimeridian crossing. + * + * @param mbr The bounding rectangle + * @param p The point + * @return true if the MBR covers the point + */ + private static boolean geodeticCoversPoint(Mbr mbr, Point p) { + double tolerance = COVERS_TOLERANCE; + double lat = p.getLat(); + double lon = p.getLon(); + + // Handle poles + if (MathUtils.doubleApprox(lat, 90.0, tolerance)) { + return coversLat(mbr, 90.0, tolerance); + } + if (MathUtils.doubleApprox(lat, -90.0, tolerance)) { + return coversLat(mbr, -90.0, tolerance); + } + + return coversLat(mbr, lat, tolerance) && + coversLonGeodetic(mbr, lon, tolerance); + } + + /** + * Tests if an MBR covers a point using cartesian coordinate system. + * + * @param mbr The bounding rectangle + * @param p The point + * @return true if the MBR covers the point + */ + private static boolean cartesianCoversPoint(Mbr mbr, Point p) { + double tolerance = COVERS_TOLERANCE; + double lat = p.getLat(); + double lon = p.getLon(); + + return coversLat(mbr, lat, tolerance) && + coversLonCartesian(mbr, lon, tolerance); + } + + /** + * Tests if an MBR covers a latitude value. + * + * @param mbr The bounding rectangle + * @param lat The latitude + * @param tolerance Tolerance for approximate equality + * @return true if south <= lat <= north (within tolerance) + */ + private static boolean coversLat(Mbr mbr, double lat, double tolerance) { + double north = mbr.getNorth() + tolerance; + double south = mbr.getSouth() - tolerance; + return lat >= south && lat <= north; + } + + /** + * Tests if an MBR covers a longitude value using geodetic rules. + * Handles antimeridian crossing and ±180° edge cases. + * + * @param mbr The bounding rectangle + * @param lon The longitude + * @param tolerance Tolerance for approximate equality + * @return true if the MBR covers the longitude + */ + private static boolean coversLonGeodetic(Mbr mbr, double lon, double tolerance) { + double west = mbr.getWest(); + double east = mbr.getEast(); + + // Check for antimeridian crossing + if (west > east) { + // MBR crosses antimeridian + return lon >= (west - tolerance) || lon <= (east + tolerance); + } else if (Math.abs(lon) == 180.0) { + double within180 = 180.0 - tolerance; + return Math.abs(west) >= within180 || Math.abs(east) >= within180; + } else { + return lon >= (west - tolerance) && lon <= (east + tolerance); + } + } + + /** + * Tests if an MBR covers a longitude value using cartesian rules. + * Simple range check without antimeridian handling. + * + * @param mbr The bounding rectangle + * @param lon The longitude + * @param tolerance Tolerance for approximate equality + * @return true if west <= lon <= east (within tolerance) + */ + private static boolean coversLonCartesian(Mbr mbr, double lon, double tolerance) { + double west = mbr.getWest() - tolerance; + double east = mbr.getEast() + tolerance; + return lon >= west && lon <= east; + } + + // ==================== Point/Segment Utilities ==================== + + /** + * Tests if a segment (Arc or LineSegment) intersects a point. + * + * @param segment Arc or LineSegment + * @param point The point to test + * @return true if the segment contains the point + */ + private static boolean intersectsPoint(Object segment, Point point) { + if (segment instanceof Arc) { + Arc arc = (Arc) segment; + // Check if point lies on arc endpoints first + Point west = arc.getWestPoint(); + Point east = arc.getEastPoint(); + + if (pointsApproxEqual(point, west) || pointsApproxEqual(point, east)) { + return true; + } + + // Use precise point-on-arc test instead of MBR approximation + return arc.pointOnArc(point); + + } else if (segment instanceof LineSegment) { + LineSegment ls = (LineSegment) segment; + return ls.pointOnSegment(point); + } else { + throw new IllegalArgumentException( + "Segment must be Arc or LineSegment, got: " + + segment.getClass().getName()); + } + } + + /** + * Checks if two points are approximately equal. + * + * @param p1 First point + * @param p2 Second point + * @return true if points are within tolerance + */ + private static boolean pointsApproxEqual(Point p1, Point p2) { + return PointIntersections.pointsApproxEqual(p1, p2, COVERS_TOLERANCE); + } + + /** + * Checks if a set contains a point within approximate tolerance. + * + * @param points Set of points + * @param target Target point to find + * @return true if set contains an approximately equal point + */ + private static boolean containsApproximatePoint(Set points, Point target) { + for (Point p : points) { + if (pointsApproxEqual(p, target)) { + return true; + } + } + return false; + } + + // ==================== LineString Helper Methods ==================== + + /** + * Converts ordinates to a list of points. + * + * @param ordinates List of coordinates [lon1, lat1, lon2, lat2, ...] + * @return List of Point objects + */ + private static List ordinatesToPoints(List ordinates) { + List points = new ArrayList<>(); + for (int i = 0; i < ordinates.size(); i += 2) { + points.add(new Point(ordinates.get(i), ordinates.get(i + 1))); + } + return points; + } + + /** + * Gets segments from a LineString. + * + * @param lineString The line string + * @param points Pre-computed points (to avoid recomputing) + * @return Array of segments (Arc for geodetic, LineSegment for cartesian) + */ + private static Object[] getLineStringSegments(cmr.spatial.shape.LineString lineString, List points) { + String coordSystem = lineString.getCoordinateSystem(); + + if ("geodetic".equals(coordSystem)) { + List arcs = new ArrayList<>(); + for (int i = 0; i < points.size() - 1; i++) { + try { + arcs.add(Arc.createArc(points.get(i), points.get(i + 1))); + } catch (IllegalArgumentException e) { + // Skip invalid arcs (duplicate or antipodal points) + // This can happen with degenerate linestrings + } + } + return arcs.toArray(); + } else { // cartesian + return LineSegment.pointsToLineSegments(points).toArray(); + } + } + + // ==================== Ring Factory Methods ==================== + + /** + * Converts a Ring data holder to its Java implementation (GeodeticRing or CartesianRing). + * This is used by PolygonIntersections to work with Ring objects from Polygon shapes. + * + * @param ring The Ring data holder with coordinate system and ordinates + * @return GeodeticRing or CartesianRing depending on coordinate system + * @throws IllegalArgumentException if ring is null or coordinate system is invalid + */ + public static Object createJavaRing(cmr.spatial.shape.Ring ring) { + if (ring == null) { + throw new IllegalArgumentException("Ring cannot be null"); + } + + List points = ordinatesToPoints(ring.getOrdinates()); + String coordSystem = ring.getCoordinateSystem(); + + if ("geodetic".equals(coordSystem)) { + return GeodeticRing.createRing(points); + } else if ("cartesian".equals(coordSystem)) { + return CartesianRing.createRing(points); + } else { + throw new IllegalArgumentException( + "Invalid coordinate system: " + coordSystem + " (expected 'geodetic' or 'cartesian')"); + } + } +} diff --git a/spatial-lib/src/java/cmr/spatial/internal/arc/Arc.java b/spatial-lib/src/java/cmr/spatial/internal/arc/Arc.java new file mode 100644 index 0000000000..044d96694e --- /dev/null +++ b/spatial-lib/src/java/cmr/spatial/internal/arc/Arc.java @@ -0,0 +1,671 @@ +package cmr.spatial.internal.arc; + +import cmr.spatial.shape.Mbr; +import cmr.spatial.shape.Point; +import cmr.spatial.math.CoordinateConversion; +import cmr.spatial.math.MathUtils; +import cmr.spatial.math.Vector; +import cmr.spatial.geometry.MbrIntersections; + +/** + * Represents an arc on a sphere between two points along a great circle. + * An arc is the shortest path between two points on a sphere, following a great circle. + * Immutable value class with cached derived information (MBRs, great circle, courses). + * + * The arc is ordered from west to east, and contains precomputed bounding rectangles + * to optimize spatial queries. + */ +public final class Arc { + + /** + * The delta in degrees to use when comparing things approximately. + */ + public static final double APPROXIMATION_DELTA = 0.0000001; + + private final Point westPoint; + private final Point eastPoint; + private final GreatCircle greatCircle; + private final double initialCourse; + private final double endingCourse; + private final Mbr mbr1; + private final Mbr mbr2; // May be null if arc only has one MBR + + /** + * Private constructor. Use createArc() factory method instead. + */ + private Arc(Point westPoint, Point eastPoint, GreatCircle greatCircle, + double initialCourse, double endingCourse, Mbr mbr1, Mbr mbr2) { + this.westPoint = westPoint; + this.eastPoint = eastPoint; + this.greatCircle = greatCircle; + this.initialCourse = initialCourse; + this.endingCourse = endingCourse; + this.mbr1 = mbr1; + this.mbr2 = mbr2; + } + + /** + * Creates an arc from two points, ordering them west-to-east. + * + * @param p1 First point + * @param p2 Second point + * @return Arc from p1 to p2 + * @throws IllegalArgumentException if points are equal or antipodal + */ + public static Arc createArc(Point p1, Point p2) { + if (p1.equals(p2)) { + throw new IllegalArgumentException("Cannot create arc from equal points"); + } + if (areAntipodal(p1, p2)) { + throw new IllegalArgumentException("Cannot create arc from antipodal points"); + } + + // Order points west to east + if (comparePoints(p1, p2) < 0) { + return createArcFromOrderedPoints(p1, p2, p1, p2); + } else { + return createArcFromOrderedPoints(p1, p2, p2, p1); + } + } + + /** + * Creates an arc from points that are already ordered west to east. + * + * @param point1 First point (original order) + * @param point2 Second point (original order) + * @param westPoint Western-most point + * @param eastPoint Eastern-most point + * @return Newly created arc + */ + private static Arc createArcFromOrderedPoints(Point point1, Point point2, + Point westPoint, Point eastPoint) { + GreatCircle gc = calculateGreatCircle(westPoint, eastPoint); + double initialCourse = calculateCourse(point1, point2); + double endingCourse = (calculateCourse(point2, point1) + 180.0) % 360.0; + + Mbr[] mbrs = calculateBoundingRectangles(westPoint, eastPoint, gc, + initialCourse, endingCourse); + + return new Arc(westPoint, eastPoint, gc, initialCourse, endingCourse, + mbrs[0], mbrs.length > 1 ? mbrs[1] : null); + } + + /** + * Calculates the great circle that contains this arc. + */ + private static GreatCircle calculateGreatCircle(Point westPoint, Point eastPoint) { + // Calculate plane vector using cross product + Vector planeVector = CoordinateConversion.lonLatCrossProduct(westPoint, eastPoint) + .normalize(); + Point pvPoint = CoordinateConversion.vectorToPoint(planeVector); + + double plon = pvPoint.getLon(); + double plat = pvPoint.getLat(); + + // Calculate northernmost and southernmost points on the great circle + double northLon = plon; + double northLat = 90.0 - Math.abs(plat); + + double southLon = (plon + 180.0 > 180.0) ? plon - 180.0 : plon + 180.0; + double southLat = Math.abs(plat) - 90.0; + + Point northernmostPoint = new Point(northLon, northLat); + Point southernmostPoint = new Point(southLon, southLat); + + return new GreatCircle(planeVector, northernmostPoint, southernmostPoint); + } + + /** + * Calculates bounding rectangles for the arc. + * May return 1 or 2 MBRs depending on whether the arc crosses a pole. + */ + private static Mbr[] calculateBoundingRectangles(Point westPoint, Point eastPoint, + GreatCircle greatCircle, + double initialCourse, + double endingCourse) { + double wLon = westPoint.getLon(); + double wLat = westPoint.getLat(); + double eLon = eastPoint.getLon(); + double eLat = eastPoint.getLat(); + + // Check if crosses north pole + if (initialCourse == 0.0 && endingCourse == 180.0) { + Mbr br1 = new Mbr(wLon, 90.0, wLon, wLat); + Mbr br2 = new Mbr(eLon, 90.0, eLon, eLat); + return new Mbr[] { br1, br2 }; + } + + // Check if crosses south pole + if (initialCourse == 180.0 && endingCourse == 0.0) { + Mbr br1 = new Mbr(wLon, wLat, wLon, -90.0); + Mbr br2 = new Mbr(eLon, eLat, eLon, -90.0); + return new Mbr[] { br1, br2 }; + } + + // Regular arc (not crossing poles) + double w = wLon; + double e = eLon; + + // If one point is at a pole, west and east longitudes should match + if (isNorthPole(westPoint) || isSouthPole(westPoint)) { + w = e; + } + if (isNorthPole(eastPoint) || isSouthPole(eastPoint)) { + e = w; + } + + // Choose north and south extents + double s = Math.min(wLat, eLat); + double n = Math.max(wLat, eLat); + + // If both on antimeridian, set west and east to same value + boolean bothAntimeridian = (Math.abs(w) == 180.0) && (Math.abs(e) == 180.0); + if (bothAntimeridian) { + w = 180.0; + e = 180.0; + } + + Mbr br = new Mbr(w, n, e, s); + Point northernmost = greatCircle.getNorthernmostPoint(); + Point southernmost = greatCircle.getSouthernmostPoint(); + + // Check if great circle extrema are within the arc's longitude range + if (coversLonGeodetic(br, northernmost.getLon())) { + // Expand north to include northernmost point + return new Mbr[] { new Mbr(br.getWest(), northernmost.getLat(), + br.getEast(), br.getSouth()) }; + } else if (coversLonGeodetic(br, southernmost.getLon())) { + // Expand south to include southernmost point + return new Mbr[] { new Mbr(br.getWest(), br.getNorth(), + br.getEast(), southernmost.getLat()) }; + } else { + return new Mbr[] { br }; + } + } + + /** + * Returns true if the MBR covers the given longitude (geodetic). + */ + private static boolean coversLonGeodetic(Mbr mbr, double lon) { + double west = mbr.getWest(); + double east = mbr.getEast(); + double tolerance = MathUtils.COVERS_TOLERANCE; + + west = west - tolerance; + east = east + tolerance; + boolean crossesAntimeridian = west > east; + + if (crossesAntimeridian) { + return lon >= west || lon <= east; + } else if (Math.abs(lon) == 180.0) { + double within180 = 180.0 - tolerance; + return Math.abs(west) >= within180 || Math.abs(east) >= within180; + } else { + return lon >= west && lon <= east; + } + } + + /** + * Compares two points from west to east. If longitudes are equal, compares south to north. + * Returns negative if p1 is west/south of p2, 0 if equal, positive if p1 is east/north of p2. + */ + private static int comparePoints(Point p1, Point p2) { + double lon1 = p1.getLon(); + double lat1 = p1.getLat(); + double lon2 = p2.getLon(); + double lat2 = p2.getLat(); + + if (lon1 == lon2) { + return Double.compare(lat1, lat2); + } + + // Compare longitudes west to east + return compareLongitudes(lon1, lon2); + } + + /** + * Compares longitudes from west to east. + * Returns -1 if l1 is west of l2, 0 if equal, 1 if l1 is east of l2. + */ + private static int compareLongitudes(double l1, double l2) { + double mod = (l1 - l2) % 360.0; + if (mod < 0) mod += 360.0; + + if (mod == 180.0) { + return Double.compare(l1, l2); + } else if (mod == 0.0) { + return (l1 == 180.0) ? -1 : 1; + } else if (mod < 180.0) { + return 1; // l1 is east of l2 + } else { + return -1; // l1 is west of l2 + } + } + + /** + * Returns true if two points are antipodal (opposite sides of the sphere). + */ + private static boolean areAntipodal(Point p1, Point p2) { + double lat1 = p1.getLat(); + double lat2 = p2.getLat(); + double lon1 = p1.getLon(); + double lon2 = p2.getLon(); + + if (isNorthPole(p1)) return isSouthPole(p2); + if (isSouthPole(p1)) return isNorthPole(p2); + + return (lat1 == -lat2) && (lon1 == antipodalLon(lon2)); + } + + /** + * Returns the longitude on the opposite side of the earth. + */ + private static double antipodalLon(double lon) { + double newLon = lon + 180.0; + return (newLon > 180.0) ? newLon - 360.0 : newLon; + } + + /** + * Calculates the initial bearing/course from p1 to p2 in degrees. + * 0 = north, 90 = east, 180 = south, 270 = west. + */ + private static double calculateCourse(Point p1, Point p2) { + double lon1 = p1.getLon(); + double lat1 = p1.getLat(); + double lon2 = p2.getLon(); + double lat2 = p2.getLat(); + + // Handle poles + if (isNorthPole(p2)) return 0.0; + if (isNorthPole(p1)) return 180.0; + if (isSouthPole(p2)) return 180.0; + if (isSouthPole(p1)) return 0.0; + + // Handle vertical lines + if (lon1 == lon2) { + return (lat1 > lat2) ? 180.0 : 0.0; + } + + // Check if crossing a pole (180 degrees longitude difference) + if (Math.abs(lon1 - lon2) == 180.0) { + double midLat = (lat1 + lat2) / 2.0; + return (midLat > 0.0) ? 0.0 : 180.0; + } + + // General case: calculate bearing + double lon1Rad = Math.toRadians(lon1); + double lat1Rad = Math.toRadians(lat1); + double lon2Rad = Math.toRadians(lon2); + double lat2Rad = Math.toRadians(lat2); + + double lonDiff = lon2Rad - lon1Rad; + double cosLat2 = Math.cos(lat2Rad); + double y = Math.sin(lonDiff) * cosLat2; + double x = Math.cos(lat1Rad) * Math.sin(lat2Rad) - + Math.sin(lat1Rad) * cosLat2 * Math.cos(lonDiff); + + double normalized = Math.toDegrees(Math.atan2(y, x)); + return (-normalized + 360.0) % 360.0; + } + + /** + * Returns true if point is at or very near the north pole. + */ + private static boolean isNorthPole(Point p) { + return MathUtils.doubleApprox(p.getLat(), 90.0, APPROXIMATION_DELTA); + } + + /** + * Returns true if point is at or very near the south pole. + */ + private static boolean isSouthPole(Point p) { + return MathUtils.doubleApprox(p.getLat(), -90.0, APPROXIMATION_DELTA); + } + + /** + * Returns true if point is at either pole. + */ + private static boolean isPole(Point p) { + return isNorthPole(p) || isSouthPole(p); + } + + // Getters + + public Point getWestPoint() { + return westPoint; + } + + public Point getEastPoint() { + return eastPoint; + } + + public GreatCircle getGreatCircle() { + return greatCircle; + } + + public double getInitialCourse() { + return initialCourse; + } + + public double getEndingCourse() { + return endingCourse; + } + + /** + * Returns the minimum bounding rectangles for this arc. + * May return 1 or 2 MBRs depending on arc geometry. + */ + public Mbr[] getMbrs() { + if (mbr2 == null) { + return new Mbr[] { mbr1 }; + } else { + return new Mbr[] { mbr1, mbr2 }; + } + } + + /** + * Returns the first minimum bounding rectangle. + */ + public Mbr getMbr1() { + return mbr1; + } + + /** + * Returns the second minimum bounding rectangle, or null if arc has only one MBR. + */ + public Mbr getMbr2() { + return mbr2; + } + + /** + * Returns true if the arc crosses the north pole. + */ + public boolean crossesNorthPole() { + return initialCourse == 0.0 && endingCourse == 180.0; + } + + /** + * Returns true if the arc crosses the south pole. + */ + public boolean crossesSouthPole() { + return initialCourse == 180.0 && endingCourse == 0.0; + } + + /** + * Returns true if the arc crosses the antimeridian (±180° longitude). + */ + public boolean crossesAntimeridian() { + return crossesAntimeridianMbr(mbr1) || + (mbr2 != null && crossesAntimeridianMbr(mbr2)); + } + + /** + * Returns true if an MBR crosses the antimeridian. + */ + private static boolean crossesAntimeridianMbr(Mbr mbr) { + return mbr.getWest() > mbr.getEast(); + } + + /** + * Returns true if the arc is vertical (north-south). + * Crossing a pole is considered vertical. + */ + public boolean isVertical() { + double wLon = westPoint.getLon(); + double eLon = eastPoint.getLon(); + + // Longitudes are equal + if (MathUtils.doubleApprox(wLon, eLon, APPROXIMATION_DELTA)) { + return true; + } + + // Crosses a pole + if (crossesNorthPole() || crossesSouthPole()) { + return true; + } + + // One point is a pole (but not both) + boolean westIsPole = isPole(westPoint); + boolean eastIsPole = isPole(eastPoint); + if (westIsPole != eastIsPole) { + return true; + } + + // Both on antimeridian + if (MathUtils.doubleApprox(180.0, Math.abs(wLon), APPROXIMATION_DELTA) && + MathUtils.doubleApprox(180.0, Math.abs(eLon), APPROXIMATION_DELTA)) { + return true; + } + + return false; + } + + /** + * Returns true if the given point lies on this arc. + * + * @param point Point to test + * @return true if point is on the arc + */ + public boolean pointOnArc(Point point) { + // Check if point is in one of the MBRs + if (!coversPointGeodetic(mbr1, point) && + (mbr2 == null || !coversPointGeodetic(mbr2, point))) { + return false; + } + + // Check if point is one of the endpoints + if (pointsApproxEqual(westPoint, point, APPROXIMATION_DELTA) || + pointsApproxEqual(eastPoint, point, APPROXIMATION_DELTA)) { + return true; + } + + // If arc is vertical and point is in MBR, it's on the arc + if (isVertical()) { + return true; + } + + // Find point on arc at the given longitude and check if it matches + Point pointAtLon = pointAtLon(point.getLon()); + return pointAtLon != null && + pointsApproxEqual(point, pointAtLon, APPROXIMATION_DELTA); + } + + /** + * Returns the point on the arc at the given longitude. + * Returns null if the longitude is outside the arc's bounds. + * Does not work for vertical arcs. + * + * @param lon Longitude in degrees + * @return Point on arc at given longitude, or null if outside bounds + */ + public Point pointAtLon(double lon) { + if (isVertical()) { + throw new IllegalStateException("Cannot call pointAtLon on vertical arc"); + } + + // Check if longitude is within arc bounds + boolean inBounds = coversLonGeodetic(mbr1, lon) || + (mbr2 != null && coversLonGeodetic(mbr2, lon)) || + crossesSouthPole() || crossesNorthPole(); + + if (!inBounds) { + return null; + } + + double lonRad = Math.toRadians(lon); + double lon1 = Math.toRadians(westPoint.getLon()); + double lon2 = Math.toRadians(eastPoint.getLon()); + double lat1 = Math.toRadians(westPoint.getLat()); + double lat2 = Math.toRadians(eastPoint.getLat()); + + double cosLat1 = Math.cos(lat1); + double cosLat2 = Math.cos(lat2); + + // From http://williams.best.vwh.net/avform.htm#Int + double top = Math.sin(lat1) * cosLat2 * Math.sin(lonRad - lon2) - + Math.sin(lat2) * cosLat1 * Math.sin(lonRad - lon1); + double bottom = cosLat1 * cosLat2 * Math.sin(lon1 - lon2); + double latRad = Math.atan(top / bottom); + + return new Point(Math.toDegrees(lonRad), Math.toDegrees(latRad)); + } + + /** + * Returns the points where the arc crosses the given latitude. + * Returns null if the arc doesn't cross that latitude. + * + * @param lat Latitude in degrees + * @return Array of points (0, 1, or 2 points), or null if no intersections + */ + public Point[] pointsAtLat(double lat) { + // Check if latitude is within arc bounds + if (!coversLatMbr(mbr1, lat) && (mbr2 == null || !coversLatMbr(mbr2, lat))) { + return null; + } + + double lat3 = Math.toRadians(lat); + double lon1 = Math.toRadians(westPoint.getLon()); + double lon2 = Math.toRadians(eastPoint.getLon()); + double lat1 = Math.toRadians(westPoint.getLat()); + double lat2 = Math.toRadians(eastPoint.getLat()); + + double lon12 = lon1 - lon2; + double sinLon12 = Math.sin(lon12); + double cosLon12 = Math.cos(lon12); + double sinLat1 = Math.sin(lat1); + double cosLat1 = Math.cos(lat1); + double sinLat2 = Math.sin(lat2); + double cosLat2 = Math.cos(lat2); + double sinLat3 = Math.sin(lat3); + double cosLat3 = Math.cos(lat3); + + // From http://williams.best.vwh.net/avform.htm#Par + double A = sinLat1 * cosLat2 * cosLat3 * sinLon12; + double B = sinLat1 * cosLat2 * cosLat3 * cosLon12 - cosLat1 * sinLat2 * cosLat3; + double C = cosLat1 * cosLat2 * sinLat3 * sinLon12; + double h = Math.sqrt(A * A + B * B); + + if (Math.abs(C) > h) { + return null; // No intersection + } + + double lonRad = Math.atan2(B, A); + double dlon = Math.acos(C / h); + + double lon31 = ((lon1 + dlon + lonRad + Math.PI) % (2 * Math.PI)) - Math.PI; + double lon32 = ((lon1 - dlon + lonRad + Math.PI) % (2 * Math.PI)) - Math.PI; + + Point p1 = new Point(Math.toDegrees(lon31), lat); + Point p2 = new Point(Math.toDegrees(lon32), lat); + + // Check if points are within MBRs + boolean p1InMbr = coversPointGeodetic(mbr1, p1) || + (mbr2 != null && coversPointGeodetic(mbr2, p1)); + boolean p2InMbr = coversPointGeodetic(mbr1, p2) || + (mbr2 != null && coversPointGeodetic(mbr2, p2)); + + if (p1InMbr && p2InMbr) { + if (pointsApproxEqual(p1, p2, APPROXIMATION_DELTA)) { + return new Point[] { p1 }; + } else { + return new Point[] { p1, p2 }; + } + } else if (p1InMbr) { + return new Point[] { p1 }; + } else if (p2InMbr) { + return new Point[] { p2 }; + } else { + return null; + } + } + + /** + * Returns true if MBR covers the given latitude. + */ + private static boolean coversLatMbr(Mbr mbr, double lat) { + double tolerance = MathUtils.COVERS_TOLERANCE; + double north = mbr.getNorth() + tolerance; + double south = mbr.getSouth() - tolerance; + return lat >= south && lat <= north; + } + + /** + * Returns true if MBR covers the given point (geodetic). + */ + private static boolean coversPointGeodetic(Mbr mbr, Point p) { + double tolerance = MathUtils.COVERS_TOLERANCE; + double lat = p.getLat(); + + // Handle poles + if (isNorthPole(p)) { + return coversLatMbr(mbr, 90.0); + } + if (isSouthPole(p)) { + return coversLatMbr(mbr, -90.0); + } + + return coversLatMbr(mbr, lat) && coversLonGeodetic(mbr, p.getLon()); + } + + /** + * Returns true if two points are approximately equal within the given delta. + */ + private static boolean pointsApproxEqual(Point p1, Point p2, double delta) { + double lat1 = p1.getLat(); + double lat2 = p2.getLat(); + double lon1 = p1.getLon(); + double lon2 = p2.getLon(); + + if (!MathUtils.doubleApprox(lat1, lat2, delta)) { + return false; + } + + if (MathUtils.doubleApprox(lon1, lon2, delta)) { + return true; + } + + // Check if both on antimeridian + if (Math.abs(lon1) == 180.0 && Math.abs(lon2) == 180.0) { + return true; + } + + // Check if both at poles + if (isNorthPole(p1) && isNorthPole(p2)) { + return true; + } + if (isSouthPole(p1) && isSouthPole(p2)) { + return true; + } + + return false; + } + + @Override + public String toString() { + return String.format("Arc{west=%s, east=%s, courses=[%.1f, %.1f]}", + westPoint, eastPoint, initialCourse, endingCourse); + } + + @Override + public boolean equals(Object obj) { + if (this == obj) return true; + if (!(obj instanceof Arc)) return false; + Arc other = (Arc) obj; + return westPoint.equals(other.westPoint) + && eastPoint.equals(other.eastPoint) + && Double.compare(initialCourse, other.initialCourse) == 0 + && Double.compare(endingCourse, other.endingCourse) == 0; + } + + @Override + public int hashCode() { + int result = 17; + result = 31 * result + westPoint.hashCode(); + result = 31 * result + eastPoint.hashCode(); + long temp = Double.doubleToLongBits(initialCourse); + result = 31 * result + (int)(temp ^ (temp >>> 32)); + temp = Double.doubleToLongBits(endingCourse); + result = 31 * result + (int)(temp ^ (temp >>> 32)); + return result; + } +} diff --git a/spatial-lib/src/java/cmr/spatial/internal/arc/ArcLineSegmentIntersections.java b/spatial-lib/src/java/cmr/spatial/internal/arc/ArcLineSegmentIntersections.java new file mode 100644 index 0000000000..8fabc22ad8 --- /dev/null +++ b/spatial-lib/src/java/cmr/spatial/internal/arc/ArcLineSegmentIntersections.java @@ -0,0 +1,723 @@ +package cmr.spatial.internal.arc; + +import cmr.spatial.shape.Mbr; +import cmr.spatial.shape.Point; +import cmr.spatial.internal.segment.LineSegment; +import cmr.spatial.math.CoordinateConversion; +import cmr.spatial.math.MathUtils; +import cmr.spatial.math.Vector; +import cmr.spatial.geometry.MbrIntersections; + +import java.util.ArrayList; +import java.util.HashSet; +import java.util.List; +import java.util.Set; + +/** + * Provides intersection functions for finding the intersection of spherical arcs and cartesian line segments. + * + * This class handles the complex problem of intersecting geodetic (spherical) arcs with + * cartesian (flat plane) line segments. The key challenge is that these geometries exist + * in different coordinate systems. + * + * Key algorithms: + * - Arc-Arc: Great circle plane intersection using cross product + * - LineSegment-LineSegment: Cartesian line equation solving + * - Arc-LineSegment (mixed): Densification strategy to bridge coordinate systems + * + * Based on cmr.spatial.arc-line-segment-intersections Clojure namespace. + */ +public class ArcLineSegmentIntersections { + + /** + * Tolerance for approximate equality comparisons. + */ + private static final double APPROXIMATION_DELTA = Arc.APPROXIMATION_DELTA; + + /** + * Default densification distance in degrees for converting line segments to arcs. + */ + private static final double DEFAULT_DENSIFICATION_DISTANCE = 0.1; + + /** + * North pole point. + */ + public static final Point NORTH_POLE = new Point(0.0, 90.0); + + /** + * South pole point. + */ + public static final Point SOUTH_POLE = new Point(0.0, -90.0); + + // ==================== Arc-Arc Intersections ==================== + + /** + * Returns the intersection points of two arcs. + * + * Algorithm: + * 1. Check if MBRs intersect (early exit if not) + * 2. Handle special case: arcs on the same great circle + * 3. Calculate plane vector cross product to find intersection line + * 4. Two antipodal points lie on this intersection line + * 5. Filter points within both arcs' bounding rectangles + * + * @param a1 First arc + * @param a2 Second arc + * @return List of intersection points (0, 1, or 2 points) + */ + public static List arcArcIntersections(Arc a1, Arc a2) { + List result = new ArrayList<>(); + + // Check if any of the arc MBRs intersect + if (!arcMbrsIntersect(a1, a2)) { + return result; + } + + GreatCircle gc1 = a1.getGreatCircle(); + GreatCircle gc2 = a2.getGreatCircle(); + + // Special case: both arcs on the same great circle + if (gc1.isEquivalent(gc2)) { + return greatCircleEquivalencyArcIntersections(a1, a2); + } + + // Default case: compute intersection using cross product + Vector pv1 = gc1.getPlaneVector(); + Vector pv2 = gc2.getPlaneVector(); + + // Compute the great circle intersection vector + // This is the cross product of the vectors defining the great circle planes + Vector intersectionVector = pv1.crossProduct(pv2); + + // Convert to points (two antipodal points on the intersection line) + Point intersectionPoint1 = CoordinateConversion.vectorToPoint(intersectionVector); + Point intersectionPoint2 = antipodal(intersectionPoint1); + + // Check if points are within bounding rectangles of both arcs + boolean point1Within = pointWithinArcBoundingRectangles(intersectionPoint1, a1, a2); + boolean point2Within = pointWithinArcBoundingRectangles(intersectionPoint2, a1, a2); + + if (point1Within) { + result.add(intersectionPoint1); + } + if (point2Within) { + result.add(intersectionPoint2); + } + + return result; + } + + /** + * Special case arc intersections for both arcs having the same great circle. + * Returns the endpoint intersections that fall within both arcs' bounds. + */ + private static List greatCircleEquivalencyArcIntersections(Arc a1, Arc a2) { + List result = new ArrayList<>(); + Set points = new HashSet<>(); + + // Collect all endpoints + points.add(a1.getWestPoint()); + points.add(a1.getEastPoint()); + points.add(a2.getWestPoint()); + points.add(a2.getEastPoint()); + + // Filter points that are within both arcs' bounding rectangles + for (Point p : points) { + if (pointWithinArcBoundingRectangles(p, a1, a2)) { + result.add(p); + } + } + + return result; + } + + /** + * Returns true if any of the arc MBRs intersect. + */ + private static boolean arcMbrsIntersect(Arc a1, Arc a2) { + Mbr[] a1Mbrs = a1.getMbrs(); + Mbr[] a2Mbrs = a2.getMbrs(); + + for (Mbr mbr1 : a1Mbrs) { + for (Mbr mbr2 : a2Mbrs) { + if (MbrIntersections.mbrsIntersect("geodetic", mbr1, mbr2)) { + return true; + } + } + } + return false; + } + + /** + * Returns true if the point is within the bounding rectangles of both arcs. + */ + private static boolean pointWithinArcBoundingRectangles(Point p, Arc a1, Arc a2) { + boolean inA1 = false; + boolean inA2 = false; + + for (Mbr mbr : a1.getMbrs()) { + if (geodeticCoversPoint(mbr, p)) { + inA1 = true; + break; + } + } + + for (Mbr mbr : a2.getMbrs()) { + if (geodeticCoversPoint(mbr, p)) { + inA2 = true; + break; + } + } + + return inA1 && inA2; + } + + // ==================== LineSegment-LineSegment Intersections ==================== + + /** + * Returns the intersection points of two line segments. + * + * Handles all special cases: + * - Both vertical: check longitude match and latitude overlap + * - Both horizontal: check latitude match and longitude overlap + * - One vertical, one horizontal: single point intersection + * - Parallel lines (same slope): check if on same line + * - Normal case: solve m1*x + b1 = m2*x + b2 + * + * @param ls1 First line segment + * @param ls2 Second line segment + * @return List of intersection points (0 or 1 point) + */ + public static List lineSegmentLineSegmentIntersections(LineSegment ls1, LineSegment ls2) { + List result = new ArrayList<>(); + Point intersection = ls1.intersection(ls2); + if (intersection != null) { + result.add(intersection); + } + return result; + } + + // ==================== Arc-LineSegment Intersections ==================== + + /** + * Returns the intersection points of an arc and a line segment. + * + * This is the most complex case as it involves geodetic (arc) and cartesian (line segment) + * coordinate systems. Strategy: + * + * 1. Pre-check: Do MBRs intersect? (early exit if not) + * 2. Handle pole intersections + * 3. Vertical line segment: convert to arc, use arc-arc intersection + * 4. Horizontal line segment: use arc.pointsAtLat() method + * 5. Vertical arc: convert to line segments at pole, use line-line intersection + * 6. General case: densify line segment, convert to arcs, find arc-arc intersections + * + * @param arc The arc + * @param ls The line segment + * @return List of intersection points + */ + public static List arcLineSegmentIntersections(Arc arc, LineSegment ls) { + List result = new ArrayList<>(); + + Mbr lsMbr = ls.getMbr(); + Mbr[] arcMbrs = arc.getMbrs(); + + // Check if MBRs intersect + boolean mbrsIntersect = false; + List intersectingMbrs = new ArrayList<>(); + + for (Mbr arcMbr : arcMbrs) { + if (MbrIntersections.mbrsIntersect("geodetic", lsMbr, arcMbr)) { + mbrsIntersect = true; + intersectingMbrs.add(arcMbr); + } + } + + if (!mbrsIntersect) { + return result; + } + + Point lsPoint1 = ls.getPoint1(); + Point lsPoint2 = ls.getPoint2(); + Point arcWest = arc.getWestPoint(); + Point arcEast = arc.getEastPoint(); + + // Handle pole intersections + if ((isNorthPole(arcWest) || isNorthPole(arcEast)) && + (isNorthPole(lsPoint1) || isNorthPole(lsPoint2))) { + result.add(NORTH_POLE); + } + + if ((isSouthPole(arcWest) || isSouthPole(arcEast)) && + (isSouthPole(lsPoint1) || isSouthPole(lsPoint2))) { + result.add(SOUTH_POLE); + } + + // Vertical line segment: treat as a vertical arc + if (ls.isVertical()) { + try { + Arc lsArc = Arc.createArc(lsPoint1, lsPoint2); + result.addAll(arcArcIntersections(arc, lsArc)); + } catch (IllegalArgumentException e) { + // Handle degenerate endpoints + if (lsPoint1.equals(lsPoint2)) { + // Line segment is a single point - check if it's on the arc + if (arc.pointOnArc(lsPoint1)) { + result.add(lsPoint1); + } + } else { + // Antipodal points: the segment spans a full meridian + // This is complex - for now, return empty (could be improved) + } + } + return removeDuplicatePoints(result); + } + + // Horizontal line segment: use arc latitude segment intersection + if (ls.isHorizontal()) { + double lat = lsPoint1.getLat(); + double lon1 = lsPoint1.getLon(); + double lon2 = lsPoint2.getLon(); + double west = Math.min(lon1, lon2); + double east = Math.max(lon1, lon2); + + result.addAll(latSegmentIntersections(arc, lat, west, east)); + return removeDuplicatePoints(result); + } + + // Vertical arc: convert to line segments + if (arc.isVertical()) { + result.addAll(verticalArcLineSegmentIntersections(arc, ls)); + return removeDuplicatePoints(result); + } + + // General case: densify line segment and convert to arcs + result.addAll(lineSegmentArcIntersectionsWithDensification(ls, arc, intersectingMbrs)); + return removeDuplicatePoints(result); + } + + /** + * Returns the points where an arc intersects a latitude segment. + * The latitude segment is defined at lat between lon-west and lon-east. + */ + private static List latSegmentIntersections(Arc arc, double lat, double lonWest, double lonEast) { + List result = new ArrayList<>(); + + // Find the points where the arc crosses that latitude (if any) + Point[] points = arc.pointsAtLat(lat); + if (points == null) { + return result; + } + + // Create MBR for the latitude segment + Mbr latSegMbr = new Mbr(lonWest, lat, lonEast, lat); + Mbr[] arcMbrs = arc.getMbrs(); + + // Filter points that are within the lat segment and arc MBRs + for (Point p : points) { + if (geodeticCoversPoint(latSegMbr, p)) { + for (Mbr arcMbr : arcMbrs) { + if (geodeticCoversPoint(arcMbr, p)) { + result.add(p); + break; + } + } + } + } + + return result; + } + + /** + * Determines the intersection points of a vertical arc and a line segment. + * Converts the arc into equivalent line segments (handling pole crossings). + */ + private static List verticalArcLineSegmentIntersections(Arc arc, LineSegment ls) { + List result = new ArrayList<>(); + + Point westPoint = arc.getWestPoint(); + Point eastPoint = arc.getEastPoint(); + List arcSegments = new ArrayList<>(); + + // A vertical arc could cross a pole. It gets divided in half at the pole in that case. + if (arc.crossesNorthPole()) { + arcSegments.add(LineSegment.createLineSegment( + westPoint, new Point(westPoint.getLon(), 90.0))); + arcSegments.add(LineSegment.createLineSegment( + eastPoint, new Point(eastPoint.getLon(), 90.0))); + } else if (arc.crossesSouthPole()) { + arcSegments.add(LineSegment.createLineSegment( + westPoint, new Point(westPoint.getLon(), -90.0))); + arcSegments.add(LineSegment.createLineSegment( + eastPoint, new Point(eastPoint.getLon(), -90.0))); + } else if (isNorthPole(eastPoint)) { + // Create a vertical line segment ignoring the original point2 lon + arcSegments.add(LineSegment.createLineSegment( + westPoint, new Point(westPoint.getLon(), 90.0))); + } else if (isNorthPole(westPoint)) { + // Create a vertical line segment ignoring the original point1 lon + arcSegments.add(LineSegment.createLineSegment( + eastPoint, new Point(eastPoint.getLon(), 90.0))); + } else if (isSouthPole(eastPoint)) { + arcSegments.add(LineSegment.createLineSegment( + westPoint, new Point(westPoint.getLon(), -90.0))); + } else if (isSouthPole(westPoint)) { + arcSegments.add(LineSegment.createLineSegment( + eastPoint, new Point(eastPoint.getLon(), -90.0))); + } else { + arcSegments.add(LineSegment.createLineSegment(westPoint, eastPoint)); + } + + // Find intersections with each segment + for (LineSegment arcSeg : arcSegments) { + Point intersection = ls.intersection(arcSeg); + if (intersection != null) { + result.add(intersection); + } + } + + return result; + } + + /** + * Performs the intersection between a line segment and an arc using densification. + * + * Densification strategy: + * 1. Subselect line segment to only portions within arc MBRs + * 2. Densify those portions into many points + * 3. Convert consecutive points into arcs + * 4. Find arc-arc intersections + */ + private static List lineSegmentArcIntersectionsWithDensification( + LineSegment ls, Arc arc, List intersectingMbrs) { + List result = new ArrayList<>(); + + // Compute the intersections of the intersecting MBRs + List mbrIntersections = new ArrayList<>(); + Mbr lsMbr = ls.getMbr(); + + for (Mbr arcMbr : intersectingMbrs) { + Mbr intersection = computeMbrIntersection(lsMbr, arcMbr); + if (intersection != null) { + mbrIntersections.add(intersection); + } + } + + // For each intersecting MBR, subselect the line segment and densify + for (Mbr mbr : mbrIntersections) { + // Simple subselection: check if endpoints are in MBR + Point p1 = ls.getPoint1(); + Point p2 = ls.getPoint2(); + + boolean p1In = cartesianCoversPoint(mbr, p1); + boolean p2In = cartesianCoversPoint(mbr, p2); + + LineSegment segmentToProcess = ls; + + // If only one point is in, we might need to clip + // For simplicity, we'll process the whole segment if any part intersects + if (!p1In && !p2In) { + // Neither point in MBR, but MBR intersects - need to find intersection points + // This is complex, so we'll just densify the whole segment + } + + // Densify the line segment + List densifiedPoints = densifyLineSegment(segmentToProcess, DEFAULT_DENSIFICATION_DISTANCE); + + // Convert to arcs + for (int i = 0; i < densifiedPoints.size() - 1; i++) { + Point arcP1 = densifiedPoints.get(i); + Point arcP2 = densifiedPoints.get(i + 1); + + // Skip if points are equal or too close + if (pointsApproxEqual(arcP1, arcP2, APPROXIMATION_DELTA)) { + continue; + } + + // Try to create arc and find intersections + try { + Arc densifiedArc = Arc.createArc(arcP1, arcP2); + List arcIntersections = arcArcIntersections(arc, densifiedArc); + result.addAll(arcIntersections); + } catch (IllegalArgumentException e) { + // Skip if arc creation fails (antipodal points, etc.) + } + } + + // Also check if any densified points lie on the arc + for (Point p : densifiedPoints) { + if (arc.pointOnArc(p)) { + result.add(p); + } + } + } + + // Remove duplicates + return removeDuplicatePoints(result); + } + + /** + * Computes the intersection of two MBRs. + * + * Uses geodetic semantics to handle the case where Arc and LineSegment MBRs + * touch at a pole but have disjoint longitude ranges. In geodetic coordinates, + * MBRs that both touch the same pole are considered to intersect at that pole. + * + * Both input MBRs are guaranteed to be non-crossing (Arc pre-splits + * antimeridian-crossing MBRs, LineSegment never crosses by design), so the + * simple Math.max/min computation is safe. + * + * Returns null if they don't intersect. + */ + private static Mbr computeMbrIntersection(Mbr mbr1, Mbr mbr2) { + if (!MbrIntersections.mbrsIntersect("geodetic", mbr1, mbr2)) { + return null; + } + + double west = Math.max(mbr1.getWest(), mbr2.getWest()); + double east = Math.min(mbr1.getEast(), mbr2.getEast()); + double north = Math.min(mbr1.getNorth(), mbr2.getNorth()); + double south = Math.max(mbr1.getSouth(), mbr2.getSouth()); + + return new Mbr(west, north, east, south); + } + + /** + * Returns points along the line segment for approximating the segment. + * Does no densification for vertical lines. + */ + private static List densifyLineSegment(LineSegment ls, double densificationDist) { + List result = new ArrayList<>(); + + if (ls.isVertical()) { + result.add(ls.getPoint1()); + result.add(ls.getPoint2()); + return result; + } + + Point p1 = ls.getPoint1(); + Point p2 = ls.getPoint2(); + double lon1 = p1.getLon(); + double lat1 = p1.getLat(); + double lon2 = p2.getLon(); + double lat2 = p2.getLat(); + + double m = ls.getM(); + + // Convert slope to angle + double angleA = Math.atan(m); + + // Calculate the difference to add for each point + double latDiff = densificationDist * Math.sin(angleA); + double lonDiff = densificationDist * Math.cos(angleA); + + // If the line is going backwards, flip the signs + if (lon1 > lon2) { + latDiff = -latDiff; + lonDiff = -lonDiff; + } + + // Calculate distance + double distance = Math.sqrt(Math.pow(lat2 - lat1, 2) + Math.pow(lon2 - lon1, 2)); + int numPoints = (int) Math.floor(distance / densificationDist); + + // Generate points + for (int i = 0; i <= numPoints; i++) { + double lon = lon1 + (lonDiff * i); + double lat = lat1 + (latDiff * i); + result.add(new Point(lon, lat)); + } + + // Ensure endpoint is included + if (!pointsApproxEqual(result.get(result.size() - 1), p2, APPROXIMATION_DELTA)) { + result.add(p2); + } + + return result; + } + + // ==================== Polymorphic Dispatcher ==================== + + /** + * Determines if line 1 and line 2 intersect. + * A line can be an arc or a line segment. + * + * @param seg1 First line (Arc or LineSegment) + * @param seg2 Second line (Arc or LineSegment) + * @return List of intersection points + */ + public static List intersections(Object seg1, Object seg2) { + if (seg1 == null || seg2 == null) { + if (seg1 == null && seg2 == null) { + throw new IllegalArgumentException("seg1 and seg2 cannot be null"); + } else if (seg1 == null) { + throw new IllegalArgumentException("seg1 cannot be null"); + } else { + throw new IllegalArgumentException("seg2 cannot be null"); + } + } + if (seg1 instanceof Arc && seg2 instanceof Arc) { + return arcArcIntersections((Arc) seg1, (Arc) seg2); + } else if (seg1 instanceof LineSegment && seg2 instanceof LineSegment) { + return lineSegmentLineSegmentIntersections((LineSegment) seg1, (LineSegment) seg2); + } else if (seg1 instanceof Arc && seg2 instanceof LineSegment) { + return arcLineSegmentIntersections((Arc) seg1, (LineSegment) seg2); + } else if (seg1 instanceof LineSegment && seg2 instanceof Arc) { + return arcLineSegmentIntersections((Arc) seg2, (LineSegment) seg1); + } else { + throw new IllegalArgumentException( + "Arguments must be Arc or LineSegment, got: " + + seg1.getClass().getName() + " and " + seg2.getClass().getName()); + } + } + + /** + * Returns true if line1 intersects line2. + */ + public static boolean intersects(Object seg1, Object seg2) { + return !intersections(seg1, seg2).isEmpty(); + } + + // ==================== Utility Methods ==================== + + /** + * Returns true if point is at or very near the north pole. + */ + private static boolean isNorthPole(Point p) { + return MathUtils.doubleApprox(p.getLat(), 90.0, APPROXIMATION_DELTA); + } + + /** + * Returns true if point is at or very near the south pole. + */ + private static boolean isSouthPole(Point p) { + return MathUtils.doubleApprox(p.getLat(), -90.0, APPROXIMATION_DELTA); + } + + /** + * Returns the point antipodal (opposite side of sphere) to the given point. + */ + private static Point antipodal(Point p) { + double lon = p.getLon(); + double lat = p.getLat(); + + // Antipodal longitude: add 180 degrees (wrapping to [-180, 180]) + double antipodalLon = lon + 180.0; + if (antipodalLon > 180.0) { + antipodalLon -= 360.0; + } + + // Antipodal latitude: negate + double antipodalLat = -lat; + + return new Point(antipodalLon, antipodalLat); + } + + /** + * Returns true if two points are approximately equal within the given delta. + */ + private static boolean pointsApproxEqual(Point p1, Point p2, double delta) { + if (!MathUtils.doubleApprox(p1.getLat(), p2.getLat(), delta)) { + return false; + } + + if (MathUtils.doubleApprox(p1.getLon(), p2.getLon(), delta)) { + return true; + } + + // Check if both on antimeridian + if (Math.abs(p1.getLon()) == 180.0 && Math.abs(p2.getLon()) == 180.0) { + return true; + } + + // Check if both at poles + if (isNorthPole(p1) && isNorthPole(p2)) { + return true; + } + if (isSouthPole(p1) && isSouthPole(p2)) { + return true; + } + + return false; + } + + /** + * Returns true if MBR covers the given point (geodetic coordinates). + */ + private static boolean geodeticCoversPoint(Mbr mbr, Point p) { + double tolerance = MathUtils.COVERS_TOLERANCE; + double lat = p.getLat(); + + // Handle poles + if (isNorthPole(p)) { + return coversLat(mbr, 90.0, tolerance); + } + if (isSouthPole(p)) { + return coversLat(mbr, -90.0, tolerance); + } + + return coversLat(mbr, lat, tolerance) && coversLonGeodetic(mbr, p.getLon(), tolerance); + } + + /** + * Returns true if MBR covers the given latitude. + */ + private static boolean coversLat(Mbr mbr, double lat, double tolerance) { + double north = mbr.getNorth() + tolerance; + double south = mbr.getSouth() - tolerance; + return lat >= south && lat <= north; + } + + /** + * Returns true if MBR covers the given longitude (geodetic). + * Handles antimeridian crossing. + */ + private static boolean coversLonGeodetic(Mbr mbr, double lon, double tolerance) { + double west = mbr.getWest(); + double east = mbr.getEast(); + + // Check for antimeridian crossing + if (west > east) { + // MBR crosses antimeridian + return lon >= (west - tolerance) || lon <= (east + tolerance); + } else if (Math.abs(lon) == 180.0) { + double within180 = 180.0 - tolerance; + return Math.abs(west) >= within180 || Math.abs(east) >= within180; + } else { + return lon >= (west - tolerance) && lon <= (east + tolerance); + } + } + + /** + * Returns true if MBR covers point in cartesian coordinates. + */ + private static boolean cartesianCoversPoint(Mbr mbr, Point point) { + double lon = point.getLon(); + double lat = point.getLat(); + + return lon >= mbr.getWest() && lon <= mbr.getEast() && + lat >= mbr.getSouth() && lat <= mbr.getNorth(); + } + + /** + * Removes duplicate points from a list based on approximate equality. + */ + private static List removeDuplicatePoints(List points) { + List result = new ArrayList<>(); + + for (Point p : points) { + boolean isDuplicate = false; + for (Point existing : result) { + if (pointsApproxEqual(p, existing, APPROXIMATION_DELTA)) { + isDuplicate = true; + break; + } + } + if (!isDuplicate) { + result.add(p); + } + } + + return result; + } +} diff --git a/spatial-lib/src/java/cmr/spatial/internal/arc/GreatCircle.java b/spatial-lib/src/java/cmr/spatial/internal/arc/GreatCircle.java new file mode 100644 index 0000000000..85043e0b43 --- /dev/null +++ b/spatial-lib/src/java/cmr/spatial/internal/arc/GreatCircle.java @@ -0,0 +1,47 @@ +package cmr.spatial.internal.arc; + +import cmr.spatial.shape.Point; +import cmr.spatial.math.Vector; + +/** + * Represents a great circle on a sphere. + * A great circle is the intersection of the sphere with a plane passing through the sphere's center. + * Immutable value class. + */ +public final class GreatCircle { + private final Vector planeVector; + private final Point northernmostPoint; + private final Point southernmostPoint; + + public GreatCircle(Vector planeVector, Point northernmostPoint, Point southernmostPoint) { + this.planeVector = planeVector; + this.northernmostPoint = northernmostPoint; + this.southernmostPoint = southernmostPoint; + } + + public Vector getPlaneVector() { + return planeVector; + } + + public Point getNorthernmostPoint() { + return northernmostPoint; + } + + public Point getSouthernmostPoint() { + return southernmostPoint; + } + + /** + * Returns true if two great circles are equivalent. + * Two great circles are equivalent if their plane vectors are parallel. + */ + public boolean isEquivalent(GreatCircle other) { + return this.planeVector.isParallel(other.planeVector); + } + + @Override + public String toString() { + return String.format("GreatCircle{north=%s, south=%s}", + northernmostPoint, southernmostPoint); + } +} diff --git a/spatial-lib/src/java/cmr/spatial/internal/ring/CartesianRing.java b/spatial-lib/src/java/cmr/spatial/internal/ring/CartesianRing.java new file mode 100644 index 0000000000..0965f90a63 --- /dev/null +++ b/spatial-lib/src/java/cmr/spatial/internal/ring/CartesianRing.java @@ -0,0 +1,301 @@ +package cmr.spatial.internal.ring; + +import cmr.spatial.shape.Point; +import cmr.spatial.shape.Mbr; +import cmr.spatial.internal.segment.LineSegment; +import cmr.spatial.math.MathUtils; + +import java.util.ArrayList; +import java.util.HashSet; +import java.util.List; +import java.util.Set; + +/** + * Represents a cartesian (flat plane) ring - a closed sequence of points + * connected by straight line segments. + * + * Used for cartesian polygons. Points must be in counter-clockwise order. + * The last point must match the first point to close the ring. + * + * This class handles: + * - Ray-casting algorithm for point-in-ring tests + * - Fixed external point (outside Earth's coordinate bounds) + * - Simple MBR calculation (no antimeridian crossing) + * + * Based on cmr.spatial.cartesian-ring Clojure namespace. + */ +public class CartesianRing { + + private static final double APPROXIMATION_DELTA = LineSegment.COVERS_TOLERANCE; + private static final double COVERS_TOLERANCE = 0.0000000001; // Matches cmr.spatial.mbr/COVERS_TOLERANCE + + /** + * Fixed external point that's outside all valid cartesian rings. + * Uses coordinates outside the Earth's bounds (lon=181, lat=91). + */ + private static final Point EXTERNAL_POINT = new Point(181.0, 91.0); + + private final List points; + private final List lineSegments; + private final Mbr mbr; + + /** + * Private constructor. Use createRing() factory method instead. + */ + private CartesianRing(List points, List lineSegments, Mbr mbr) { + this.points = points; + this.lineSegments = lineSegments; + this.mbr = mbr; + } + + /** + * Creates a cartesian ring from a list of points. + * Calculates derived fields: line segments and MBR. + * + * @param points List of points forming the ring (last must equal first) + * Can be either cmr.spatial.shape.Point or cmr.spatial.point.Point (Clojure records) + * @return A new CartesianRing + */ + public static CartesianRing createRing(List points) { + // Convert Clojure Points to Java Points if necessary + List javaPoints = convertToJavaPoints(points); + List lineSegments = pointsToLineSegments(javaPoints); + Mbr mbr = calculateMbr(lineSegments); + + return new CartesianRing(javaPoints, lineSegments, mbr); + } + + /** + * Converts a list of objects (potentially Clojure Point records) to Java Points. + */ + private static List convertToJavaPoints(List points) { + List result = new ArrayList<>(); + for (Object p : points) { + if (p instanceof Point) { + result.add((Point) p); + } else { + // Assume it's a Clojure Point record with lon/lat fields + try { + java.lang.reflect.Field lonField = p.getClass().getField("lon"); + java.lang.reflect.Field latField = p.getClass().getField("lat"); + double lon = ((Number) lonField.get(p)).doubleValue(); + double lat = ((Number) latField.get(p)).doubleValue(); + result.add(new Point(lon, lat)); + } catch (Exception e) { + throw new IllegalArgumentException( + "Could not convert point to Java Point: " + p.getClass().getName(), e); + } + } + } + return result; + } + + /** + * Converts a list of points to a list of line segments. + */ + private static List pointsToLineSegments(List points) { + List segments = new ArrayList<>(); + for (int i = 0; i < points.size() - 1; i++) { + Point p1 = points.get(i); + Point p2 = points.get(i + 1); + + // Skip if points are equal + if (!pointsEqual(p1, p2)) { + segments.add(LineSegment.createLineSegment(p1, p2)); + } + } + return segments; + } + + /** + * Calculates the minimum bounding rectangle for the ring. + * For cartesian rings, the MBR does not cross the antimeridian. + */ + private static Mbr calculateMbr(List segments) { + if (segments.isEmpty()) { + return new Mbr(0.0, 0.0, 0.0, 0.0); + } + + Mbr result = segments.get(0).getMbr(); + + for (int i = 1; i < segments.size(); i++) { + Mbr segmentMbr = segments.get(i).getMbr(); + result = unionMbrsCartesian(result, segmentMbr); + } + + return result; + } + + /** + * Unions two MBRs (cartesian - no antimeridian crossing). + */ + private static Mbr unionMbrsCartesian(Mbr mbr1, Mbr mbr2) { + double west = Math.min(mbr1.getWest(), mbr2.getWest()); + double north = Math.max(mbr1.getNorth(), mbr2.getNorth()); + double east = Math.max(mbr1.getEast(), mbr2.getEast()); + double south = Math.min(mbr1.getSouth(), mbr2.getSouth()); + + return new Mbr(west, north, east, south); + } + + /** + * Tests if a point is covered by (inside) this ring using ray-casting algorithm. + * + * Algorithm (matching Clojure implementation): + * 1. First check if MBR covers the point (optimization) + * 2. Check if point is in the point set + * 3. Cast a ray from the test point to the external point + * 4. Count unique intersections (rounded to 5 decimal places) + * 5. Odd count = inside, even count = outside + * + * @param point The point to test + * @return true if the ring covers the point + */ + public boolean coversPoint(Point point) { + // Only do real intersection if the MBR covers the point (matches Clojure) + if (!cartesianMbrCoversPoint(mbr, point)) { + return false; + } + + // Check if point is one of the ring's points (matches Clojure point-set check) + for (Point p : points) { + if (pointsEqual(p, point)) { + return true; + } + } + + // Create the test ray from point to external point (matches Clojure direction) + LineSegment crossingLine = LineSegment.createLineSegment(point, EXTERNAL_POINT); + + // Find all intersections with ring segments and round them + Set intersections = new HashSet<>(); + for (LineSegment segment : lineSegments) { + Point intersection = crossingLine.intersection(segment); + if (intersection != null) { + // Round to 5 decimal places to match Clojure INTERSECTION_POINT_PRECISION + Point rounded = roundPoint(intersection); + intersections.add(rounded); + } + } + + // Check if point itself is one of the intersections (matches Clojure) + if (intersections.contains(roundPoint(point))) { + return true; + } + + // Odd number of intersections = inside, even = outside + return intersections.size() % 2 == 1; + } + + /** + * Checks if an MBR covers a point in cartesian space. + * Matches cmr.spatial.mbr/cartesian-covers-point? with default tolerance. + */ + private static boolean cartesianMbrCoversPoint(Mbr mbr, Point point) { + double lon = point.getLon(); + double lat = point.getLat(); + double delta = COVERS_TOLERANCE; + + // Check latitude with tolerance (covers-lat?) + double north = mbr.getNorth() + delta; + double south = mbr.getSouth() - delta; + if (lat < south || lat > north) { + return false; + } + + // Check longitude with tolerance (cartesian-lon-range-covers-lon?) + double west = mbr.getWest() - delta; + double east = mbr.getEast() + delta; + boolean crossesAntimeridian = west > east; + + if (crossesAntimeridian) { + return lon >= west || lon <= east; + } else { + return lon >= west && lon <= east; + } + } + + /** + * Rounds a point to a fixed precision to handle near-duplicates. + * Uses precision = 5 to match Clojure INTERSECTION_POINT_PRECISION. + */ + private static Point roundPoint(Point p) { + double scale = 100000.0; // 5 decimal places (10^5) + double lon = Math.round(p.getLon() * scale) / scale; + double lat = Math.round(p.getLat() * scale) / scale; + return new Point(lon, lat); + } + + /** + * Returns true if two points are equal (with tolerance). + */ + private static boolean pointsEqual(Point p1, Point p2) { + return MathUtils.doubleApprox(p1.getLon(), p2.getLon(), APPROXIMATION_DELTA) && + MathUtils.doubleApprox(p1.getLat(), p2.getLat(), APPROXIMATION_DELTA); + } + + /** + * Determines the winding order of the ring. + * Uses the sum-over-area-under-edges algorithm. + * + * @return CLOCKWISE or COUNTER_CLOCKWISE + */ + public WindingOrder calculateWinding() { + double sum = 0.0; + + for (LineSegment segment : lineSegments) { + Point p1 = segment.getPoint1(); + Point p2 = segment.getPoint2(); + + double x1 = p1.getLon(); + double y1 = p1.getLat(); + double x2 = p2.getLon(); + double y2 = p2.getLat(); + + // Sum (x2 - x1) * (y2 + y1) + sum += (x2 - x1) * (y2 + y1); + } + + return (sum >= 0.0) ? WindingOrder.CLOCKWISE : WindingOrder.COUNTER_CLOCKWISE; + } + + // Getters + + public List getPointsList() { + return points; + } + + public Point[] getPoints() { + return points.toArray(new Point[0]); + } + + public List getLineSegmentsList() { + return lineSegments; + } + + public LineSegment[] getLineSegments() { + return lineSegments.toArray(new LineSegment[0]); + } + + public Mbr getMbr() { + return mbr; + } + + public static Point getExternalPoint() { + return EXTERNAL_POINT; + } + + /** + * Enum for winding order. + */ + public enum WindingOrder { + CLOCKWISE, + COUNTER_CLOCKWISE + } + + @Override + public String toString() { + return String.format("CartesianRing{points=%d, segments=%d}", + points.size(), lineSegments.size()); + } +} diff --git a/spatial-lib/src/java/cmr/spatial/internal/ring/GeodeticRing.java b/spatial-lib/src/java/cmr/spatial/internal/ring/GeodeticRing.java new file mode 100644 index 0000000000..8111e207e0 --- /dev/null +++ b/spatial-lib/src/java/cmr/spatial/internal/ring/GeodeticRing.java @@ -0,0 +1,750 @@ +package cmr.spatial.internal.ring; + +import cmr.spatial.shape.Point; +import cmr.spatial.shape.Mbr; +import cmr.spatial.internal.arc.Arc; +import cmr.spatial.internal.arc.ArcLineSegmentIntersections; +import cmr.spatial.math.MathUtils; + +import java.util.ArrayList; +import java.util.HashSet; +import java.util.List; +import java.util.Set; + +/** + * Represents a geodetic (spherical) ring - a closed sequence of points on a sphere + * connected by great circle arcs. + * + * Used for geodetic polygons. Points must be in counter-clockwise order when viewed + * from above the surface. The last point must match the first point to close the ring. + * + * This class handles: + * - Ray-casting algorithm for point-in-ring tests + * - Pole containment (rings can contain north/south poles) + * - Multiple external points for antipodal safety + * - Antimeridian crossing + * + * Based on cmr.spatial.geodetic-ring Clojure namespace. + */ +public class GeodeticRing { + + private static final double EXTERNAL_POINT_PRECISION = 4.0; + private static final double INTERSECTION_POINT_PRECISION = 5.0; // Matches Clojure + private static final double APPROXIMATION_DELTA = Arc.APPROXIMATION_DELTA; + private static final double MBR_COVERS_TOLERANCE = 0.0000000001; // Matches Clojure mbr/COVERS_TOLERANCE + + private final List points; + private final List arcs; + private final Mbr mbr; + private final List externalPoints; + private final boolean containsNorthPole; + private final boolean containsSouthPole; + + /** + * Private constructor. Use createRing() factory method instead. + */ + private GeodeticRing(List points, List arcs, Mbr mbr, + List externalPoints, boolean containsNorthPole, + boolean containsSouthPole) { + this.points = points; + this.arcs = arcs; + this.mbr = mbr; + this.externalPoints = externalPoints; + this.containsNorthPole = containsNorthPole; + this.containsSouthPole = containsSouthPole; + } + + /** + * Creates a geodetic ring from a list of points. + * Calculates derived fields: arcs, MBR, external points, pole containment. + * + * @param points List of points forming the ring (last must equal first) + * Can be either cmr.spatial.shape.Point or cmr.spatial.point.Point (Clojure records) + * @return A new GeodeticRing + */ + public static GeodeticRing createRing(List points) { + // Convert Clojure Points to Java Points if necessary + List javaPoints = convertToJavaPoints(points); + List arcs = pointsToArcs(javaPoints); + + // Calculate course rotation direction and pole containment + RotationDirection courseRotation = calculateRotationDirection(arcs); + boolean containsNorthPole = determineNorthPoleContainment(javaPoints, arcs, courseRotation); + boolean containsSouthPole = determineSouthPoleContainment(javaPoints, arcs, courseRotation); + + // Calculate MBR + Mbr mbr = calculateMbr(arcs, containsNorthPole, containsSouthPole, javaPoints); + + // Calculate external points + List externalPoints = calculateExternalPoints(mbr, containsNorthPole, containsSouthPole); + + return new GeodeticRing(javaPoints, arcs, mbr, externalPoints, + containsNorthPole, containsSouthPole); + } + + /** + * Converts a list of objects (potentially Clojure Point records) to Java Points. + */ + private static List convertToJavaPoints(List points) { + List result = new ArrayList<>(); + for (Object p : points) { + if (p instanceof Point) { + result.add((Point) p); + } else { + // Assume it's a Clojure Point record with lon/lat fields + try { + java.lang.reflect.Field lonField = p.getClass().getField("lon"); + java.lang.reflect.Field latField = p.getClass().getField("lat"); + double lon = ((Number) lonField.get(p)).doubleValue(); + double lat = ((Number) latField.get(p)).doubleValue(); + result.add(new Point(lon, lat)); + } catch (Exception e) { + throw new IllegalArgumentException( + "Could not convert point to Java Point: " + p.getClass().getName(), e); + } + } + } + return result; + } + + /** + * Converts a list of points to a list of arcs. + */ + private static List pointsToArcs(List points) { + List arcs = new ArrayList<>(); + for (int i = 0; i < points.size() - 1; i++) { + Point p1 = points.get(i); + Point p2 = points.get(i + 1); + + // Skip creating arc if points are equal (would throw exception) + if (!pointsEqual(p1, p2)) { + try { + arcs.add(Arc.createArc(p1, p2)); + } catch (IllegalArgumentException e) { + // Skip antipodal points or other invalid arcs + } + } + } + return arcs; + } + + /** + * Calculates the rotation direction of the ring. + */ + private static RotationDirection calculateRotationDirection(List arcs) { + if (arcs.isEmpty()) { + return RotationDirection.NONE; + } + + // Collect all initial and ending courses from each arc + List courses = new ArrayList<>(); + for (Arc arc : arcs) { + courses.add(arc.getInitialCourse()); + courses.add(arc.getEndingCourse()); + } + + // Add the first course again to complete the turn + courses.add(courses.get(0)); + + // Calculate net bearing change by summing angle deltas + double netTurn = 0.0; + for (int i = 0; i < courses.size() - 1; i++) { + netTurn += angleDelta(courses.get(i), courses.get(i + 1)); + } + + // Counter-clockwise: ~360 degrees, Clockwise: ~-360 degrees, None: ~0 degrees + if (Math.abs(netTurn) < 0.01) { + return RotationDirection.NONE; + } else if (netTurn > 0.0) { + return RotationDirection.COUNTER_CLOCKWISE; + } else { + return RotationDirection.CLOCKWISE; + } + } + + /** + * Find the difference between a pair of angles (in degrees). + * Matches the Clojure angle-delta function. + */ + private static double angleDelta(double a1, double a2) { + double a2Shifted = a2; + + // Shift angle 2 so it is always greater than angle 1 + if (a2 < a1) { + a2Shifted = a2 + 360.0; + } + + // Calculate left turn amount + double leftTurn = a2Shifted - a1; + + // Determine which is smaller: turning left or turning right + if (leftTurn == 180.0) { + // Can't determine which is smaller, return 0 + return 0.0; + } else if (leftTurn > 180.0) { + // Turning right is smaller (returns negative value) + return leftTurn - 360.0; + } else { + return leftTurn; + } + } + + /** + * Determines if the ring contains the north pole. + */ + private static boolean determineNorthPoleContainment(List points, List arcs, + RotationDirection rotation) { + // Check if any point is the north pole + for (Point p : points) { + if (isNorthPole(p)) { + return true; + } + } + + // Check if any arc crosses the north pole + for (Arc arc : arcs) { + if (arc.crossesNorthPole()) { + return true; + } + } + + // If rotation is NONE, check longitude rotation direction + if (rotation == RotationDirection.NONE) { + RotationDirection lonRotation = calculateLongitudeRotation(points); + return lonRotation == RotationDirection.COUNTER_CLOCKWISE; + } + + return false; + } + + /** + * Determines if the ring contains the south pole. + */ + private static boolean determineSouthPoleContainment(List points, List arcs, + RotationDirection rotation) { + // Check if any point is the south pole + for (Point p : points) { + if (isSouthPole(p)) { + return true; + } + } + + // Check if any arc crosses the south pole + for (Arc arc : arcs) { + if (arc.crossesSouthPole()) { + return true; + } + } + + // If rotation is NONE, check longitude rotation direction + if (rotation == RotationDirection.NONE) { + RotationDirection lonRotation = calculateLongitudeRotation(points); + return lonRotation == RotationDirection.CLOCKWISE; + } + + return false; + } + + /** + * Calculates the rotation direction based on longitudes only. + */ + private static RotationDirection calculateLongitudeRotation(List points) { + if (points.size() < 2) { + return RotationDirection.NONE; + } + + // Extract longitudes + List longitudes = new ArrayList<>(); + for (Point p : points) { + longitudes.add(p.getLon()); + } + + // Calculate net rotation using angle deltas + double netTurn = 0.0; + for (int i = 0; i < longitudes.size() - 1; i++) { + netTurn += angleDelta(longitudes.get(i), longitudes.get(i + 1)); + } + + // Determine direction + if (Math.abs(netTurn) < 0.01) { + return RotationDirection.NONE; + } else if (netTurn > 0.0) { + return RotationDirection.COUNTER_CLOCKWISE; + } else { + return RotationDirection.CLOCKWISE; + } + } + + /** + * Calculates the minimum bounding rectangle for the ring. + */ + private static Mbr calculateMbr(List arcs, boolean containsNorthPole, + boolean containsSouthPole, List points) { + if (arcs.isEmpty()) { + return new Mbr(-180.0, 90.0, 180.0, -90.0); + } + + Mbr result = null; + + // Union all arc MBRs + for (Arc arc : arcs) { + Mbr[] arcMbrs = arc.getMbrs(); + for (Mbr arcMbr : arcMbrs) { + if (result == null) { + result = arcMbr; + } else { + result = unionMbrs(result, arcMbr); + } + } + } + + // Expand to include north pole if contained but not explicitly crossed + if (containsNorthPole && !anyPointIsNorthPole(points) && !anyArcCrossesNorthPole(arcs)) { + result = new Mbr(-180.0, 90.0, 180.0, result.getSouth()); + } + + // Expand to include south pole if contained but not explicitly crossed + if (containsSouthPole && !anyPointIsSouthPole(points) && !anyArcCrossesSouthPole(arcs)) { + result = new Mbr(-180.0, result.getNorth(), 180.0, -90.0); + } + + return result; + } + + /** + * Unions two MBRs handling geodetic wrapping. + * Faithful translation of cmr.spatial.mbr/union. + */ + private static Mbr unionMbrs(Mbr m1, Mbr m2) { + double north = Math.max(m1.getNorth(), m2.getNorth()); + double south = Math.min(m1.getSouth(), m2.getSouth()); + + boolean m1Crosses = m1.getWest() > m1.getEast(); + boolean m2Crosses = m2.getWest() > m2.getEast(); + + double west, east; + + if (m1Crosses && m2Crosses) { + // Both cross antimeridian + west = Math.min(m1.getWest(), m2.getWest()); + east = Math.max(m1.getEast(), m2.getEast()); + + // If result would cover whole world, set it to that + if (west <= east) { + west = -180.0; + east = 180.0; + } + } else if (m1Crosses || m2Crosses) { + // One crosses the antimeridian + // Make m1 be the one that crosses + Mbr crossing = m1Crosses ? m1 : m2; + Mbr notCrossing = m1Crosses ? m2 : m1; + + double w1 = crossing.getWest(); + double e1 = crossing.getEast(); + double w2 = notCrossing.getWest(); + double e2 = notCrossing.getEast(); + + // We could expand m1 to the east or to the west. Pick the shorter of the two. + double westDist = w1 - w2; + double eastDist = e2 - e1; + + if (westDist <= 0.0 || eastDist <= 0.0) { + // Non-crossing is already contained + west = w1; + east = e1; + } else if (eastDist < westDist) { + // Expand east + west = w1; + east = e2; + } else { + // Expand west + west = w2; + east = e1; + } + + // If result would cover whole world, set it to that + if (west <= east) { + west = -180.0; + east = 180.0; + } + } else { + // Neither crosses the antimeridian + // Ensure m1.west <= m2.west + if (m1.getWest() > m2.getWest()) { + Mbr temp = m1; + m1 = m2; + m2 = temp; + } + + double w1 = m1.getWest(); + double e1 = m1.getEast(); + double w2 = m2.getWest(); + double e2 = m2.getEast(); + + west = Math.min(w1, w2); + east = Math.max(e1, e2); + + // Check if it's shorter to cross the antimeridian + double dist = east - west; + double altWest = w2; + double altEast = e1; + double altDist = (180.0 - altWest) + (altEast - (-180.0)); + + if (altDist < dist) { + west = altWest; + east = altEast; + } + } + + return new Mbr(west, north, east, south); + } + + /** + * Calculates external points (points guaranteed to be outside the ring). + * Used for ray-casting algorithm. + * Matches the logic in cmr.spatial.mbr/external-points. + */ + private static List calculateExternalPoints(Mbr mbr, boolean containsNorthPole, + boolean containsSouthPole) { + List result = new ArrayList<>(); + + // Cannot determine external points if ring contains both poles + if (containsNorthPole && containsSouthPole) { + return result; + } + + // Find the biggest area around the MBR to place external points + boolean crossesAntimeridian = mbr.getWest() > mbr.getEast(); + + double northDist = 90.0 - mbr.getNorth(); + double southDist = mbr.getSouth() - (-90.0); + double westDist = crossesAntimeridian ? 0.0 : mbr.getWest() - (-180.0); + double eastDist = crossesAntimeridian ? 0.0 : 180.0 - mbr.getEast(); + + double biggestDist = Math.max(Math.max(northDist, southDist), + Math.max(westDist, eastDist)); + + double w, n, e, s; + + if (biggestDist == northDist) { + // Place points in area north of MBR + w = -180.0; + n = 90.0; + e = 180.0; + s = mbr.getNorth(); + } else if (biggestDist == southDist) { + // Place points in area south of MBR + w = -180.0; + n = mbr.getSouth(); + e = 180.0; + s = -90.0; + } else if (!crossesAntimeridian && biggestDist == westDist) { + // Place points in area west of MBR + w = -180.0; + n = 90.0; + e = mbr.getWest(); + s = -90.0; + } else if (!crossesAntimeridian && biggestDist == eastDist) { + // Place points in area east of MBR + w = mbr.getEast(); + n = 90.0; + e = 180.0; + s = -90.0; + } else if (crossesAntimeridian) { + // Place points in area between east and west (the gap when crossing antimeridian) + w = mbr.getEast(); + n = 90.0; + e = mbr.getWest(); + s = -90.0; + } else { + throw new RuntimeException("Logic error in calculateExternalPoints: " + + "Could not determine biggest distance area"); + } + + // Find 3 points within the area: left, middle, right + // These are distributed along the longitude range at the mid-latitude + double midLon = mid(w, e); + double rightLon = mid(w, midLon); + double leftLon = mid(midLon, e); + double midLat = mid(n, s); + + result.add(roundPointForExternalPoints(new Point(leftLon, midLat))); + result.add(roundPointForExternalPoints(new Point(midLon, midLat))); + result.add(roundPointForExternalPoints(new Point(rightLon, midLat))); + + return result; + } + + /** + * Calculates the midpoint between two values. + */ + private static double mid(double a, double b) { + return (a + b) / 2.0; + } + + /** + * Tests if a point is covered by (inside) this ring using ray-casting algorithm. + * Matches the Clojure implementation in geodetic_ring.clj exactly. + * + * @param point The point to test + * @return true if the ring covers the point + */ + public boolean coversPoint(Point point) { + // Check if ring contains north or south pole (matches Clojure) + if (containsNorthPole && isNorthPole(point)) { + return true; + } + if (containsSouthPole && isSouthPole(point)) { + return true; + } + + // Only do real intersection if the MBR covers the point (matches Clojure) + if (!geodeticMbrCoversPoint(mbr, point)) { + return false; + } + + // Check if point is in the point set (matches Clojure) + for (Point p : points) { + if (pointsEqual(p, point)) { + return true; + } + } + + // Choose an external point for ray-casting + Point externalPoint = chooseExternalPoint(point); + if (externalPoint == null) { + throw new RuntimeException("Could not find suitable external point for ring coverage test. " + + "This can happen when the ring contains both poles."); + } + + // Create the test arc from point to external point (matches Clojure) + Arc crossingArc; + try { + crossingArc = Arc.createArc(point, externalPoint); + } catch (IllegalArgumentException e) { + // Points are equal or antipodal - shouldn't happen with proper external point + return false; + } + + // Find all intersections with ring arcs and round them + Set intersections = new HashSet<>(); + for (Arc arc : arcs) { + List arcIntersections = ArcLineSegmentIntersections.arcArcIntersections(crossingArc, arc); + for (Point p : arcIntersections) { + // Round to 5 decimal places to match Clojure INTERSECTION_POINT_PRECISION + Point rounded = roundPoint(p); + intersections.add(rounded); + } + } + + // Check if point itself is one of the intersections (matches Clojure) + if (intersections.contains(roundPoint(point))) { + return true; + } + + // Odd number of intersections = inside, even = outside + return intersections.size() % 2 == 1; + } + + /** + * Checks if an MBR covers a point in geodetic space. + * Matches cmr.spatial.mbr/geodetic-covers-point? and geodetic-lon-range-covers-lon? + */ + private static boolean geodeticMbrCoversPoint(Mbr mbr, Point point) { + double lon = point.getLon(); + double lat = point.getLat(); + + // Check latitude bounds + if (lat < mbr.getSouth() || lat > mbr.getNorth()) { + return false; + } + + // Check longitude - matches geodetic-lon-range-covers-lon? exactly + double west = mbr.getWest() - MBR_COVERS_TOLERANCE; + double east = mbr.getEast() + MBR_COVERS_TOLERANCE; + boolean crossesAntimeridian = west > east; + + if (crossesAntimeridian) { + // Crosses antimeridian: longitude is in [west, 180] or [-180, east] + return lon >= west || lon <= east; + } else if (Math.abs(lon) == 180.0) { + // Special case: point is on antimeridian (±180) + // Check if west or east is within tolerance of ±180 + double within180 = 180.0 - MBR_COVERS_TOLERANCE; + return Math.abs(west) >= within180 || Math.abs(east) >= within180; + } else { + // Normal case: doesn't cross antimeridian + return lon >= west && lon <= east; + } + } + + /** + * Checks if a point is at the North Pole. + */ + private static boolean isNorthPole(Point p) { + return Math.abs(p.getLat() - 90.0) < APPROXIMATION_DELTA; + } + + /** + * Checks if a point is at the South Pole. + */ + private static boolean isSouthPole(Point p) { + return Math.abs(p.getLat() + 90.0) < APPROXIMATION_DELTA; + } + + /** + * Chooses an external point that is not equal or antipodal to the test point. + * Matches choose-external-point in Clojure. + */ + private Point chooseExternalPoint(Point point) { + // Round to EXTERNAL_POINT_PRECISION (4) to match Clojure + Point rounded = roundPointForExternalPoints(point); + Point antipodal = antipodalPoint(rounded); + + for (Point ep : externalPoints) { + // Use exact equality since both are rounded to same precision + if (!pointsExactlyEqual(ep, rounded) && !pointsExactlyEqual(ep, antipodal)) { + return ep; + } + } + + return null; + } + + /** + * Checks if two points are exactly equal (no approximation). + * Used when comparing rounded points. + */ + private static boolean pointsExactlyEqual(Point p1, Point p2) { + return p1.getLon() == p2.getLon() && p1.getLat() == p2.getLat(); + } + + /** + * Returns the antipodal point (opposite side of the sphere). + */ + private static Point antipodalPoint(Point p) { + double lon = p.getLon(); + double lat = p.getLat(); + + double antiLon = lon + 180.0; + if (antiLon > 180.0) antiLon -= 360.0; + + return new Point(antiLon, -lat); + } + + /** + * Rounds a point for external point generation (precision 4). + * Matches EXTERNAL_POINT_PRECISION in Clojure. + */ + private static Point roundPointForExternalPoints(Point p) { + double scale = Math.pow(10, EXTERNAL_POINT_PRECISION); + double lon = Math.round(p.getLon() * scale) / scale; + double lat = Math.round(p.getLat() * scale) / scale; + return new Point(lon, lat); + } + + /** + * Rounds a point for intersection calculations (precision 5). + * Matches INTERSECTION_POINT_PRECISION in Clojure. + */ + private static Point roundPoint(Point p) { + double scale = Math.pow(10, INTERSECTION_POINT_PRECISION); // 100000.0 + double lon = Math.round(p.getLon() * scale) / scale; + double lat = Math.round(p.getLat() * scale) / scale; + return new Point(lon, lat); + } + + /** + * Returns true if any point in the list is at the north pole. + */ + private static boolean anyPointIsNorthPole(List points) { + for (Point p : points) { + if (isNorthPole(p)) return true; + } + return false; + } + + /** + * Returns true if any point in the list is at the south pole. + */ + private static boolean anyPointIsSouthPole(List points) { + for (Point p : points) { + if (isSouthPole(p)) return true; + } + return false; + } + + /** + * Returns true if any arc crosses the north pole. + */ + private static boolean anyArcCrossesNorthPole(List arcs) { + for (Arc arc : arcs) { + if (arc.crossesNorthPole()) return true; + } + return false; + } + + /** + * Returns true if any arc crosses the south pole. + */ + private static boolean anyArcCrossesSouthPole(List arcs) { + for (Arc arc : arcs) { + if (arc.crossesSouthPole()) return true; + } + return false; + } + + /** + * Returns true if two points are equal (with tolerance). + */ + private static boolean pointsEqual(Point p1, Point p2) { + return MathUtils.doubleApprox(p1.getLon(), p2.getLon(), APPROXIMATION_DELTA) && + MathUtils.doubleApprox(p1.getLat(), p2.getLat(), APPROXIMATION_DELTA); + } + + // Getters + + public List getPointsList() { + return points; + } + + public Point[] getPoints() { + return points.toArray(new Point[0]); + } + + public List getArcsList() { + return arcs; + } + + public Arc[] getArcs() { + return arcs.toArray(new Arc[0]); + } + + public Mbr getMbr() { + return mbr; + } + + public boolean containsNorthPole() { + return containsNorthPole; + } + + public boolean containsSouthPole() { + return containsSouthPole; + } + + /** + * Enum for rotation direction. + */ + public enum RotationDirection { + CLOCKWISE, + COUNTER_CLOCKWISE, + NONE + } + + @Override + public String toString() { + return String.format("GeodeticRing{points=%d, arcs=%d, northPole=%b, southPole=%b}", + points.size(), arcs.size(), containsNorthPole, containsSouthPole); + } +} diff --git a/spatial-lib/src/java/cmr/spatial/internal/segment/LineSegment.java b/spatial-lib/src/java/cmr/spatial/internal/segment/LineSegment.java new file mode 100644 index 0000000000..1ff8cce425 --- /dev/null +++ b/spatial-lib/src/java/cmr/spatial/internal/segment/LineSegment.java @@ -0,0 +1,641 @@ +package cmr.spatial.internal.segment; + +import cmr.spatial.shape.Mbr; +import cmr.spatial.shape.Point; +import cmr.spatial.math.MathUtils; +import cmr.spatial.geometry.MbrIntersections; + +import java.util.ArrayList; +import java.util.List; + +/** + * Represents a cartesian line segment (flat plane geometry). + * A line segment is defined by two endpoints and has properties like slope and y-intercept. + * + * This class handles special cases: + * - Vertical lines (undefined slope) + * - Horizontal lines (zero slope) + * - Parallel lines (same slope) + * + * For intersections, the intersection point must be within both line segments' MBRs. + * + * Based on cmr.spatial.line-segment Clojure namespace. + */ +public class LineSegment { + + /** + * Tolerance used for determining if points are on the line segment. + */ + public static final double COVERS_TOLERANCE = 0.00001; + + /** + * Tolerance used in MBR coverage checks (matches cmr.spatial.mbr/COVERS_TOLERANCE). + */ + public static final double MBR_COVERS_TOLERANCE = 0.0000000001; + + /** + * Tolerance used for the covers method during point intersections. + * Longitudes and latitudes technically outside the bounding rectangle but within + * this tolerance will be considered covered by the bounding rectangle. + */ + public static final double INTERSECTION_COVERS_TOLERANCE = 0.0000001; + + private final Point point1; + private final Point point2; + private final boolean vertical; + private final boolean horizontal; + private final Double m; // Slope (can be null for vertical lines) + private final Double b; // Y-intercept (can be null for vertical lines) + private final Mbr mbr; // Bounding rectangle + + /** + * Private constructor. Use createLineSegment factory method instead. + */ + private LineSegment(Point point1, Point point2, boolean vertical, boolean horizontal, + Double m, Double b, Mbr mbr) { + this.point1 = point1; + this.point2 = point2; + this.vertical = vertical; + this.horizontal = horizontal; + this.m = m; + this.b = b; + this.mbr = mbr; + } + + /** + * Creates a new line segment from two points. + * Calculates slope, y-intercept, and bounding rectangle. + * + * @param p1 First endpoint + * @param p2 Second endpoint + * @return A new LineSegment + */ + public static LineSegment createLineSegment(Point p1, Point p2) { + double lon1 = p1.getLon(); + double lat1 = p1.getLat(); + double lon2 = p2.getLon(); + double lat2 = p2.getLat(); + + boolean vertical = (lon1 == lon2); + boolean horizontal = (lat1 == lat2); + + Double m; + Double b; + + if (vertical) { + m = Double.POSITIVE_INFINITY; + b = null; // Undefined for vertical lines + } else if (horizontal) { + m = 0.0; + b = lat1; + } else { + m = (lat2 - lat1) / (lon2 - lon1); + b = lat1 - (m * lon1); + } + + // Create MBR (does not cross antimeridian for cartesian polygons) + Mbr mbr = pointsToMbr(lon1, lat1, lon2, lat2); + + return new LineSegment(p1, p2, vertical, horizontal, m, b, mbr); + } + + /** + * Creates an MBR that covers both points but does not cross the antimeridian. + */ + private static Mbr pointsToMbr(double lon1, double lat1, double lon2, double lat2) { + double north = Math.max(lat1, lat2); + double south = Math.min(lat1, lat2); + double west, east; + + if (lon2 > lon1) { + west = lon1; + east = lon2; + } else { + west = lon2; + east = lon1; + } + + return new Mbr(west, north, east, south); + } + + /** + * Checks if a point lies on this line segment. + * Uses tolerance for approximate equality. + * + * @param point The point to test + * @return true if the point is on the segment (within tolerance) + */ + public boolean pointOnSegment(Point point) { + // Check if point equals either endpoint + if (MathUtils.doubleApprox(point.getLon(), point1.getLon(), COVERS_TOLERANCE) && + MathUtils.doubleApprox(point.getLat(), point1.getLat(), COVERS_TOLERANCE)) { + return true; + } + if (MathUtils.doubleApprox(point.getLon(), point2.getLon(), COVERS_TOLERANCE) && + MathUtils.doubleApprox(point.getLat(), point2.getLat(), COVERS_TOLERANCE)) { + return true; + } + + // Check if point is within MBR + if (!cartesianCoversPoint(mbr, point, COVERS_TOLERANCE)) { + return false; + } + + // For horizontal lines, check latitude + if (horizontal) { + return MathUtils.doubleApprox(point1.getLat(), point.getLat(), COVERS_TOLERANCE); + } + + // For other lines, calculate expected longitude at the point's latitude + Double expectedLon = latToLon(point.getLat()); + if (expectedLon == null) { + return false; + } + + return MathUtils.doubleApprox(expectedLon, point.getLon(), COVERS_TOLERANCE); + } + + /** + * Finds the intersection point between this line segment and another. + * Returns null if they don't intersect. + * + * @param other The other line segment + * @return The intersection point, or null if no intersection + */ + public Point intersection(LineSegment other) { + boolean ls1Vert = this.vertical; + boolean ls2Vert = other.vertical; + boolean ls1Horz = this.horizontal; + boolean ls2Horz = other.horizontal; + + // Handle vertical line cases + if (ls1Vert) { + if (ls2Vert) { + return intersectionBothVertical(this, other); + } else if (ls2Horz) { + return intersectionHorizontalAndVertical(other, this); + } else { + return intersectionOneVertical(this, other); + } + } + + if (ls2Vert) { + if (ls1Horz) { + return intersectionHorizontalAndVertical(this, other); + } else { + return intersectionOneVertical(other, this); + } + } + + // Handle horizontal line cases + if (ls1Horz && ls2Horz) { + return intersectionBothHorizontal(this, other); + } + + // Handle parallel lines (same slope) + if (this.m.equals(other.m)) { + return intersectionParallel(this, other); + } + + // Normal intersection + return intersectionNormal(this, other); + } + + /** + * Creates intermediate points along the line segment for densification. + * + * @param numPoints Number of intermediate points to create + * @return List of points including endpoints and intermediate points + */ + public List densifyLineSegment(int numPoints) { + if (numPoints <= 0) { + throw new IllegalArgumentException("numPoints must be > 0, got: " + numPoints); + } + + List points = new ArrayList<>(); + + // For vertical lines, just return endpoints + if (vertical) { + points.add(point1); + points.add(point2); + return points; + } + + double lon1 = point1.getLon(); + double lat1 = point1.getLat(); + double lon2 = point2.getLon(); + double lat2 = point2.getLat(); + + // Calculate differences per step + double lonDiff = (lon2 - lon1) / numPoints; + double latDiff = (lat2 - lat1) / numPoints; + + // Create points + for (int i = 0; i <= numPoints; i++) { + double lon = lon1 + (lonDiff * i); + double lat = lat1 + (latDiff * i); + points.add(new Point(lon, lat)); + } + + return points; + } + + /** + * Solves the line equation for longitude at a given latitude. + * For line equation y = mx + b, this returns x = (y - b) / m + * + * @param lat The latitude (y-coordinate) + * @return The longitude (x-coordinate) at the given latitude, or null if out of bounds + */ + public Double latToLon(double lat) { + if (horizontal) { + throw new IllegalStateException( + "Cannot determine longitude of points at a given latitude in a horizontal line"); + } + + if (vertical) { + // For vertical lines, check if lat is within range + if (lat >= mbr.getSouth() && lat <= mbr.getNorth()) { + return point1.getLon(); + } + return null; + } + + // Check if latitude is within MBR bounds + if (lat < mbr.getSouth() || lat > mbr.getNorth()) { + return null; + } + + // Solve for x: x = (y - b) / m + return (lat - b) / m; + } + + /** + * Converts an MBR to a list of line segments representing its edges. + * + * @param mbr The minimum bounding rectangle + * @return List of line segments forming the MBR edges + */ + public static List mbrToLineSegments(Mbr mbr) { + List segments = new ArrayList<>(); + + double west = mbr.getWest(); + double east = mbr.getEast(); + double north = mbr.getNorth(); + double south = mbr.getSouth(); + + // Check for degenerate cases + if (west == east && north == south) { + throw new IllegalArgumentException( + "This function doesn't work for an MBR that's a single point."); + } + + if (west == east) { + // Zero width MBR - single vertical line + segments.add(createLineSegment( + new Point(west, north), + new Point(east, south) + )); + return segments; + } + + if (north == south) { + // Zero height MBR - single or two horizontal lines + if (MbrIntersections.crossesAntimeridian(mbr)) { + segments.add(createLineSegment( + new Point(west, north), + new Point(180.0, north) + )); + segments.add(createLineSegment( + new Point(-180.0, north), + new Point(east, north) + )); + } else { + segments.add(createLineSegment( + new Point(west, north), + new Point(east, south) + )); + } + return segments; + } + + // Normal rectangle + if (MbrIntersections.crossesAntimeridian(mbr)) { + // Crosses antimeridian - need 6 segments + Point ul = new Point(west, north); + Point ur = new Point(east, north); + Point lr = new Point(east, south); + Point ll = new Point(west, south); + Point rightTop = new Point(180.0, north); + Point rightBot = new Point(180.0, south); + Point leftTop = new Point(-180.0, north); + Point leftBot = new Point(-180.0, south); + + segments.add(createLineSegment(ul, rightTop)); + segments.add(createLineSegment(leftTop, ur)); + segments.add(createLineSegment(ur, lr)); + segments.add(createLineSegment(lr, leftBot)); + segments.add(createLineSegment(rightBot, ll)); + segments.add(createLineSegment(ll, ul)); + } else { + // Normal rectangle - 4 segments + Point ul = new Point(west, north); + Point ur = new Point(east, north); + Point lr = new Point(east, south); + Point ll = new Point(west, south); + + segments.add(createLineSegment(ul, ur)); // Top + segments.add(createLineSegment(ur, lr)); // Right + segments.add(createLineSegment(lr, ll)); // Bottom + segments.add(createLineSegment(ll, ul)); // Left + } + + return segments; + } + + /** + * Creates line segments connecting consecutive points in a list. + * + * @param points List of points to connect + * @return List of line segments + */ + public static List pointsToLineSegments(List points) { + List segments = new ArrayList<>(); + + for (int i = 0; i < points.size() - 1; i++) { + segments.add(createLineSegment(points.get(i), points.get(i + 1))); + } + + return segments; + } + + // ========== Private intersection helper methods ========== + + /** + * Returns the intersection point of two vertical line segments. + */ + private static Point intersectionBothVertical(LineSegment ls1, LineSegment ls2) { + Mbr mbr1 = ls1.mbr; + Mbr mbr2 = ls2.mbr; + double lon1 = ls1.point1.getLon(); + double lon2 = ls2.point1.getLon(); + + if (lon1 != lon2) { + return null; + } + + double ls1North = mbr1.getNorth(); + double ls1South = mbr1.getSouth(); + double ls2North = mbr2.getNorth(); + double ls2South = mbr2.getSouth(); + + // Check which latitude ranges overlap + if (MathUtils.withinRange(ls2North, ls1South, ls1North)) { + return new Point(lon1, ls2North); + } else if (MathUtils.withinRange(ls2South, ls1South, ls1North)) { + return new Point(lon1, ls2South); + } else if (MathUtils.withinRange(ls1South, ls2South, ls2North)) { + return new Point(lon1, ls1South); + } + + return null; + } + + /** + * Returns the intersection point of one vertical and one non-vertical line. + */ + /** + * Returns the intersection of one vertical line and another non-vertical line. + * Matches Clojure intersection-one-vertical. + */ + private static Point intersectionOneVertical(LineSegment vertLs, LineSegment ls) { + double lon = vertLs.point1.getLon(); + Mbr mbr = ls.mbr; + Mbr vertMbr = vertLs.mbr; + + // Calculate latitude at the vertical line's longitude + Double lat = ls.lonToLat(lon); + if (lat == null) { + return null; + } + + Point point = new Point(lon, lat); + + // Check if point is within both MBRs (uses default tolerance, not INTERSECTION_COVERS_TOLERANCE) + if (cartesianCoversPoint(mbr, point, MBR_COVERS_TOLERANCE) && + cartesianCoversPoint(vertMbr, point, MBR_COVERS_TOLERANCE)) { + return point; + } + + return null; + } + + /** + * Returns the intersection point of two horizontal line segments. + */ + private static Point intersectionBothHorizontal(LineSegment ls1, LineSegment ls2) { + Mbr mbr1 = ls1.mbr; + Mbr mbr2 = ls2.mbr; + double lat1 = ls1.point1.getLat(); + double lat2 = ls2.point1.getLat(); + + if (lat1 != lat2) { + return null; + } + + double ls1East = mbr1.getEast(); + double ls1West = mbr1.getWest(); + double ls2East = mbr2.getEast(); + double ls2West = mbr2.getWest(); + + // Check which longitude ranges overlap + if (MathUtils.withinRange(ls2East, ls1West, ls1East)) { + return new Point(ls2East, lat1); + } else if (MathUtils.withinRange(ls2West, ls1West, ls1East)) { + return new Point(ls2West, lat1); + } else if (MathUtils.withinRange(ls1West, ls2West, ls2East)) { + return new Point(ls1West, lat1); + } + + return null; + } + + /** + * Returns the intersection of one horizontal and one vertical line. + */ + private static Point intersectionHorizontalAndVertical(LineSegment horizLs, LineSegment vertLs) { + Mbr horizMbr = horizLs.mbr; + double horizLat = horizMbr.getNorth(); + double horizWest = horizMbr.getWest(); + double horizEast = horizMbr.getEast(); + + Mbr vertMbr = vertLs.mbr; + double vertLon = vertMbr.getWest(); + double vertSouth = vertMbr.getSouth(); + double vertNorth = vertMbr.getNorth(); + + if (MathUtils.withinRange(horizLat, vertSouth, vertNorth) && + MathUtils.withinRange(vertLon, horizWest, horizEast)) { + return new Point(vertLon, horizLat); + } + + return null; + } + + /** + * Returns the intersection of two parallel line segments. + */ + private static Point intersectionParallel(LineSegment ls1, LineSegment ls2) { + // They only intersect if y-intercepts are the same + if (!ls1.b.equals(ls2.b)) { + return null; + } + + // Find the common intersecting MBR + Mbr intersectionMbr = mbrIntersection(ls1.mbr, ls2.mbr); + if (intersectionMbr == null) { + return null; + } + + // Use the longitude to find a point + double lon = intersectionMbr.getWest(); + Double lat = ls1.lonToLat(lon); + + if (lat != null) { + return new Point(lon, lat); + } + + return null; + } + + /** + * Returns the intersection of two normal (non-vertical, non-horizontal, non-parallel) line segments. + */ + private static Point intersectionNormal(LineSegment ls1, LineSegment ls2) { + double m1 = ls1.m; + double b1 = ls1.b; + double m2 = ls2.m; + double b2 = ls2.b; + Mbr mbr1 = ls1.mbr; + Mbr mbr2 = ls2.mbr; + + // Solve for intersection: m1*x + b1 = m2*x + b2 + // x = (b2 - b1) / (m1 - m2) + double lon = (b2 - b1) / (m1 - m2); + double lat = m1 * lon + b1; + + Point point = new Point(lon, lat); + + // Check if intersection point is within both MBRs + if (cartesianCoversPoint(mbr1, point, INTERSECTION_COVERS_TOLERANCE) && + cartesianCoversPoint(mbr2, point, INTERSECTION_COVERS_TOLERANCE)) { + return point; + } + + return null; + } + + /** + * Returns the latitude at the given longitude for this line segment. + */ + private Double lonToLat(double lon) { + if (vertical) { + throw new IllegalStateException( + "Cannot determine latitude of points at a given longitude in a vertical line"); + } + + // Check if longitude is within MBR bounds (with tolerance) + if (!cartesianLonRangeCoversLon(mbr.getWest(), mbr.getEast(), lon, COVERS_TOLERANCE)) { + return null; + } + + // Calculate: y = mx + b + return m * lon + b; + } + + /** + * Checks if a cartesian longitude range covers a specific longitude. + */ + private static boolean cartesianLonRangeCoversLon(double west, double east, double lon, double tolerance) { + return lon >= (west - tolerance) && lon <= (east + tolerance); + } + + /** + * Checks if an MBR covers a point in cartesian coordinates. + * Matches Clojure cartesian-covers-point? with tolerance. + */ + private static boolean cartesianCoversPoint(Mbr mbr, Point point, double tolerance) { + double lon = point.getLon(); + double lat = point.getLat(); + + // Check latitude (covers-lat?) + double north = mbr.getNorth() + tolerance; + double south = mbr.getSouth() - tolerance; + if (lat < south || lat > north) { + return false; + } + + // Check longitude (cartesian-lon-range-covers-lon?) + double west = mbr.getWest() - tolerance; + double east = mbr.getEast() + tolerance; + boolean crossesAntimeridian = west > east; // Check AFTER adjusting + + if (crossesAntimeridian) { + return lon >= west || lon <= east; + } else { + return lon >= west && lon <= east; + } + } + + /** + * Returns the intersection of two MBRs, or null if they don't intersect. + */ + private static Mbr mbrIntersection(Mbr mbr1, Mbr mbr2) { + if (!MbrIntersections.mbrsIntersect("cartesian", mbr1, mbr2)) { + return null; + } + + double west = Math.max(mbr1.getWest(), mbr2.getWest()); + double east = Math.min(mbr1.getEast(), mbr2.getEast()); + double north = Math.min(mbr1.getNorth(), mbr2.getNorth()); + double south = Math.max(mbr1.getSouth(), mbr2.getSouth()); + + return new Mbr(west, north, east, south); + } + + // ========== Getters ========== + + public Point getPoint1() { + return point1; + } + + public Point getPoint2() { + return point2; + } + + public boolean isVertical() { + return vertical; + } + + public boolean isHorizontal() { + return horizontal; + } + + public Double getM() { + return m; + } + + public Double getB() { + return b; + } + + public Mbr getMbr() { + return mbr; + } + + @Override + public String toString() { + return String.format("LineSegment{p1=%s, p2=%s, vertical=%b, horizontal=%b, m=%s, b=%s}", + point1, point2, vertical, horizontal, m, b); + } +} diff --git a/spatial-lib/src/java/cmr/spatial/math/CoordinateConversion.java b/spatial-lib/src/java/cmr/spatial/math/CoordinateConversion.java new file mode 100644 index 0000000000..fc9741d151 --- /dev/null +++ b/spatial-lib/src/java/cmr/spatial/math/CoordinateConversion.java @@ -0,0 +1,71 @@ +package cmr.spatial.math; + +import cmr.spatial.shape.Point; + +/** + * Coordinate conversion utilities for spherical geometry. + * Converts between lon/lat points and 3D vectors on unit sphere. + */ +public class CoordinateConversion { + + /** + * A safer version of the cross product for lon/lat points. + * Based on formula from http://williams.best.vwh.net/intersect.htm + * + * Original formulas (simplified to avoid duplicate calculations): + * x = sin(lat1-lat2) * sin((lon1+lon2)/2) * cos((lon1-lon2)/2) + * - sin(lat1+lat2) * cos((lon1+lon2)/2) * sin((lon1-lon2)/2) + * y = sin(lat1-lat2) * cos((lon1+lon2)/2) * cos((lon1-lon2)/2) + * + sin(lat1+lat2) * sin((lon1+lon2)/2) * sin((lon1-lon2)/2) + * z = cos(lat1) * cos(lat2) * sin(lon1-lon2) + */ + public static Vector lonLatCrossProduct(Point p1, Point p2) { + double lon1Rad = Math.toRadians(p1.getLon()); + double lat1Rad = Math.toRadians(p1.getLat()); + double lon2Rad = Math.toRadians(p2.getLon()); + double lat2Rad = Math.toRadians(p2.getLat()); + + double b = (lon1Rad + lon2Rad) / 2.0; + double sinB = Math.sin(b); + double cosB = Math.cos(b); + + double c = (lon1Rad - lon2Rad) / 2.0; + double e = Math.sin(lat1Rad - lat2Rad) * Math.cos(c); + double f = Math.sin(lat1Rad + lat2Rad) * Math.sin(c); + + double x = e * sinB - f * cosB; + double y = e * cosB + f * sinB; + double z = Math.cos(lat1Rad) * Math.cos(lat2Rad) * Math.sin(lon1Rad - lon2Rad); + + return new Vector(x, y, z); + } + + /** + * Converts a point (lon/lat) to a 3D vector on the unit sphere. + */ + public static Vector pointToVector(Point p) { + double lonRad = Math.toRadians(p.getLon()); + double latRad = Math.toRadians(p.getLat()); + + double cosLat = Math.cos(latRad); + double x = cosLat * Math.cos(lonRad); + double y = -1.0 * cosLat * Math.sin(lonRad); + double z = Math.sin(latRad); + + return new Vector(x, y, z); + } + + /** + * Converts a 3D vector to a point (lon/lat in degrees). + */ + public static Point vectorToPoint(Vector v) { + double x = v.getX(); + double y = v.getY(); + double z = v.getZ(); + + double lonRad = Math.atan2(-1.0 * y, x); + double latRad = Math.atan2(z, Math.sqrt(x * x + y * y)); + + return new Point(Math.toDegrees(lonRad), Math.toDegrees(latRad)); + } +} diff --git a/spatial-lib/src/java/cmr/spatial/math/MathUtils.java b/spatial-lib/src/java/cmr/spatial/math/MathUtils.java new file mode 100644 index 0000000000..9f99544da5 --- /dev/null +++ b/spatial-lib/src/java/cmr/spatial/math/MathUtils.java @@ -0,0 +1,121 @@ +package cmr.spatial.math; + +/** + * Utility functions for spatial math operations. + * Translates Clojure spatial math functions to Java for use in pure Java intersection testing. + */ +public class MathUtils { + + /** + * Delta tolerance for floating point approximate equality. + * Matches CMR spatial-lib DELTA constant. + */ + public static final double DELTA = 0.000000001; + + /** + * Tolerance for covers method. Longitudes and latitudes technically outside the + * bounding rectangle but within this tolerance will be considered covered. + */ + public static final double COVERS_TOLERANCE = 0.0000000001; + + /** + * Returns true if v is within min and max (inclusive). + * Implements the within-range? macro from Clojure. + */ + public static boolean withinRange(double v, double min, double max) { + return v >= min && v <= max; + } + + /** + * Returns true if two ranges intersect. + * Implements the range-intersects? macro from Clojure. + * Returns true if range2 intersects range1. + * + * @param r1min Lower bound of range 1 + * @param r1max Upper bound of range 1 + * @param r2min Lower bound of range 2 + * @param r2max Upper bound of range 2 + * @return true if the ranges intersect + */ + public static boolean rangeIntersects(double r1min, double r1max, + double r2min, double r2max) { + return withinRange(r2min, r1min, r1max) + || withinRange(r2max, r1min, r1max) + || withinRange(r1min, r2min, r2max); + } + + /** + * Returns true if a and b are approximately equal within DELTA tolerance. + * Implements approx= protocol for doubles. + */ + public static boolean doubleApprox(double a, double b) { + return doubleApprox(a, b, DELTA); + } + + /** + * Returns true if a and b are approximately equal within the given delta tolerance. + */ + public static boolean doubleApprox(double a, double b, double delta) { + return Math.abs(a - b) <= delta; + } + + /** + * Returns the absolute value of a double. + */ + public static double abs(double v) { + return Math.abs(v); + } + + /** + * Radius of Earth in meters for distance calculations. + * Matches EARTH_RADIUS_METERS in math.clj. + */ + public static final double EARTH_RADIUS_METERS = 6367435.0; + + /** + * Calculates the angular distance between two points in radians. + * Uses the Haversine formula. + * From: http://williams.best.vwh.net/avform.htm#Dist + * + * @param p1 First point + * @param p2 Second point + * @return Angular distance in radians + */ + public static double angularDistance(cmr.spatial.shape.Point p1, cmr.spatial.shape.Point p2) { + double lon1Rad = Math.toRadians(p1.getLon()); + double lat1Rad = Math.toRadians(p1.getLat()); + double lon2Rad = Math.toRadians(p2.getLon()); + double lat2Rad = Math.toRadians(p2.getLat()); + + // Haversine formula: a = sin²(Δlat/2) + cos(lat1) * cos(lat2) * sin²(Δlon/2) + double sinSqLat = sinSquared((lat1Rad - lat2Rad) / 2.0); + double sinSqLon = sinSquared((lon1Rad - lon2Rad) / 2.0); + double part1 = sinSqLat; + double part2 = Math.cos(lat1Rad) * Math.cos(lat2Rad) * sinSqLon; + + // Clamp to [0, 1] to prevent NaN from floating-point error + double x = part1 + part2; + x = Math.min(1.0, Math.max(0.0, x)); + + return 2.0 * Math.asin(Math.sqrt(x)); + } + + /** + * Calculates the distance between two points in meters. + * + * @param p1 First point + * @param p2 Second point + * @return Distance in meters + */ + public static double distance(cmr.spatial.shape.Point p1, cmr.spatial.shape.Point p2) { + return angularDistance(p1, p2) * EARTH_RADIUS_METERS; + } + + /** + * Helper for computing angular distance: sin²(v) = sin(v)² + */ + private static double sinSquared(double v) { + double s = Math.sin(v); + return s * s; + } +} diff --git a/spatial-lib/src/java/cmr/spatial/math/Vector.java b/spatial-lib/src/java/cmr/spatial/math/Vector.java new file mode 100644 index 0000000000..8171028228 --- /dev/null +++ b/spatial-lib/src/java/cmr/spatial/math/Vector.java @@ -0,0 +1,122 @@ +package cmr.spatial.math; + +/** + * Represents a 3D vector for spherical geometry calculations. + * Immutable value class. + */ +public final class Vector { + private final double x; + private final double y; + private final double z; + + public Vector(double x, double y, double z) { + this.x = x; + this.y = y; + this.z = z; + } + + public double getX() { + return x; + } + + public double getY() { + return y; + } + + public double getZ() { + return z; + } + + /** + * Calculates the length (magnitude) of this vector. + */ + public double length() { + return Math.sqrt(x * x + y * y + z * z); + } + + /** + * Returns a normalized version of this vector (length = 1). + * Throws IllegalArgumentException if vector length is 0. + */ + public Vector normalize() { + double len = length(); + if (len == 0.0) { + throw new IllegalArgumentException("Cannot normalize a vector with length 0"); + } + return new Vector(x / len, y / len, z / len); + } + + /** + * Calculates the cross product of this vector with another. + * Returns a vector perpendicular to both input vectors. + */ + public Vector crossProduct(Vector other) { + double px = this.y * other.z - other.y * this.z; + double py = this.z * other.x - other.z * this.x; + double pz = this.x * other.y - this.y * other.x; + return new Vector(px, py, pz); + } + + /** + * Returns a vector pointing in the opposite direction. + */ + public Vector opposite() { + return new Vector(-x, -y, -z); + } + + /** + * Vector equality tolerance: roughly equivalent to 0.9 meters on Earth's surface. + */ + public static final double VECTOR_EQUAL_DELTA = 0.0000001; + + /** + * Tests if two vectors are approximately equal within default tolerance. + */ + public boolean approxEquals(Vector other) { + return approxEquals(other, VECTOR_EQUAL_DELTA); + } + + /** + * Tests if two vectors are approximately equal within specified tolerance. + */ + public boolean approxEquals(Vector other, double delta) { + return MathUtils.doubleApprox(this.x, other.x, delta) + && MathUtils.doubleApprox(this.y, other.y, delta) + && MathUtils.doubleApprox(this.z, other.z, delta); + } + + /** + * Returns true if two vectors are parallel (point in same or opposite directions). + * Assumes both vectors are normalized. + */ + public boolean isParallel(Vector other) { + return this.approxEquals(other) || this.opposite().approxEquals(other); + } + + @Override + public String toString() { + return String.format("Vector{x=%.6f, y=%.6f, z=%.6f}", x, y, z); + } + + @Override + public boolean equals(Object obj) { + if (this == obj) return true; + if (!(obj instanceof Vector)) return false; + Vector other = (Vector) obj; + return Double.compare(x, other.x) == 0 + && Double.compare(y, other.y) == 0 + && Double.compare(z, other.z) == 0; + } + + @Override + public int hashCode() { + int result = 17; + long temp = Double.doubleToLongBits(x); + result = 31 * result + (int)(temp ^ (temp >>> 32)); + temp = Double.doubleToLongBits(y); + result = 31 * result + (int)(temp ^ (temp >>> 32)); + temp = Double.doubleToLongBits(z); + result = 31 * result + (int)(temp ^ (temp >>> 32)); + return result; + } +} diff --git a/spatial-lib/src/java/cmr/spatial/relations/ShapeIntersections.java b/spatial-lib/src/java/cmr/spatial/relations/ShapeIntersections.java new file mode 100644 index 0000000000..be76635c4d --- /dev/null +++ b/spatial-lib/src/java/cmr/spatial/relations/ShapeIntersections.java @@ -0,0 +1,250 @@ +package cmr.spatial.relations; + +import cmr.spatial.shape.*; +import cmr.spatial.geometry.PointIntersections; +import cmr.spatial.geometry.PointMbrIntersections; +import cmr.spatial.geometry.MbrIntersections; +import cmr.spatial.geometry.LineStringIntersections; +import cmr.spatial.geometry.RingIntersections; +import cmr.spatial.geometry.PolygonIntersections; +import cmr.spatial.geometry.CircleIntersections; + +/** + * Creates intersection testing functions for spatial shapes. + * Routes shape intersection tests to appropriate Java or Clojure implementations. + * + * Implemented in Java: + * - Point-to-Point, Point-to-Mbr + * - Mbr-to-Mbr + * - LineString, Ring, Polygon, and Circle intersections + * + * Delegated to Clojure: + * - Arc and LineSegment intersection tests (complex geometry algorithms) + * - Winding number algorithm for point-in-polygon + */ +public class ShapeIntersections { + + /** + * Creates a ShapePredicate that tests if shapes intersect with the given shape. + * Dispatches to appropriate intersection logic based on shape types. + * + * Supports: Point, Mbr, LineString, Ring, Polygon, Circle + * + * For Clojure shape records (GeodeticRing, CartesianRing, etc.), use + * the Clojure API (relations/shape->intersects-fn) instead. + * + * @param shape The target shape to test intersections against + * @return A ShapePredicate that returns true if another shape intersects this one + * @throws UnsupportedOperationException if shape type is not supported + */ + public static ShapePredicate createIntersectsFn(Object shape) { + if (shape instanceof Point) { + return createPointIntersectsFn((Point) shape); + } else if (shape instanceof Mbr) { + return createMbrIntersectsFn((Mbr) shape); + } else if (shape instanceof LineString) { + return createLineStringIntersectsFn((LineString) shape); + } else if (shape instanceof Ring) { + return createRingIntersectsFn((Ring) shape); + } else if (shape instanceof Polygon) { + return createPolygonIntersectsFn((Polygon) shape); + } else if (shape instanceof Circle) { + return createCircleIntersectsFn((Circle) shape); + } else { + throw new UnsupportedOperationException( + "createIntersectsFn not implemented for shape type: " + shape.getClass().getName() + + ". Use Clojure API (relations/shape->intersects-fn) for Clojure shape records."); + } + } + + /** + * Creates a predicate for Point intersection testing. + * Handles Point-to-Point, Point-to-Mbr, Point-to-LineString, Point-to-Ring, and Point-to-Polygon intersections. + */ + private static ShapePredicate createPointIntersectsFn(Point point) { + return otherShape -> { + if (otherShape instanceof Point) { + return PointIntersections.pointsIntersect(point, (Point) otherShape); + } else if (otherShape instanceof Mbr) { + // Point intersects Mbr if Mbr covers the point + return PointMbrIntersections.pointIntersectsMbr(point, (Mbr) otherShape); + } else if (otherShape instanceof LineString) { + // Point intersects LineString if LineString covers the point + return LineStringIntersections.coversPoint((LineString) otherShape, point); + } else if (otherShape instanceof Ring) { + // Point intersects Ring if Ring covers the point + Object javaRing = RingIntersections.createJavaRing((Ring) otherShape); + return RingIntersections.coversPoint(javaRing, point); + } else if (otherShape instanceof Polygon) { + // Point intersects Polygon if Polygon covers the point + return PolygonIntersections.coversPoint((Polygon) otherShape, point); + } else if (otherShape instanceof Circle) { + // Point intersects Circle if Circle covers the point + return CircleIntersections.coversPoint((Circle) otherShape, point); + } else { + throw new UnsupportedOperationException( + "Point intersection not implemented for shape type: " + otherShape.getClass().getName()); + } + }; + } + + /** + * Creates a predicate for Mbr intersection testing. + * Handles Mbr-to-Point, Mbr-to-Mbr, Mbr-to-LineString, Mbr-to-Ring, and Mbr-to-Polygon intersections. + */ + private static ShapePredicate createMbrIntersectsFn(Mbr mbr) { + return otherShape -> { + if (otherShape instanceof Point) { + // Mbr intersects Point if Mbr covers the point + return PointMbrIntersections.pointIntersectsMbr((Point) otherShape, mbr); + } else if (otherShape instanceof Mbr) { + return MbrIntersections.mbrsIntersect(mbr, (Mbr) otherShape); + } else if (otherShape instanceof LineString) { + // Mbr intersects LineString if LineString intersects the Mbr + return LineStringIntersections.intersectsMbr((LineString) otherShape, mbr); + } else if (otherShape instanceof Ring) { + // Mbr intersects Ring if Ring intersects the Mbr + Object javaRing = RingIntersections.createJavaRing((Ring) otherShape); + return RingIntersections.intersectsBr(javaRing, mbr); + } else if (otherShape instanceof Polygon) { + // Mbr intersects Polygon if Polygon intersects the Mbr + return PolygonIntersections.intersectsBr((Polygon) otherShape, mbr); + } else if (otherShape instanceof Circle) { + // Mbr intersects Circle if Circle intersects the Mbr + return CircleIntersections.intersectsBr((Circle) otherShape, mbr); + } else { + // Delegate to Clojure for other types (should not happen with ES plugin) + throw new UnsupportedOperationException( + "Intersection not implemented for shape type: " + otherShape.getClass().getName()); + } + }; + } + + /** + * Creates a predicate for LineString intersection testing. + * Handles LineString-to-Point, LineString-to-Mbr, LineString-to-LineString, + * LineString-to-Ring, and LineString-to-Polygon intersections. + */ + private static ShapePredicate createLineStringIntersectsFn(LineString lineString) { + return otherShape -> { + if (otherShape instanceof Point) { + return LineStringIntersections.coversPoint(lineString, (Point) otherShape); + } else if (otherShape instanceof Mbr) { + return LineStringIntersections.intersectsMbr(lineString, (Mbr) otherShape); + } else if (otherShape instanceof LineString) { + return LineStringIntersections.intersectsLineString(lineString, (LineString) otherShape); + } else if (otherShape instanceof Ring) { + // LineString intersects Ring if Ring intersects the LineString + Object javaRing = RingIntersections.createJavaRing((Ring) otherShape); + return RingIntersections.intersectsLineString(javaRing, lineString); + } else if (otherShape instanceof Polygon) { + // LineString intersects Polygon if Polygon intersects the LineString + return PolygonIntersections.intersectsLineString((Polygon) otherShape, lineString); + } else if (otherShape instanceof Circle) { + // LineString intersects Circle if Circle intersects the LineString + return CircleIntersections.intersectsLineString((Circle) otherShape, lineString); + } else { + // Delegate to Clojure for other types (should not happen with ES plugin) + throw new UnsupportedOperationException( + "Intersection not implemented for shape type: " + otherShape.getClass().getName()); + } + }; + } + + /** + * Creates a predicate for Ring intersection testing. + * Handles Ring-to-Point, Ring-to-Mbr, Ring-to-LineString, Ring-to-Ring, + * and Ring-to-Polygon intersections. + */ + private static ShapePredicate createRingIntersectsFn(Ring ring) { + Object javaRing = RingIntersections.createJavaRing(ring); + return otherShape -> { + if (otherShape instanceof Point) { + return RingIntersections.coversPoint(javaRing, (Point) otherShape); + } else if (otherShape instanceof Mbr) { + return RingIntersections.intersectsBr(javaRing, (Mbr) otherShape); + } else if (otherShape instanceof LineString) { + return RingIntersections.intersectsLineString(javaRing, (LineString) otherShape); + } else if (otherShape instanceof Ring) { + Object otherJavaRing = RingIntersections.createJavaRing((Ring) otherShape); + return RingIntersections.intersectsRing(javaRing, otherJavaRing); + } else if (otherShape instanceof Polygon) { + // Ring intersects Polygon if Polygon intersects the Ring + return PolygonIntersections.intersectsRing((Polygon) otherShape, ring); + } else if (otherShape instanceof Circle) { + // Ring intersects Circle if Circle intersects the Ring + return CircleIntersections.intersectsRing((Circle) otherShape, ring); + } else { + // Delegate to Clojure for other types (should not happen with ES plugin) + throw new UnsupportedOperationException( + "Intersection not implemented for shape type: " + otherShape.getClass().getName()); + } + }; + } + + /** + * Creates a predicate for Polygon intersection testing. + * Handles Polygon-to-Point, Polygon-to-Mbr, Polygon-to-LineString, Polygon-to-Ring, + * Polygon-to-Polygon, and Polygon-to-Circle intersections. + */ + private static ShapePredicate createPolygonIntersectsFn(Polygon polygon) { + return otherShape -> { + if (otherShape instanceof Point) { + return PolygonIntersections.coversPoint(polygon, (Point) otherShape); + } else if (otherShape instanceof Mbr) { + return PolygonIntersections.intersectsBr(polygon, (Mbr) otherShape); + } else if (otherShape instanceof LineString) { + return PolygonIntersections.intersectsLineString(polygon, (LineString) otherShape); + } else if (otherShape instanceof Ring) { + return PolygonIntersections.intersectsRing(polygon, (Ring) otherShape); + } else if (otherShape instanceof Polygon) { + return PolygonIntersections.intersectsPolygon(polygon, (Polygon) otherShape); + } else if (otherShape instanceof Circle) { + // Polygon intersects Circle if Circle intersects the Polygon + return CircleIntersections.intersectsPolygon((Circle) otherShape, polygon); + } else { + // Delegate to Clojure for other types (should not happen with ES plugin) + throw new UnsupportedOperationException( + "Intersection not implemented for shape type: " + otherShape.getClass().getName()); + } + }; + } + + /** + * Creates a predicate for Circle intersection testing. + * Handles Circle-to-Point, Circle-to-Mbr, Circle-to-LineString, Circle-to-Ring, + * and Circle-to-Polygon intersections. + */ + private static ShapePredicate createCircleIntersectsFn(Circle circle) { + return otherShape -> { + if (otherShape == null) { + throw new IllegalArgumentException("Intersection not implemented for null shape"); + } + + if (otherShape instanceof Point) { + return CircleIntersections.coversPoint(circle, (Point) otherShape); + } else if (otherShape instanceof Mbr) { + return CircleIntersections.intersectsBr(circle, (Mbr) otherShape); + } else if (otherShape instanceof LineString) { + return CircleIntersections.intersectsLineString(circle, (LineString) otherShape); + } else if (otherShape instanceof Ring) { + return CircleIntersections.intersectsRing(circle, (Ring) otherShape); + } else if (otherShape instanceof Polygon) { + return CircleIntersections.intersectsPolygon(circle, (Polygon) otherShape); + } else if (otherShape instanceof Circle) { + // Circle-to-circle: convert both to polygons + Polygon circlePoly = CircleIntersections.circleToPolygon(circle); + return CircleIntersections.intersectsPolygon((Circle) otherShape, circlePoly); + } else { + // Delegate to Clojure for other types + throw new UnsupportedOperationException( + "Intersection not implemented for shape type: " + otherShape.getClass().getName()); + } + }; + } + + // Clojure code uses protocol dispatch (relations.clj:shape->intersects-fn) + // and never calls this Java API with Clojure shape records. + // ES plugin only uses Java shapes (Point, Mbr, LineString, Ring, Polygon, Circle). + // Therefore, no Clojure delegation is needed. +} diff --git a/spatial-lib/src/java/cmr/spatial/relations/ShapePredicate.java b/spatial-lib/src/java/cmr/spatial/relations/ShapePredicate.java new file mode 100644 index 0000000000..414f8a65e5 --- /dev/null +++ b/spatial-lib/src/java/cmr/spatial/relations/ShapePredicate.java @@ -0,0 +1,16 @@ +package cmr.spatial.relations; + +/** + * A functional interface for testing if a shape intersects another shape. + * Used by the intersection testing logic in the spatial plugin. + */ +@FunctionalInterface +public interface ShapePredicate { + /** + * Tests if the given shape intersects the target shape. + * + * @param shape The shape to test + * @return true if the shape intersects the target shape + */ + boolean intersects(Object shape); +} diff --git a/spatial-lib/src/java/cmr/spatial/serialize/OrdsInfoShapes.java b/spatial-lib/src/java/cmr/spatial/serialize/OrdsInfoShapes.java new file mode 100644 index 0000000000..8a7463482d --- /dev/null +++ b/spatial-lib/src/java/cmr/spatial/serialize/OrdsInfoShapes.java @@ -0,0 +1,156 @@ +package cmr.spatial.serialize; + +import cmr.spatial.shape.*; +import java.util.*; + +/** + * Deserialization of ordinates and ords-info into spatial shapes. + * + * This converts the serialized format (ords-info + ords) back into shape objects. + * The ords-info is a sequence of [type, size, type, size, ...] pairs where: + * type - integer mapping to a shape type (1=geodetic-polygon, 2=geodetic-hole, etc) + * size - number of ordinates this shape uses + * The ords is a flat sequence of all ordinates. + */ +public class OrdsInfoShapes { + + private static final double MULTIPLICATION_FACTOR = 10000000.0; + + // Type mappings - must match cmr.spatial.serialize/shape-type->integer + private static final Map INT_TO_TYPE = new HashMap<>(); + static { + INT_TO_TYPE.put(1, "geodetic-polygon"); + INT_TO_TYPE.put(2, "geodetic-hole"); + INT_TO_TYPE.put(3, "br"); + INT_TO_TYPE.put(4, "point"); + INT_TO_TYPE.put(5, "geodetic-line-string"); + INT_TO_TYPE.put(6, "cartesian-polygon"); + INT_TO_TYPE.put(7, "cartesian-hole"); + INT_TO_TYPE.put(8, "cartesian-line-string"); + } + + /** + * Converts stored ordinate value to double for spatial calculations. + */ + public static double storedToOrdinate(long v) { + return ((double) v) / MULTIPLICATION_FACTOR; + } + + /** + * Converts ords-info and ords into a list of shapes. + * + * @param ordsInfo List of type/size pairs describing each shape (may be null for granules without spatial data) + * @param ords Flat list of all ordinates (may be null) + * @return List of SpatialShape objects (empty list if ordsInfo is null) + */ + public static List ordsInfoToShapes(List ordsInfo, List ords) { + // Handle nil case - granules without spatial data + if (ordsInfo == null) { + return new ArrayList<>(); + } + if (ordsInfo.size() % 2 != 0) { + throw new IllegalArgumentException( + String.format("ordsInfo must contain pairs of (type,size), got odd length: %d. Contents: %s", + ordsInfo.size(), ordsInfo)); + } + + List shapes = new ArrayList<>(); + List ordinates = convertToDoubleList(ords); + int ordsIndex = 0; + + // Process pairs of (type, size) from ordsInfo + for (int i = 0; i < ordsInfo.size(); i += 2) { + int intType = ((Number) ordsInfo.get(i)).intValue(); + int size = ((Number) ordsInfo.get(i + 1)).intValue(); + + String type = INT_TO_TYPE.get(intType); + if (type == null) { + throw new IllegalArgumentException( + "Could not get a shape type from integer [" + intType + "]. Ords info: " + ordsInfo); + } + + // Extract ordinates for this shape + if (ordsIndex + size > ordinates.size()) { + throw new IllegalArgumentException( + String.format("Ords info entry at index %d (type=%s, size=%d) exceeds ordinates length %d (ordsIndex=%d).", + i / 2, type, size, ordinates.size(), ordsIndex)); + } + List shapeOrds = new ArrayList<>(); + for (int j = 0; j < size; j++) { + shapeOrds.add(ordinates.get(ordsIndex + j)); + } + ordsIndex += size; + + // Convert ordinates to shape and add to results + SpatialShape shape = storedOrdsToShape(type, shapeOrds); + + // If this is a hole, add it to the rings of the last shape + if ("geodetic-hole".equals(type) || "cartesian-hole".equals(type)) { + if (!shapes.isEmpty() && shapes.get(shapes.size() - 1) instanceof Polygon) { + Polygon lastPolygon = (Polygon) shapes.get(shapes.size() - 1); + lastPolygon.getMutableRings().add((Ring) shape); + } else { + throw new IllegalStateException( + String.format("Orphaned hole shape without preceding polygon (type=%s, ordsInfo index=%d).", + type, i / 2)); + } + } else { + shapes.add(shape); + } + } + + return shapes; + } + + /** + * Converts a type and ordinates into a spatial shape. + */ + private static SpatialShape storedOrdsToShape(String type, List ords) { + switch (type) { + case "geodetic-polygon": + return new Polygon("geodetic", new ArrayList<>(Arrays.asList(new Ring("geodetic", ords)))); + + case "cartesian-polygon": + return new Polygon("cartesian", new ArrayList<>(Arrays.asList(new Ring("cartesian", ords)))); + + case "geodetic-hole": + return new Ring("geodetic", ords, true); + + case "cartesian-hole": + return new Ring("cartesian", ords, true); + + case "br": + if (ords.size() != 4) { + throw new IllegalArgumentException("MBR requires exactly 4 ordinates, got " + ords.size()); + } + return new Mbr(ords.get(0), ords.get(1), ords.get(2), ords.get(3)); + + case "point": + if (ords.size() != 2) { + throw new IllegalArgumentException("Point requires exactly 2 ordinates, got " + ords.size()); + } + return new Point(ords.get(0), ords.get(1)); + + case "geodetic-line-string": + return new LineString("geodetic", ords); + + case "cartesian-line-string": + return new LineString("cartesian", ords); + + default: + throw new IllegalArgumentException("Unknown ords shape type: " + type); + } + } + + /** + * Converts a list of numbers to a list of doubles, converting stored ordinates to actual values. + */ + private static List convertToDoubleList(List input) { + List result = new ArrayList<>(); + for (Object obj : input) { + long val = ((Number) obj).longValue(); + result.add(storedToOrdinate(val)); + } + return result; + } +} diff --git a/spatial-lib/src/java/cmr/spatial/shape/Circle.java b/spatial-lib/src/java/cmr/spatial/shape/Circle.java new file mode 100644 index 0000000000..61ad063c3b --- /dev/null +++ b/spatial-lib/src/java/cmr/spatial/shape/Circle.java @@ -0,0 +1,33 @@ +package cmr.spatial.shape; + +/** + * Java representation of a Circle shape. + * A circle is defined by a center point and a radius in meters. + */ +public class Circle implements SpatialShape { + private final Point center; + private final double radius; + + public Circle(Point center, double radius) { + this.center = center; + this.radius = radius; + } + + public Point getCenter() { + return center; + } + + public double getRadius() { + return radius; + } + + @Override + public String getType() { + return "circle"; + } + + @Override + public String toString() { + return String.format("Circle{center=%s, radius=%.2f}", center, radius); + } +} diff --git a/spatial-lib/src/java/cmr/spatial/shape/LineString.java b/spatial-lib/src/java/cmr/spatial/shape/LineString.java new file mode 100644 index 0000000000..92902c8776 --- /dev/null +++ b/spatial-lib/src/java/cmr/spatial/shape/LineString.java @@ -0,0 +1,34 @@ +package cmr.spatial.shape; + +import java.util.*; + +/** + * Java representation of a LineString shape. + */ +public class LineString implements SpatialShape { + private final String coordinateSystem; + private final List ordinates; + + public LineString(String coordinateSystem, List ordinates) { + this.coordinateSystem = coordinateSystem; + this.ordinates = new ArrayList<>(ordinates); + } + + public String getCoordinateSystem() { + return coordinateSystem; + } + + public List getOrdinates() { + return Collections.unmodifiableList(ordinates); + } + + @Override + public String getType() { + return coordinateSystem + "-line-string"; + } + + @Override + public String toString() { + return String.format("LineString{cs=%s, ords=%d}", coordinateSystem, ordinates.size()); + } +} diff --git a/spatial-lib/src/java/cmr/spatial/shape/Mbr.java b/spatial-lib/src/java/cmr/spatial/shape/Mbr.java new file mode 100644 index 0000000000..f627e001b4 --- /dev/null +++ b/spatial-lib/src/java/cmr/spatial/shape/Mbr.java @@ -0,0 +1,45 @@ +package cmr.spatial.shape; + +/** + * Java representation of a Minimum Bounding Rectangle (MBR). + * Stores west, north, east, south coordinates. + */ +public class Mbr implements SpatialShape { + private final double west; + private final double north; + private final double east; + private final double south; + + public Mbr(double west, double north, double east, double south) { + this.west = west; + this.north = north; + this.east = east; + this.south = south; + } + + public double getWest() { + return west; + } + + public double getNorth() { + return north; + } + + public double getEast() { + return east; + } + + public double getSouth() { + return south; + } + + @Override + public String getType() { + return "br"; + } + + @Override + public String toString() { + return String.format("Mbr{w=%f, n=%f, e=%f, s=%f}", west, north, east, south); + } +} diff --git a/spatial-lib/src/java/cmr/spatial/shape/Point.java b/spatial-lib/src/java/cmr/spatial/shape/Point.java new file mode 100644 index 0000000000..a7a8c4a017 --- /dev/null +++ b/spatial-lib/src/java/cmr/spatial/shape/Point.java @@ -0,0 +1,133 @@ +package cmr.spatial.shape; + +/** + * Java representation of a Point shape. + * Implements geodetic equality semantics matching the Clojure Point implementation. + */ +public class Point implements SpatialShape { + private static final double APPROXIMATION_DELTA = 0.000000001; + private static final int PRIME = 31; + + // Precomputed hash codes matching Clojure Point implementation + private static final int NORTH_POLE_HASH; + private static final int SOUTH_POLE_HASH; + private static final int INITIAL_AM_HASH; + + static { + // Matches Clojure: (int (+ (* PRIME (+ PRIME lon-hash)) (.hashCode (Double. 90.0)))) + int lonHash = Double.hashCode(0.0); + NORTH_POLE_HASH = PRIME * (PRIME + lonHash) + Double.hashCode(90.0); + SOUTH_POLE_HASH = PRIME * (PRIME + lonHash) + Double.hashCode(-90.0); + // Matches Clojure: (* PRIME (+ PRIME (.hashCode (Double. 180.0)))) + INITIAL_AM_HASH = PRIME * (PRIME + Double.hashCode(180.0)); + } + + private final double lon; + private final double lat; + + public Point(double lon, double lat) { + this.lon = lon; + this.lat = lat; + } + + public double getLon() { + return lon; + } + + public double getLat() { + return lat; + } + + @Override + public String getType() { + return "point"; + } + + @Override + public String toString() { + return String.format("Point{lon=%f, lat=%f}", lon, lat); + } + + /** + * Checks if this point is at the north pole (latitude ~90). + */ + private boolean isNorthPole() { + return Math.abs(lat - 90.0) < APPROXIMATION_DELTA; + } + + /** + * Checks if this point is at the south pole (latitude ~-90). + */ + private boolean isSouthPole() { + return Math.abs(lat + 90.0) < APPROXIMATION_DELTA; + } + + /** + * Checks if this point is on the antimeridian (longitude ~180 or ~-180). + */ + private boolean onAntimeridian() { + return Math.abs(Math.abs(lon) - 180.0) < APPROXIMATION_DELTA; + } + + /** + * Normalizes a coordinate value to a grid aligned with APPROXIMATION_DELTA. + * This ensures values within epsilon of each other get normalized to the same value. + */ + private static double normalize(double value) { + return Math.round(value / APPROXIMATION_DELTA) * APPROXIMATION_DELTA; + } + + @Override + public boolean equals(Object other) { + if (this == other) return true; + if (other == null || getClass() != other.getClass()) return false; + + Point otherPoint = (Point) other; + + // Geodetic special cases + // Any longitude value at north pole is considered equivalent + if (isNorthPole() && otherPoint.isNorthPole()) { + return true; + } + + // Any longitude value at south pole is considered equivalent + if (isSouthPole() && otherPoint.isSouthPole()) { + return true; + } + + // -180 and 180 are considered equivalent longitudes + if (onAntimeridian() && otherPoint.onAntimeridian()) { + return Math.abs(lat - otherPoint.lat) < APPROXIMATION_DELTA; + } + + // Normal equality with normalized comparison + double normalizedLon = normalize(lon); + double normalizedOtherLon = normalize(otherPoint.lon); + double normalizedLat = normalize(lat); + double normalizedOtherLat = normalize(otherPoint.lat); + + return normalizedLon == normalizedOtherLon && normalizedLat == normalizedOtherLat; + } + + @Override + public int hashCode() { + // Geodetic special cases for consistent hashing + if (isNorthPole()) { + return NORTH_POLE_HASH; + } + if (isSouthPole()) { + return SOUTH_POLE_HASH; + } + if (Math.abs(Math.abs(lon) - 180.0) < APPROXIMATION_DELTA) { + // Points on antimeridian hash to same value regardless of sign + // Use normalized latitude for consistent hashing + int latHash = Double.hashCode(normalize(lat)); + return INITIAL_AM_HASH + latHash; + } + + // Normal case: use normalized values for consistent hashing + int lonHash = Double.hashCode(normalize(lon)); + int latHash = Double.hashCode(normalize(lat)); + return PRIME * (PRIME + lonHash) + latHash; + } +} diff --git a/spatial-lib/src/java/cmr/spatial/shape/Polygon.java b/spatial-lib/src/java/cmr/spatial/shape/Polygon.java new file mode 100644 index 0000000000..8e920c9aeb --- /dev/null +++ b/spatial-lib/src/java/cmr/spatial/shape/Polygon.java @@ -0,0 +1,52 @@ +package cmr.spatial.shape; + +import java.util.*; + +/** + * Java representation of a Polygon shape. + * Consists of an outer ring (boundary) and optional inner rings (holes). + */ +public class Polygon implements SpatialShape { + private final String coordinateSystem; + private final List rings; + + public Polygon(String coordinateSystem, List rings) { + this.coordinateSystem = coordinateSystem; + this.rings = new ArrayList<>(rings); + } + + public String getCoordinateSystem() { + return coordinateSystem; + } + + public List getRings() { + return Collections.unmodifiableList(rings); + } + + /** + * Internal method for shape deserialization only. + * Provides mutable access to rings list for OrdsInfoShapes to add holes during construction. + * Do not use in application code - use getRings() instead. + */ + public List getMutableRings() { + return rings; + } + + public Ring getBoundary() { + return rings.isEmpty() ? null : rings.get(0); + } + + public List getHoles() { + return rings.size() <= 1 ? Collections.emptyList() : rings.subList(1, rings.size()); + } + + @Override + public String getType() { + return coordinateSystem + "-polygon"; + } + + @Override + public String toString() { + return String.format("Polygon{cs=%s, rings=%d}", coordinateSystem, rings.size()); + } +} diff --git a/spatial-lib/src/java/cmr/spatial/shape/Ring.java b/spatial-lib/src/java/cmr/spatial/shape/Ring.java new file mode 100644 index 0000000000..eaef9c95ea --- /dev/null +++ b/spatial-lib/src/java/cmr/spatial/shape/Ring.java @@ -0,0 +1,45 @@ +package cmr.spatial.shape; + +import java.util.*; + +/** + * Java representation of a Ring (closed sequence of points). + * Used as a component of Polygon shapes, and as a hole in Polygon objects. + */ +public class Ring implements SpatialShape { + private final String coordinateSystem; + private final List ordinates; + private final boolean isHole; + + public Ring(String coordinateSystem, List ordinates) { + this(coordinateSystem, ordinates, false); + } + + public Ring(String coordinateSystem, List ordinates, boolean isHole) { + this.coordinateSystem = coordinateSystem; + this.ordinates = new ArrayList<>(ordinates); + this.isHole = isHole; + } + + public String getCoordinateSystem() { + return coordinateSystem; + } + + public List getOrdinates() { + return Collections.unmodifiableList(ordinates); + } + + public boolean isHole() { + return isHole; + } + + @Override + public String getType() { + return coordinateSystem + (isHole ? "-hole" : "-ring"); + } + + @Override + public String toString() { + return String.format("Ring{cs=%s, ords=%d}", coordinateSystem, ordinates.size()); + } +} diff --git a/spatial-lib/src/java/cmr/spatial/shape/SpatialShape.java b/spatial-lib/src/java/cmr/spatial/shape/SpatialShape.java new file mode 100644 index 0000000000..1863bb7649 --- /dev/null +++ b/spatial-lib/src/java/cmr/spatial/shape/SpatialShape.java @@ -0,0 +1,9 @@ +package cmr.spatial.shape; + +/** + * Base interface for spatial shapes used by the Java implementation. + * These are minimal POJOs used internally for shape deserialization and intersection testing. + */ +public interface SpatialShape { + String getType(); +} diff --git a/spatial-lib/src/test/java/cmr/spatial/geometry/IntersectionTestsTest.java b/spatial-lib/src/test/java/cmr/spatial/geometry/IntersectionTestsTest.java new file mode 100644 index 0000000000..54fafde6d1 --- /dev/null +++ b/spatial-lib/src/test/java/cmr/spatial/geometry/IntersectionTestsTest.java @@ -0,0 +1,204 @@ +package cmr.spatial.geometry; + +import cmr.spatial.shape.Point; +import cmr.spatial.shape.Mbr; +import org.junit.Test; + +import static org.junit.Assert.*; + +public class IntersectionTestsTest { + + // ==================== Point-to-Point Intersection Tests ==================== + + @Test + public void testPointsIntersect_SamePoint() { + Point p1 = new Point(10.0, 20.0); + Point p2 = new Point(10.0, 20.0); + assertTrue("Points at same location should intersect", PointIntersections.pointsIntersect(p1, p2)); + } + + @Test + public void testPointsIntersect_DifferentPoints() { + Point p1 = new Point(10.0, 20.0); + Point p2 = new Point(11.0, 21.0); + assertFalse("Different points should not intersect", PointIntersections.pointsIntersect(p1, p2)); + } + + @Test + public void testPointsIntersect_ApproximatelyEqual() { + Point p1 = new Point(10.0, 20.0); + Point p2 = new Point(10.0 + 1e-10, 20.0 + 1e-10); // Within DELTA + assertTrue("Points within DELTA should be approximately equal", PointIntersections.pointsIntersect(p1, p2)); + } + + @Test + public void testPointsIntersect_NorthPole() { + Point p1 = new Point(0.0, 90.0); + Point p2 = new Point(180.0, 90.0); // Different longitude at north pole + assertTrue("Any longitude at north pole should be equivalent", PointIntersections.pointsIntersect(p1, p2)); + } + + @Test + public void testPointsIntersect_SouthPole() { + Point p1 = new Point(0.0, -90.0); + Point p2 = new Point(180.0, -90.0); // Different longitude at south pole + assertTrue("Any longitude at south pole should be equivalent", PointIntersections.pointsIntersect(p1, p2)); + } + + @Test + public void testPointsIntersect_Antimeridian() { + Point p1 = new Point(180.0, 45.0); + Point p2 = new Point(-180.0, 45.0); // Both on antimeridian + assertTrue("-180 and 180 should be equivalent at same latitude", PointIntersections.pointsIntersect(p1, p2)); + } + + @Test + public void testPointsIntersect_DifferentLatitudes() { + Point p1 = new Point(10.0, 20.0); + Point p2 = new Point(10.0, 21.0); // Same longitude, different latitude + assertFalse("Points at different latitudes should not intersect", PointIntersections.pointsIntersect(p1, p2)); + } + + // ==================== Point-to-MBR Intersection Tests ==================== + + @Test + public void testPointIntersectsMbr_PointInside() { + Point point = new Point(10.0, 20.0); + Mbr mbr = new Mbr(-30.0, 50.0, 30.0, -10.0); // w, n, e, s + assertTrue("Point inside MBR should intersect", PointMbrIntersections.pointIntersectsMbr(point, mbr)); + } + + @Test + public void testPointIntersectsMbr_PointOutside() { + Point point = new Point(50.0, 20.0); + Mbr mbr = new Mbr(-30.0, 50.0, 30.0, -10.0); // w, n, e, s (east bound is 30) + assertFalse("Point outside MBR should not intersect", PointMbrIntersections.pointIntersectsMbr(point, mbr)); + } + + @Test + public void testPointIntersectsMbr_PointOnEdge() { + Point point = new Point(30.0, 20.0); + Mbr mbr = new Mbr(-30.0, 50.0, 30.0, -10.0); // w, n, e, s (east bound is 30) + assertTrue("Point on edge of MBR should intersect", PointMbrIntersections.pointIntersectsMbr(point, mbr)); + } + + @Test + public void testPointIntersectsMbr_NorthPoleTouches() { + Point point = new Point(0.0, 90.0); + Mbr mbr = new Mbr(-30.0, 90.0, 30.0, -10.0); // w, n, e, s (north at 90) + assertTrue("North pole point should intersect MBR touching north pole", + PointMbrIntersections.pointIntersectsMbr(point, mbr)); + } + + @Test + public void testPointIntersectsMbr_SouthPoleTouches() { + Point point = new Point(0.0, -90.0); + Mbr mbr = new Mbr(-30.0, 50.0, 30.0, -90.0); // w, n, e, s (south at -90) + assertTrue("South pole point should intersect MBR touching south pole", + PointMbrIntersections.pointIntersectsMbr(point, mbr)); + } + + @Test + public void testPointIntersectsMbr_AntimeridianCrossingMbr() { + Point point = new Point(170.0, 20.0); // Near antimeridian on positive side + Mbr mbr = new Mbr(160.0, 50.0, -160.0, -10.0); // Crosses antimeridian + assertTrue("Point near antimeridian should intersect MBR crossing antimeridian", + PointMbrIntersections.pointIntersectsMbr(point, mbr)); + } + + // ==================== MBR-to-MBR Intersection Tests ==================== + + @Test + public void testMbrsIntersect_SimpleOverlap() { + Mbr mbr1 = new Mbr(-30.0, 50.0, 30.0, -10.0); // w, n, e, s + Mbr mbr2 = new Mbr(0.0, 40.0, 60.0, 10.0); // Overlaps with mbr1 + assertTrue("Overlapping MBRs should intersect", MbrIntersections.mbrsIntersect(mbr1, mbr2)); + } + + @Test + public void testMbrsIntersect_NoOverlap() { + Mbr mbr1 = new Mbr(-30.0, 50.0, -10.0, -10.0); + Mbr mbr2 = new Mbr(40.0, 50.0, 60.0, 10.0); // Completely to the east + assertFalse("Non-overlapping MBRs should not intersect", MbrIntersections.mbrsIntersect(mbr1, mbr2)); + } + + @Test + public void testMbrsIntersect_EdgeTouch() { + Mbr mbr1 = new Mbr(-30.0, 50.0, 0.0, -10.0); + Mbr mbr2 = new Mbr(0.0, 40.0, 30.0, 10.0); // Touch at edge (lon=0) + assertTrue("MBRs touching at edge should intersect", MbrIntersections.mbrsIntersect(mbr1, mbr2)); + } + + @Test + public void testMbrsIntersect_PolesTouch() { + Mbr mbr1 = new Mbr(-30.0, 90.0, 30.0, 50.0); // Touches north pole + Mbr mbr2 = new Mbr(0.0, 90.0, 60.0, 40.0); // Also touches north pole + assertTrue("MBRs touching at north pole should intersect in geodetic", + MbrIntersections.mbrsIntersect("geodetic", mbr1, mbr2)); + } + + @Test + public void testMbrsIntersect_ContainedMbr() { + Mbr mbr1 = new Mbr(-30.0, 50.0, 30.0, -10.0); + Mbr mbr2 = new Mbr(-10.0, 30.0, 10.0, 0.0); // Completely inside mbr1 + assertTrue("Contained MBR should intersect", MbrIntersections.mbrsIntersect(mbr1, mbr2)); + } + + @Test + public void testMbrsIntersect_AntimeridianNonCrossing() { + Mbr mbr1 = new Mbr(150.0, 50.0, 170.0, -10.0); + Mbr mbr2 = new Mbr(160.0, 40.0, 180.0, 0.0); // Overlaps, neither crosses + assertTrue("Non-crossing MBRs near antimeridian should intersect if they overlap", + MbrIntersections.mbrsIntersect(mbr1, mbr2)); + } + + @Test + public void testMbrsIntersect_OneCrossesAntimeridian() { + Mbr mbr1 = new Mbr(160.0, 50.0, -160.0, -10.0); // Crosses antimeridian + Mbr mbr2 = new Mbr(-170.0, 40.0, -150.0, 0.0); // On the other side + assertTrue("MBR crossing antimeridian should intersect with MBR on far side", + MbrIntersections.mbrsIntersect(mbr1, mbr2)); + } + + @Test + public void testMbrsIntersect_BothCrossAntimeridian() { + Mbr mbr1 = new Mbr(160.0, 50.0, -170.0, -10.0); // Crosses antimeridian + Mbr mbr2 = new Mbr(170.0, 40.0, -150.0, 0.0); // Also crosses antimeridian + assertTrue("Both MBRs crossing antimeridian should intersect", + MbrIntersections.mbrsIntersect(mbr1, mbr2)); + } + + // ==================== Edge Case Tests ==================== + + @Test + public void testCrossesAntimeridian() { + Mbr mbr_crosses = new Mbr(160.0, 50.0, -160.0, -10.0); + Mbr mbr_not_crosses = new Mbr(-30.0, 50.0, 30.0, -10.0); + + assertTrue("MBR with west > east should cross antimeridian", + MbrIntersections.crossesAntimeridian(mbr_crosses)); + assertFalse("MBR with west < east should not cross antimeridian", + MbrIntersections.crossesAntimeridian(mbr_not_crosses)); + } + + @Test + public void testSplitAcrossAntimeridian() { + Mbr mbr_crosses = new Mbr(160.0, 50.0, -160.0, -10.0); + var parts = MbrIntersections.splitAcrossAntimeridian(mbr_crosses); + + assertEquals("Crossing MBR should split into 2 parts", 2, parts.size()); + assertEquals("First part should go from west to 180", 160.0, parts.get(0).getWest(), 0.0); + assertEquals("First part should go from west to 180", 180.0, parts.get(0).getEast(), 0.0); + assertEquals("Second part should go from -180 to east", -180.0, parts.get(1).getWest(), 0.0); + assertEquals("Second part should go from -180 to east", -160.0, parts.get(1).getEast(), 0.0); + } + + @Test + public void testSplitAcrossAntimeridian_NoCross() { + Mbr mbr_not_crosses = new Mbr(-30.0, 50.0, 30.0, -10.0); + var parts = MbrIntersections.splitAcrossAntimeridian(mbr_not_crosses); + + assertEquals("Non-crossing MBR should return 1 part", 1, parts.size()); + assertEquals("Part should be the same as original", mbr_not_crosses, parts.get(0)); + } +} diff --git a/spatial-lib/test/cmr/spatial/test/ESPluginJavaIntersectionsTest.java b/spatial-lib/test/cmr/spatial/test/ESPluginJavaIntersectionsTest.java new file mode 100644 index 0000000000..952373996a --- /dev/null +++ b/spatial-lib/test/cmr/spatial/test/ESPluginJavaIntersectionsTest.java @@ -0,0 +1,189 @@ +package cmr.spatial.test; + +import cmr.spatial.geometry.*; +import cmr.spatial.shape.*; +import cmr.spatial.relations.ShapeIntersections; +import cmr.spatial.relations.ShapePredicate; +import cmr.spatial.serialize.OrdsInfoShapes; +import org.junit.Test; +import static org.junit.Assert.*; + +import java.util.Arrays; +import java.util.List; + +/** + * Tests to verify that all ES plugin shape intersections work in pure Java + * without any Clojure delegation. + * + * The ES plugin uses OrdsInfoShapes to deserialize shapes from stored format, + * which creates Java shapes (Point, Mbr, LineString, Ring, Polygon). + * All shape-to-shape intersections must work without calling Clojure. + */ +public class ESPluginJavaIntersectionsTest { + + @Test + public void testPointIntersections() { + Point p1 = new Point(0, 0); + Point p2 = new Point(0, 0); + Point p3 = new Point(10, 10); + + // Point-to-Point + ShapePredicate pred = ShapeIntersections.createIntersectsFn(p1); + assertTrue("Point intersects itself", pred.test(p2)); + assertFalse("Point does not intersect distant point", pred.test(p3)); + + // Point-to-Mbr + Mbr mbr = new Mbr(-1, 1, 1, -1); + assertTrue("Point intersects containing MBR", pred.test(mbr)); + + // Point-to-LineString + LineString line = new LineString("geodetic", Arrays.asList(-1.0, 0.0, 1.0, 0.0)); + assertTrue("Point intersects line through it", pred.test(line)); + + // Point-to-Ring + Ring ring = new Ring("geodetic", Arrays.asList(-2.0, -2.0, 2.0, -2.0, 2.0, 2.0, -2.0, 2.0, -2.0, -2.0)); + assertTrue("Point intersects containing ring", pred.test(ring)); + + // Point-to-Polygon + Polygon polygon = new Polygon("geodetic", Arrays.asList(ring)); + assertTrue("Point intersects containing polygon", pred.test(polygon)); + } + + @Test + public void testMbrIntersections() { + Mbr mbr1 = new Mbr(-1, 1, 1, -1); + Mbr mbr2 = new Mbr(0, 2, 2, 0); + + ShapePredicate pred = ShapeIntersections.createIntersectsFn(mbr1); + + // Mbr-to-Point + Point p = new Point(0, 0); + assertTrue("MBR intersects contained point", pred.test(p)); + + // Mbr-to-Mbr + assertTrue("MBR intersects overlapping MBR", pred.test(mbr2)); + + // Mbr-to-LineString + LineString line = new LineString("geodetic", Arrays.asList(-0.5, 0.0, 0.5, 0.0)); + assertTrue("MBR intersects line through it", pred.test(line)); + + // Mbr-to-Ring + Ring ring = new Ring("geodetic", Arrays.asList(-0.5, -0.5, 0.5, -0.5, 0.5, 0.5, -0.5, 0.5, -0.5, -0.5)); + assertTrue("MBR intersects ring inside it", pred.test(ring)); + + // Mbr-to-Polygon + Polygon polygon = new Polygon("geodetic", Arrays.asList(ring)); + assertTrue("MBR intersects polygon inside it", pred.test(polygon)); + } + + @Test + public void testLineStringIntersections() { + LineString line1 = new LineString("geodetic", Arrays.asList(-1.0, 0.0, 1.0, 0.0)); + + ShapePredicate pred = ShapeIntersections.createIntersectsFn(line1); + + // LineString-to-Point + Point p = new Point(0, 0); + assertTrue("LineString intersects point on it", pred.test(p)); + + // LineString-to-Mbr + Mbr mbr = new Mbr(-0.5, 0.5, 0.5, -0.5); + assertTrue("LineString intersects MBR it passes through", pred.test(mbr)); + + // LineString-to-LineString + LineString line2 = new LineString("geodetic", Arrays.asList(0.0, -1.0, 0.0, 1.0)); + assertTrue("LineString intersects crossing line", pred.test(line2)); + + // LineString-to-Ring + Ring ring = new Ring("geodetic", Arrays.asList(-2.0, -2.0, 2.0, -2.0, 2.0, 2.0, -2.0, 2.0, -2.0, -2.0)); + assertTrue("LineString intersects containing ring", pred.test(ring)); + + // LineString-to-Polygon + Polygon polygon = new Polygon("geodetic", Arrays.asList(ring)); + assertTrue("LineString intersects containing polygon", pred.test(polygon)); + } + + @Test + public void testRingIntersections() { + Ring ring1 = new Ring("geodetic", Arrays.asList(-2.0, -2.0, 2.0, -2.0, 2.0, 2.0, -2.0, 2.0, -2.0, -2.0)); + + ShapePredicate pred = ShapeIntersections.createIntersectsFn(ring1); + + // Ring-to-Point + Point p = new Point(0, 0); + assertTrue("Ring intersects point inside it", pred.test(p)); + + // Ring-to-Mbr + Mbr mbr = new Mbr(-1, 1, 1, -1); + assertTrue("Ring intersects MBR inside it", pred.test(mbr)); + + // Ring-to-LineString + LineString line = new LineString("geodetic", Arrays.asList(-1.0, 0.0, 1.0, 0.0)); + assertTrue("Ring intersects line inside it", pred.test(line)); + + // Ring-to-Ring + Ring ring2 = new Ring("geodetic", Arrays.asList(-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0)); + assertTrue("Ring intersects ring inside it", pred.test(ring2)); + + // Ring-to-Polygon + Polygon polygon = new Polygon("geodetic", Arrays.asList(ring2)); + assertTrue("Ring intersects polygon inside it", pred.test(polygon)); + } + + @Test + public void testPolygonIntersections() { + Ring outerRing = new Ring("geodetic", Arrays.asList(-2.0, -2.0, 2.0, -2.0, 2.0, 2.0, -2.0, 2.0, -2.0, -2.0)); + Polygon polygon1 = new Polygon("geodetic", Arrays.asList(outerRing)); + + ShapePredicate pred = ShapeIntersections.createIntersectsFn(polygon1); + + // Polygon-to-Point + Point p = new Point(0, 0); + assertTrue("Polygon intersects point inside it", pred.test(p)); + + // Polygon-to-Mbr + Mbr mbr = new Mbr(-1, 1, 1, -1); + assertTrue("Polygon intersects MBR inside it", pred.test(mbr)); + + // Polygon-to-LineString + LineString line = new LineString("geodetic", Arrays.asList(-1.0, 0.0, 1.0, 0.0)); + assertTrue("Polygon intersects line inside it", pred.test(line)); + + // Polygon-to-Ring + Ring ring = new Ring("geodetic", Arrays.asList(-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0)); + assertTrue("Polygon intersects ring inside it", pred.test(ring)); + + // Polygon-to-Polygon + Polygon polygon2 = new Polygon("geodetic", Arrays.asList(ring)); + assertTrue("Polygon intersects polygon inside it", pred.test(polygon2)); + } + + @Test + public void testOrdsInfoShapesDeserializationWorkflow() { + // Simulate ES plugin workflow: deserialize shapes from stored format + + // Point shape + List ordsInfo1 = Arrays.asList(4, 2); // type 4 = point, size 2 + List ords1 = Arrays.asList(0L, 0L); // lon=0, lat=0 + List shapes1 = OrdsInfoShapes.ordsInfoToShapes(ordsInfo1, ords1); + assertEquals(1, shapes1.size()); + assertTrue(shapes1.get(0) instanceof Point); + + // Geodetic polygon shape + List ordsInfo2 = Arrays.asList(1, 10); // type 1 = geodetic-polygon, size 10 + List ords2 = Arrays.asList( + -20000000L, -20000000L, // -2, -2 + 20000000L, -20000000L, // 2, -2 + 20000000L, 20000000L, // 2, 2 + -20000000L, 20000000L, // -2, 2 + -20000000L, -20000000L // -2, -2 (closing) + ); + List shapes2 = OrdsInfoShapes.ordsInfoToShapes(ordsInfo2, ords2); + assertEquals(1, shapes2.size()); + assertTrue(shapes2.get(0) instanceof Polygon); + + // Verify intersection works end-to-end + ShapePredicate pred = ShapeIntersections.createIntersectsFn(shapes1.get(0)); + assertTrue("Deserialized shapes can intersect", pred.test(shapes2.get(0))); + } +} diff --git a/spatial-lib/test/cmr/spatial/test/arc_segment_intersections.clj b/spatial-lib/test/cmr/spatial/test/arc_segment_intersections.clj index 317bd7efb9..5db85369a2 100644 --- a/spatial-lib/test/cmr/spatial/test/arc_segment_intersections.clj +++ b/spatial-lib/test/cmr/spatial/test/arc_segment_intersections.clj @@ -69,6 +69,13 @@ ;; line segment at south pole [10 -90 30 -90] [0 -90 10 -85] [0 -90] + ;; Geodetic pole-touching: disjoint longitude ranges but both touch north pole + ;; Tests that geodetic semantics correctly identify pole intersection + [0 90 90 85] [180 90 -170 85] [0 90] + + ;; Geodetic pole-touching: disjoint longitude ranges at south pole + [45 -90 135 -85] [-135 -90 -45 -85] [0 -90] + ;; line segment along antimeridian [180 10 180 20] [175 15 -175 15] [180.0 15.0547] diff --git a/spatial-lib/test/cmr/spatial/test/serialize.clj b/spatial-lib/test/cmr/spatial/test/serialize.clj index 43859e26ea..f5dba8aa69 100644 --- a/spatial-lib/test/cmr/spatial/test/serialize.clj +++ b/spatial-lib/test/cmr/spatial/test/serialize.clj @@ -75,3 +75,8 @@ parsed-shapes (srl/ords-info->shapes ords-info ords)] (= shapes parsed-shapes)))) +(deftest ords-info->shapes-nil-test + "Test that nil ords-info (granules without spatial data) returns empty list" + (is (= [] (srl/ords-info->shapes nil nil))) + (is (= [] (srl/ords-info->shapes nil [])))) +