text:
To find islands of 1 and 0 in matrix
Given a two dimensional matrix of 0 and 1s. Find the number of island for 1s and 0s where neighbours are only in the horizontal and vertical.
Given input:

1 1 1 0
1 1 1 0

output = 1 1
Number of 1s island = 1

xxx-
xxx-

Number of 0s island = 1

---x
---x

------------------------------

Given input:

0 0 0 0
1 1 1 1
0 0 0 0
1 1 1 1

output = 2 2
Number of 1s island = 2

----
xxxx  <-- an island of 1s
----
xxxx  <-- another island of 1s

Number of 0s island = 2

xxxx  <-- an island
----
xxxx  <-- another island
----

------------------------------

Given input:

1 0 0
0 0 0
0 0 1
output = 2 1
Number for 1's island = 2:

x--  <-- an island of 1s
---
--x  <-- an island of 1s

Number of 0's island = 1:

-xx  \
xxx   > 1 big island of 0s
xx-  /

------------------------------

Given input:

1 1 0
1 0 0
output = 1 1
Number for 1's island =1 and number of 0's island = 1

------------------------------

Given input:

1 1
1 1
output = 1 0
Number for 1's island =1 and number of 0's island = 0

Can we assume the output is a 2-item list? e.g. 1 1 in either order?

yes but in order for 1's and then 0's island

You should add a testcase like [[1,0];[0,1]] to make sure diagonal connectivity is not included

I'd suggest the output can be in either order as long as the order is specified - it doesn't add any value to force an order

Welcome to the site!

What was answered in the comments should be clarified in the body of the challenge. And more specifically, if you really want us to return 1's before 0's, it should be clearly stated.

I'd like to see a PMA/Snails answer https://codegolf.stackexchange.com/questions/47311/language-design-2-d-pattern-matching/47495#47495

Suggested test case: 11111 / 10001 / 10101 / 10001 / 11111 → 2 1

Should we output only 2 numbers, or we output the ascii art?

@tsh It is clear fom the last 2 test cases that only the counts are needed.

J, 57 bytes
,&([:(0#@-.~~.@,)](*@[*[:>./((,-)#:i.3)|.!.0])^:_ i.@$)-.

Try it online!
This is one of those where the idea is incredibly simple (and I think fun), but executing it had some mechanical lengthiness which masks the simplicity... eg, shifting the original matrix in all directions with 0 fill is the verbose ((,-)#:i.3) |.!.0.
It's likely this mechanical lengthiness can be golfed further, and I may try tomorrow evening, but I'll post the crux of it now.
Say our input is:
0 0 0 0
1 1 1 1
0 0 0 0
1 1 1 1

We start with a matrix of unique integers the same size:
 0  1  2  3
 4  5  6  7
 8  9 10 11
12 13 14 15

Then for each cell we find the max of all its neighbors, and multiply by the input mask:
 0  0  0  0
 8  9 10 11
 0  0  0  0
13 14 15 15

We iterate this process until the matrix stops changing:
 0  0  0  0
11 11 11 11
 0  0  0  0
15 15 15 15

And then count the number of unique, non-zero elements.  That tells us the number of 1-islands.
We apply the same process to "1 minus the input" to get the number of 0-islands.

This looks very much like a "flood-fill" mechanism, really neat.

APL (Dyalog Unicode), 29 28 bytesSBCS
-1 thanks to @Adám
{≢∪∨.∧⍨⍣≡2>+/↑|∘.-⍨⍸⍵}¨⊂,~∘⊂

Try it online!
⊂,~∘⊂ the matrix and its negation
{ }¨ for each of them do
⍸⍵ list of pairs of coords of 1s
+/↑|∘.-⍨ matrix of manhattan distances
2> neighbour matrix
∨.∧⍨⍣≡ transitive closure
≢∪ number of unique rows

this is really clever.  could you elaborate on why the final line is guaranteed to work -- ie, why unique rows is equivalent to the answer.  also, is "transitive closure" is like J's ^:_?

@Jonah see chat

MATL, 14 12 bytes
,G@-K&1ZIugs

Try it online! Or verify all test cases.
Explanation
,        % Do twice
  G      %   Push input
  @      %   Push iteration index: first 0, then 1
  -      %   Subtract. This converts 0 and 1 into -1 and 0 in the second iteration
  K      %   Push 4
  &1ZI   %   Label connected components of matrix using 4-connectedness. Zeros in the
         %   matrix are background. This replaces the nonzeros by 1, 2, 3, ..., where
         %   each number defines a connected component
  u      %   Unique values. This gives [0; 1; 2; ..., L], where L is the number of
         %   connected components.
  g      %   Convert nonzeros to 1
  s      %   Sum. This gives L, to be output
         % End (implicit).
         % Display stack (implicit)

JavaScript (ES7),  138 ... 107  106 bytes
Returns an array [ones, zeros].
f=(m,X,Y,V=.5,c=[0,0])=>m.map((r,y)=>r.map((v,x)=>V-v|(x-X)**2+(y-Y)**2>1||f(m,x,y,v,r[c[v^1]++,x]=2)))&&c

Try it online!
How?
We use a recursive function. During the initial call, we look for \$0\$'s and \$1\$'s. Whenever we find such a starting point, we increment the corresponding island counter (\$c[0]\$ or \$c[1]\$) and enter the recursive phase to flood-fill the area of similar adjacent cells with \$2\$'s.
To save bytes, the exact same code is used for both the root iteration and the recursive iterations, but it behaves a bit differently.
During the first iteration:

We start with \$V=0.5\$ so that \$V-v\$ is rounded to \$0\$ for both \$v=0\$ and \$v=1\$, thanks to the bitwise OR.
\$X\$ and \$Y\$ are undefined and the computed quadrance \$(x-X)^2+(y-Y)^2\$ evaluates to NaN. It forces the inequality test to fail, no matter the value of \$(x,y)\$.

During the recursive iterations:

The parameter \$c\$ is set to an integer (\$2\$) instead of an array. Because numbers are objects in JS, the expression c[v ^ 1]++ is valid if \$c\$ is a number -- although it has no effect at all. It means that we can execute this statement unconditionally, without knowing if we're currently looking for starting points or flood-filling.

Commented
f = (                 // f is a recursive function taking:
  m,                  //   m[]  = input binary matrix
  X, Y,               //   X, Y = coordinates of the previous cell, initially undefined
  V = .5,             //   V    = value of the previous cell, initially set to 0.5
                      //          so that the integer part of V - v is 0 for v = 0 or 1
  c = [0, 0]          //   c[]  = array of counters of 1's and 0's islands
) =>                  //          (or an integer when called recursively)
  m.map((r, y) =>     // for each row r[] at position y in m[]:
    r.map((v, x) =>   //   for each value v at position x in r[]:
      V - v |         //     abort if |V - v| ≥ 1
      (x - X) ** 2 +  //     or X and Y are defined and the quadrance between
      (y - Y) ** 2    //     (X, Y) and (x, y)
      > 1 ||          //     is greater than 1
      f(              //     otherwise, do a recursive call to f:
        m,            //       leave m[] unchanged
        x, y,         //       pass the new coordinates
        v,            //       pass the new reference value
        r[c[v ^ 1]++, //       increment c[v ^ 1] (ineffective if c is an integer)
          x           //       and set the current cell ...
        ] = 2         //       ... to 2
      )               //     end of recursive call
    )                 //   end of inner map()
  ) && c              // end of outer map(); return c

This code does not work for large matrices like 100*100, with just 1s or 0s due to stack overflow.

@KBjoy Unless explicitly specified otherwise in the challenge, our default rule is that we don't care about implementation limits as long as the underlying algorithm work in theory for any input. (Here is a meta post about that, but there's probably a more relevant one somewhere.)

Wolfram Language (Mathematica), 64 62 bytes
Max@MorphologicalComponents[#,CornerNeighbors->1<0]&/@{#,1-#}&

Try it online!
Thanks to attinat: we can write 1<0 instead of False and save two bytes.
un-golfed version:
F[M_] := {Max[MorphologicalComponents[M,   CornerNeighbors -> False]],
          Max[MorphologicalComponents[1-M, CornerNeighbors -> False]]}

There is, of course, a Mathematica builtin MorphologicalComponents that takes an array (or an image) and returns the same with the pixels of each morphologically connected island replaced by the island index. Taking the Max of this result gives the number of islands (the background zeros are left at zero, and the island index starts at 1). We need to do this separately for the array (giving the number of 1-islands) and one minus the array (giving the number of 0-islands). To make sure diagonal neighbors don't count as neighbors, the option CornerNeighbors->False needs to be given.

-2 bytes since inequalities have higher precedence than Rule

K (ngn/k), 60 55 51 50 46 bytes
{#?{|/'x*\:x}/2>+/x*x:x-\:'x:(0,#*x)\&,/x}'~:\

Try it online!
~:\ a pair of the input and its negation (literally: negate iterate-converge)
{ }' for each
,/x flatten the arg
& where are the 1s? - list of indices
(0,#*x)\ divmod width(input) to get two separate lists for ys and xs
x-\:'x: per-axis distances ∆x and ∆y
x*x: square them
+/ add ∆x² and ∆y²
2> neighbour matrix
{|/'x*\:x}/ transitive closure
#? count unique rows

After seeing your answer I'm glad I didn't try to tackle this one in K :)

@streetster haha, thanks! that's not the effect i intended :) i'd actually like to encourage people to learn (any dialect of) k and golf in it

Python 3, 144 127 bytes
This solution uses cv2's awesome image processing power. Despite cv's less awesome, super long and readable method names, it beats both other Python answers!
Golfed:
import cv2,numpy as n
f=lambda b:n.amax(cv2.connectedComponents(b*255,0,4)[1])
def g(a):b=n.array(a,n.uint8);print(f(1-b),f(b))

Expanded:
import cv2
import numpy as np

# Finds the number of connected 1 regions
def get_components(binary_map):
    _, labels = cv2.connectedComponents(binary_map*255, connectivity=4) # default connectivity is 8
    # labels is a 2d array of the binary map but with 0, 1, 2, etc. marking the connected regions
    components = np.amax(labels)
    return components

# Takes a 2d array of 0s and 1s and returns the number of connected regions
def solve(array):
    binary_map = np.array(input_map, dtype=np.uint8)
    black_regions = get_components(1 - binary_map) # 0s
    white_regions = get_components(binary_map) # 1s
    return (black_regions, white_regions)

I'm not too familiar with Python, but why do you need the explicit argument names? Isn't just 4 instead of connectivity=4 and n.uint8 instead of dtype=n.uint8 possible?

@KevinCruijssen, you need the argument names if you skip optional arguments. Taking a peek at the docs, I actually don't have to skip, which saves me a lot of bytes. Thanks!

Ah ok, I thought it was something like that, but when I just looked at the docs I could only find a single cv2.connectedComponents method, so I was confused and thought there might have been a different reason for needing the argument names. As I said, I'm not too familiar with Python. All I've learned from it is from here on CCGC. ;) But it makes sense to use the variable names to skip other optional arguments.

Very nice! I found an online compiler which includes the cv2 module here.

J, 46 44 43 bytes
-1 byte thanks to @miles
,&#&~.&([:+./ .*~^:_:2>1#.[:|@-"1/~4$.$.)-.

Try it online!
tests and the ,& -. wrapper stolen from @jonah's answer
,& -. for the input and its negation do:
4$.$. (y,x) coordinates of the 1s as an n×2 matrix
1#.[:|@-"1/~ manhattan distances: abs(∆x)+abs(∆y)
2> neighbour matrix
[:+./ .*~^:_: transitive closure
#&~.&( ) number of unique rows

you can compose the length and unique to save another byte, ie ,&#&~. to avoid the cap [:

Jelly, 20 bytes
,CŒṪ€µạþ`§Ịæ*LṠQL«L)

Try it online!
Adapted from my Jelly answer to caird's shape grouping challenge. Unfortunately, handling empty lists had multiple consequences that made me use § instead of complex numbers and add «L. The outer structure was shamelessly stolen from caird's answer.
,CŒṪ€µạþ`§Ịæ*LṠQL«L)    Main link. Input = binary matrix
,CŒṪ€                   Classify coords of zeros and ones (Ġ doesn't work)
     µ             )    Operate on each list...
      ạþ`§Ị             Proximity matrix; m[i][j]==1 if two points are adjacent
           æ*L          m^(length of input list)
              ṠQL       Sign (convert positives to 1) and count unique rows
                 «L     Handle empty list to make it 0

Retina 0.8.2, 155 bytes
s`1(.*)
;$1a
}+`(?<=(.)*)(1|;)(.*¶(?<-1>.)*(?(1)$))?(?!\2)[1;]
;$3;
s`0(.*)
:$1b
}+`(?<=(.)*)(0|:)(.*¶(?<-1>.)*(?(1)$))?(?!\2)[0:]
:$3:
\W+(a*)(b*)
$.1 $.2

Try it online! Link includes test case. Explanation:
s`1(.*)
;$1a

If there's a 1, change it to ; and append an a to the end of the input so that it's out of the way.
}+`(?<=(.)*)(1|;)(.*¶(?<-1>.)*(?(1)$))?(?!\2)[1;]
;$3;

Flood fill any more adjacent 1s with ;s.
}

Repeat until all of the islands of 1s have been turned into ;s.
s`0(.*)
:$1b

If there's a 0, change it to : and append a b to the end of the input so that it's out of the way.
+`(?<=(.)*)(0|:)(.*¶(?<-1>.)*(?(1)$))?(?!\2)[0:]
:$3:

Flood fill any more adjacent 0s with :s.
}

Repeat until all of the islands of 0s have been turned into :s.
\W+(a*)(b*)
$.1 $.2

Separately count the numbers of islands of 1s and 0s.

Jelly, 44 42 bytes
ŒJfⱮ+€¥Ø.,UŻ¤œịḢ¥Ƈ⁹œịƇ€ɗⱮ,¬$fƇⱮ`ẎQ$€QƊÐL€Ẉ

Try it online!
A monadic link accepting a list of lists of integers as its argument and returning a list of the number of 1 and 0 islands in that order.
Thanks to @JonathanAllan for pointing out a bug in my code when there were islands that were diagonally adjoined.
Explanation (needs updating)
Step 1
Generate list of all matrix indices each with the indices of its neighbour to the right (unless on the right side) and down (unless at the bottom)
ŒJ            | Multi-dimensional indices (e.g. [1,1],[1,2],[1,3],[2,1],[2,2],[2,3])
      ¥       | Following as as a dyad:
  fⱮ          | - Filter the indices by each of:
    +€      ¤ |   - The indices added to the following
       Ø.     |     - 0,1
         ,U   |     - Paired with itself reversed [0,1],[1,0]
           Ż  |     - Prepended with zero 0,[0,1],[1,0]

Step 2
Split these indices by whether there was 1 or 0 in input. Returns one list of indices with neighbours for 1s and another for 0s.
  Ƈþ   | Filter each member of the output of stage 1 using the following criteria:
œị   $ | - Corresponding value for the multi-dimensional indices in each of the following as a monad:
   ,¬  |   - The input paired with its inverse

Step 3
Merge lists with members in common and output counts
           ƲÐL€  | For each of the outputs from stage 2, do the following as a monad and repeat until no changes
¹Ƈ               | - Filter out empty lists (only needed on first pass through but included here to save a byte)
  fƇⱮ`           | - Take each list of indices and filter the list of indices for those containing a match for any of them
        $€       | - For each resulting list of lists:
      Ẏ          |   - Tighten (concatenate top level of lists)
       Q         |   - Uniquify
          Q      | - Uniquify
               Ẉ | Finally output the lengths of the final lists

@JonathanAllan thanks. Should be fixed now, albeit at a cost of 6 bytes

Haskell, 228 227 225 224 bytes
import Data.List
z=zipWith
a!b=div(max(a*a)(a*b))a
l x=z(!)(z(!)x(0:x))$tail x++[0]
s=(\x->length.($x).filter<$>[(>0),(<0)]).nub.(>>=id).(until=<<((==)=<<))((.)>>=id$transpose.map l).z(\i->z(\j x->2^i*j*(2*x-1))[1,3..])[1..]

Try it online!
Explanation:
The idea for this solution is as follows: Initialise the matrix with unique values in each cell, positive for 1 and negative for 0. Then repeatedly compare each cell with its neighbours and, if the neighbour has the same sign but a number with a larger absolute value, replace the cell's number wit the neighbour's number. Once this hits a fixed point, count the number of distinct positive numbers for the number of 1 regions and the distinct negative numbers for the number of 0 regions.
In code:
s=(\x->length.($x).filter<$>[(>0),(<0)]).nub.(>>=id).(until=<<((==)=<<))((.)>>=id$transpose.map l).z(\i->z(\j x->2^i*j*(2*x-1))[1,3..])[1..]

can be separated into the preprocessing (assigning numbers to cells), the iteration, and the postprocessing (counting cells)
Preprocessing
The preprocessing part is the function
z(\i->z(\j x->2^i*j*(2*x-1))[1,3..])[1..]

Which uses z as abbreviation for zipWith to shave off a few bytes.
What we do here is zip the two dimensional array with integer indices at the rows and odd integer indices at the columns. We do this since we can build a unique integer from a pair of integers (i,j) using the formula (2^i)*(2j+1). If we generate only odd integers for j, we can skip calculating the 2*j+1, saving three bytes.
With the unique number, we now only have to multiply in a sign based on the value in the matrix, which is obtained as 2*x-1
Iteration
The iteration is done by
(until=<<((==)=<<))((.)>>=id$transpose.map l)

Since the input is in the form of a list of lists, we perform the neighbour comparison on each row, transpose the matrix, perform the comparison on each row again (which due to the transpose is what was the columns before) and transpose again. The code that accomplishes one of these steps is
((.)>>=id$transpose.map l)
where l is the comparison function (detailed below) and transpose.map l performs one half of the comparison and transposition steps. (.)>>=id performs its argument twice, being the pointfree form of \f -> f.f and by one byte shorter in this case due to operator precedence rules.
l is defined in the row above as l x=z(!)(z(!)x(0:x))$tail x++[0]. This code performs a comparison operator (!) (see below) on every cell with first its left neighbour, and then with its right neighbour, by zipping the list x with the right shifted list 0:x and the left shifted list tail x++[0] in turn. We use zeroes to pad the shifted lists, since they can never occur in the preprocessed matrix.
a!b is defined in the row above this as a!b=div(max(a*a)(a*b))a. What we want to do here is the following case distinction:

If sgn(a) = -sgn(b), we have two opposing areas in the matrix and do not wish to unify them, so a stays unchanged
If sgn(b) = 0, we have the corner case where b is padding and therefore a stays unchanged
If sgn(a) = sgn(b), we wish to unify the two areas and take the one with the bigger absolute value (for convenience's sake).

Note that sgn(a) can never be 0. We accomplish this with the formula given. If the signs of a and b differ, a*b is less or equal to zero, while a*a is always greater than zero, so we pick it as the maximum and divide with a to get back a. Otherwise, max(a*a)(a*b) is abs(a)*max(abs(a),(abs(b)), and by dividing this by a, we get sgn(a)*max(abs(a),abs(b)), which is the number with the larger absolute value.
To iterate the function ((.)>>=id$transpose.map l) until it reached a fixed point, we use (until=<<((==)=<<)), which is taken from this stackoverflow answer.
Postprocessing
For postprocessing, we use the part
(\x->length.($x).filter<$>[(>0),(<0)]).nub.(>>=id)

which is just a collection of steps.
(>>=id) squashes the list of lists into a single list,
nub gets rid of doubles,
(\x->length.($x).filter<$>[(>0),(<0)]) partitions the list into a pair of lists, one for positive and one for negative numbers, and calculates their lengths.

Jelly, 23 bytes
,CµŒṪŒ!ŒṖ€ẎạƝ€§€ỊȦƲƇẈṂ)

Try it online!
Somewhat of a brute force approach, times out if there are more than \$7\$ \$0\$s or \$7\$ \$1\$s
How it works
,CµŒṪŒ!ŒṖ€ẎạƝ€§€ỊȦƲƇẈṂ) - Main link. Takes a matrix B on the left
 C                      - Complement of each element
,                       - Pair; [B, 1-B]
  µ                   ) - Over each matrix M:
   ŒṪ                   -   Get the [x, y] coords of each element of 1 in M
     Œ!                 -   All permutations
         €              -   Over each permutation:
       ŒṖ               -     Get all partitions
          Ẏ             -   Tighten into a list of partitions
                  ƲƇ    -   Keep the partitions P for which the following is true:
             €          -     Over each partition:
            Ɲ           -       Over each overlapping pair:
           ạ            -         Absolute differences
               €        -     Over each difference:
              §         -       Get the sums
                Ị       -     Insignificant; Equal to 1 or 0
                 Ȧ      -     There are no zeros anywhere here?
                     Ẉ  -   Lengths of each
                      Ṃ -   Minimum

Java 10, 359 355 281 280 261 246 240 236 bytes
int[][]M;m->{int c[]={0,0},j=m[0].length,i=m.length*j,t;for(M=m;i-->0;c[~t&1]+=t<2?f(t,i/j,i%j):0)t=M[i/j][i%j];return c;}int f(int v,int x,int y){try{if(M[x][y]==v){M[x][y]|=2;f(v,x,y-f(v,x-f(v,x,y+f(v,x+1,y)),y));}}finally{return 1;}}

-74 bytes thanks to @NahuelFouilleul.
-10 bytes thanks to @ceilingcat.
Try it online.
Explanation:
int[][]M;              // Integer-matrix on class-level, uninitialized

m->{                   // Method with integer-matrix parameter and integer-array return-type
  int c[]={0,0}        //  Counters for the islands of 1s/0s, starting both at 0
      j=m[0].length,   //  Amount of columns
      i=m.length*j,    //  Index of the cells
      t;               //  Temp-value to decrease the byte-count
  for(M=m;             //  Set the class-level matrix to the input-matrix
      i-->0            //  Loop over the rows
      ;                //    After every iteration:
       c[~t&1]+=       //     Increase the corresponding counter by:
        t<2?           //      If the cell's value `t` is 0 or 1:
         f(t,i/j,i%j)  //       Call the recursive flood-fill method with value `t`
                       //       And increase the counter by 1
        :              //      Else:
         0;            //       Keep the counter the same by increasing by 0
    t=M[i/j][i%j];     //     Set `t` to the value of the current cell
  return c;}           //  After the nested loops: return the counters-array as result

// Recursive method with value and cell-coordinate as parameters,
// This method will flood-fill the matrix, where 0 becomes 2 and 1 becomes 3
int f(int v,int x,int y){
  try{if(M[x][y]==v){  //   If the cell contains the given value:
    M[x][y]|=2;        //    Fill the cell with 0→2 or 1→3 depending on the value
       f(v,x+1,y)      //    Do a recursive call downwards
      f(v,x,y+...)     //    Do a recursive call towards the right
     f(v,x-...,y)      //    Do a recursive call upwards
    f(v,x,y-...);}     //    Do a recursive call towards the left
                       //    (method `f` always result in 1, so we use it to our advantage to
                       //     save bytes here by nesting, instead of having four loose calls)
  }finally{return 1;}} //  Ignore any ArrayIndexOutOfBoundsExceptions with a finally-return,
                       //  which is shorter than manual checks
                       //  And return 1 to increase the counter

-74 bytes, removing clone and using |=2 : 0 -> 2 and 1 -> 3, however >0 was changed to ==1

sorry i had to remove the tests so that the tio link fit in comments

@NahuelFouilleul Thanks, smart using |=2! And I could still use <2 instead of ==1 for -1 byte by first checking for 0 (and thus they are changed to 2, and then using the <2 to check for 1 (which are changed to 3).

Python 3, 167 bytes
def f(m):
 n=[0,0];i=-2
 for r in m:
  j=0;i+=1
  for c in r:n[c^1]+=1-((i>=0)*(m[i][j]==c)*(1+({*r[:j]}=={c})*({*m[i][:j]}=={c^1}))or(j>0)*(r[j-1]==c));j+=1
 print(n)

Try it online!

Python 2, 168 bytes
def f(m):
 n=[0,0];i=-2
 for r in m:
	j=0;i+=1
	for c in r:n[c^1]+=1-((i>=0)*(m[i][j]==c)*(1+(set(r[:j])=={c})*(set(m[i][:j])=={c^1}))or(j>0)*(r[j-1]==c));j+=1
 print n

Try it online!
-2 bytes thanks to Kevin Cruijssen
+2 bytes formatting fix
Explanation
A counter is kept for 0s and 1s.
For each entry in the matrix, the following actions are performed:

Incrase counter for current value by 1
If the same value exists directly above or to the left, decrease by 1

This results in a false positive for left-aligned cases like
0 0 1
1 1 1

or
0 1
1 1

If such a situation arises, the counter is decreased by 1.
Return value is [#1, #0]

I'm afraid the OP mentioned in the second comment the order should be [#1, #0]. Bit pointless imo to enforce this, but it is what it is for now. Anyway, you can golf the {not c} to {c^1}, and fix the issue I mentioned by changing n[c]+= to n[c^1]+= in a similar matter. Nice answer though, +1 from me. :)

Ah, you're right. Thanks!

Perl 5 (-0777p), 110 bytes
May be improved, uses a regex to replace 1 with 3, then 0 with 2.
/
/;$m="(.{@-})?";sub f{($a,$b,$c)=@_;1while s/$b$m\K$a|$a(?=$m$b)/$b/s||s/$a/$b/&&++$c;$c}$_=f(1,3).$".f(0,2)

TIO

T-SQL 2008, 178 bytes
Input is a table variable.

x and y are the coordinates
v is the values 0 and 1(could also handle other numeric values)

The test data used in this example:
100
000
001

DECLARE @ table(x int, y int, v int)

INSERT @ values
(1,1,1),(1,2,0),(1,3,0),
(2,1,0),(2,2,0),(2,3,0),
(3,1,0),(3,2,0),(3,3,1)
SELECT*,y-x*99r INTO # FROM @
WHILE @@rowcount>0UPDATE #
SET r=b.r
FROM #,# b
WHERE abs(#.x-b.x)+abs(#.y-b.y)=1and #.v=b.v and #.r>b.r
SELECT v,count(distinct r)FROM #
GROUP BY v

Try it online

R, 194 172 bytes
function(m,u=!1:2){for(i in 1:2){w=which(m==i-1,T)
N=1:nrow(w)
A=!!N
for(s in N){u[i]=u[i]+A[s]
while(any(s)){A[s]=F
s=c(N[as.matrix(dist(w))[s[1],]==1&A],s[-1])}}}
rev(u)}

Try it online!
Perform a depth-first search starting in each cell of the matrix that is equal to 1 (or zero).

-2 bytes thanks to @Giuseppe