Density for a layer

okguy

It would be nice to have a module to be added in the menus to help us to calculate the density of a layer.
Thanks to this module : http://www.klayout.de/useful_scripts.html#cell_bbox.rbm
and this routine : klayout.de/forum/comments.php?DiscussionID=164

But, I don't know how to make a window to select a layer to calculate its density :(

Thanks,
OkGuy

Matthias

Hi OkGuy,

here is some derived script that adds some "user interface". It computes the density on the selected layer and asks for the tile dimensions in a input dialog. It will also use the cell's bounding box to derive the tiling area.

This is hopefully some acceptable compromise. "Real" user interface are possible using qtruby or the Qt integration which will come with the next version.

Best regards,

Matthias

$density_map = RBA::Action.new
$density_map.title = "Density Map"
$density_map.on_triggered do 

  app = RBA::Application.instance
  mw = app.main_window

  lv = mw.current_view
  if !lv || !lv.active_cellview
    raise "Multiclip: No view or no layout active"
  end

  cl = lv.current_layer
  if cl.is_null? || cl.current.layer_index < 0
    raise "No layer selected to create the density map from"
  end
  input_layer = cl.current.layer_index

  value = RBA::InputDialog.get_string("Tile Size", "Enter the tile size in micron (x,y)", "50,50");
  if !value.has_value?
    return
  end
  value = value.value
  if value !~ /^\d*(\.\d*)?,\d*(\.\d*)?$/
    raise "Expected a tile size (x, y) in micron"
  end

  va = value.split(/,/)
  xw = va[0].to_f
  yw = va[1].to_f
  if xw < 1.0e-6 || yw < 1.0e-6
    raise "Invalid tile size"
  end

  # Let's start ...

  # Set up the environment:
  app = RBA::Application.instance
  mw = app.main_window

  lv = mw.current_view
  if !lv || !lv.active_cellview
    raise "Density map: No view or no layout active"
  end

  # obtain the pointers to the layout and cell 
  lay = lv.active_cellview.layout
  top = lv.active_cellview.cell_index
  top_name = lay.cell_name(top)
  top_bbox = lay.cell(top).bbox
  dbu = lay.dbu
  bbox = lv.active_cellview.cell.bbox
  xmin = bbox.left * dbu
  xmax = bbox.right * dbu
  ymin = bbox.bottom * dbu
  ymax = bbox.top * dbu

  # Start to collect the data:
  data = []

  # proceed row by row
  # (for each row we use the multi_clip_into which is more efficient than single clips per window)
  nrows = 0
  y = ymin
  while y < ymax-1e-6

    # Prepare a new layout to receive the clips
    # Hint: we need to clip all layers (clip does not support clipping of one layer alone currently).
    cl = RBA::Layout.new
    cl.dbu = dbu
    lay.layer_indices.each do |li|
      cl.insert_layer_at(li, lay.get_info(li))
    end

    # Prepare the clip boxes for this row.
    # Note: because clip only works with boxes that overlap the cell's bounding box currently, we
    # have to operate with a subset of fields to support the general case.
    boxes = []
    x = xmin
    ncolumns = 0
    colstart = nil
    while x < xmax-1e-6
      b = RBA::Box.new((0.5 + x / dbu).floor, (0.5 + y / dbu).floor, (0.5 + (x + xw) / dbu).floor, (0.5 + (y + yw) / dbu).floor);
      if b.overlaps?(top_bbox)
        colstart ||= ncolumns 
        boxes.push(b)
      end
      x += xw
      ncolumns += 1
    end

    columns = []
    ncolumns.times { columns.push(0.0) }

    $stdout.write "Running row y=#{y}..#{y+yw} "

    if colstart

      col = colstart

      # Actually do the clip
      cells = lay.multi_clip_into(top, cl, boxes)

      ep = RBA::ShapeProcessor.new
      merged = RBA::Shapes.new

      # Compute the area within area box
      cells.each do |c|

        # merge the shapes to polygons and compute the area
        ep.merge(cl, cl.cell(c), input_layer, merged, true, 0, false, false)
        a = 0
        merged.each do |m|
          a += m.polygon.area
        end

        # compute and store the density
        a = a * dbu * dbu / (xw * xw)
        columns[col] = a

        col += 1

        $stdout.write((0.5 + a * 100).floor.to_s)
        $stdout.write " "
        $stdout.flush

      end

    end

    puts ""

    columns.each { |d| data.push(d) }

    y += yw
    nrows += 1

  end

  if nrows * ncolumns > 0
    img = RBA::Image.new(ncolumns, nrows, data)
    img.pixel_width = xw
    img.pixel_height = yw
    img.trans = RBA::DCplxTrans.new(1.0, 0.0, false, RBA::DPoint.new(xmin, ymin))
    if $density_map_image_id
      lv.replace_image($density_map_image_id, img) 
    else
      lv.insert_image(img) 
      $density_map_image_id = img.id
    end
  end

end


app = RBA::Application.instance
mw = app.main_window

menu = mw.menu
menu.insert_separator("@toolbar.end", "name")
menu.insert_item("@toolbar.end", "density_map", $density_map)

okguy

Matthias,

I modified it to include it in the tools menu.
However, I did expected a figure:
Layer density = ..% of the selected area

How to calculate the area of the selected layer / the total area ?

Thanks,
OkGuy

#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# DESCRIPTION: Compute area of selected shapes
#
# Run the script with
# klayout -rm density_map.rbm ...
# or put the script as "density_map.rbm" into the installation path (on Unix for version <=0.21:
# set $KLAYOUTPATH to the installation folder).
#

class MenuAction < RBA::Action
def initialize( title, shortcut, &action )
self.title = title
self.shortcut = shortcut
@action = action
end
def triggered
@action.call( self )
end
private
@action
end

# $density_map = RBA::Action.new
# $density_map.title = "Density Map"

$density_map = MenuAction.new( "Layer density map", "" ) do

app = RBA::Application.instance
mw = app.main_window

lv = mw.current_view
if !lv || !lv.active_cellview
raise "Multiclip: No view or no layout active"
end

cl = lv.current_layer
if cl.is_null? || cl.current.layer_index < 0
raise "No layer selected to create the density map from"
end
input_layer = cl.current.layer_index

value = RBA::InputDialog.get_string("Tile Size", "Enter the tile size in micron (x,y)", "50,50");
if !value.has_value?
return
end
value = value.value
if value !~ /^\d*(\.\d*)?,\d*(\.\d*)?$/
raise "Expected a tile size (x, y) in micron"
end

va = value.split(/,/)
xw = va[0].to_f
yw = va[1].to_f
if xw < 1.0e-6 || yw < 1.0e-6
raise "Invalid tile size"
end

# Let's start ...

# Set up the environment:
app = RBA::Application.instance
mw = app.main_window

lv = mw.current_view
if !lv || !lv.active_cellview
raise "Density map: No view or no layout active"
end

# obtain the pointers to the layout and cell
lay = lv.active_cellview.layout
top = lv.active_cellview.cell_index
top_name = lay.cell_name(top)
top_bbox = lay.cell(top).bbox
dbu = lay.dbu
bbox = lv.active_cellview.cell.bbox
xmin = bbox.left * dbu
xmax = bbox.right * dbu
ymin = bbox.bottom * dbu
ymax = bbox.top * dbu

# Start to collect the data:
data = []

# proceed row by row
# (for each row we use the multi_clip_into which is more efficient than single clips per window)
nrows = 0
y = ymin
while y < ymax-1e-6

# Prepare a new layout to receive the clips
# Hint: we need to clip all layers (clip does not support clipping of one layer alone currently).
cl = RBA::Layout.new
cl.dbu = dbu
lay.layer_indices.each do |li|
cl.insert_layer_at(li, lay.get_info(li))
end

# Prepare the clip boxes for this row.
# Note: because clip only works with boxes that overlap the cell's bounding box currently, we
# have to operate with a subset of fields to support the general case.
boxes = []
x = xmin
ncolumns = 0
colstart = nil
while x < xmax-1e-6
b = RBA::Box.new((0.5 + x / dbu).floor, (0.5 + y / dbu).floor, (0.5 + (x + xw) / dbu).floor, (0.5 + (y + yw) / dbu).floor);
if b.overlaps?(top_bbox)
colstart ||= ncolumns
boxes.push(b)
end
x += xw
ncolumns += 1
end

columns = []
ncolumns.times { columns.push(0.0) }

$stdout.write "Running row y=#{y}..#{y+yw} "

if colstart

col = colstart

# Actually do the clip
cells = lay.multi_clip_into(top, cl, boxes)

ep = RBA::ShapeProcessor.new
merged = RBA::Shapes.new

# Compute the area within area box
cells.each do |c|

# merge the shapes to polygons and compute the area
ep.merge(cl, cl.cell(c), input_layer, merged, true, 0, false, false)
a = 0
merged.each do |m|
a += m.polygon.area
end

# compute and store the density
a = a * dbu * dbu / (xw * xw)
columns[col] = a

col += 1

$stdout.write((0.5 + a * 100).floor.to_s)
$stdout.write " "
$stdout.flush

end

end

puts ""

columns.each { |d| data.push(d) }

y += yw
nrows += 1

end

if nrows * ncolumns > 0
img = RBA::Image.new(ncolumns, nrows, data)
img.pixel_width = xw
img.pixel_height = yw
img.trans = RBA::DCplxTrans.new(1.0, 0.0, false, RBA::DPoint.new(xmin, ymin))
if $density_map_image_id
lv.replace_image($density_map_image_id, img)
else
lv.insert_image(img)
$density_map_image_id = img.id
end
end

end


app = RBA::Application.instance
mw = app.main_window

menu = mw.menu
menu.insert_item("tools_menu.end", "density_map", $density_map)

Matthias

Hi OkGuy,

from the explanation you gave initially I thought you might be interested in the density map. At least that is what the forum entry is about. Could you be somewhat specific what exactly you require? Do you need the total density? Do you need the density map at all?

Even if you require the overall density it may be a good idea to compute a density map first and then average over the individual tiles. That avoids memory allocation issues with large layouts because the algorithm works on a tile-by-tile basis. You still need to enter a tile size which should be reasonably small (i.e 500x500 micron max). Please note that the chip dimensions are rounded to the tile size, to the results are most accurate when the tile size is an integer fraction of the chip size.

The code change would be to replace that piece:

if nrows * ncolumns > 0
  img = RBA::Image.new(ncolumns, nrows, data)
  img.pixel_width = xw
  img.pixel_height = yw
  img.trans = RBA::DCplxTrans.new(1.0, 0.0, false, RBA::DPoint.new(xmin, ymin))
  if $density_map_image_id
    lv.replace_image($density_map_image_id, img) 
  else
    lv.insert_image(img) 
    $density_map_image_id = img.id
  end
end

by something like:

if nrows * ncolumns > 0
  density = 0
  data.each { |d| density += d }
  density /= (nrows * ncolumns)
  RBA::MessageBox::info("Density", "The total density is #{'%.1f' % (density*100)}%", RBA::MessageBox::b_ok) 
end

Regards,

Matthias

okguy

Matthias,

Sorry, for my bad explanation.
It works fine, even on a very large dimension.

Thank you very much!
OkGuy

mikamar

Hi-

I tried using it.
I get completely black map in all my attempts

When I modify macro to only report the total density, I get a result different than zero which seems to make sense.
So i don't understand why the density map is all black.

Itamar

Matthias

Hi Itamar,

the script should print out the densities in percent on the console (if you run it in the macro editor for example). What are those numbers? Do they make sense?

If yes, then it's somehow related to the data mapping inside the image.

Matthias

mikamar

Hi Matthias-

I think the problem is when doing density for via layers which have inherent low density.
When I ran it for poly i got a valid map.
Maybe it's possible to do the color scaling based on relative layer density and not based on fixed 0-100% density

This is the printout for VIA2:
Running row y=0.0..500.0 1 1 1 1 1 1 2 0
Running row y=500.0..1000.0 1 1 0 0 0 0 1 0
Running row y=1000.0..1500.0 1 0 0 1 0 0 0 0
Running row y=1500.0..2000.0 0 0 0 1 0 0 0 0
Running row y=2000.0..2500.0 1 0 0 1 0 0 0 0
Running row y=2500.0..3000.0 1 0 0 1 0 0 0 0
Running row y=3000.0..3500.0 1 0 0 1 0 0 1 0
Running row y=3500.0..4000.0 0 0 0 0 0 0 0 0

This is the printout for poly:
Running row y=0.0..500.0 19 19 31 33 29 34 36 0
Running row y=500.0..1000.0 19 24 20 40 27 45 42 0
Running row y=1000.0..1500.0 41 41 40 26 42 48 41 0
Running row y=1500.0..2000.0 47 51 51 34 51 52 49 0
Running row y=2000.0..2500.0 45 49 49 32 48 49 47 0
Running row y=2500.0..3000.0 46 49 49 33 49 50 47 0
Running row y=3000.0..3500.0 40 40 42 34 47 40 39 0
Running row y=3500.0..4000.0 1 0 0 2 2 0 0 0

Itamar

Matthias

Sure that's possible. Just exploit the power of Ruby:

# normalize to 0..1
max = data.max
max < 1e-10 && raise("Overall density too low")
data = data.collect { |d| d / max }

# now the original code follows:
if nrows * ncolumns > 0
  img = RBA::Image.new(ncolumns, nrows, data)
  img.pixel_width = xw
  img.pixel_height = yw
  img.trans = RBA::DCplxTrans.new(1.0, 0.0, false, RBA::DPoint.new(xmin, ymin))
  if $density_map_image_id
    lv.replace_image($density_map_image_id, img) 
  else
    lv.insert_image(img) 
    $density_map_image_id = img.id
  end
end

Matthias

dmohata

Hi Matthias,

Thank you for the free platform and the discussion forum.
I have a problem where I need to compute density of different layers and then scale them and merge them into one image. Is this something you think the above code can be modified to?
Sorry I am a beginner with Ruby coding so asking for your help with the development.

Also, when I used the above code, the image is in grey scale, is it possible to convert to visible color scale.

Thanks,
Dheeraj

davidnhutch

Something to help you get started.

The line that makes the image is

img = RBA::Image.new(ncolumns, nrows, data)

Look at the class Image under Public Constructors. The above is using the fourth version of the constructor listed there, because it has three arguments: ncolumns, nrows, data.

When you see "double[]" in front of one of the Public Constructor's arguments, it means that value should be an array, which "data" is.

Anyway instead of the fourth version of the Public Constructor, you want one of the other versions where you put red, green, and blue separately. For example for the result to be red you need first to make an array of zeros the same size as your data array, then to create the image with r, g, b arguments:

zeros = data.map{|i| i.map{|j| 0 }} # There's probably a better way to do this... Basically just make an array of zeros
img = RBA::Image.new(ncolumns, nrows, data, zeros, zeros)

If you want it to be purple you can do:

halfs = data.map{|i| i.map{|j| j/2.0 }}
zeros = data.map{|i| i.map{|j| 0 }}
img = RBA::Image.new(ncolumns, nrows, halfs, zeros, halfs)

To get a red image overlaid on a purple one, just call the first one img1 and the second one img2 or something. I haven't tried this but I guess it should work.

HTH

P.S. This code has been improved slightly and given a GUI: See TRT. After downloading it per the instructions, look under TRT menu > Layout > Density calculation.

Matthias

Hi David,

thanks a lot - perfect intro :-)

I realize, that the image cannot be saved to PNG directly. If you want that you'll need to traverse the image pixel by pixel and generate a QImage object. This can be saved. I am overhauling the Ruby/Python API currently and I guess that interfaces to QImage and/or PNG/JPG files make some sense.

Thanks, Matthias

dmohata

Thank you so much David. I was able to change the color:)
TRT code did not work for density calculations though:(

Mathias,
I have one more issue with using the code (just cumbersome).
Currently I am running the code by clicking macros and then execute. However, I would like the "layer density" tab to appear right when I open any layout. Could you suggest how to modify the code to achieve this?

Thanks,
Dheeraj

garry

Hi, I recently used the calculation of the density map and i find that Region#area() calculate faster with this problem. This is my test code(python):

#compare two method to calculate the density map
import pya
import time

pixel = 9
pattern_layer = 1000

gds_file = "C:/Users/Aeolus/Desktop/testcase1.gds"
layout = pya.Layout.new()
layout.read(gds_file)
top_cell = layout.cell('top')

bbox = top_cell.bbox()
clip_boxes = [pya.Box.new(x, y, x+pixel, y+pixel) for x in range(bbox.p1.x, bbox.p2.x, pixel) for y in range(bbox.p1.y, bbox.p2.y, pixel)]

#method one
start1 = time.clock()

clip_layout = pya.Layout.new()
cell_indexes = layout.multi_clip_into(top_cell.cell_index(), clip_layout, clip_boxes)
#end1 = time.clock()
for cell_index in cell_indexes:
    cell = clip_layout.cell(cell_index)
    area = 0
    for shape in cell.each_shape(clip_layout.layer(pattern_layer, 0)):
        area += shape.area()
    density = area/(pixel^2)

end1 = time.clock()

#method two
start2 = time.clock()

region = pya.Region.new(top_cell.shapes(layout.layer(pattern_layer, 0)))
bbox = region.bbox()

for x in range(bbox.p1.x, bbox.p2.x, pixel):
    for y in range(bbox.p1.y, bbox.p2.y, pixel):
        clip_box = pya.Box.new(x, y, x+pixel, y+pixel)
        area = region.area(clip_box)
        density = area/(pixel^2)

end2 = time.clock()

print("method one time cost:%f"%(end1-start1))
print("method two time cost:%f"%(end2-start2))

the result of my gds file is:

method one time cost:4.806324
method two time cost:1.133012

you can test the code with yourself by replacing the "gds_file", "pattern_layer" and top cell name.

Matthias

Hi Garry,

Thank you for your inputs.

Eventually the best method probably depends on the nature of the layout and the pixel size. The most efficient method should be the tiling processor which provides an engine for doing operations on a regular tile grid.

The advantage of this approach is that it not only picks the best method but also supports distribution the computation on multiple CPUs.

I have code for this somewhere. I just need to dig this out.

Matthias

garry

Hi, Matthias,

Thank you for your reply. I really happy and look forward to getting your code.

Thanks,
Garry

laurent_c

The main difference I have foreseen between the 2 methods is that the first one (multi_clip) will take in account the shapes of the subcells, while the second one (region) will only consider the shapes of the topcell.
By the way, Matthias, Is there a way with region to check the shapes from the subcells ?
Regards,
Laurent

Matthias

@garry: Here is a sample script using the TilingProcessor for density computation. It employs an Image object to store the density values:

import pya

pixel_size = 10

threads = 4

ly = pya.CellView.active().layout()
# takes density from layer 100
li = ly.layer(100, 0)
si = ly.top_cell().begin_shapes_rec(li)

tp = pya.TilingProcessor()

# this example uses an image for output
img = pya.Image()

tp.tile_size(pixel_size, pixel_size)
tp.output("res", img)
tp.input("input", si)
tp.dbu = ly.dbu
tp.threads = threads
# inside tiling processor compute the density and output to the image
tp.queue("_tile && (var d = to_f(input.area(_tile.bbox)) / to_f(_tile.bbox.area);_output(res, d))")
tp.execute("Density map")

# place the image into the view - it's already properly configured by the 
# tiling processor
# you could use the "data mapping" feature of the images for creating a heat map
# for example.
pya.LayoutView.current().insert_image(img)

@laurent: when you construct a region from a RecursiveShapeIterator, it will work an all shapes of a layer, including the ones from subcells.

Matthias

14darcia

Hi Matthias,

I'm using the above scrip to get the density map which works great! However, as there's no way to export the image, I exported it as a string to a text file. It actually works, for the most part but I'm a little confused by the metadata within the "img" variable. Specifically:

mono:matrix=(100,0,1585.58) (0,100,-904.316) (0,0,1);min_value=0;max_value=1;is_visible=true;z_position=0;brightness=0;contrast=0;gamma=1;red_gain=1;green_gain=1;blue_gain=1;color_mapping=[0,'#000000';1,'#ffffff';];width=30;height=51

A lot of this is self explanatory, but I'm just looking for the positional data. I can tell the (100,0,1585.58) (0,100,-904.316) contains the size of the unit cell, as well as some other positional data. But what exactly does 1585 and -904 correspond to? It doesn't seem to exactly match the position of any of the squares.

Thanks,
Adam

Matthias

Hi Adam,

"matrix" is a transformation matrix which transforms pixel coordinates (x=0,1,2..., y=0,1,2...) into micrometer values. The format lists the rows of the matrix. It's a 3x3 matrix for a general 2d transformation with homogeneous coordinates (affine plus perspective). "Homogeneous" coordinates are formed by adding a 1 as the third coordinates (x, y, 1) before feeding this into the multiplication with the 3x3 matrix. The concept is usually explained in terms of 4x4 matrices for 3d point transformations - for 2d, the z component can be thought to be 0 and the matrix is reduced to 3x3.

In your case, there is no perspective transformation and the third row is (0,0,1). The scale and rotate part of the matrix is (100,0)/(0,100) which means no rotation but a scaling with 100. This means the pixel dimension is 100µm. Finally the position of the center of the image is 1585.58,-904.316µm. This should correspond to the center your tiling arena. A little rounding may be involved because the arena will be fit an integer number of pixels.

In general, the string is rather an internal format used for serializing the image information in session files. It's likely that there are undocumented details and the format will be modified if necessary. If you need to export this information to another program, a better idea might be to use the Python API of the Image object and write the pixel values to a file in the format of your choice.

Matthias

14darcia

Thank you so much for the detailed explanation!

shijie

Hello Matthias,
This density calculation is very covenient for global density check in Win, but there is something strange about the results, when I do the density map calculation with insufficient RAM, the results is not correct compared with the results from Calibre nmDRC, unless I reboot my pc and try it again.

Best Regards,
Shijie

Matthias

Hi,

maybe you need to check the log for "bad_alloc" errors (File/Log Viewer).

Basically errors during the computation should make the density computation fail, but I cannot rule out entirely that errors get unnoticed. In this case, a tile may get skipped.

But this means your pixel size is pretty big and the number of polygons per pixel gets kind of high.

If possible, you could try to reduce the pixel size. Or you need to reduce the number of threads to avoid parallel computation. Or to plug in more memory and switch to 64bit if you're not already there yet.

Matthias

tomas2004

Hello Matthias,

I'm looking into the "heat map part" and added the code below to create a kind of rainbow map ranging from blue (0) <> cyan <> green <> yellow <> red (1):

datamapping = ImageDataMapping.new

datamapping.add_colormap_entry(0.0, 0x0000ff) # blue

datamapping.add_colormap_entry(0.25, 0x00ffff) # cyan

datamapping.add_colormap_entry(0.5, 0x00ff00) # green

datamapping.add_colormap_entry(0.75, 0xffff00) # yellow

datamapping.add_colormap_entry(1.0, 0xff0000) # red

img.data_mapping = datamapping

Is this the way to do it?

I noticed that the color map entry for 0.0 does not show up as blue, but it is black instead. Is that intended?

Is there some more information about color, lcolor and rcolor besides https://www.klayout.de/doc-qt5/code/class_ImageDataMapping.html#method15?

Cheers and enjoy the x-mas holidays!

Tomas

Matthias

Hi @tomas2004,

I'm afraid there isn't much more documentation than the provided one - apart from the code, which of course is also documentation. That's Open Source

Your data mapping is basically correct and 0 should show up a blue. Maybe you can save the image in KLayout's image format (XML) using "img.write(...)" and share a test case?

I can't confirm there is something wrong with the color mapping, but I see some issues when I enable to color corrections in that case (brightness, contrast, gamma). Maybe there is actually some bug. I have created a ticket: https://github.com/KLayout/klayout/issues/1959

Another side note: the values for the color map entries are normalized - 0.0 for low end and 1.0 for high end. Low and high end need to be se by "min_value" and "max_value".

Thanks,

Matthias