Search Results: "Dima Kogan"

28 June 2022

Dima Kogan: vnlog 1.33 released

This is a minor release to the vnlog toolkit that adds a few convenience options to the vnl-filter tool. The new options are

vnl-filter -l
Prints out the existing columns, and exits. I've been low-level wanting this for years, but never acutely-enough to actually write it. Today I finally did it.

vnl-filter --sub-abs
Defines an absolute-value abs() function in the default awk mode. I've been low-level wanting this for years as well. Previously I'd use --perl just to get abs(), or I'd explicitly define it: = sub 'abs(x) return x>0?x:-x; '=. Typing all that out was becoming tiresome, and now I don't need to anymore.

vnl-filter --begin ... and vnl-filter --end ...
Theses add BEGIN and END clauses. They're useful to, for instance, use a perl module in BEGIN, or to print out some final output in END. Previously you'd add these inside the --eval block, but that was awkward because BEGIN and END would then appear inside the while(<>) loop. And there was no clear was to do it in the normal -p mode (no --eval). Clearly these are all minor, since the toolkit is now mature. It does everything I want it to, that doesn't require lots of work to implement. The big missing features that I want would patch the underlying GNU coreutils instead of vnlog:
  • The sort tool can select different sorting modes, but join works only with alphanumeric sorting. join should have similarly selectable sorting modes. In the vnlog wrappe I can currently do something like vnl-join --vnl-sort n. This would pre-sort the input alphanumerically, and then post-sort it numerically. That is slow for big datasets. If join could handle numerically-sorted data directly, neither the pre- or post-sorts would be needed
  • When joining on a numerical field, join should be able to do some sort of interpolation when given fields that don't match exactly.
Both of these probably wouldn't take a ton of work to implement, and I'll look into it someday.

16 June 2022

Dima Kogan: Ricoh GR IIIx 802.11 reverse engineering

I just got a fancy new camera: Ricoh GR IIIx. It's pretty great, and I strongly recommend it to anyone that wants a truly pocketable camera with fantastic image quality and full manual controls. One annoyance is the connectivity. It does have both Bluetooth and 802.11, but the only official method of using them is some dinky closed phone app. This is silly. I just did some reverse-engineering, and I now have a functional shell script to download the last few images via 802.11. This is more convenient than plugging in a wire or pulling out the memory card. Fortunately, Ricoh didn't bend over backwards to make the reversing difficult, so to figure it out I didn't even need to download the phone app, and sniff the traffic. When you turn on the 802.11 on the camera, it says stuff about essid and password, so clearly the camera runs its own access point. Not ideal, but it's good-enough. I connected, and ran nmap to find hosts and open ports: only port 80 on is open. Pointing curl at it yields some error, so I need to figure out the valid endpoints. I downloaded the firmware binary, and tried to figure out what's in it:
dima@shorty:/tmp$ binwalk fwdc243b.bin
3036150       0x2E53F6        Cisco IOS microcode, for "8"
3164652       0x3049EC        Certificate in DER format (x509 v3), header length: 4, sequence length: 5412
5472143       0x537F8F        Copyright string: "Copyright ("
6128763       0x5D847B        PARity archive data - file number 90
10711634      0xA37252        gzip compressed data, maximum compression, from Unix, last modified: 2022-02-15 05:47:23
13959724      0xD5022C        MySQL ISAM compressed data file Version 11
24829873      0x17ADFB1       MySQL MISAM compressed data file Version 4
24917663      0x17C369F       MySQL MISAM compressed data file Version 4
24918526      0x17C39FE       MySQL MISAM compressed data file Version 4
24921612      0x17C460C       MySQL MISAM compressed data file Version 4
24948153      0x17CADB9       MySQL MISAM compressed data file Version 4
25221672      0x180DA28       MySQL MISAM compressed data file Version 4
25784158      0x1896F5E       Cisco IOS microcode, for "\"
26173589      0x18F6095       MySQL MISAM compressed data file Version 4
28297588      0x1AFC974       MySQL ISAM compressed data file Version 6
28988307      0x1BA5393       MySQL ISAM compressed data file Version 3
28990184      0x1BA5AE8       MySQL MISAM index file Version 3
29118867      0x1BC5193       MySQL MISAM index file Version 3
29449193      0x1C15BE9       JPEG image data, JFIF standard 1.01
29522133      0x1C278D5       JPEG image data, JFIF standard 1.08
29522412      0x1C279EC       Copyright string: "Copyright ("
29632931      0x1C429A3       JPEG image data, JFIF standard 1.01
29724094      0x1C58DBE       JPEG image data, JFIF standard 1.01
The gzip chunk looks like what I want:
dima@shorty:/tmp$ tail -c+10711635 fwdc243b.bin> /tmp/tst.gz
dima@shorty:/tmp$ < /tmp/tst.gz gunzip   file -
/dev/stdin: ASCII cpio archive (SVR4 with no CRC)
dima@shorty:/tmp$ < /tmp/tst.gz gunzip > tst.cpio
OK, we have some .cpio thing. It's plain-text. I grep around it in, looking for GET and POST and such, and I see various URI-looking things at /v1/..... Grepping for that I see
dima@shorty:/tmp$ strings tst.cpio   grep /v1/
GET /v1/debug/revisions
GET /v1/ping
GET /v1/photos
GET /v1/props
PUT /v1/params/device
PUT /v1/params/lens
PUT /v1/params/camera
GET /v1/liveview
GET /v1/transfers
POST /v1/device/finish
POST /v1/device/wlan/finish
POST /v1/lens/focus
POST /v1/camera/shoot
POST /v1/camera/shoot/compose
POST /v1/camera/shoot/cancel
GET /v1/photos/ / 
GET /v1/photos/ / /info
PUT /v1/photos/ / /transfer
/v1/changes message received.
/v1/changes issue event.
/v1/changes new websocket connection.
/v1/changes websocket connection closed. reason( )
/v1/transfers, transferState( ), afterIndex( ), limit( )
Jackpot. I pointed curl at most of these, and they do interesting things. Generally they all spit out JSON. /v1/liveview sends out a sequence of JPEG images. The thing I care about is /v1/photos/DIRECTORY/FILE and /v1/photos/DIRECTORY/FILE/info. The result is a script I just wrote to connect to the camera, download N images, and connect back to the original access point: Kinda crude, but works for now. I'll improve it with time. After I did this I found an old thread from 2015 where somebody was using an apparently-compatible camera, and wrote a fancier tool:

29 November 2021

Dima Kogan: GL_image_display

I just spent an unspeakable number of days typing to produce something that sounds very un-impressive: an FLTK widget that can display an image. The docs and code live here. The big difference from the usual image-drawing widget is that this one uses OpenGL internally, so after the initial image load, the common operations (drawing, redrawing, panning and zooming) are very fast. I have high-resolution images in my projects, and this will make my tools much nicer. Three separate interfaces are available: The FLTK widgets have built-in interactive panning/zooming, and the library can draw line overlays. So nice applications can be built quickly. I already added some early disabled-by-default support into the mrcal-stereo tool to visualize the rectification and report sensitivities:

8 November 2021

Dima Kogan: mrcal 2.0: triangulation and stereo

mrcal is my big toolkit for geometric computer vision: making models (camera calibration) and using models (mapping, ranging, etc). Since the release of mrcal 1.0 back in February I've been busy using the tools in the field, fixing things and improving things. Today I'm happy to finally be able to announce the release of mrcal 2.0. A big part of this release is maintenance and cleanup that resulted from me heavily using the tools over the course of this past year, and improving whatever was bugging me. The most notable result of that effort, is that splined models are no longer "experimental". They work well and they're awesome. Go try them. And there're a number of new features, most notably nice dense stereo support and nice sparse triangulation support (with uncertainty propagation!) These are awesome. Go try them. As before, the tour of mrcal provides a good overview of the capabilities of the toolkit, and is a good place to start reading the documentation. Reading these docs would be very illuminating for anybody that calibrates cameras, even for those that have no intent to actually use the mrcal tools. Let me know if you try it out! The most list of most notable improvements, from the release notes:

13 March 2021

Dima Kogan: Making the Supernova E3 tail light brighter

I got a dynamo hub for my bike and a fancy headlight. It's sweet, but I'm discovering that there're no standards for making tail lights work, so I just had to do some light reverse-engineering and soldering. And this is the findings. Different manufacturers do tail lights differently. Most tail lights are not connected directly to the dynamo, but to the headlight instead. Busch+M ller tail lights take an AC signal that looks very similar to what the hub is producing. You're not supposed to hook it up to the hub directly, but it does appear to work, and it's not clear how the headlight's tail-light output is different from the hub input. I haven't scoped it. Supernova tail lights work differently. Some guy on the internet reverse-engineered the headlight circuit showing an LM317 regulator producing 5.9V for the tail light. I have a Supernova E3 tail light (original one; model E161). The case says "6V", which is close to the 5.9V they give it. It wants its 6V, but I don't have a Supernova headlight, so I don't have 6V to give it. I do have a USB charger, so I have 5V instead. Giving it 5V does appear to work, but that results in less brightness than I would like. Presumably the voltage difference is to blame? I took apart the tail light. The circuitry is encased in hot glue. Scraping that off, we can see the action:
The circuit looks like this:
More or less, this is as expected. The circuit was designed to receive 6V, so the resistances (180 ohms) were selected to produce a certain amount of current. Given 5V, there's less current and less brightness. I can fix this by reducing the amount of resistance to bring the current levels up. I hooked up a power supply to produce 5V and 5.9V, and I measured the voltages in the circuit in those two states. Assume little accuracy in all of this
5V 5.9V
V across input diode 0.8V 0.8V
V across the resistors 2.37V 3.25V
V across the LEDs 1.83V 1.85V
As expected, the voltages across the diodes are stable, and the resistors see the bulk of the voltage difference. The circuit designers wanted 3.25/180 ~ 18mA. With my reduced voltage I was getting 2.37/180 ~ 13mA instead. So I was 28% less bright than I should have been. That doesn't sound like a lot, but it's hard to tell by just looking at the thing. In any case, I can reduce the resistances to get the higher current at 5V: I want a resistance R of 180/3.25*2.37 ~ 131 ohms. In the interest of doing less work, I simply added more resistors in parallel instead of replacing the existing ones. So I need to add in parallel a resistance R such that 180R/(180+R) = 131. So I want a parallel R ~ 485 ohms. Looking through my box of parts, I don't have any nice surface-mount resistors with anywhere near the right resistance. But I do have through-hole ones at 470 ohms. Close-enough. I did some soldering gymnastics:
And I'm done. Haven't done any night rides outside yet with this setup, but in theory this should be bright-enough. And since the resistors are just burning off the excess voltage, and I'm giving it less voltage to burn off, I'm being more efficient than I would have been with 6V and the stock resistors.

28 February 2021

Dima Kogan: mrcal: principled camera calibrations

This is a big deal. In my day job I work with images captured by cameras, using those images to infer something about the geometry of the scene being observed. Naturally, to get good results you need to have a good estimate of the behavior of the lens (the "intrinsics"), and of the relative geometry of the cameras (the "extrinsics"; if there's more than one camera). The usual way to do this is to perform a "calibration" procedure to compute the intrinsics and extrinsics, and then to use the resulting "camera model" to process the subsequent images. Wikipedia has an article. And from experience, the most common current toolkit to do this appears to be OpenCV. People have been doing this for a while, but for whatever reason the existing tools all suck. They make basic questions like "how much data should I gather for a calibration?" and "how good is this calibration I just computed?" and "how different are these two models?" unanswerable. This is clearly seen from the links above. The wikipedia article talks about fitting a pinhole model to lenses, even though no real lenses follow this model (telephoto lenses do somewhat; wider lenses don't at all). And the OpenCV tutorial cheerfully says that
Re-projection error gives a good estimation of just how exact the found
parameters are. The closer the re-projection error is to zero, the more accurate
the parameters we found are.
This statement is trivially proven false: throw away most of your calibration data, and your reprojection error becomes very low. But we can all agree that a calibration computed from less data is actually worse. Right? All the various assumptions and hacks in the existing tooling are fine as long as you don't need a whole lot of accuracy out of your results. I need a lot of accuracy, however, so all the existing tools don't work for my applications. So I built a new set of tools, and have been using them with great delight. I just got the blessing to do a public release, so I'm announcing it here. The tools are mrcal does a whole lot to produce calibrations that are as good as possible, and it will tell you just how good they are, and it includes visualization capabilities for extensive user feedback. An overview of the capabilities of the toolkit (with lots of pretty pictures!) is at the tour of mrcal. There's a lot of documentation and examples, but up to now I have been the primary user of the tools. So I expect this to be somewhat rough when others look at it. Bug reports and patches are welcome. mrcal is an excellent base, but it's nowhere near "done". The documentation has some notes about the planned features and improvements, and I'm always reachable by email. Let me know if you try it out!

27 February 2021

Dima Kogan: horizonator: terrain renderer based on SRTM DEMs

Check this out:
I just resurrected and cleaned up an old tool I had lying around. It's now nice and usable by others. This tool loads terrain data, and renders it from the ground, simulating what a human or a camera would see. This is useful for armchair exploring or for identifying peaks. This was relatively novel when I wrote it >10 years ago, but there are a number of similar tools in existence now. This implementation is still useful in that it's freely licensed and contains APIs, so fancier processing can be performed on its output. Sources and (barely-complete-enough) documentation live here:

22 February 2021

Dima Kogan: feedgnuplot: labelled bar charts and a guide

I just released feedgnuplot 1.57, which includes two new pieces that I've long thought about adding:

Labelled bar charts
I've thought about adding these for a while, but had no specific need for them. Finally, somebody asked for it, and I wrote the code. Now that I can, I will probably use these all the time. The new capability can override the usual numerical tic labels on the x axis, and instead use text from a column in the data stream. The most obvious use case is labelled bar graphs:
echo "# label value
      aaa     2
      bbb     3
      ccc     5
      ddd     2"   \
feedgnuplot --vnl \
            --xticlabels \
            --with 'boxes fill solid border lt -1' \
            --ymin 0 --unset grid
But the usage is completely generic. All --xticlabels does, is to accept a data column as labels for the x-axis tics. Everything else that's supported by feedgnuplot and gnuplot works as before. For instance, I can give a domain, and use a style that takes y values and a color:
echo "# x label y color
        5 aaa   2 1
        6 bbb   3 2
       10 ccc   5 4
       11 ddd   2 1"   \
feedgnuplot --vnl --domain \
            --xticlabels \
            --tuplesizeall 3 \
            --with 'points pt 7 ps 2 palette' \
            --xmin 4 --xmax 12 \
            --ymin 0 --ymax 6 \
            --unset grid
And we can use gnuplot's support for clustered histograms:
echo "# x label a b
        5 aaa   2 1
        6 bbb   3 2
       10 ccc   5 4
       11 ddd   2 1"   \
vnl-filter -p label,a,b   \
feedgnuplot --vnl \
            --xticlabels \
            --set 'style data histogram' \
            --set 'style histogram cluster gap 2' \
            --set 'style fill solid border lt -1' \
            --autolegend \
            --ymin 0 --unset grid
Or we can stack the bars on top of one another:
echo "# x label a b
        5 aaa   2 1
        6 bbb   3 2
       10 ccc   5 4
       11 ddd   2 1"   \
vnl-filter -p label,a,b   \
feedgnuplot --vnl \
            --xticlabels \
            --set 'style data histogram' \
            --set 'style histogram rowstacked' \
            --set 'boxwidth 0.8' \
            --set 'style fill solid border lt -1' \
            --autolegend \
            --ymin 0 --unset grid
This is gnuplot's "row stacking". It also supports "column stacking", which effectively transposes the data, and it's not obvious to me that makes sense in the context of feedgnuplot. Similarly, it can label y and/or z axes; I can't think of a specific use case, so I don't have a realistic usage in mind, and I don't support that yet. If anybody can think of a use case, email me. Notes and limitations:
  • Since with --domain you can pass in both an x value and a tic label, it is possible to give it conflicting tic labels for the same x value. gnuplot itself has this problem too, and it just takes the last label it has for a given x. This is probably good-enough.
  • feedgnuplot uses whitespace-separated columns with no escape mechanism, so the field labels cannot have whitespace in it. Fixing this is probably not worth the effort.
  • These tic labels do not count towards the tuplesize
  • I really need to add a similar feature to gnuplotlib. This will happen when I need it or when somebody asks for it, whichever comes first.

A feedgnuplot guide
This fills in a sorely needed missing part of the documentation: the main feedgnuplot website now has a page containing examples and corresponding graphical output. This serves as a tutorial and a gallery demonstrating some usages. It's somewhat incomplete, since it can't show streaming plots, or real-world interfacing with stuff that produces data: some of those usages remain the the README. It's a million times better than what I had before though, which was nothing. Internally this is done just like the gnuplotlib guide: the thing is an org-mode document with org-babel snippets that are evaluated by emacs to make the images. There's some fancy emacs lisp to tie it all together. Works great!

29 July 2020

Dima Kogan: An awk corner case?

So even after years and years of experience, core tools still find ways to surprise me. Today I tried to do some timestamp comparisons with mawk (vnl-filter, to be more precise), and ran into a detail of the language that made it not work. Not a bug, I guess, since both mawk and gawk are affected. I'll claim "language design flaw", however. Let's say I'm processing data with unix timestamps in it (seconds since the epoch). gawk and recent versions of mawk have strftime() for that:
$ date
Wed Jul 29 15:31:13 PDT 2020
$ date +"%s"
$ date +"%s"   mawk ' print strftime("%H",$1) '
And let's say I want to do something conditional on them. I want only data after 9:00 each day:
$ date +"%s"   mawk 'strftime("%H",$1) >= 9  print "Yep. After 9:00" '
That's right. No output. But it is 15:31 now, and I confirmed above that strftime() reports the right time, so it should know that it's after 9:00, but it doesn't. What gives? As we know, awk (and perl after it) treat numbers and strings containing numbers similarly: 5+5 and ="5"+5= both work the same, which is really convenient. This can only work if it can be inferred from context whether we want a number or a string; it knows that addition takes two numbers, so it knows to convert ="5"= into a number in the example above. But what if an operator is ambiguous? Then it picks a meaning based on some internal logic that I don't want to be familiar with. And apparently awk implements string comparisons with the same < and > operators, as numerical comparisons, creating the ambiguity I hit today. strftime returns strings, and you get silent, incorrect behavior that then demands debugging. How to fix? By telling awk to treat the output of strftime() as a number:
$ date +"%s"   mawk '0+strftime("%H",$1) >= 9  print "Yep. After 9:00" '
Yep. After 9:00
With the benefit of hindsight, they really should not have reused any operators for both number and string operations. Then these ambiguities wouldn't occur, and people wouldn't be grumbling into their blogs decades after these decisions were made.

23 July 2020

Dima Kogan: Finding long runs of "notable" data in a log

Here's yet another instance where the data processing I needed done could be acomplished entirely in the shell, with vnlog tools. I have some time-series data in a text table. Via some join and filter operations, I have boiled down this table to a sequence of time indices where something interesting happened. For instance let's say it looks like this: t.vnl
# time
I'd like to find the longest contiguous chunk of time where the interesting thing kept happening. How? Like this!
$ < t.vnl vnl-filter -p 'time,d=diff(time)'  
          vnl-uniq -c -f -1  
          vnl-filter 'd==1' -p 'count=count+1,time=time-1'  
          vnl-sort -nrk count  
# count time
9       4679
5       2011
5       1976
4       1986
3       7291
3       7281
Bam! So the longest run was 9-frames-long, starting at time = 4679.

18 July 2020

Dima Kogan: Converting images while extracting a tarball

Last week at the lab I received a data dump: a gzip-compressed tarball with lots of images in it. The images are all uncompressed .pgm, with the whole tarball weighing in at ~ 1TB. I tried to extract it, and after chugging all day, it ran out of disk space. Added more disk, tried again: out of space again. Just getting a listing of the archive contents (tar tvfz) took something like 8 hours. Clearly this is unreasonable. I made an executive decision to use .jpg files instead: I'd take the small image quality hit for the massive gains in storage efficiency. But the tarball has .pgm and just extracting the thing is challenging. So I'm now extracting the archive, and converting all the .pgm images to .jpg as soon as they hit disk. How? Glad you asked! I'm running two parallel terminal sessions (I'm using screen, but you can do whatever you like).

Session 1
< archive.tar.gz unpigz -p20   tar xv
Here I'm just extracting the archive to disk normally. Using unpigz instead of plain, old tar to get parallelization.

Session 2
inotifywait -r PATH -e close_write -m   mawk -Winteractive '/pgm$/   print $1$3  '   parallel -v -j10 'convert   -quality 96  . .jpg && rm  '
This is the secret sauce. I'm using inotifywait to tell me when any file is closed for writing in a subdirectory of PATH. Then I mawk it to only tell me when .pgm files are done being written, then I convert them to .jpg, and delete the .pgm when that's done. I'm using GNU Parallel to parallelize the image conversion. Otherwise the image conversion doesn't keep up. This is going to take at least all day, but I'm reasonably confident that it will actually finish successfully, and I can then actually do stuff with the data.

16 July 2020

Dima Kogan: Visualizing a broken internet connection

We've all seen our ISPs crap out periodically. Things feel slow, we then go ping some server, and see lost packets or slow response times. When bitching to the ISP it's useful to have evidence so you can clearly explain exactly how they are fucking up. This happened to me again last week, so I wrote a quick oneliner to do the visualization. Here it is:
( echo '# i dt'; ping   perl -ne 'BEGIN   $ =1  next if /PING/; s/.*seq=([0-9]+).*time=([0-9]+).*/\1 \2/ && print')   tee /tmp/ping.vnl   vnl-filter --stream -p i,dt,di='diff(i)'   vnl-filter --stream --has di   feedgnuplot  --domain --stream --with 'linespoints pt 7 palette' --tuplesizeall 3 --ymin 0 --set 'cbrange [1:]' --xlabel "ping index" --ylabel "Response time (ms)" --title "ping response times. Consecutive indices shown as colors"
You run that, and you get a realtime-updating plot of ping times and missed packets. The one-liner also logs the data to a file, so the saved data can be re-visualized. Showing what happens on my currently-working internet connection:
< /tmp/ping.vnl vnl-filter -p i,dt,di='diff(i)'   vnl-filter --has di   feedgnuplot  --domain --with 'linespoints pt 7 palette' --tuplesizeall 3 --ymin 0 --set 'cbrange [1:]' --xlabel "ping index" --ylabel "Response time (ms)" --title "ping response times. Consecutive indices shown as colors" --hardcopy /tmp/isp-good.svg
On a broken network it looks like this (I edited the log to show the kind of broken-ness I was seeing earlier):
The x axis is the ping index. The y axis is the response time. Misses responses are shows as gaps in the points and also as colors, for legibility. As usual, this uses the vnlog and feedgnuplot tools, so
sudo apt install feedgnuplot vnlog
if you want to try this

20 June 2020

Dima Kogan: OpenCV C API transition. A rant.

I just went through a debugging exercise that was so ridiculous, I just had to write it up. Some of this probably should go into a bug report instead of a rant, but I'm tired. And clearly I don't care anymore. OK, so I'm doing computer vision work. OpenCV has been providing basic functions in this area, so I have been using them for a while. Just for really, really basic stuff, like projection. The C API was kinda weird, and their error handling is a bit ridiculous (if you give it arguments it doesn't like, it asserts!), but it has been working fine for a while. At some point (around OpenCV 3.0) somebody over there decided that they didn't like their C API, and that this was now a C++ library. Except the docs still documented the C API, and the website said it supported C, and the code wasn't actually removed. They just kinda stopped testing it and thinking about it. So it would mostly continue to work, except some poor saps would see weird failures; like this and this, for instance. OpenCV 3.2 was the last version where it was mostly possible to keep using the old C code, even when compiling without optimizations. So I was doing that for years. So now, in 2020, Debian is finally shipping a version of OpenCV that definitively does not work with the old code, so I had to do something. Over time I stopped using everything about OpenCV, except a few cvProjectPoints2() calls. So I decided to just write a small C++ shim to call the new version of that function, expose that with =extern "C"= to the rest of my world, and I'd be done. And normally I would be, but this is OpenCV we're talking about. I wrote the shim, and it didn't work. The code built and ran, but the results were wrong. After some pointless debugging, I boiled the problem down to this test program:
#include <opencv2/calib3d.hpp>
#include <stdio.h>
int main(void)
    double fx = 1000.0;
    double fy = 1000.0;
    double cx = 1000.0;
    double cy = 1000.0;
    double _camera_matrix[] =
          fx,  0, cx,
          0,  fy, cy,
          0,   0,  1  ;
    cv::Mat camera_matrix(3,3, CV_64FC1, _camera_matrix);
    double pp[3] =  1., 2., 10. ;
    double qq[2] =  444, 555 ;
    int N=1;
    cv::Mat object_points(N,3, CV_64FC1, pp);
    cv::Mat image_points (N,2, CV_64FC1, qq);
    // rvec,tvec
    double _zero3[3] =  ;
    cv::Mat zero3(1,3,CV_64FC1, _zero3);
    cv::projectPoints( object_points,
                       cv::noArray(), 0.0);
    fprintf(stderr, "manually-projected no-distortion: %f %f\n",
            pp[0]/pp[2] * fx + cx,
            pp[1]/pp[2] * fy + cy);
    fprintf(stderr, "opencv says: %f %f\n", qq[0], qq[1]);
    return 0;
This is as trivial as it gets. I project one point through a pinhole camera, and print out the right answer (that I can easily compute, since this is trivial), and what OpenCV reports:
$ g++ -I/usr/include/opencv4 -o tst -lopencv_calib3d -lopencv_core && ./tst
manually-projected no-distortion: 1100.000000 1200.000000
opencv says: 444.000000 555.000000
Well that's no good. The answer is wrong, but it looks like it didn't even write anything into the output array. Since this is supposed to be a thin shim to C code, I want this thing to be filling in C arrays, which is what I'm doing here:
double qq[2] =  444, 555 ;
int N=1;
cv::Mat image_points (N,2, CV_64FC1, qq);
This is how the C API has worked forever, and their C++ API works the same way, I thought. Nothing barfed, not at build time, or run time. Fine. So I went to figure this out. In the true spirit of C++, the new API is inscrutable. I'm passing in cv::Mat, but the API wants cv::InputArray for some arguments and cv::OutputArray for others. Clearly cv::Mat can be coerced into either of those types (and that's what you're supposed to do), but the details are not meant to be understood. You can read the snazzy C++-style documentation. Clicking on "OutputArray" in the doxygen gets you here. Then I guess you're supposed to click on "_OutputArray", and you get here. Understand what's going on now? Me neither. Stepping through the code revealed the problem. cv::projectPoints() looks like this:
void cv::projectPoints( InputArray _opoints,
                        InputArray _rvec,
                        InputArray _tvec,
                        InputArray _cameraMatrix,
                        InputArray _distCoeffs,
                        OutputArray _ipoints,
                        OutputArray _jacobian,
                        double aspectRatio )
    _ipoints.create(npoints, 1, CV_MAKETYPE(depth, 2), -1, true);
I.e. they're allocating a new data buffer for the output, and giving it back to me via the OutputArray object. This object already had a buffer, and that's where I was expecting the output to go. Instead it went to the brand-new buffer I didn't want. Issues: Well that's just super. I can call the C++ function, copy the data into the place it's supposed to go to, and then deallocate the extra buffer. Or I can pull out the meat of the function I want into my project, and then I can drop the OpenCV dependency entirely. Clearly that's the way to go. So I go poking back into their code to grab what I need, and here's what I see:
static void cvProjectPoints2Internal( const CvMat* objectPoints,
                  const CvMat* r_vec,
                  const CvMat* t_vec,
                  const CvMat* A,
                  const CvMat* distCoeffs,
                  CvMat* imagePoints, CvMat* dpdr CV_DEFAULT(NULL),
                  CvMat* dpdt CV_DEFAULT(NULL), CvMat* dpdf CV_DEFAULT(NULL),
                  CvMat* dpdc CV_DEFAULT(NULL), CvMat* dpdk CV_DEFAULT(NULL),
                  CvMat* dpdo CV_DEFAULT(NULL),
                  double aspectRatio CV_DEFAULT(0) )
Looks familiar? It should. Because this is the original C-API function they replaced. So in their quest to move to C++, they left the original code intact, C API and everything, un-exposed it so you couldn't call it anymore, and made a new, shitty C++ wrapper for people to call instead. CvMat is still there. I have no words. Yes, this is a massive library, and maybe other parts of it indeed did make some sort of non-token transition, but this thing is ridiculous. In the end, here's the function I ended up with (licensed as OpenCV; see the comment)
// The implementation of project_opencv is based on opencv. The sources have
// been heavily modified, but the opencv logic remains. This function is a
// cut-down cvProjectPoints2Internal() to keep only the functionality I want and
// to use my interfaces. Putting this here allows me to drop the C dependency on
// opencv. Which is a good thing, since opencv dropped their C API
// from opencv-4.2.0+dfsg/modules/calib3d/src/calibration.cpp
// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
// Third party copyrights are property of their respective owners.
// Redistribution and use in source and binary forms, with or without modification,
// are permitted provided that the following conditions are met:
//   * Redistribution's of source code must retain the above copyright notice,
//     this list of conditions and the following disclaimer.
//   * Redistribution's in binary form must reproduce the above copyright notice,
//     this list of conditions and the following disclaimer in the documentation
//     and/or other materials provided with the distribution.
//   * The name of the copyright holders may not be used to endorse or promote products
//     derived from this software without specific prior written permission.
// This software is provided by the copyright holders and contributors "as is" and
// any express or implied warranties, including, but not limited to, the implied
// warranties of merchantability and fitness for a particular purpose are disclaimed.
// In no event shall the Intel Corporation or contributors be liable for any direct,
// indirect, incidental, special, exemplary, or consequential damages
// (including, but not limited to, procurement of substitute goods or services;
// loss of use, data, or profits; or business interruption) however caused
// and on any theory of liability, whether in contract, strict liability,
// or tort (including negligence or otherwise) arising in any way out of
typedef union
        double x,y;
    double xy[2];
typedef union
        double x,y,z;
    double xyz[3];
void project_opencv( // outputs
                     point2_t* q,
                     point3_t* dq_dp,               // may be NULL
                     double* dq_dintrinsics_nocore, // may be NULL
                     // inputs
                     const point3_t* p,
                     int N,
                     const double* intrinsics,
                     int Nintrinsics)
    const double fx = intrinsics[0];
    const double fy = intrinsics[1];
    const double cx = intrinsics[2];
    const double cy = intrinsics[3];
    double k[12] =  ;
    for(int i=0; i<Nintrinsics-4; i++)
        k[i] = intrinsics[i+4];
    for( int i = 0; i < N; i++ )
        double z_recip = 1./p[i].z;
        double x = p[i].x * z_recip;
        double y = p[i].y * z_recip;
        double r2      = x*x + y*y;
        double r4      = r2*r2;
        double r6      = r4*r2;
        double a1      = 2*x*y;
        double a2      = r2 + 2*x*x;
        double a3      = r2 + 2*y*y;
        double cdist   = 1 + k[0]*r2 + k[1]*r4 + k[4]*r6;
        double icdist2 = 1./(1 + k[5]*r2 + k[6]*r4 + k[7]*r6);
        double xd      = x*cdist*icdist2 + k[2]*a1 + k[3]*a2 + k[8]*r2+k[9]*r4;
        double yd      = y*cdist*icdist2 + k[2]*a3 + k[3]*a1 + k[10]*r2+k[11]*r4;
        q[i].x = xd*fx + cx;
        q[i].y = yd*fy + cy;
        if( dq_dp )
            double dx_dp[] =   z_recip, 0,       -x*z_recip  ;
            double dy_dp[] =   0,       z_recip, -y*z_recip  ;
            for( int j = 0; j < 3; j++ )
                double dr2_dp = 2*x*dx_dp[j] + 2*y*dy_dp[j];
                double dcdist_dp = k[0]*dr2_dp + 2*k[1]*r2*dr2_dp + 3*k[4]*r4*dr2_dp;
                double dicdist2_dp = -icdist2*icdist2*(k[5]*dr2_dp + 2*k[6]*r2*dr2_dp + 3*k[7]*r4*dr2_dp);
                double da1_dp = 2*(x*dy_dp[j] + y*dx_dp[j]);
                double dmx_dp = (dx_dp[j]*cdist*icdist2 + x*dcdist_dp*icdist2 + x*cdist*dicdist2_dp +
                                k[2]*da1_dp + k[3]*(dr2_dp + 4*x*dx_dp[j]) + k[8]*dr2_dp + 2*r2*k[9]*dr2_dp);
                double dmy_dp = (dy_dp[j]*cdist*icdist2 + y*dcdist_dp*icdist2 + y*cdist*dicdist2_dp +
                                k[2]*(dr2_dp + 4*y*dy_dp[j]) + k[3]*da1_dp + k[10]*dr2_dp + 2*r2*k[11]*dr2_dp);
                dq_dp[i*2 + 0].xyz[j] = fx*dmx_dp;
                dq_dp[i*2 + 1].xyz[j] = fy*dmy_dp;
        if( dq_dintrinsics_nocore )
            dq_dintrinsics_nocore[(Nintrinsics-4)*(2*i + 0) + 0] = fx*x*icdist2*r2;
            dq_dintrinsics_nocore[(Nintrinsics-4)*(2*i + 1) + 0] = fy*(y*icdist2*r2);
            dq_dintrinsics_nocore[(Nintrinsics-4)*(2*i + 0) + 1] = fx*x*icdist2*r4;
            dq_dintrinsics_nocore[(Nintrinsics-4)*(2*i + 1) + 1] = fy*y*icdist2*r4;
            if( Nintrinsics-4 > 2 )
                dq_dintrinsics_nocore[(Nintrinsics-4)*(2*i + 0) + 2] = fx*a1;
                dq_dintrinsics_nocore[(Nintrinsics-4)*(2*i + 1) + 2] = fy*a3;
                dq_dintrinsics_nocore[(Nintrinsics-4)*(2*i + 0) + 3] = fx*a2;
                dq_dintrinsics_nocore[(Nintrinsics-4)*(2*i + 1) + 3] = fy*a1;
                if( Nintrinsics-4 > 4 )
                    dq_dintrinsics_nocore[(Nintrinsics-4)*(2*i + 0) + 4] = fx*x*icdist2*r6;
                    dq_dintrinsics_nocore[(Nintrinsics-4)*(2*i + 1) + 4] = fy*y*icdist2*r6;
                    if( Nintrinsics-4 > 5 )
                        dq_dintrinsics_nocore[(Nintrinsics-4)*(2*i + 0) + 5] = fx*x*cdist*(-icdist2)*icdist2*r2;
                        dq_dintrinsics_nocore[(Nintrinsics-4)*(2*i + 1) + 5] = fy*y*cdist*(-icdist2)*icdist2*r2;
                        dq_dintrinsics_nocore[(Nintrinsics-4)*(2*i + 0) + 6] = fx*x*cdist*(-icdist2)*icdist2*r4;
                        dq_dintrinsics_nocore[(Nintrinsics-4)*(2*i + 1) + 6] = fy*y*cdist*(-icdist2)*icdist2*r4;
                        dq_dintrinsics_nocore[(Nintrinsics-4)*(2*i + 0) + 7] = fx*x*cdist*(-icdist2)*icdist2*r6;
                        dq_dintrinsics_nocore[(Nintrinsics-4)*(2*i + 1) + 7] = fy*y*cdist*(-icdist2)*icdist2*r6;
                        if( Nintrinsics-4 > 8 )
                            dq_dintrinsics_nocore[(Nintrinsics-4)*(2*i + 0) + 8] = fx*r2; //s1
                            dq_dintrinsics_nocore[(Nintrinsics-4)*(2*i + 1) + 8] = fy*0; //s1
                            dq_dintrinsics_nocore[(Nintrinsics-4)*(2*i + 0) + 9] = fx*r4; //s2
                            dq_dintrinsics_nocore[(Nintrinsics-4)*(2*i + 1) + 9] = fy*0; //s2
                            dq_dintrinsics_nocore[(Nintrinsics-4)*(2*i + 0) + 10] = fx*0;//s3
                            dq_dintrinsics_nocore[(Nintrinsics-4)*(2*i + 1) + 10] = fy*r2; //s3
                            dq_dintrinsics_nocore[(Nintrinsics-4)*(2*i + 0) + 11] = fx*0;//s4
                            dq_dintrinsics_nocore[(Nintrinsics-4)*(2*i + 1) + 11] = fy*r4; //s4
This does only the stuff I need: projection only (no geometric transformation), and gradients in respect to the point coordinates and distortions only. Gradients in respect to fxy and cxy are trivial, and I don't bother reporting them. So now I don't compile or link against OpenCV, my code builds and runs on Debian/sid and (surprisingly) it runs much faster than before. Apparently there was a lot of pointless overhead happening. Alright. Rant over.

3 June 2020

Dima Kogan: vnlog now functional on *BSD and OSX

So somebody finally bugged me about supporting vnlog tools on OSX. I was pretty sure that between all the redirection, process communication, and file descriptor management something was Linux-specific, but apparently not: everything just works. There were a few uninteresting issues with tool paths, core tool and linker flags and so on, but it was all really minor. I have a report that the test suite passes on OSX, and I verified it on FreeBSD. I made a new 1.28 release tag, but it exists mostly for the benefit of any OSX or *BSD people who'd want to make a package for their system. Progress!

4 May 2020

Dima Kogan: Redirection and /dev/stderr

I just hit some unexpected unix-ism (or maybe Linux-ism?) that I'd like to mention here. I found a workaround, but if anybody knows what's actually happening and can enlighten me, please send me an email, and I'll update the post. OK, so I have a program that does stuff, and outputs its progress to stderr. And I'm invoking this from a shell script that also outputs stuff to stderr. A minimized sample:
echo "Outer print 1111" > /dev/stderr
perl -E 'say STDERR ("1111 1111 1111 1111 1111 1111 1111 1111\n")'
echo "Outer print 2222" > /dev/stderr
perl -E 'say STDERR ("2222 2222 2222 2222 2222 2222 2222 2222\n")'
Let's call this This works just fine. Now let's say I want to log the output to a file. Naturally I do this:
./ 2> /tmp/log
Surprisingly, this does not work:
$ bash 2> /tmp/log && cat /tmp/log
Outer print 2222
2222 2222 2222 2222 2222 2222 2222 2222
$ bash 2> /tmp/log && xxd /tmp/log
0000000: 4f75 7465 7220 7072 696e 7420 3232 3232  Outer print 2222
0000010: 0a00 0000 0000 0000 0000 0000 0000 0000  ................
0000020: 0000 0000 0000 0000 0032 3232 3220 3232  .........2222 22
0000030: 3232 2032 3232 3220 3232 3232 2032 3232  22 2222 2222 222
0000040: 3220 3232 3232 2032 3232 3220 3232 3232  2 2222 2222 2222
0000050: 0a0a                                     ..
The 1111 prints are gone, and we now have a block of binary 0 bytes in there for some reason. I had been assuming that writes to stderr block until they're finished. So I can, without thinking about it, mix writing to fd 2 and opening, and then writing to /dev/stderr. Apparently the right way to do this is to write to fd 2 in the script directly, instead of opening /dev/stderr:
echo "Outer print 1111" >&2
perl -E 'say STDERR ("1111 1111 1111 1111 1111 1111 1111 1111\n")'
echo "Outer print 2222" >&2
perl -E 'say STDERR ("2222 2222 2222 2222 2222 2222 2222 2222\n")'

Update So I poked at it for a bit more, and it turns out that the obvious scapegoat (buffering) is actually not the issue. The problem is this:

$ ls -l /dev/stderr(:A)
-rw-r--r-- 1 dima dima 36 May  4 17:02 /tmp/log
:A is the zsh syntax for recursive link-following. /dev/stderr is thus ultimately a link to /tmp/log. So the file on disk /tmp/log is being opened for writing and written-to twice at the same time
  • once in the outer 2>/tmp/log redirection
  • again in the inner >/dev/stderr redirection
Now we know. Thanks to Larry Doolittle for pointers.

23 March 2020

Dima Kogan: org-babel for documentation

So I just gave a talk at SCaLE 18x about numpysane and gnuplotlib, two libraries I wrote to make using numpy bearable. With these two, it's actually quite nice! Prior to the talk I overhauled the documentation for both these projects. The gnuplotlib docs now have a tutorial/gallery page, which is interesting-enough to write about. Check it out! Mostly it is a sequence of Clearly you want the plots in the documentation to correspond to the code, so you want something to actually run each code snippet to produce each plot. Automatically. I don't want to maintain these manually, and periodically discover that the code doesn't make the plot I claim it does or worse: that the code barfs. This is vaguely what Jupyter notebooks do, but they're ridiculous, so I'm doing something better: That's it. The git repo is hosted by github, which has a rudimentary renderer for .org documents. I'm committing the .svg files, so that's enough to get rendered documentation that looks nice. Note that the usual workflow is to use org to export to html, but here I'm outsourcing that job to github; I just make the .svg files, and that's enough. Look at the link again: gnuplotlib tutorial/gallery. This is just a .org file committed to the git repo. github is doing its normal org->html thing to display this file. This has drawbacks too: github is ignoring the :noexport: tag on the init section at the end of the file, so it's actually showing all the emacs lisp goop that makes this work (described below!). It's at the end, so I guess this is good-enough. Those of us that use org-mode would be completely unsurprised to hear that the talk is also written as .org document. And the slides that show gnuplotlib plots use the same org-babel system to render the plots. It's all oh-so-nice. As with anything as flexible as org-babel, it's easy to get into a situation where you're bending it to serve a not-quite-intended purpose. But since this all lives in emacs, you can make it do whatever you want with a bit of emacs lisp. I ended up advising a few things (mailing list post here). And I stumbled on an (arguable) bug in emacs that needed working around (mailing list post here). I'll summarize both here.

Handling large Local Variables blocks
The advises I ended up with ended up longer than emacs expected, which made emacs not evaluate them when loading the buffer. As I discovered (see the mailing list post) the loading code looks for the string Local Variables in the last 3000 bytes of the buffer only, and I exceeded that. Stefan Monnier suggested a workaround in this post. Instead of the normal Local Variables block at the end:
Local Variables:
eval: (progn ... ...
             ... ...
             LONG chunk of emacs-lisp
I do this:
(progn ;;local-config
   lisp lisp lisp
   as long as I want
Local Variables:
eval: (progn (re-search-backward "^(progn ;;local-config") (eval (read (current-buffer))))
So emacs sees a small chunk of code that searches backwards through the buffer (as far back as needed) for the real lisp to evaluate. As an aside, this blog is also an .org document, and the lisp snippets above are org-babel blocks that I'm not evaluating. The exporter knows to respect the emacs-lisp syntax highlighting, however.

OK, so what was all the stuff I needed to tell org-babel to do specially here? First off, org needed to be able to communicate to the Python session the name of the file to write the plot to. I do this by making the whole plist for this org-babel snippet available to python:
;; THIS advice makes all the org-babel parameters available to python in the
;; _org_babel_params dict. I care about _org_babel_params['_file'] specifically,
;; but everything is available
(defun dima-org-babel-python-var-to-python (var)
  "Convert an elisp value to a python variable.
  Like the original, but supports (a . b) cells and symbols
  (if (listp var)
      (if (listp (cdr var))
          (concat "[" (mapconcat #'org-babel-python-var-to-python var ", ") "]")
        (format "\"\"\"%s\"\"\"" var))
    (if (symbolp var)
        (format "\"\"\"%s\"\"\"" var)
      (if (eq var 'hline)
         (if (and (stringp var) (string-match "[\n\r]" var)) "\"\"%S\"\"" "%S")
         (if (stringp var) (substring-no-properties var) var))))))
(defun dima-alist-to-python-dict (alist)
  "Generates a string defining a python dict from the given alist"
  (let ((keyvalue-list
         (mapcar (lambda (x)
                   (format "%s = %s, "
                            "[^a-zA-Z0-9_]" "_"
                            (symbol-name (car x)))
                           (dima-org-babel-python-var-to-python (cdr x))))
     "dict( "
     (apply 'concat keyvalue-list)
(defun dima-org-babel-python-pass-all-params (f params)
    "_org_babel_params = "
    (dima-alist-to-python-dict params))
   (funcall f params)))
   :around #'dima-org-babel-python-pass-all-params))
So if there's a :file plist key, the python code can grab that, and write the plot to that filename. But I don't really want to specify an output file for every single org-babel snippet. All I really care about is that each plot gets a unique filename. So I omit the :file key entirely, and use this advice to generate one for me:
;; This sets a default :file tag, set to a unique filename. I want each demo to
;; produce an image, but I don't care what it is called. I omit the :file tag
;; completely, and this advice takes care of it
(defun dima-org-babel-python-unique-plot-filename
    (f &optional arg info params)
  (funcall f arg info
           (cons (cons ':file
                       (format "guide-%d.svg"
                               (condition-case nil
                                   (setq dima-unique-plot-number (1+ dima-unique-plot-number))
                                 (error (setq dima-unique-plot-number 0)))))
   :around #'dima-org-babel-python-unique-plot-filename))
This uses the dima-unique-plot-number integer to keep track of each plot. I increment this with each plot. Getting closer. It isn't strictly required, but it'd be nice if each plot had the same output filename each time I generated it. So I want to reset the plot number to 0 each time:
;; If I'm regenerating ALL the plots, I start counting the plots from 0
(defun dima-reset-unique-plot-number
    (&rest args)
    (setq dima-unique-plot-number 0))
   :after #'dima-reset-unique-plot-number))
Finally, I want to lie to the user a little bit. The code I'm actually executing writes each plot to an .svg. But the code I'd like the user to see should use the default output: an interactive, graphical window. I do that by tweaking the python session to tell the gnuplotlib object to write to .svg files from org by default, instead of using the graphical terminal:
;; I'm using github to display, so I'm not using the "normal" org
;; exporter. I want the demo text to not contain the hardcopy= tags, but clearly
;; I need the hardcopy tag when generating the plots. I add some python to
;; override gnuplotlib.plot() to add the hardcopy tag somewhere where the reader
;; won't see it. But where to put this python override code? If I put it into an
;; org-babel block, it will be rendered, and the :export tags will be ignored,
;; since github doesn't respect those (probably). So I put the extra stuff into
;; an advice. Whew.
(defun dima-org-babel-python-set-demo-output (f body params)
    (insert body)
    (when (search-forward "import gnuplotlib as gp" nil t)
       "if not hasattr(gp.gnuplotlib, 'orig_init'):\n"
       "    gp.gnuplotlib.orig_init = gp.gnuplotlib.__init__\n"
       "gp.gnuplotlib.__init__ = lambda self, *args, **kwargs: gp.gnuplotlib.orig_init(self, *args, hardcopy=_org_babel_params['_file'] if 'file' in _org_babel_params['_result_params'] else None, **kwargs)\n"))
    (setq body (buffer-substring-no-properties (point-min) (point-max))))
  (funcall f body params))
   :around #'dima-org-babel-python-set-demo-output))
And that's it. The advises in the talk are slightly different, in uninteresting ways. Some of this should be upstreamed to org-babel somehow. Now entirely clear which part, but I'll cross that bridge when I get to it.

15 March 2020

Dima Kogan: numpysane and broadcasting in C

Since the beginning, the numpysane library provided a broadcast_define() function to decorate existing Python routines to give them broadcasting awareness. This was very useful, but slow. I just did lots of typing, and now I have a flavor of this in C (the numpysane_pywrap module; new in numpysane 0.22). As expected, you get fast C loops! And similar to the rest of this library, this is a port of something in PDL: PDL::PP. Full documentation lives here: After writing this I realized that there was something similar available in numpy this whole time: I haven't looked too deeply into this yet, but 2 things are clear: There's a design difference: the numpy implementation uses function callbacks, while I generate C code. Code generation is what PDL::PP does, and when I thought about it earlier, it seemed like doing this with function pointers would be too painful. I guess it's doable, though. And at least in one case, the gufuncs aren't doing the right broadcasting thing:
>>> a = np.arange(5).reshape(5,1)
>>> b = np.arange(3)
>>> np.matmul(a,b)
ValueError: matmul: Input operand 1 has a mismatch in
   its core dimension 0, with gufunc signature
   (n?,k),(k,m?)->(n?,m?) (size 3 is different from 1)
This should work. And if you do this with numpysane.broadcast_define() or with numpysane_pywrap, it does work. I'll look at it later to figure out what it's doing.

10 July 2016

Bits from Debian: New Debian Developers and Maintainers (May and June 2016)

The following contributors got their Debian Developer accounts in the last two months: The following contributors were added as Debian Maintainers in the last two months: Congratulations!

4 October 2015

Johannes Schauer: new sbuild release 0.66.0

I just released sbuild 0.66.0-1 into unstable. It fixes a whopping 30 bugs! Thus, I'd like to use this platform to: And a super big thank you to Roger Leigh who, despite having resigned from Debian, was always available to give extremely helpful hints, tips, opinion and guidance with respect to sbuild development. Thank you! Here is a list of the major changes since the last release: