[OTHER] HDR Imaging

[Bra1nsen]

The problem occurs because exposure and analogue gain changes take a few frames to happen (the sensor can't apply them instantly, and there is also some pipe-lining going on the software stack, so it's "a few frames").

The AGC algorithm unfortunately only remembers its current state, not the state it was in "a few frames ago". This means it may apply the wrong digital gain to the image. As I implied earlier, digital gain affects only the processed image, not the raw Bayer image which is why those would be OK.

As I also said earlier, there are patches under review for this in libcamera, I would hope they'll get merged in a week or two, but I can't give any guarantees.

REFERENCE ISSUE#399 (12th comment or so)

Any Updates here?

---

Moreover I saw a hdr function got implemented, how satisfied are you with the results? Maybe I can help further:

https://www.youtube.com/watch?v=rdib2wja9W8&t=2s&ab_channel=ColumbiaComputerVision

Cheers, Paul (RaspberrySkyCam)

[Bra1nsen]

It would be nice to take a really good HDR sky image or constellation. Maybe to promote Raspberry PiCam2. It would be impressive to get the maximum out of the Sony sensor.

The new Raspberry All Sky Camera prototype is really decent!

[davidplowman]

Hi, I've just posted a new release 0.3.7 which fixes the digital gain issue that we had in libcamera. So the code that I sent to change the exposure on every frame should give you correct processed images too now, as well as correct raw images.

One slight caveat is that we have a sensor/driver bug related to the imx477 which means that for this to work correctly all the time, you need to use a fixed framerate, but if I recall correctly, I think you may be doing this anyway.

But glad to hear you seem to have things working well!

[Bra1nsen]

Okay great. I will test it soon and give feedback!

In the meantime, one quick other question. Where can I find information about the Bayering process?

I would like to understand what formulas are used for the transformation of the raw pixel array into the r,g,b image file. Do you have a link for me?

[davidplowman]

De-Bayering (also "demosaicing", "demosaicking", "CFA interpolation") algorithms tend to get quite complex.

If you don't mind reducing the resolution by 2 in each direction, then the "use the red and blue, average the greens" is actually not terrible.

If you want to maintain the full resolution, that's where it gets tricky. The simplest algorithm is probably "bilinear interpolation" (for example, here) though the results are generally not very good.

After that it gets more difficult. Many algorithms try to interpolate along directions where there seems to be least change (i.e. interpolating along and not across edges). Some algorithms start by interpolating the missing greens everywhere and then interpolate R - G and B - G as the colours are often correlated and you get less bad aliasing. If you check the Wikipedia page it lists quite a few papers and algorithms... be warned that it's a deep hole!

The hardware on the Pi uses a relatively complex algorithm that is using a number of these techniques.

[Bra1nsen]

Yes my current dataset consists of processed r,g,b images (.tga files).

It would be great to know exactly which demosaicing algorithmic (math) Raspberry uses, it has to be somewhere?

I am currently reconstructing the camera response function from the exposure time series. (including channelwise) With the known response function, the algorithm can fuse the multiple photographs into a single, high dynamic range radiance map whose pixel values are proportional to the true radiance values in the scene. (debevec)

[davidplowman]

The demosaic algorithm on the Raspberry Pi is proprietary and owned by Broadcom, so I can't say too much about it. I would say it's at the more complex end of the various algorithms I've seen in use.

[Bra1nsen]

So as a Raspberry employee I would have access to the algorithms? Interesting.

Is that the formula metadata "lux" beeing calculated? (lux.cpp)

        double estimatedLux = shutterSpeedRatio * gainRatio *
                      apertureRatio * apertureRatio *
                      yRatio * referenceLux_;

Where does this formula come from?

---

Maybe it would gain benefit to fight a map() function for W/m², so the camera could measure that accurately and returning it. Kind Regards .. Paul

[Bra1nsen]

I imagine to capture images like that:

working wiht numpy and matplotlib

[davidplowman]

It's well-known that Raspberry Pis use Broadcom silicon. As part of that agreement, Raspberry Pi gets to use Broadcom's firmware, but this all happens under an NDA, so the details cannot be discussed with anyone else.

When the sensor is calibrated, the tools are given a scene of known lux and known exposure time and gain, which is stored in the tuning file. After that, this calculation is simply assuming that a scene where the pixels are about the same brightness, but the exposure time is halved, must be twice as bright (i.e. just assuming a linear relationship).

There are numerous techniques for generating HDR images and tone-mapping them. Some work in the raw domain but others can even use processed RGB images (for example, OpenCV's "Mertens Merge" function).

[Bra1nsen]

Wouldn't it be cool if the HQ Cam (IMX477) could measure solar radiation. I think the raw Bayer arrays approach (stored as .npy) is pretty good.

In principle, I want to achieve exactly the same thing as you did back then. Only I don't measure Lux, but W/m².

The Metadatalux was probably determined using a lux meter, wasn't it?

[davidplowman]

Yes, exactly, an instrument like one of those. The tuning tool has to be fed images where each one is labelled with the measured lux level and also the colour temperature of the illumination (in our case, the same device produces both numbers).

[cpixip]

Just a short comment: the raw image data of a camera is (mostly) correlated linearily with scene radiance. So having exact exposure times and raw image data, it is only a matter of calibration to "convert" a camera sensor to a measurement instrument.

Another (not based on raw data) approach is the Debevec approach (http://www.pauldebevec.com/Research/HDR/debevec-siggraph97.pdf) - this assumes that you are working not on the raw image data, but on the usual output format (say: jpg). In this case, you need to recover the gain curves of the processing pipeline (in our case: libcamera). Here's an example of some gain curves, left is for the Raspberry PI v1 camera, center for the See3CAM_CU135 camera, and right, for the HQ camera sensor. To obtain these gain curves, both Raspberry Pi cameras were operated under the old Broadcom-ISP (not: libcamera):

This estimate of the gain curve would have to be re-done with a different ISP, say, libcamera. Sadly, there are a bunch of automatic algorithms in the standard libcamera approach and it is not that easy to turn them off. So the results are yet less precise than they could be. There is a first attempt to approach this issue, have a look here (https://forums.raspberrypi.com/viewtopic.php?t=343449&sid=f99105e9ce7e5b383c3f96dc1eadd987). Basically, one needs to modify the tuning file which libcamera uses to get from raw image data to normal image output.

While the "Mertens Merge" approach creates pleasing low dynamic range images from a stack of differently exposed input images, it should not be considered a real HDR. There is no linear relationship between scene radiance and image data and also it is usually sufficient to store the result of the Mertens algorithm in a standard 8-bit/channel format. Contrary, a real HDR requires floating point data, and if viewed without further processing, looks rather dull. Only if a tone-mapping step is performed on the HDR (in effect converting it back to a LDR (low dynamic range) image, a normal, viewable image is obtained.

In effect, a Mertens merge short-cuts HDR-creation and tone-mapping, going directly from a stack of differently exposed LDR images to a nicely viewable output LDR.

[Bra1nsen]

@cpixip cpixip, Thanks for pariticipating. I linked the exact same Paper 4 hours ago (debevec) ;D If you like add me in LinkedIn.

As you can see, I already fused the images according to their "camera_response_curve" characteristics.

So the resultign hdr pixel values are proportional to the true radiance values in the scene.

BUT

the image array has 3 channels. Iam wondering how to convert the r,g,b values to a grayscale while keeping radiometric properties true.

Since the 3 colours are way different in terms of radiance intensity. Maybe using a weight function for grayscaling. Or just using single channels, which is my current approach. Examining channelwise.

[cpixip]

Sorry, overlooked your link to the paper. With respect on how to derive a meaningful number from three RGB-values - for me this was always a rabbit hole I tried to avoid (there's a nice overview table of a lot of the different quantities involved at the bottom of the German Wikipedia "Strahldichte" page https://de.wikipedia.org/wiki/Strahldichte).

Well, here's another overview paper (https://www.mikewoodconsulting.com/articles/Protocol%20Spring%202014%20-%20Photopic%20Curves.pdf) which might be useful. Technically, you map the response fct of the sensor (which you have displayed above from the IMX477 data sheet) to the desired response fct. you want your measurements to be done with. Some idea of how that can be done might be further found in this application note (section 2.5/Calibration with Matrix Based Algorithm https://ams.com/documents/20143/36005/AS7341_AN000633_2-00.pdf. But again, I would not suggest to go down that rabbit hole...

Some further comments: your above figure just shows the CFA-filter fcts. Be aware that in the Raspberry HQ camera, there is also an IR-blocking filter. So the real response curves of the HQ camera do look more like this:

Also, if you are using the Debevec-approach, you are always using the libcamera stack. This has advantages, like noticably smaller image size, faster frame rates and the like. It has however the disadvantage that you rely on libcamera - which features by default some automatic adjustment algorithms. That will change locally your transfer (gain) fct. Whether this is a concern for your application depends. If this is of no concern, you actually might even fare well with simply assuming that your RGB-values are sRGB-encoded (they are not sRGB within the context of libcamera) and calculate the Luminance value in Lab or YUV space by the usual formulas. An even simpler approach would be to use the lux-value reported in the metadata of libcamera. The above formula you quoted is simply describing the relationship between the various variables which affect the raw image values; once you include the appropriate calibration constant (I think that constant would be "referenceLux_"), you have a value which should be not too far away from the real thing. Advantage: the people which created the IMX477 tuning have already done this work for you. And, as the lux computation is one of the first things in the processing pipeline, this value is also not affected by any of the automatic algorithms implemented further down the libcamera stack.

[Bra1nsen]

Thanks.

Thanks for the tip about the ir filter. these are very interesting plots. I would like to stay in touch with you, perhaps via email? (paulmatteschk@skycam.info)

We will start with linear regression model between lux values (currently average from timeseries) and global horizontal irradiance (ghi)

basically a mapping function() from Lux to W/m²

Also I'm considering to integrate Environment- and Sensortemperature to the model.

@davidplowman Are there more metadata values which could be useful?

What does the noise..cpp? https://git.linuxtv.org/libcamera.git/tree/src/ipa/raspberrypi/controller/rpi/noise.cpp

[davidplowman]

Are there more metadata values which could be useful?

What does the noise..cpp? https://git.linuxtv.org/libcamera.git/tree/src/ipa/raspberrypi/controller/rpi/noise.cpp

It's hard to say exactly what you might find useful. These is quite a bit of documentation here (check out chapter 5) which would help you understand what's in the tuning file.

As regards those noise parameters, they would let you increase or decrease the amount of denoising applied to the processed images. They're described in section 5.4 of the documentation, though it would be more usual to change the "rpi.sdn" settings (section 5.6).

[davidplowman]

Hi, is there anything else to look at here, or can we consider closing this one? Thanks!

[Bra1nsen]

Is Raspberry interested in measuring Solar Irradiance?

my linear regression modell for estimating solar irradiance shows high accuracy... I would like to add that feature to the library with your help.

Cheers, Paul

[Bra1nsen]

When the sensor is calibrated, the tools are given a scene of known lux and known exposure time and gain, which is stored in the tuning file. After that, this calculation is simply assuming that a scene where the pixels are about the same brightness, but the exposure time is halved, must be twice as bright (i.e. just assuming a linear relationship).

Hey David, how are you? I'm at the documentation right now and I have to ask you again about the calculation of the "Lux" metadata.

Specifically, this parameter:

        double yRatio = currentY * (65536 / numBins) / referenceY_;

In order to output the lux value, the image brightness is taken into account, right? Is the image brightness of the raw bayer array used? And where does the 65536 come from, they look familiar to me, but I can't place them right now.

[Bra1nsen]

Furthermore, maybe the role of the lens plays an important role, where can I define the reference aperture?

For example, the image sensor is utilized more (in terms of area) with one lens and than with the other.

[davidplowman]

Hi

Specifically, this parameter:

      double yRatio = currentY * (65536 / numBins) / referenceY_;

In order to output the lux value, the image brightness is taken into account, right? Is the image brightness of the raw bayer array used? And where does the 65536 come from, they look familiar to me, but I can't place them right now.

Yes, the value currentY is the average value calculated from the image histogram, so it lies between zero and numBins. However, referenceY_ is read from the tuning file and is in the range 0 to 65536 (which is 2^16). So this really is just calculating the ratio of the current brightness to the reference image brightness, but currentY has to be scaled up to the same range as referenceY_, hence multiplying by (65536 / numBins).

As regards aperture, this is a problem because the system has no idea what the current aperture is, and there is absolutely no way for it to find out. So we normally put values in the tuning file for an "in-between" kind of aperture value, but in general, you would need to scale the final lux value by the current aperture... if you know what it is. I think currentAperture and referenceAperture_ are always just 1.0 because we have no good values to use instead!

[Bra1nsen]

The dynamic range is 12 bits, so why choose 2^16?

        config = camera.create_still_configuration(main={"size": (1014, 760), "format": "RGB888"},
                                                   raw={"format": "SRGGB12", "size": (2032, 1520)})

I would like to understand exactly how the mean of the histograms is determined. So in my example, the referenceY (Why symbol Y?)

2032 x 1520 pixels each pixel 2^16 mean average of all pixels right?

The aperture would primarily play a role in the radiance (area) of the imx477. maybe you can determine the experimental

[davidplowman]

We use 16 bits in the tuning file because we prefer to normalise everything in them so that the person generating or changing it doesn't have to worry about the bit depth of the sensor.

The symbol Y is widely used to refer to the luminance (or brightness) of pixels, which is what we're interested in here.

The reference Y is normally obtained by using the camera tuning tool. I've never calculated it in any other way, but you would take the mean of all the raw green pixels after subtracting the pixel black level and applying any lens shading correction, and then scale it up to 16 bits (so for 12-bit raw modes, that would be x16).

[Bra1nsen]

Is there a way to get referenceY as metadata?

For doing it manually:

picam2 = Picamera2()
config = picam2.create_still_configuration(raw={"format": "SRGGB12", "size" : (2032, 1520)}, buffer_count=2)
picam2.configure(config)

picam2.set_controls({"ExposureTime": 3000, "AnalogueGain": 1, "ColourGains": (1.0,1.0)})
picam2.start()
array = picam2.capture_array("raw")
picam2.stop()

bayer = array.view(np.uint16)

bggr = np.zeros((2*760, 2*1016), dtype=np.uint16)
bggr[1::2, ::2]=bayer[1::2, 0::2] #green
bggr[::2, 1::2]=bayer[0::2, 1::2] #green

referenceY = np.average(bggr)

How can I subtract all black levels pixels? & is lenshading correction really necessary?

[davidplowman]

Unfortunately that reference Y doesn't come out as metadata. For 12-bit raw data, the black level is 256 so that's what you need to subtract. I wasn't quite sure about the bggr array above, it looked like half the pixels (the B and R ones) are left as zero so you'd get half the average that you want? Does this work:

referenceY= (np.average(bayer[1::2, 0::2]) + np.average(bayer[0::2, 1::2])) / 2 - 256

Finally there's a scaling of x16 if you want a number to put into the tuning file.

[Bra1nsen]

I get negative values. Not sure how to validate/check if numbers are correct.

Basically I would like to get referenceY from every "raw" frame. The list images contains 9 raw image frames (exposure timeseries)

from picamera2 import Picamera2

images = [] #list elements (index, array) with array = (1520, 2032) for each raw bayer frame

referenceY = (np.average(images[0][1][1::2, 0::2]) + np.average(images[0][1][0::2, 1::2])) / 2 - 256

[Bra1nsen]

[Bra1nsen]

I shortend the code to focus on the essentiell problem. Could you please provide further recommendations for calculating referenceY from an raw bayer image frame? Moreover it would be great to get some kind of validation.

[Bra1nsen]

I still dont fully understand why we are dividing by 2 and subtracting 256 (black level), could please explain that further?

[cpixip]

Well, there are two green channels in the raw image file. Normally, you add them together. To get the values of this combined green channel back into the range of the red and blue channel (12 bit), you need to divide by two.

In pseudo-code: green = 0.5*(green1_raw + green2_raw).

The reason you are getting negative values after subtraction of the blacklevel is not unusual. Blacklevel values are only approximate values for the whole camera array, but single pixels vary. You just can ignore these - image intensities near the blacklevel of a sensor are basically noise, there is no usefull information hidden in such data.

Of course, you need to subtract the correct blacklevel for the sensor you are using. The blacklevel given in the tuning file of the IMX477 is 4096 - which translates into a blacklevel = 256 for the 12 bit image data you are working with. To make it even more explicit: a blacklevel of 4096 in reference to 16 bit data (that is the tuning file convention) corresponds to a real 256-value in reference to the 12 bit data your sensor is delivering and which you are processing.

[davidplowman]

@cpixip Thank you very much for replying to the question, I've been away for a few days! @Bra1nsen The explanation is all exactly right, please let us know if there's anything else!

[Bra1nsen]

Dark current can have a negative impact on image quality. Some image sensors can compensate for this by automatically adjusting the black level. Understood.

#CAMERA SENSOR SONY IMX477 SRGGB12
#FISHEYE LENS 185° ALL SKY IMAGES
#INPUT EXPOSURE SERIES/STACK WITH 9 FRAMES
#EXPOSURETIMES = exp_time = [60,106,166,273,470,773,1258,2032,3291] µs (1e-6s)

How to interpret the y-axis values from referenceY (average brightness greens pixels, correct?)

def referenceY(timestamp, images):
    try:
        referenceY0 = (np.average(images[0][1][1::2, 0::2]) + np.average(images[0][1][0::2, 1::2])) / 2 - 256
        referenceY1 = (np.average(images[1][1][1::2, 0::2]) + np.average(images[1][1][0::2, 1::2])) / 2 - 256
        referenceY2 = (np.average(images[2][1][1::2, 0::2]) + np.average(images[2][1][0::2, 1::2])) / 2 - 256
        referenceY3 = (np.average(images[3][1][1::2, 0::2]) + np.average(images[3][1][0::2, 1::2])) / 2 - 256
        referenceY4 = (np.average(images[4][1][1::2, 0::2]) + np.average(images[4][1][0::2, 1::2])) / 2 - 256
        referenceY5 = (np.average(images[5][1][1::2, 0::2]) + np.average(images[5][1][0::2, 1::2])) / 2 - 256
        referenceY6 = (np.average(images[6][1][1::2, 0::2]) + np.average(images[6][1][0::2, 1::2])) / 2 - 256
        referenceY7 = (np.average(images[7][1][1::2, 0::2]) + np.average(images[7][1][0::2, 1::2])) / 2 - 256
        referenceY8 = (np.average(images[8][1][1::2, 0::2]) + np.average(images[8][1][0::2, 1::2])) / 2 - 256

        refY = [referenceY0, referenceY1, referenceY2, referenceY3, referenceY4, referenceY5, referenceY6, referenceY7,
                referenceY8]
    except:
        pass

    return refY

[Bra1nsen]

interesting ref8 (longest exposure) mapst best, maybe I have to integrate longer exposure times for the raw images.

[cpixip]

The referenceY-values should not look like the plot above, assuming that they are created with the exposures you listed above. So check your code. Most important, check the exposure time and digital gain returned for the specifc capture. Most probably the exposure time used is not the exposure time you have been requesting.

Another point: some exposure times in the set of exposure times you have selected are rather short. It is a little challenging for the hardware. A neutral density filter could be used to shift the set if exposure times to a more suitable range.

[Bra1nsen]

Thanks for your help cpixip.

The script provided by davidplowman ensures that the exposure times are correct. Furthermore, the gain should be constant at 1.0. But thanks for the tip, I haven't dealt with gain very much until now.

picam2.set_controls({"ExposureTime": exp, "AnalogueGain": 1.0, "ColourGains": (1.0, 1.0), "FrameDurationLimits": (25000, 50000)})

Since I am calculating the luminance from the whole image array, the value may be washed out by the pixels not exposed by the lens.

I would try two things now: -Integrate longer exposure times -refY refer to the relevant pixels of the sky area

Any other ideas?

[Bra1nsen]

I will adjust the exposure time for better ranging. Iam also thinking about night sky/star images. Anyway Iam still confused by the y-axis values...

referenceY0 = (np.average(images[0][1][1::2, 0::2]) + np.average(images[0][1][0::2, 1::2])) / 2 - 256

[Bra1nsen]

[cgobat]

@cpixip could you share where you sourced the data for those transmission curves that account for the IR filter?

[cpixip]

If you mean these curves

These curves are calculated. For that, I measured the response of the IR-filter with a NIR-spectrometer and took the curves of the raw CFA directly from the IMX477 data sheet.

[cgobat]

Cool, thanks! Would you mind sharing the data of your IR filter response measurements?

[cpixip]

Well, here we go. This here

HQ Cam IR_Filter_Transmission_UV und VIS.csv

should be the measured IR-filter response. In fact, I think the Raspberry Pi Foundation published at one point in time the IR-filter response as well. So you could search for this and digitize it. As far as I remember, my measurements and the published curve mostly agree with each other.

Two further comments. This here

IMX477.json

contains the digitized filter responses of the CFA of the IMX477, according to the data sheet. I think they are pretty close to ground truth. From these two data sets you should be able to calculate the above plot.

Also, it is not published very widely, but the Raspberry Pi foundation changed or is going to change the IR-filter on the HQ cameras they sell. The new filter is only slightly different than the old one, here's

HQ Cam New IR.csv

this data set. Color science is only mildly affected; up to my knowledge, there is no way to figure out which IR-filter a specifc unit is using - older ones should of course feature the old "Hoya CM500", but newer units might already use the "CU700 UVIR" filter.

The CFA and IR filter curves are combined in my Jupyter notebook via the following code segment:

# fct for combining two filters
def applyFilter(multiSpectra,singleSpectra,name='filter applied'):
    '''
    folds a multi-spectral distribution with a 
    single spectral distribution
    '''
    newData = []
    for channels in multiSpectra.labels:
        newData.append(multiSpectra.signals[channels]*singleSpectra)
    return colour.MultiSpectralDistributions(newData,name=name,labels = multiSpectra.labels)

using the colour package referenced here.

[cgobat]

Very helpful—thank you for sharing!

[Bra1nsen]

@cpixip if no ir filter is used, what do you think, how far will the red channel be spectral sensitive ~ 1100nm?

here ist the csv, which ends by 700nm: https://github.com/bluegreen-labs/raspberry_pi_camera_responses

[Bra1nsen]

Do you also have absolut spectral responses from some camera sensors? Such data seems hard to find.

[cgobat]

@Bra1nsen note also that so far we were talking about the HQ camera (which has the IMX477 sensor). The link and plot you shared are for the IMX219, and the bluegreen-labs/raspberry_pi_camera_responses repo doesn't appear to have data for the HQ.

[cpixip]

if no ir filter is used, what do you think, how far will the red channel be spectral sensitive ~ 1100nm?

Probably not. The general behaviour of all these camera sensors is pretty standard: the CFA filter curves approach each other (that is, red sensitivity similar to green and blue sensitivity, thus no longer any "color vision") and the response dies out anyway already at 900 to 1000 nm. Your milage might vary, of course.

[Bra1nsen]

@cgobat Haha yes true, I know, thanks for mention it :)

@cpixip Thanks for the response, helped!

I took the these curves and plotte them within the solar spectrum:

It would be nice if the quantum efficiency curves of the hq cam (imx477) would also be available...

[cgobat]

The JSON file that @cpixip already linked has the IMX477 response curves, or here is a standardized NumPy binary format data product derived from that JSON.

[Bra1nsen]

As far as I understood, these are only the relative responses of color channels to each other...

Which is particularly useful in color imaging to ensure accurate color reproduction...

The Quantum Efficiency (QE), on the other hand, is a measure of the absolute efficiency of a sensor in converting incoming photons to electrons

[cpixip]

Yes, these are relative responses. Some people did mesure QE for the RP v1 and v2 cams using an integrating sphere and a monochromator. For the sensor's color science, QE is irrelevant. To get a usable QE, you would need to measure the full system, including the specific lens you are using.