

function isobj(filepath)
    base, ext = splitext(filepath)
    ext == ".obj" ? true : false
end


"""
    loadobj(shapepath::String; scale=1, message=true) -> nodes, faces
"""
function loadobj(shapepath::String; scale=1, static=true, message=true)

    nodes = Vector{Float64}[]
    faces = Vector{Int64}[]

    open(shapepath, "r") do f
        for line in eachline(f)
            line == "" && continue  # Skip blank lines
            data = split(line)

            if data[1] == "v"
                node = parse.(Float64, data[2:4])  # (x, y, z) coordinates [km]
                push!(nodes, node)
            elseif data[1] == "f"
                face = parse.(Int64, data[2:4])  # indices of three vertices
                push!(faces, face)
            end
        end
    end

    nodes *= scale  # if scale is 1000, converted [km] to [m]
    
    if static == true
        nodes = [SVector{3,Float64}(node) for node in nodes]
        faces = [SVector{3,Int64}(face) for face in faces]
    end
    
    if message == true
        println("+-----------------------------+")
        println("|        Load OBJ file        |")
        println("+-----------------------------+")
        println(" Nodes: ", length(nodes))
        println(" Faces: ", length(faces))
    end

    return nodes, faces
end


"""
    saveobj(path::String, nodes, faces)
"""
function saveobj(path::String, nodes, faces)
    open(path, "w") do f
        println(f, "# Generated by obj.jl written by M. Kanamaru")
        println(f, "# Vertices: ", length(nodes))
        println(f, "# Faces:    ", length(faces))
        # println(f)
        
        for node in nodes
            println(f, "v  ", node[1], "  ", node[2], "  ", node[3])
        end
        
        println(f)
        
        for face in faces
            println(f, "f  ", face[1], "  ", face[2], "  ", face[3])
        end
    end
end


"""
Make a perturbed shape model
"""
function perturbnodes!(nodes, noise_weight)

    for node in nodes
        noise = rand()*2 - 1  # random number in [-1, 1]
        noise *= noise_weight

        dr = normalize(node) * noise / 1000  # [km]
        node .+= dr
    end
    nodes
end

perturbnodes(nodes, noise_weight) = perturbnodes!(deepcopy(nodes), noise_weight)


# ****************************************************************

getFaceCenters(nodes, faces) = [getcenter(nodes[face]) for face in faces]
getFaceNormals(nodes, faces) = [getnormal(nodes[face]) for face in faces]
getFaceAreas(nodes, faces) = [getarea(nodes[face]) for face in faces]

