A spiral sequence

JavaScript (ES6),  267 .. 206  199 bytes

Returns an array containing the \$N\$ first terms of the sequence.

n=>(F=v=>++i<n?F([...v,(A=N[i]=[1,!j++||d+1,j%L?d:(j%=L*6)?++d:L++&&d++].map(k=>N[k=i-k].push(i)&&k),g=k=>v[E='every']((V,x)=>V-k|N[x][E](y=>A[E](z=>v[y]-v[z])))?k:g(-~k))()]):v)([L=1],N=[[i=j=d=0]])

Try it online!

How?

Definitions

By convention, we will call corner-cell a cell that has only one edge in common with the previous layer of the spiral and side-cell a cell that has two edges in common with the previous layer. As suggested by Ourous, we can also think of them as order-1 cells and order-2 cells, respectively.

Corner-cells are shown in yellow below. All other cells are side-cells (except the center cell which is a special case).

About cell neighbors

We don't really need to keep track of the coordinates of the cells on the grid. The only thing that we need to know is the list of neighbor cells for each cell in the spiral at any given time.

In the following diagrams, neighbors in the previous layer are shown in light shade and additional neighbors in the current layer in darker shade.

A cell has 2 neighbors among the previous cells if:

• it's the first side-cell of a new layer (like 8)
• or it's a corner-cell, but not the last one of the layer (like 9)

A cell has 3 neighbors among the previous cells if:

• it's a side-cell, but not the first one of the layer (like 10)
• or it's the last corner-cell of the current layer (like 19)

Implementation of cell neighbors

The index of the current cell (starting at \$1\$) is stored in \$i\$. The list of neighbors of a cell \$n\$ is stored in \$A[n]\$.

For obscure golfing reasons, we first compute the opposite offsets of the neighbors. For instance, \$1\$ means \$-1\$, which means the previous cell.

[                    //
  1,                 // the previous cell is always a neighbor of the current cell
  !j++ || d + 1,     // if this is not the first cell of the layer, the cell at -(d + 1)
                     // is a neighbor (otherwise, we insert 1 twice; doing it that way
                     // saves bytes and having duplicate neighbors is not a problem)
  j % L ?            // if this is a side-cell:
    d                //   the cell at -d is a neighbor
  :                  // else (corner-cell):
    (j %= L * 6) ?   //   if this is not the last cell:
      ++d            //     insert the dummy duplicate neighbor at -(d + 1); increment d
    :                //   else (last cell):
      L++ && d++     //     the cell at -d is a neighbor; increment L; increment d
]                    //

In the above code:

• \$L\$ is the index of the current layer, starting at \$1\$ (and not counting the center cell).
• \$j\$ is the index of the current cell within the current layer, going from \$1\$ to \$6\times L\$.
• \$d\$ is the distance of the current cell to those of the previous layer. It is incremented each time we go through a corner cell.

We then process a map() loop which converts the offsets \$k\$ into indices (\$i-k\$) and push the current cell as a new neighbor for all neighbor cells, to enforce neighborhood symmetry.

.map(k =>
  N[k = i - k].push(i) && k
)

Finding the next term of the sequence

Now that we have an up-to-date list of neighbors for all cells, we can finally compute the next term \$k\$ of the sequence, using another recursive helper function.

The value of a cell \$n\$ is stored in \$v[n]\$.

( g =                 // g = recursive function taking
  k =>                // the candidate value k
    v.every((V, x) => // for each previous cell of value V at position x, make sure that:
      V - k           //   V is not equal to k
      |               //   OR
      N[x].every(y => //   for each neighbor y of x:
        A.every(z =>  //     for each neighbor z of the current cell:
          v[y] - v[z] //       the value of y is not equal to the value of z
        )             //     end
      )               //   end
    )                 // end
    ?                 // if the above conditions are fulfilled:
      k               //   stop recursion and return k
    :                 // else:
      g(-~k)          //   try again with k + 1
)()                   // initial call to g with k undefined (this will cause V - k to be
                      // evaluated as NaN and force the 1st iteration to fail)

Clean, 284 279 272 262 bytes

import StdEnv
l=[0,-1,-1,0,1,1]
c(u,v)(p,q)=(u-p)^2+(v-q)^2<2||(u-p)*(q-v)==1
$[h:t]m=hd[[e: $t[(h,e):m]]\\e<-[1..]|and[e<>j\\(u,v)<-m|c h u,(p,q)<-m|q==v,(i,j)<-m|c p i]]

$(scan(\(a,b)(u,v)=(a-u,b-v))(0,0)[(i,j)\\n<-[1..],i<-[1,1:l]&j<-l,_<-[max(~j<<i)1..n]])[]

Try it online!

Generates the sequence forever.

Hexagon Mapping

Most of the code goes into mapping hexagons uniquely to (x,y) coordinates so that there's a single, simple function to determine adjacency which holds for all point mappings.

The mapped points look like this:

              ---
        --- < 2,-2> ---       x-axis ___.X'
  --- < 1,-2> === < 2,-1> ---  /__.X'
< 0,-2> === < 1,-1> === < 2, 0>'
  === < 0,-1> === < 1, 0> ===
<-1,-1> === < 0, 0> === < 1, 1>
  === <-1, 0> === < 0, 1> ===
<-2, 0> === <-1, 1> === < 0, 2>.__
  --- <-2, 1> === <-1, 2> ---  \  'Y.___
        --- <-2, 2> ---       y-axis    'Y.
              ---

From there, determining adjacency is trivial, and occurs when one of:

• x1 == x2 and abs(y1-y2) == 1
• y1 == y2 and abs(x1-x2) == 1
• y1 == y2 - 1 and x2 == x1 - 1
• y1 == y2 + 1 and x2 == x1 + 1
• x1 == x2 and y1 == y2

Point Generation

Notice that when traversing the hexagon in a spiral the differences recur for each layer n:

1. n steps of (1,0)
2. n-1 steps of (1,-1)
3. n steps of (0,-1)
4. n steps of (-1,0)
5. n steps of (-1,1)
6. n steps of (0,1)

This generates the points in the right order by taking sums of prefixes of this sequence:

scan(\(a,b)(u,v)=(a-u,b-v))(0,0)[(i,j)\\n<-[1..],i<-[1,1:l]&j<-l,_<-[max(~j<<i)1..n]]

Bringing it Together

The code that actually finds the sequence from the question is just:

$[h:t]m=hd[[e: $t[(h,e):m]]\\e<-[1..]|and[e<>j\\(u,v)<-m|c h u,(p,q)<-m|q==v,(i,j)<-m|c p i]]

Which in turn is mostly filtering by and[r<>j\\(u,v)<-m|c h u,(p,q)<-m|q==v,(i,j)<-m|c p i]

This filter takes points from m (the list of already-mapped points) by:

• Ignoring natural numbers that are equal to any j
• For every (i,j) where i is adjacent to p
• For every (p,q) where the value q is equal to v
• For every (u,v) where u is adjacent to the current point