I guess you can solve it easily by taking multiple samples. E.g. divide your texel into a 4x4 grid, then do lens correction and one texture fetch for each subpixel sample individually, take average result of all samples.
You could also use N random subpixel positions, where N is larger in areas of more distortion to optimize.

It depends on the lens correction distortion how much bicubic filter can help. If one pixel maps to an area of about 5x5 pixels in the photo, a cubic filter covering 3x3 pixels will still under sample. And you want a average of an area, not a better point sample, so that's why the better filter might not be good enough. To improve this, mip maps would help ofc as they make an average of an area.

I'll draw a picture to illustrate how multiple samples should fix it:

On the left is a destination pixel with 6 random samples. You map each of them with your distortion fixing math to the source image, which might be a larger area shown on the right.
A bilinear lookup should suffice for each sample, and the average of all samples gives a good estimate if sample count is high enough. Simple Monte Carlo integration.

So all you need is a way to generate random sub pixel positions, usually done using hash functions. Pseudo code would be like this:

const int sampleCount = 8;
const int dimensions = 2;
int seed = (int(curPixel.x) * screenWidth + int(curPixel.y)) * sampleCount * dimensions; // we ensure each sub sample gets it's own unique seed
vec3 sum (0);
for (int i=0; i<sampleCount; i++)
{
	float subOffsetU = randF (seed + i * dimensions + 0) - .5f; // assuming hash returns values between 0 and 1
	float subOffsetV = randF (seed + i * dimensions + 1) - .5f; // assuming hash returns values between 0 and 1
	vec2 subSampleCoords = pixelUV + vec2(subOffsetU, subOffsetV);
	vec3 sample = TextureFetch (ProjectionUnDistort(subSampleCoords));
	sum += sample;
}
vec3 averagedResult = sum / sampleCount;

Some example c++ hash function i'm using would be this:

	inline uint32_t randI6 (uint32_t v) // pcg
	{
		uint32_t state = v * 747796405u + 2891336453u;
		uint32_t word = ((state >> ((state >> 28u) + 4u)) ^ state) * 277803737u;
		return (word >> 22u) ^ word;
	}
	
	inline float randF (uint32_t time)
	{
		uint32_t s = randI6(time);
		return float(s) / float(0x100000000LL);
	}

I forgot to mention that such hash functions are no perfect random number generators. It can happen, e.g. if image width has a certain value, that the generated random samples show patterns, and then the method no longer works properly. Without visualizing sub sample positions that's hard to detect - quite an annoying problem. That's why you can find hundreds of different hash functions in shadertoy, but none is perfect. To fix it, we can try adding some constant arbitrary numbers like so:

int seed = (int(curPixel.x) * (screenWidth + 111) + int(curPixel.y)) * sampleCount * (dimensions + 1);

AhmedSaleh said:
How would I use the above approach without random samples ?

A regular grid, like subdividing our square into a NxN grid for N^2 samples. Often works better anyway.

AhmedSaleh said:
How would I use the c++ hash functions inside the shader ?

They should just work after adopting syntax issues. Otherwise shadertoy has many examples.

1. Is there a need for the mask that I created ?

no.

AhmedSaleh said:
There is no need to do the filtering over the whole image, that's why I created the mask..

I would make it work first, and after quality is fine you could use the mask for an adaptive sample count for better performance.

AhmedSaleh said:
How about this solution that I did ?

I see it's sampling a regular grid per texel, so yes this should give you high quality multi sampling.

But i miss the projection mapping in your inner loop?
I assumed you have some fancy analytical mapping from planar to fisheye projection, and you would need to do this for each sample before the texture lookup.